Mercurial > repos > pointsource_algs / file comparison
comparison: src/sliding_pdps.rs
src/sliding_pdps.rs
- changeset 95
- d9318317ccac
- parent 93
- 8f1e57821d76
- child 96
- 5ba25d740dfc
equal
deleted
inserted
replaced
| 94:88d994f2d63d | 95:d9318317ccac |
|---|---|
| 136 H: Conjugable<Y, F, Codomain = F>, | 136 H: Conjugable<Y, F, Codomain = F>, |
| 137 for<'b> H::Conjugate<'b>: Prox<Y>, | 137 for<'b> H::Conjugate<'b>: Prox<Y>, |
| 138 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>, | 138 Plot: Plotter<P::ReturnMapping, S, DiscreteMeasure<U, F>>, |
| 139 { | 139 { |
| 140 // Check parameters | 140 // Check parameters |
| 141 /*ensure!( | 141 ensure!( |
| 142 config.τ0 > 0.0 | 142 config.τ0 > 0.0 && config.σp0 > 0.0, |
| 143 && config.τ0 < 1.0 | 143 //&& config.τ0 < 1.0 && config.σp0 < 1.0 |
| 144 && config.σp0 > 0.0 | |
| 145 && config.σp0 < 1.0 | |
| 146 && config.σd0 > 0.0 | |
| 147 && config.σp0 * config.σd0 <= 1.0, | |
| 148 "Invalid step length parameters" | 144 "Invalid step length parameters" |
| 149 );*/ | 145 ); |
| 150 config.transport.check()?; | 146 config.transport.check()?; |
| 151 | 147 |
| 152 // Initialise iterates | 148 // Initialise iterates |
| 153 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); | 149 let mut μ = μ0.unwrap_or_else(|| DiscreteMeasure::new()); |
| 154 let mut γ = Transport::new(); | 150 let mut γ = Transport::new(); |
| 163 let is_fb = nKz == 0.0; | 159 let is_fb = nKz == 0.0; |
| 164 let idOpZ = IdOp::new(); | 160 let idOpZ = IdOp::new(); |
| 165 let opKz_adj = opKz.adjoint(); | 161 let opKz_adj = opKz.adjoint(); |
| 166 let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?; | 162 let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?; |
| 167 | 163 |
| 168 // We need to satisfy | 164 let (τ, σ_p, σ_d, κ) = if is_fb { |
| 169 // | 165 (config.τ0 / l, config.σp0 / l_z, 0.0, 0.0) |
| 170 // τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1 | 166 } else { |
| 171 // ^^^^^^^^^^^^^^^^^^^^^^^^^ | 167 ensure!( |
| 172 // with 1 > σ_p L_z and 1 > τ L. | 168 //config.τ0 < 1.0 && |
| 173 // | 169 config.σp0 < 1.0 && config.σd0 > 0.0 && config.σp0 * config.σd0 <= 1.0, |
| 174 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition | 170 "Invalid step length parameters" |
| 175 // ^^^^^ < 1. then we solve τ from the rest. | 171 ); |
| 176 // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below. | 172 // We need to satisfy |
| 177 let σ_d = if is_fb { 0.0 } else { config.σd0 / nKz }; | 173 // |
| 178 let σ_p = config.σp0 / (l_z + config.σd0 * nKz); | 174 // τσ_dM(1-σ_p L_z)/(1 - τ L) + [σ_p L_z + σ_pσ_d‖K_z‖^2] < 1 |
| 179 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0} | 175 // ^^^^^^^^^^^^^^^^^^^^^^^^^ |
| 180 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L) | 176 // with 1 > σ_p L_z and 1 > τ L. |
| 181 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0}) | 177 // |
| 182 let φ = 1.0 - config.σp0; | 178 // To do so, we first solve σ_p and σ_d from standard PDPS step length condition |
| 183 let a = 1.0 - σ_p * l_z; | 179 // ^^^^^ < 1. then we solve τ from the rest. |
| 184 let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l); | 180 // If opKZ is the zero operator, then we set σ_d = 0 for τ to be calculated correctly below. |
| 185 let ψ = 1.0 - τ * l; | 181 let σ_d = config.σd0 / nKz; |
| 186 let β = σ_p * config.σd0 * nKz / a; // σ_p * σ_d * (nKz * nK_z) / a; | 182 let σ_p = config.σp0 / (l_z + config.σd0 * nKz); |
| 187 ensure!(β < 1.0); | 183 // Observe that = 1 - ^^^^^^^^^^^^^^^^^^^^^ = 1 - σ_{p,0} |
| 188 // Now we need κ‖K_μ(π_♯^1 - π_♯^0)γ‖^2 ≤ (1/θ - τ[ℓ_F + ℓ]) ∫ c_2 dγ for κ defined as: | 184 // We get the condition τσ_d M (1-σ_p L_z) < (1-σ_{p,0})*(1-τ L) |
| 189 let κ = τ * σ_d * ψ / ((1.0 - β) * ψ - τ * σ_d * bigM); | 185 // ⟺ τ [ σ_d M (1-σ_p L_z) + (1-σ_{p,0}) L ] < (1-σ_{p,0}) |
| 190 // The factor two in the manuscript disappears due to the definition of 𝚹 being | 186 let φ = 1.0 - config.σp0; |
| 191 // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2. | 187 let a = 1.0 - σ_p * l_z; |
| 188 let τ = config.τ0 * φ / (σ_d * bigM * a + φ * l); | |
| 189 let ψ = 1.0 - τ * l; | |
| 190 let β = σ_p * config.σd0 * nKz / a; // σ_p * σ_d * (nKz * nK_z) / a; | |
| 191 ensure!(β < 1.0); | |
| 192 // Now we need κ‖K_μ(π_♯^1 - π_♯^0)γ‖^2 ≤ (1/θ - τ[ℓ_F + ℓ]) ∫ c_2 dγ for κ defined as: | |
| 193 // The factor two in the manuscript disappears due to the definition of 𝚹 being | |
| 194 // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2. | |
| 195 let κ = τ * σ_d * ψ / ((1.0 - β) * ψ - τ * σ_d * bigM); | |
| 196 (τ, σ_p, σ_d, κ) | |
| 197 }; | |
| 192 | 198 |
| 193 let mut τθ_or_adaptive = prox_penalty.get_transport_steplength( | 199 let mut τθ_or_adaptive = prox_penalty.get_transport_steplength( |
| 194 f.curvature_bound_components(config.guess), | 200 f.curvature_bound_components(config.guess), |
| 195 &config.transport, | 201 &config.transport, |
| 196 0.0, | 202 0.0, |
