# HG changeset patch # User Tuomo Valkonen # Date 1784439006 -7200 # Node ID 677a5fd1b0144ec7f70b5e5f69ff2d64960b1d27 # Parent df92e78cc3f4bbe6832f192670e51523d01d25b4 Added an extra step heuristic diff -r df92e78cc3f4 -r 677a5fd1b014 src/fb.rs --- a/src/fb.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/fb.rs Sun Jul 19 07:30:06 2026 +0200 @@ -188,17 +188,14 @@ // Run the algorithm for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) { + let maybe_μ_base = config.merge_now(&state).then(|| μ.clone()); + let μ_base_len = μ.len(); + // Calculate smooth part of surrogate model. - // TODO: optimise τ to be applied to residual. let mut τv = f.differential(&μ) * τ; - // Save current base point for merge - let μ_base_len = μ.len(); - let maybe_μ_base = config.merge_now(&state).then(|| μ.clone()); - - // Insert and reweigh - let (maybe_d, _within_tolerances) = prox_penalty - .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + // Do spike insertiona nd finite-dimensional weight optimisation + prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; stats.inserted += μ.len() - μ_base_len; @@ -218,12 +215,22 @@ stats.pruned += prune_with_stats(&mut μ); + // Do extra weight optimisation step heuristic + for _ in 0..config.extra_weight_optimisation_steps { + τv = f.differential(&μ) * τ; + prox_penalty.reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + } + + if config.extra_weight_optimisation_steps > 0 { + stats.pruned += prune_with_stats(&mut μ); + } + let iter = state.iteration(); stats.this_iters += 1; // Give statistics if needed state.if_verbose(|| { - plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ); + plotter.plot_spikes(iter, None, Some(&τv), &μ); full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new())) }); @@ -292,14 +299,12 @@ // Run the algorithm for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) { + let μ_base_len = μ.len(); // Calculate smooth part of surrogate model. let mut τv = f.differential(&μ) * τ; - let μ_base_len = μ.len(); - - // Insert new spikes and reweigh - let (maybe_d, _within_tolerances) = prox_penalty - .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + // Do spike insertiona nd finite-dimensional weight optimisation + prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; stats.inserted += μ.len() - μ_base_len; @@ -310,6 +315,12 @@ warned_merging = true; } + // Do extra weight optimisation step heuristic + for _ in 0..config.extra_weight_optimisation_steps { + τv = f.differential(&μ) * τ; + prox_penalty.reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + } + // Update inertial prameters let λ_prev = λ; λ = 2.0 * λ_prev / (λ_prev + (4.0 + λ_prev * λ_prev).sqrt()); @@ -334,7 +345,7 @@ // Give statistics if needed state.if_verbose(|| { - plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ_prev); + plotter.plot_spikes(iter, None, Some(&τv), &μ_prev); full_stats(&μ_prev, ε, std::mem::replace(&mut stats, IterInfo::new())) }); diff -r df92e78cc3f4 -r 677a5fd1b014 src/forward_pdps.rs --- a/src/forward_pdps.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/forward_pdps.rs Sun Jul 19 07:30:06 2026 +0200 @@ -197,15 +197,29 @@ }; let mut stats = IterInfo::new(); + let zy_step = |mut z, mut y, mut z_tmp| { + // Do z variable primal update + opKz_adj.apply_add(&mut z_tmp, &y); + z_tmp.axpy(1.0, &z, -σ_p); + let z_new = fnR.prox(σ_p, z_tmp); + //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z); + // Do dual update + z.axpy(1.0 + ω, &z_new, -ω); + //z = (z - &z_new) * (-ω) + &z_new; + opKz.gemv(&mut y, σ_d, z, 1.0); + let y = starH.prox(σ_d, y); + (z_new, y) + }; + // Run the algorithm for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { // Calculate initial transport - let Pair(v, mut z_tmp) = f.differential(Pair(&μ, &z)); + let Pair(v, z_tmp) = f.differential(Pair(&μ, &z)); let mut τv = v * τ; let μ_base = μ.clone(); // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. - let (maybe_d, _within_tolerances) = prox_penalty.insert_and_reweigh( + prox_penalty.insert_and_reweigh( &mut μ, &mut τv, τ, @@ -216,6 +230,8 @@ &mut stats, )?; + (z, y) = zy_step(z, y, z_tmp); + stats.inserted += μ.len() - μ_base.len(); // Merge spikes. @@ -241,17 +257,27 @@ // Prune spikes with zero weight. stats.pruned += prune_with_stats(&mut μ); - // Do z variable primal update - opKz_adj.apply_add(&mut z_tmp, &y); - z_tmp.axpy(1.0, &z, -σ_p); - let z_new = fnR.prox(σ_p, z_tmp); - //let z_new = fnR.prox(σ_p, z_tmp * (-σ_p) + &z); - // Do dual update - z.axpy(1.0 + ω, &z_new, -ω); - //z = (z - &z_new) * (-ω) + &z_new; - opKz.gemv(&mut y, σ_d, z, 1.0); - y = starH.prox(σ_d, y); - z = z_new; + for _ in 0..config.insertion.extra_weight_optimisation_steps { + // Do extra weight optimisation step heuristic + let Pair(v, _z_tmp) = f.differential(Pair(&μ, &z)); + τv = v * τ; + prox_penalty.reweigh( + &mut μ, + &mut τv, + τ, + ε, + &config.insertion, + ®, + &state, + &mut stats, + )?; + //(z, y) = zy_step(z, y, z_tmp); + } + + if config.insertion.extra_weight_optimisation_steps > 0 { + // Prune spikes with zero weight. + stats.pruned += prune_with_stats(&mut μ); + } // Update step length parameters // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); @@ -261,7 +287,7 @@ stats.this_iters += 1; state.if_verbose(|| { - plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ); + plotter.plot_spikes(iter, None, Some(&τv), &μ); full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) }); diff -r df92e78cc3f4 -r 677a5fd1b014 src/lib.rs --- a/src/lib.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/lib.rs Sun Jul 19 07:30:06 2026 +0200 @@ -238,6 +238,10 @@ /// Scaling heuristic factor pub scaling_heuristic: Option, + #[arg(long)] + /// Extra finite dimensional steps to take + pub extra_weight_optimisation_steps: Option, + #[arg(long, value_names = &["ε", "θ", "p"])] /// Set the tolerance to ε_k = ε/(1+θk)^p pub tolerance: Option>, diff -r df92e78cc3f4 -r 677a5fd1b014 src/pdps.rs --- a/src/pdps.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/pdps.rs Sun Jul 19 07:30:06 2026 +0200 @@ -206,8 +206,7 @@ let μ_base = μ.clone(); // Insert and reweigh - let (maybe_d, _within_tolerances) = prox_penalty - .insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; + prox_penalty.insert_and_reweigh(&mut μ, &mut τv, τ, ε, config, ®, &state, &mut stats)?; // Prune and possibly merge spikes if config.merge_now(&state) { @@ -234,7 +233,7 @@ stats.this_iters += 1; state.if_verbose(|| { - plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv), &μ); + plotter.plot_spikes(iter, None, Some(&τv), &μ); full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new())) }); diff -r df92e78cc3f4 -r 677a5fd1b014 src/prox_penalty.rs --- a/src/prox_penalty.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/prox_penalty.rs Sun Jul 19 07:30:06 2026 +0200 @@ -65,8 +65,9 @@ /// Iterations between merging heuristic tries pub merge_every: usize, - // /// Save $μ$ for postprocessing optimisation - // pub postprocessing : bool + + /// Additional weight optimisation steps + pub extra_weight_optimisation_steps: usize, } #[replace_float_literals(F::cast_from(literal))] @@ -85,7 +86,7 @@ fitness_merging: false, merge_every: 10, merge_tolerance_mult: 2.0, - // postprocessing : false, + extra_weight_optimisation_steps: 0, } } } @@ -118,6 +119,7 @@ Reg: RegTerm, Domain: Space + Clone, { + /// Unused, but required as a plotter parametrisation. type ReturnMapping: Mapping; /// Returns the type of this proximality penalty @@ -126,8 +128,7 @@ /// Insert new spikes into `μ` to approximately satisfy optimality conditions /// with the forward step term fixed to `τv`. /// - /// May return `τv + w` for `w` a subdifferential of the regularisation term `reg`, - /// as well as an indication of whether the tolerance bounds `ε` are satisfied. + /// Returns an indication of whether the tolerance bounds `ε` are satisfied. /// /// `τv` is mutable to allow [`alg_tools::bounds::MinMaxMapping`] optimisation to /// refine data. Actual values of `τv` are not supposed to be mutated. @@ -144,7 +145,23 @@ reg: &Reg, state: &AlgIteratorIteration, stats: &mut IterInfo, - ) -> DynResult<(Option, bool)> + ) -> DynResult + where + I: AlgIterator; + + /// A variant of [`insert_and_reweigh`] that only does finite-dimensional weight optimisation, + /// without inserting spikes. + fn reweigh( + &self, + μ: &mut DiscreteMeasure, + τv: &mut PreadjointCodomain, + τ: F, + ε: F, + config: &InsertionConfig, + reg: &Reg, + state: &AlgIteratorIteration, + stats: &mut IterInfo, + ) -> DynResult<()> where I: AlgIterator; diff -r df92e78cc3f4 -r 677a5fd1b014 src/prox_penalty/radon_squared.rs --- a/src/prox_penalty/radon_squared.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/prox_penalty/radon_squared.rs Sun Jul 19 07:30:06 2026 +0200 @@ -52,7 +52,7 @@ reg: &Reg, _state: &AlgIteratorIteration, stats: &mut IterInfo, - ) -> DynResult<(Option, bool)> + ) -> DynResult where I: AlgIterator, { @@ -63,7 +63,25 @@ let violation = reg.find_tolerance_violation(τv, τ, ε, skip_by_rough_check, config); reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats); - Ok((None, true)) + Ok(true) + } + + fn reweigh( + &self, + μ: &mut DiscreteMeasure, + τv: &mut M, + τ: F, + ε: F, + config: &InsertionConfig, + reg: &Reg, + _state: &AlgIteratorIteration, + stats: &mut IterInfo, + ) -> DynResult<()> + where + I: AlgIterator, + { + reg.solve_oc_radonsq(μ, τv, τ, ε, None, config, stats); + Ok(()) } fn merge_spikes( diff -r df92e78cc3f4 -r 677a5fd1b014 src/prox_penalty/wave.rs --- a/src/prox_penalty/wave.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/prox_penalty/wave.rs Sun Jul 19 07:30:06 2026 +0200 @@ -52,7 +52,7 @@ reg: &Reg, state: &AlgIteratorIteration, stats: &mut IterInfo, - ) -> DynResult<(Option, bool)> + ) -> DynResult where I: AlgIterator, { @@ -70,7 +70,7 @@ // Add points to support until within error tolerance or maximum insertion count reached. let mut count = 0; - let (within_tolerances, d) = 'insertion: loop { + let within_tolerances = 'insertion: loop { if μ.len() > 0 { // Form finite-dimensional subproblem. The subproblem references to the original μ^k // from the beginning of the iteration are all contained in the immutable c and g. @@ -115,13 +115,13 @@ // Find a spike to insert, if needed let (ξ, _v_ξ, in_bounds) = match reg.find_tolerance_violation(&mut d, τ, ε, skip_by_rough_check, config) { - None => break 'insertion (true, d), + None => break 'insertion true, Some(res) => res, }; // Break if maximum insertion count reached if count >= max_insertions { - break 'insertion (in_bounds, d); + break 'insertion in_bounds; } // No point in optimising the weight here; the finite-dimensional algorithm is fast. @@ -139,7 +139,55 @@ println!("{}", err.red()); } - Ok((Some(d), within_tolerances)) + Ok(within_tolerances) + } + + fn reweigh( + &self, + μ: &mut DiscreteMeasure, + τv: &mut M, + τ: F, + ε: F, + config: &InsertionConfig, + reg: &Reg, + _state: &AlgIteratorIteration, + stats: &mut IterInfo, + ) -> DynResult<()> + where + I: AlgIterator, + { + if μ.len() > 0 { + let op𝒟norm = self.opnorm_bound(Radon, Linfinity)?; + let ω0 = self.apply(&*μ); + + // Form finite-dimensional subproblem. The subproblem references to the original μ^k + // from the beginning of the iteration are all contained in the immutable c and g. + // TODO: observe negation of -τv after switch from minus_τv: finite-dimensional + // problems have not yet been updated to sign change. + let à = self.findim_matrix(μ.iter_locations()); + let g̃ = DVector::from_iterator( + μ.len(), + μ.iter_locations() + .map(|ζ| ω0.apply(ζ) - τv.apply(ζ)) + .map(F::to_nalgebra_mixed), + ); + let mut x = μ.masses_dvector(); + + // The gradient of the forward component of the inner objective is C^*𝒟Cx - g̃. + // We have |C^*𝒟Cx|_2 = sup_{|z|_2 ≤ 1} ⟨z, C^*𝒟Cx⟩ = sup_{|z|_2 ≤ 1} ⟨Cz|𝒟Cx⟩ + // ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟Cx|_∞ ≤ sup_{|z|_2 ≤ 1} |Cz|_ℳ |𝒟| |Cx|_ℳ + // ≤ sup_{|z|_2 ≤ 1} |z|_1 |𝒟| |x|_1 ≤ sup_{|z|_2 ≤ 1} n |z|_2 |𝒟| |x|_2 + // = n |𝒟| |x|_2, where n is the number of points. Therefore + let Ã_normest = op𝒟norm * F::cast_from(μ.len()); + + // Solve finite-dimensional subproblem. + stats.inner_iters += reg.solve_findim(&Ã, &g̃, τ, &mut x, Ã_normest, ε, config); + + // Update masses of μ based on solution of finite-dimensional subproblem. + μ.set_masses_dvector(&x); + } + + Ok(()) } fn merge_spikes( diff -r df92e78cc3f4 -r 677a5fd1b014 src/run.rs --- a/src/run.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/run.rs Sun Jul 19 07:30:06 2026 +0200 @@ -123,6 +123,9 @@ final_merging: cli.final_merging.unwrap_or(g.final_merging), fitness_merging: cli.fitness_merging.unwrap_or(g.fitness_merging), inner: override_inner(g.inner), + extra_weight_optimisation_steps: cli + .extra_weight_optimisation_steps + .unwrap_or(g.extra_weight_optimisation_steps), tolerance: cli .tolerance .as_ref() @@ -797,7 +800,6 @@ DefaultAlgorithm::SlidingFB, DefaultAlgorithm::RadonSlidingFB, DefaultAlgorithm::FW, - DefaultAlgorithm::FWRelax, DefaultAlgorithm::PDPS, ], DataTermType::L1 => &[DefaultAlgorithm::PDPS], diff -r df92e78cc3f4 -r 677a5fd1b014 src/sliding_fb.rs --- a/src/sliding_fb.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/sliding_fb.rs Sun Jul 19 07:30:06 2026 +0200 @@ -1041,7 +1041,7 @@ // Solve finite-dimensional subproblem several times until the dual variable for the // regularisation term conforms to the assumptions made for the transport above. - let (maybe_d, _within_tolerances, mut τv̆, μ̆) = 'adapt_transport: loop { + let (mut τv̆, μ̆) = 'adapt_transport: loop { // Set initial guess for μ=μ^{k+1}. γ.μ̆_into(&mut μ); let μ̆ = μ.clone(); @@ -1054,7 +1054,7 @@ let mut τv̆ = f.differential(&μ̆) * τ; // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. - let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh( + prox_penalty.insert_and_reweigh( &mut μ, &mut τv̆, τ, @@ -1081,7 +1081,7 @@ &config.insertion.refinement, &mut attempts, ) { - break 'adapt_transport (maybe_d, within_tolerances, τv̆, μ̆); + break 'adapt_transport (τv̆, μ̆); } stats.get_transport_mut().readjustment_iters += 1; @@ -1103,20 +1103,39 @@ ®, Some(|μ̃: &RNDM| f.apply(μ̃)), ); - if m > 0 { - stats.merged += m; - v = f.differential(&μ); - } + //if m > 0 { + stats.merged += m; + //v = f.differential(&μ); + //} } γ.prune_compat(&mut μ, &mut stats); + // Do extra weight optimisation step heuristic + for _ in 1..config.insertion.extra_weight_optimisation_steps { + τv̆ = f.differential(&μ) * τ; + prox_penalty.reweigh( + &mut μ, + &mut τv̆, + τ, + ε, + &config.insertion, + ®, + &state, + &mut stats, + )?; + } + + if config.insertion.extra_weight_optimisation_steps > 0 { + γ.prune_compat(&mut μ, &mut stats); + } + let iter = state.iteration(); stats.this_iters += 1; // Give statistics if requested state.if_verbose(|| { - plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv̆), &μ); + plotter.plot_spikes(iter, None, Some(&τv̆), &μ); full_stats(&μ, ε, std::mem::replace(&mut stats, IterInfo::new())) }); diff -r df92e78cc3f4 -r 677a5fd1b014 src/sliding_pdps.rs --- a/src/sliding_pdps.rs Thu Jul 16 19:51:51 2026 +0300 +++ b/src/sliding_pdps.rs Sun Jul 19 07:30:06 2026 +0200 @@ -242,7 +242,7 @@ // Solve finite-dimensional subproblem several times until the dual variable for the // regularisation term conforms to the assumptions made for the transport above. - let (maybe_d, _within_tolerances, mut τv̆, z_new, μ̆) = 'adapt_transport: loop { + let (mut τv̆, z_new, μ̆) = 'adapt_transport: loop { // Set initial guess for μ=μ^{k+1}. γ.μ̆_into(&mut μ); let μ̆ = μ.clone(); @@ -253,7 +253,7 @@ // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. - let (maybe_d, within_tolerances) = prox_penalty.insert_and_reweigh( + prox_penalty.insert_and_reweigh( &mut μ, &mut τv̆, τ, @@ -286,12 +286,20 @@ &config.insertion.refinement, &mut attempts, ) { - break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new, μ̆); + break 'adapt_transport (τv̆, z_new, μ̆); } }; γ.get_transport_stats(&mut stats, &μ); + // Do dual update + // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] + z.axpy(1.0 + ω, &z_new, -ω); + //z = (z - &z_new) * (-ω) + &z_new; + opKz.gemv(&mut y, σ_d, z, 1.0); + y = starH.prox(σ_d, y); + z = z_new; + // Merge spikes. // This crucially expects the merge routine to be stable with respect to spike locations, // and not to performing any pruning. That is be to done below simultaneously for γ. @@ -313,13 +321,34 @@ γ.prune_compat(&mut μ, &mut stats); - // Do dual update - // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] - z.axpy(1.0 + ω, &z_new, -ω); - //z = (z - &z_new) * (-ω) + &z_new; - opKz.gemv(&mut y, σ_d, z, 1.0); - y = starH.prox(σ_d, y); - z = z_new; + for _ in 0..config.insertion.extra_weight_optimisation_steps { + // Do extra weight optimisation step heuristic + let Pair(v, mut _z_tmp) = f.differential(Pair(&μ, &z)); + τv̆ = v * τ; + prox_penalty.reweigh( + &mut μ, + &mut τv̆, + τ, + ε, + &config.insertion, + ®, + &state, + &mut stats, + )?; + + // opKz_adj.apply_add(&mut z_tmp, &y); + // z_tmp.axpy(1.0, &z, -σ_p); + // let z_new = fnR.prox(σ_p, z_tmp); + // z.axpy(1.0 + ω, &z_new, -ω); + // //z = (z - &z_new) * (-ω) + &z_new; + // opKz.gemv(&mut y, σ_d, z, 1.0); + // y = starH.prox(σ_d, y); + // z = z_new; + } + + if config.insertion.extra_weight_optimisation_steps > 0 { + γ.prune_compat(&mut μ, &mut stats); + } // Update step length parameters // let ω = pdpsconfig.acceleration.accelerate(&mut τ, &mut σ, γ); @@ -329,7 +358,7 @@ stats.this_iters += 1; state.if_verbose(|| { - plotter.plot_spikes(iter, maybe_d.as_ref(), Some(&τv̆), &μ); + plotter.plot_spikes(iter, None, Some(&τv̆), &μ); full_stats(&μ, &z, ε, std::mem::replace(&mut stats, IterInfo::new())) });