Mercurial > repos > pointsource_algs / changeset
changeset
New simplified sliding a posteriori rule. Enable scaling heuristic from fixed `measure` crate. “Fast” spread Lipschitz factor fixes, etc.
Sun, 19 Jul 2026 07:34:39 +0200
- author
- Tuomo Valkonen <tuomov@iki.fi>
- date
- Sun, 19 Jul 2026 07:34:39 +0200
- changeset 72
- e9a460a0e638
- parent 71
- e2953ffd4e0b
- child 73
- 9c6432200aba
New simplified sliding a posteriori rule. Enable scaling heuristic from fixed `measure` crate. “Fast” spread Lipschitz factor fixes, etc.
--- a/.gitignore Fri May 15 14:40:02 2026 -0500 +++ b/.gitignore Sun Jul 19 07:34:39 2026 +0200 @@ -1,1 +1,12 @@ -.hgignore \ No newline at end of file +syntax:glob +out/ +test/ +target/ +debug_out/ +**/pointsource??_*.txt +flamegraph.svg +DEADJOE +**/*.orig +results_new/ +results_4000/ +results/
--- a/Cargo.lock Fri May 15 14:40:02 2026 -0500 +++ b/Cargo.lock Sun Jul 19 07:34:39 2026 +0200 @@ -33,7 +33,7 @@ [[package]] name = "alg_tools" -version = "0.4.1-dev" +version = "0.5.2-dev" dependencies = [ "anyhow", "colored", @@ -647,12 +647,11 @@ [[package]] name = "measures" -version = "0.1.0" +version = "0.2.0" dependencies = [ "alg_tools", "nalgebra", "numeric_literals", - "regex", "serde", ] @@ -817,7 +816,7 @@ [[package]] name = "pointsource_algs" -version = "3.0.2-dev" +version = "3.1.0-dev" dependencies = [ "GSL", "alg_tools", @@ -944,9 +943,9 @@ [[package]] name = "regex" -version = "1.11.1" +version = "1.12.4" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "b544ef1b4eac5dc2db33ea63606ae9ffcfac26c1416a2806ae0bf5f56b201191" +checksum = "f1292b7759ae1cb9ec195452d1390a074f0cd8541ab7a5a8c31cd6db45d4a6ba" dependencies = [ "aho-corasick", "memchr", @@ -956,9 +955,9 @@ [[package]] name = "regex-automata" -version = "0.4.9" +version = "0.4.14" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "809e8dc61f6de73b46c85f4c96486310fe304c434cfa43669d7b40f711150908" +checksum = "6e1dd4122fc1595e8162618945476892eefca7b88c52820e74af6262213cae8f" dependencies = [ "aho-corasick", "memchr", @@ -967,9 +966,9 @@ [[package]] name = "regex-syntax" -version = "0.8.5" +version = "0.8.11" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c" +checksum = "d6f6ff9a378485b298a5286656da665ba74413d36db0979633275d2e708145d4" [[package]] name = "rustc_version"
--- a/Cargo.toml Fri May 15 14:40:02 2026 -0500 +++ b/Cargo.toml Sun Jul 19 07:34:39 2026 +0200 @@ -1,6 +1,6 @@ [package] name = "pointsource_algs" -version = "3.0.2-dev" +version = "3.1.0-dev" edition = "2021" rust-version = "1.85" authors = ["Tuomo Valkonen <tuomov@iki.fi>"] @@ -22,13 +22,13 @@ categories = ["mathematics", "science", "computer-vision"] [dependencies.alg_tools] -version = "~0.4.1-dev" +version = "~0.5.2-dev" path = "../alg_tools" default-features = false features = ["nightly"] [dependencies.measures] -version = "~0.1.0" +version = "~0.2.0" path = "../measures" [dependencies] @@ -54,7 +54,7 @@ default = [] [build-dependencies] -regex = "~1.11.0" +regex = "~1.12.0" [profile.release] debug = true
--- a/README.md Fri May 15 14:40:02 2026 -0500 +++ b/README.md Sun Jul 19 07:34:39 2026 +0200 @@ -27,14 +27,14 @@ 1. Install the [Rust] infrastructure (including Cargo) with [rustup]. 2. Install a “nightly” release of the Rust compiler. With rustup, installed in the previous step, this can be done with - ```console - rustup toolchain install nightly - ``` + + rustup toolchain install nightly + 3. Install [GNU Scientific Library]. On a Mac with [Homebrew] installed, this can be done with - ```console - brew install gsl - ``` + + brew install gsl + For other operating systems, suggestions are available in the [rust-GSL] crate documentation. If not correctly installed, you may need to pass extra `RUSTFLAGS` options to Cargo in the following steps to locate the @@ -56,21 +56,21 @@ ### Building and running the experiments To compile and install the program, use -```console -cargo install --path=. -``` + + cargo install --path=. + When doing this for the first time, several dependencies will be downloaded. Now you can run the default set of experiments with -``` -pointsource_experiments -o results -``` + + pointsource_experiments -o results --max-iter 4000 1d_fast 2d_fast 1d_tv_fast 2d_tv_fast + The `-o results` option tells `pointsource_algs` to write results in the `results` directory. The option is required. Alternatively, you may build and run the program without installing with -```console -cargo run --release -- -o results -``` + + cargo run --release -- -o results --maxiter 4000 1d_fast 2d_fast 1d_tv_fast 2d_tv_fast + The double-dash separates the options for the Cargo build system and `pointsource_experiments`. @@ -83,10 +83,10 @@ If you are interested in the program internals, the integrated source code documentation may be built and opened with -```console -cargo doc # build dependency docs -misc/cargo-d --open # build and open KaTeX-aware docs for this crate -``` + + cargo doc # build dependency docs + misc/cargo-d --open # build and open KaTeX-aware docs for this crate + The `cargo-d` script ensures that KaTeX mathematics is rendered in the generated documentation through an ugly workaround. Unfortunately, `rustdoc`, akin to Rust largely itself, is stuck in 80's 7-bit gringo ASCII
--- a/src/experiments.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/experiments.rs Sun Jul 19 07:34:39 2026 +0200 @@ -156,18 +156,6 @@ //#[replace_float_literals(F::cast_from(literal))] impl DefaultExperiment { - // fn default_list() -> Vec<Self> { - // use DefaultExperiment::*; - // [ - // Experiment1D, - // Experiment1DFast, - // Experiment2D, - // Experiment2DFast, - // Experiment1D_L1, - // ] - // .into() - // } - /// Convert the experiment shorthand into a runnable experiment configuration. fn get_experiment( &self, @@ -193,58 +181,41 @@ make_float_constant!(Hat1 = 0.16); make_float_constant!(HatBias = 0.05); - // We use a different step length for PDPS in 2D experiments - // let pdps_2d = (DefaultAlgorithm::PDPS, - // AlgorithmOverrides { - // tau0 : Some(3.0), - // sigma0 : Some(0.99 / 3.0), - // .. Default::default() - // } - // ); - // let radon_pdps_2d = (DefaultAlgorithm::RadonPDPS, - // AlgorithmOverrides { - // tau0 : Some(3.0), - // sigma0 : Some(0.99 / 3.0), - // .. Default::default() - // } - // ); - let sliding_fb_cut_gaussian = (DefaultAlgorithm::SlidingFB, AlgorithmOverrides { - theta0: Some(0.3), - ..Default::default() - }); - // let higher_cpos = |alg| (alg, - // AlgorithmOverrides { - // transport_tolerance_pos : Some(1000.0), - // .. Default::default() - // } - // ); - let higher_cpos_merging = |alg| { + let basic = |alg| { (alg, AlgorithmOverrides { - transport_tolerance_pos: Some(1000.0), - merge: Some(true), - fitness_merging: Some(true), + gradv_lipest_mult: Some(0.04), + transport_tolerance: Some(100.0), + ..Default::default() + }) + }; + let radon = |alg| { + (alg, AlgorithmOverrides { + gradv_lipest_mult: Some(0.04), + transport_tolerance: Some(1000.0), ..Default::default() }) }; - let higher_cpos_merging_steptune = |alg| { + let radon_2dtv = |alg| { (alg, AlgorithmOverrides { - transport_tolerance_pos: Some(1000.0), - theta0: Some(0.3), - merge: Some(true), - fitness_merging: Some(true), + gradv_lipest_mult: Some(0.04), + transport_tolerance: Some(10000.0), + sigma0: Some(0.15), ..Default::default() }) }; - let much_higher_cpos_merging_steptune = |alg| { - (alg, AlgorithmOverrides { - transport_tolerance_pos: Some(10000.0), - sigma0: Some(0.15), - theta0: Some(0.3), - merge: Some(true), - fitness_merging: Some(true), - ..Default::default() - }) - }; + macro_rules! overrides { + ($a:ident, $b:ident) => { + HashMap::from([ + $a(DefaultAlgorithm::SlidingFB), + $a(DefaultAlgorithm::SlidingPDPS), + $a(DefaultAlgorithm::ForwardPDPS), + $b(DefaultAlgorithm::RadonFB), + $b(DefaultAlgorithm::RadonSlidingFB), + $b(DefaultAlgorithm::RadonSlidingPDPS), + ]) + }; + } + // We add a hash of the experiment name to the configured // noise seed to not use the same noise for different experiments. let mut h = DefaultHasher::new(); @@ -273,11 +244,7 @@ kernel_plot_width, noise_seed, default_merge_radius, - algorithm_overrides: HashMap::from([ - sliding_fb_cut_gaussian, - higher_cpos_merging(DefaultAlgorithm::RadonFB), - higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), - ]), + algorithm_overrides: overrides!(basic, radon), }, }) } @@ -298,10 +265,7 @@ kernel_plot_width, noise_seed, default_merge_radius, - algorithm_overrides: HashMap::from([ - higher_cpos_merging(DefaultAlgorithm::RadonFB), - higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), - ]), + algorithm_overrides: overrides!(basic, radon), }, }) } @@ -323,11 +287,7 @@ kernel_plot_width, noise_seed, default_merge_radius, - algorithm_overrides: HashMap::from([ - sliding_fb_cut_gaussian, - higher_cpos_merging(DefaultAlgorithm::RadonFB), - higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), - ]), + algorithm_overrides: overrides!(basic, radon), }, }) } @@ -348,10 +308,7 @@ kernel_plot_width, noise_seed, default_merge_radius, - algorithm_overrides: HashMap::from([ - higher_cpos_merging(DefaultAlgorithm::RadonFB), - higher_cpos_merging(DefaultAlgorithm::RadonSlidingFB), - ]), + algorithm_overrides: overrides!(basic, radon), }, }) } @@ -482,10 +439,7 @@ kernel_plot_width, noise_seed, default_merge_radius, - algorithm_overrides: HashMap::from([ - higher_cpos_merging_steptune(DefaultAlgorithm::RadonForwardPDPS), - higher_cpos_merging_steptune(DefaultAlgorithm::RadonSlidingPDPS), - ]), + algorithm_overrides: overrides!(basic, radon), }, }, }) @@ -522,14 +476,7 @@ kernel_plot_width, noise_seed, default_merge_radius, - algorithm_overrides: HashMap::from([ - much_higher_cpos_merging_steptune( - DefaultAlgorithm::RadonForwardPDPS, - ), - much_higher_cpos_merging_steptune( - DefaultAlgorithm::RadonSlidingPDPS, - ), - ]), + algorithm_overrides: overrides!(basic, radon_2dtv), }, }, })
--- a/src/forward_model.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/forward_model.rs Sun Jul 19 07:34:39 2026 +0200 @@ -43,40 +43,20 @@ BetterThanZero, } -/// Curvature error control: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a). -/// -/// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8). -/// Thus, subject to `guess` being correct, returns factor $(ℓ_F, Θ²)$ such that -/// $B_{F'(μ)} dγ ≤ ℓ_F c_2$ and $⟨F'(μ+Δ)-F'(μ)|Δ⟩ ≤ Θ²|γ|(c_2)$, where $Δ=(π_♯^1-π_♯^0)γ$. -/// The latter is the firm transport Lipshitz property of (3.8b) and Lemma 5.11. -/// -/// This trait is supposed to be implemented by the data term $F$, in the basic case a -/// [`Mapping`] from [`RNDM<N, F>`] to a [`Float`] `F`. -/// The generic implementation for operators that satisfy [`BasicCurvatureBoundEstimates`] -/// uses Remark 5.15 and Example 5.16 for (4.2d) and (5.2a), and (5.2b); -/// and Lemma 3.8 for (3.8). +/// Curvature error control. pub trait BoundedCurvature<F: Float = f64> { - /// Returns $(ℓ_F, Θ²)$ or individual errors for each. + /// Returns an estimate of $(ℓ_F, ℓ_∇v, Θ²)$ or individual errors for each. fn curvature_bound_components( &self, guess: BoundedCurvatureGuess, - ) -> (DynResult<F>, DynResult<F>); + ) -> (DynResult<F>, DynResult<F>, DynResult<F>); } -/// Curvature error control: helper bounds for (4.2d), (5.2a), (5.2b), (5.15a), and (5.16a) -/// for quadratic dataterms $F(μ) = \frac{1}{2}\|Aμ-b\|^2$. -/// -/// This trait is to be implemented by the [`Linear`] operator $A$, in the basic from -/// [`RNDM<N, F>`] to a an Euclidean space. -/// It is used by implementations of [`BoundedCurvature`] for $F$. -/// -/// Based on Lemma 5.11 and Example 5.12, the helper bound for (5.15a) and (5.16a) is (3.8). -/// This trait provides the factor $θ²$ of (3.8) as determined by Lemma 3.8. -/// To aid in calculating (4.2d), (5.2a), (5.2b), motivated by Example 5.16, it also provides -/// $ℓ_F^0$ such that $∇v^k$ $ℓ_F^0 \|Aμ-b\|$-Lipschitz. Here $v^k := F'(∪^k)$. +/// Curvature error control: helper Lipschitz-like bounds for the quadratic dataterms +/// $F(μ) = \frac{1}{2}\|Aμ-b\|^2$. pub trait BasicCurvatureBoundEstimates<F: Float = f64> { - /// Returns $(ℓ_F^0, Θ²)$ or individual errors for each. - fn basic_curvature_bound_components(&self) -> (DynResult<F>, DynResult<F>); + /// Returns $(ℓ_F, ℓ_{∇v}^0, Θ²)$ or individual errors for each. + fn basic_curvature_bound_components(&self) -> (DynResult<F>, DynResult<F>, DynResult<F>); } impl<F, A, Z, const N: usize> BoundedCurvature<F> for QuadraticDataTerm<F, RNDM<N, F>, A> @@ -89,13 +69,13 @@ fn curvature_bound_components( &self, guess: BoundedCurvatureGuess, - ) -> (DynResult<F>, DynResult<F>) { + ) -> (DynResult<F>, DynResult<F>, DynResult<F>) { match guess { BoundedCurvatureGuess::BetterThanZero => { let opA = self.operator(); let b = self.data(); - let (ℓ_F0, θ2) = opA.basic_curvature_bound_components(); - (ℓ_F0.map(|l| l * b.norm2()), θ2) + let (ℓ_F, ℓ_gradv_0, θ2) = opA.basic_curvature_bound_components(); + (ℓ_F, ℓ_gradv_0.map(|l| l * b.norm2()), θ2) } } }
--- a/src/forward_model/bias.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/forward_model/bias.rs Sun Jul 19 07:34:39 2026 +0200 @@ -138,13 +138,13 @@ fn curvature_bound_components( &self, guess: BoundedCurvatureGuess, - ) -> (DynResult<F>, DynResult<F>) { + ) -> (DynResult<F>, DynResult<F>, DynResult<F>) { match guess { BetterThanZero => { let opA = &self.operator().0; let b = self.data(); - let (ℓ_F0, θ2) = opA.basic_curvature_bound_components(); - (ℓ_F0.map(|l| l * b.norm2()), θ2) + let (ℓ_F, ℓ_gradv_0, θ2) = opA.basic_curvature_bound_components(); + (ℓ_F, ℓ_gradv_0.map(|l| l * b.norm2()), θ2) } } }
--- a/src/forward_model/sensor_grid.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/forward_model/sensor_grid.rs Sun Jul 19 07:34:39 2026 +0200 @@ -118,14 +118,7 @@ ) -> Self { let base_sensor = Convolution(sensor.clone(), spread.clone()); let bt = BT::new(domain, depth); - let mut sensorgrid = SensorGrid { - domain, - sensor_count, - sensor, - spread, - base_sensor, - bt, - }; + let mut sensorgrid = SensorGrid { domain, sensor_count, sensor, spread, base_sensor, bt }; for (x, id) in sensorgrid.grid().into_iter().zip(0usize..) { let s = sensorgrid.shifted_sensor(x); @@ -173,7 +166,7 @@ }); w.iter() .zip(d.iter()) - .map(|(&wi, &di)| (wi / di).ceil()) + .map(|(&wi, &di)| (wi / di).ceil() + F::ONE) .reduce(F::mul) .unwrap() } @@ -349,14 +342,15 @@ for<'b> <Convolution<S, P> as DifferentiableMapping<Loc<N, F>>>::Differential<'b>: Lipschitz<L2, FloatType = F>, { - fn basic_curvature_bound_components(&self) -> (DynResult<F>, DynResult<F>) { + fn basic_curvature_bound_components(&self) -> (DynResult<F>, DynResult<F>, DynResult<F>) { let n_ψ = self.max_overlapping(); let ψ_diff_lip = self.base_sensor.diff_ref().lipschitz_factor(L2); let ψ_lip = self.base_sensor.lipschitz_factor(L2); - let ℓ_F0 = ψ_diff_lip.map(|l| (2.0 * n_ψ).sqrt() * l); + let ℓ_gradv_0 = ψ_diff_lip.map(|l| n_ψ * l); let Θ2 = ψ_lip.map(|l| 4.0 * n_ψ * l.powi(2)); + let ℓ_F = Ok(0.0); // convex problem - (ℓ_F0, Θ2) + (ℓ_F, ℓ_gradv_0, Θ2) } }
--- a/src/forward_pdps.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/forward_pdps.rs Sun Jul 19 07:34:39 2026 +0200 @@ -128,7 +128,7 @@ KOpZ::SimpleAdjoint: GEMV<F, Y, Z>, Y: ClosedEuclidean<F>, for<'b> &'b Y: Instance<Y>, - Z: ClosedEuclidean<F>, + Z: ClosedEuclidean<F> + AXPY, for<'b> &'b Z: Instance<Z>, R: Prox<Z, Codomain = F>, H: Conjugable<Y, F, Codomain = F>, @@ -200,7 +200,8 @@ // Run the algorithm for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { // Calculate initial transport - let Pair(mut τv, τz) = f.differential(Pair(&μ, &z)) * τ; + let Pair(v, mut 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. @@ -233,7 +234,7 @@ ε, ins, ®, - is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))), + Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))), ); } @@ -241,14 +242,14 @@ stats.pruned += prune_with_stats(&mut μ); // Do z variable primal update - let mut z_new = τz; - opKz_adj.gemv(&mut z_new, -σ_p, &y, -σ_p / τ); - z_new = fnR.prox(σ_p, z_new + &z); + 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 - // opKμ.gemv(&mut y, σ_d*(1.0 + ω), &μ, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1}] - opKz.gemv(&mut y, σ_d * (1.0 + ω), &z_new, 1.0); - // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b - opKz.gemv(&mut y, -σ_d * ω, z, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b + 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;
--- a/src/kernels/ball_indicator.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/kernels/ball_indicator.rs Sun Jul 19 07:34:39 2026 +0200 @@ -176,7 +176,7 @@ if N == 1 { 2.0 * r } else if N == 2 { - r * r + 2.0 * r * r } else { (2.0 * r).powi(N as i32) * F::cast_from(factorial(N)) }
--- a/src/kernels/base.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/kernels/base.rs Sun Jul 19 07:34:39 2026 +0200 @@ -437,24 +437,28 @@ dlip: F, ) -> DynResult<F> { // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have - // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j ≠ i} ψ_j(x_j) + // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j < i} ψ_j(y_j) ∏_{j > i} ψ_j(x_j) // by a simple recursive argument. In particular, if ψ_i=g for all i, j, we have - // |ψ(x) - ψ(y)| ≤ ∑_i L_g M_g^{n-1}|x-y|, where L_g is the Lipschitz factor of g, and - // M_g a bound on it. + // |ψ(x) - ψ(y)| ≤ ∑_i L_g|x_i-y_i| M_g^{n-1} ≤ √n L_g M_g^{n-1} |x-y|₂ + // where L_g is the Lipschitz factor of g, and M_g a bound on it. // - // We also have in the general case ∇ψ(x) = ∑_i ∇ψ_i(x_i) ∏_{j ≠ i} ψ_j(x_j), whence - // using the previous formula for each i with f_i=∇ψ_i and f_j=ψ_j for j ≠ i, we get - // ∇ψ(x) - ∇ψ(y) = ∑_i[ ∇ψ_i(x_i)∏_{j ≠ i} ψ_j(x_j) - ∇ψ_i(y_i)∏_{j ≠ i} ψ_j(y_j)] - // = ∑_i[ [∇ψ_i(x_i) - ∇ψ_j(x_j)] ∏_{j ≠ i}ψ_j(x_j) - // + [∑_{k ≠ i} [ψ_k(x_k) - ∇ψ_k(x_k)] ∏_{j ≠ i, k}ψ_j(x_j)]∇ψ_i(x_i)]. - // With $ψ_i=g for all i, j, it follows that - // |∇ψ(x) - ∇ψ(y)| ≤ ∑_i L_{∇g} M_g^{n-1} + ∑_{k ≠ i} L_g M_g^{n-2} M_{∇g} - // = n [L_{∇g} M_g^{n-1} + (n-1) L_g M_g^{n-2} M_{∇g}]. - // = n M_g^{n-2}[L_{∇g} M_g + (n-1) L_g M_{∇g}]. - if N >= 2 { + // We also have in the general case [∇ψ(x)]_i = ψ_i'(x_i) ∏_{j ≠ i} ψ_j(x_j), whence + // from above, the Lipschitz factor of [∇ψ(x)]_i is + // L' = √(L_{g'}^2M_g^{2(n-1)} + (n-1)L_g^2M_{g'}^2M_g^{2(n-2)}) if n > 2, + // in particular + // L' = √(L_{g'}^2M_g^{2(n-1)} + L_g^2M_{g'}^2) if n > 2, and + // Now the Lipschitz factor of ∇ψ is √n L' from + // ‖∇ψ(x)-∇ψ(y)‖ = √(∑_{i=1}^n [∇ψ(x)-∇ψ(y)]_i^2) ≤ √(∑_{i=1}^n (L')^2) ≤ √n L' + if N > 2 { Ok(F::cast_from(N) - * bound.powi((N - 2) as i32) - * (dlip * bound + F::cast_from(N - 1) * lip * dbound)) + * (dlip.powi(2) * bound.powi(2 * (N - 1) as i32) + + F::cast_from(N - 1) + * lip.powi(2) + * dbound.powi(2) + * bound.powi(2 * (N - 2) as i32)) + .sqrt()) + } else if N == 2 { + Ok((F::TWO * ((dlip * bound).powi(2) + (lip * dbound).powi(2))).sqrt()) } else if N == 1 { Ok(dlip) } else {
--- a/src/kernels/gaussian.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/kernels/gaussian.rs Sun Jul 19 07:34:39 2026 +0200 @@ -421,6 +421,7 @@ } */ +#[replace_float_literals(F::cast_from(literal))] impl<F: Float, R, C, S, const N: usize> Convolution<CubeIndicator<R, N>, BasicCutGaussian<C, S, N>> where R: Constant<Type = F>,
--- a/src/kernels/hat_convolution.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/kernels/hat_convolution.rs Sun Jul 19 07:34:39 2026 +0200 @@ -97,42 +97,30 @@ #[inline] fn apply<I: Instance<Loc<N, S::Type>>>(&self, y: I) -> Self::Codomain { let σ = self.radius(); - y.decompose().product_map(|x| self.value_1d_σ1(x / σ) / σ) + y.decompose().product_map(|x| self.value_1d_σ1(x / σ)) / σ.powi(N as i32) } } #[replace_float_literals(S::Type::cast_from(literal))] -impl<S, const N: usize> Lipschitz<L1> for HatConv<S, N> -where - S: Constant, -{ - type FloatType = S::Type; - #[inline] - fn lipschitz_factor(&self, L1: L1) -> DynResult<Self::FloatType> { - // For any ψ_i, we have - // ∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i) - // = [ψ_1(x_1)-ψ_1(y_1)] ∏_{i=2}^N ψ_i(x_i) - // + ψ_1(y_1)[ ∏_{i=2}^N ψ_i(x_i) - ∏_{i=2}^N ψ_i(y_i)] - // = ∑_{j=1}^N [ψ_j(x_j)-ψ_j(y_j)]∏_{i > j} ψ_i(x_i) ∏_{i < j} ψ_i(y_i) - // Thus - // |∏_{i=1}^N ψ_i(x_i) - ∏_{i=1}^N ψ_i(y_i)| - // ≤ ∑_{j=1}^N |ψ_j(x_j)-ψ_j(y_j)| ∏_{j ≠ i} \max_j |ψ_j| - let σ = self.radius(); - let l1d = self.lipschitz_1d_σ1() / (σ * σ); - let m1d = self.value_1d_σ1(0.0) / σ; - Ok(l1d * m1d.powi(N as i32 - 1)) - } -} - impl<S, const N: usize> Lipschitz<L2> for HatConv<S, N> where S: Constant, { type FloatType = S::Type; #[inline] - fn lipschitz_factor(&self, L2: L2) -> DynResult<Self::FloatType> { - self.lipschitz_factor(L1) - .map(|l1| l1 * <S::Type>::cast_from(N).sqrt()) + fn lipschitz_factor(&self, _L2: L2) -> DynResult<Self::FloatType> { + // For arbitrary ψ(x) = ∏_{i=1}^n ψ_i(x_i), we have + // ψ(x) - ψ(y) = ∑_i [ψ_i(x_i)-ψ_i(y_i)] ∏_{j < i} ψ_j(y_j) ∏_{j > i} ψ_j(x_j) + // by a simple recursive argument. In particular, if ψ_i=g for all i, j, we have + // |ψ(x) - ψ(y)| ≤ ∑_i L_g|x_i-y_i| M_g^{n-1} ≤ √n L_g M_g^{n-1} |x-y|₂ + // where L_g is the Lipschitz factor of g, and M_g a bound on it. + let σ = self.radius(); + let l1d = self.lipschitz_1d_σ1(); + let m1d = self.value_1d_σ1(0.0); + Ok( + S::Type::cast_from(N).sqrt() * l1d * m1d.powi((N - 1) as i32) + / σ.powi((N - 1 + 2) as i32), + ) } } @@ -146,9 +134,9 @@ fn differential_impl<I: Instance<Loc<N, S::Type>>>(&self, y0: I) -> Self::Derivative { let y = y0.decompose(); let σ = self.radius(); - let σ2 = σ * σ; - let vs = y.map(|x| self.value_1d_σ1(x / σ) / σ); - product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ) / σ2) + let tmp = y.map(|x| x / σ); + let vs = tmp.map(|x| self.value_1d_σ1(x)); + product_differential(&tmp, &vs, |x| self.diff_1d_σ1(x)) / σ.powi((N + 1) as i32) } } @@ -164,11 +152,14 @@ fn diff_lipschitz_factor(&self, _l2: L2) -> DynResult<F> { let σ = self.radius(); product_differential_lipschitz_factor::<F, N>( - self.value_1d_σ1(0.0) / σ, - self.lipschitz_1d_σ1() / (σ * σ), - self.maxabsdiff_1d_σ1() / (σ * σ), - self.lipschitz_diff_1d_σ1() / (σ * σ), + self.value_1d_σ1(0.0), + self.lipschitz_1d_σ1(), + self.maxabsdiff_1d_σ1(), + self.lipschitz_diff_1d_σ1(), ) + .map(|l| l / σ.powi((N + 2) as i32)) + // We have f(x)=f₀(x/σ)/σ^N, where f₀(z)=∏_{i=1}^N f₁(z_i). + // Thus |f'(x)-f'(y)| = |f₀'(x/σ)-f₀'(y/σ)|/σ^{N+1} ≤ L_{f₀'}/σ^{N+2}. } } @@ -195,15 +186,15 @@ /// Computes the differential of the kernel for $n=1$ with $σ=1$. #[inline] fn diff_1d_σ1(&self, x: F) -> F { - let y = x.abs(); - if y >= 1.0 { + if x >= 1.0 || x <= -1.0 { 0.0 - } else if y > 0.5 { - -8.0 * (y - 1.0).powi(2) - } else - /* 0 ≤ y ≤ 0.5 */ - { - (24.0 * y - 16.0) * y + } else if x > 0.5 { + -8.0 * (x - 1.0).powi(2) + } else if x < -0.5 { + 8.0 * (x + 1.0).powi(2) + } else { + /* 0 ≤ y ≤ 0.5 */ + (24.0 * x.abs() - 16.0) * x } } @@ -221,22 +212,6 @@ 2.0 } - /// Computes the second differential of the kernel for $n=1$ with $σ=1$. - #[inline] - #[allow(dead_code)] - fn diff2_1d_σ1(&self, x: F) -> F { - let y = x.abs(); - if y >= 1.0 { - 0.0 - } else if y > 0.5 { - -16.0 * (y - 1.0) - } else - /* 0 ≤ y ≤ 0.5 */ - { - 48.0 * y - 16.0 - } - } - /// Computes the differential of the kernel for $n=1$ with $σ=1$. #[inline] fn lipschitz_diff_1d_σ1(&self) -> F { @@ -332,7 +307,7 @@ // = ∫_{x-β}^{x+β} u_σ(z) d z // = (1/σ)∫_{x-β}^{x+β} u_1(z/σ) d z // = ∫_{(x-β)/σ}^{(x+β)/σ} u_1(z) d z - // = [χ_{-β/σ, β/σ} * u_1](x/σ) + // = (χ_{[-β/σ, β/σ]} * u_1)(x/σ) // $$ self.value_1d_σ1(x / σ, β / σ) }) @@ -354,14 +329,14 @@ let Convolution(ref ind, ref hatconv) = self; let β = ind.r.value(); let σ = hatconv.radius(); - let σ2 = σ * σ; let vs = y.map(|x| self.value_1d_σ1(x / σ, β / σ)); - product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ) / σ2) + product_differential(&*y, &vs, |x| self.diff_1d_σ1(x / σ, β / σ)) / σ } } -/// Integrate $f$, whose support is $[c, d]$, on $[a, b]$. +/// Integrate $f'$, whose support is $[c, d]$, on $[a, b]$. +/// The value $f$ is given, not the derivative. /// If $b > d$, add $g()$ to the result. #[inline] #[replace_float_literals(F::cast_from(literal))] @@ -403,12 +378,6 @@ let a = x - β; let b = x + β; - #[inline] - fn pow4<F: Float>(x: F) -> F { - let y = x * x; - y * y - } - // Observe the factor 1/6 at the front from the antiderivatives below. // The factor 4 is from normalisation of the original function. (4.0 / 6.0) @@ -419,37 +388,28 @@ -0.5, // (2/3) (y+1)^3 on -1 < y ≤ -1/2 // The antiderivative is (2/12)(y+1)^4 = (1/6)(y+1)^4 - |y| pow4(y + 1.0), + |y| (y + 1.0).powi(4), || { i( a, b, -0.5, - 0.0, + 0.5, // -2 y^3 - 2 y^2 + 1/3 on -1/2 < y ≤ 0 // The antiderivative is -1/2 y^4 - 2/3 y^3 + 1/3 y - |y| y * (-y * y * (y * 3.0 + 4.0) + 2.0), + // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2 + // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y + |y| y * (y * y * (y.abs() * 3.0 - 4.0) + 2.0), || { i( a, b, - 0.0, 0.5, - // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2 - // The antiderivative is 1/2 y^4 - 2/3 y^3 + 1/3 y - |y| y * (y * y * (y * 3.0 - 4.0) + 2.0), - || { - i( - a, - b, - 0.5, - 1.0, - // -(2/3) (y-1)^3 on 1/2 < y ≤ 1 - // The antiderivative is -(2/12)(y-1)^4 = -(1/6)(y-1)^4 - |y| -pow4(y - 1.0), - || 0.0, - ) - }, + 1.0, + // -(2/3) (y-1)^3 on 1/2 < y ≤ 1 + // The antiderivative is -(2/12)(y-1)^4 = -(1/6)(y-1)^4 + |y| -(y - 1.0).powi(4), + || 0.0, ) }, ) @@ -466,41 +426,30 @@ let a = x - β; let b = x + β; - // The factor 4 is from normalisation of the original function. - 4.0 * i( + i( a, b, -1.0, -0.5, // (2/3) (y+1)^3 on -1 < y ≤ -1/2 - |y| (2.0 / 3.0) * (y + 1.0).powi(3), + |y| (8.0 / 3.0) * (y + 1.0).powi(3), || { i( a, b, -0.5, - 0.0, - // -2 y^3 - 2 y^2 + 1/3 on -1/2 < y ≤ 0 - |y| -2.0 * (y + 1.0) * y * y + (1.0 / 3.0), + 0.5, + // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2 + |y| 8.0 * (y.abs() - 1.0) * y * y + (4.0 / 3.0), || { i( a, b, - 0.0, 0.5, - // 2 y^3 - 2 y^2 + 1/3 on 0 < y < 1/2 - |y| 2.0 * (y - 1.0) * y * y + (1.0 / 3.0), - || { - i( - a, - b, - 0.5, - 1.0, - // -(2/3) (y-1)^3 on 1/2 < y ≤ 1 - |y| -(2.0 / 3.0) * (y - 1.0).powi(3), - || 0.0, - ) - }, + 1.0, + // -(2/3) (y-1)^3 on 1/2 < y ≤ 1 + |y| -(8.0 / 3.0) * (y - 1.0).powi(3), + || 0.0, ) }, ) @@ -525,6 +474,7 @@ } */ +#[replace_float_literals(F::cast_from(literal))] impl<F: Float, R, C, const N: usize> Convolution<CubeIndicator<R, N>, HatConv<C, N>> where R: Constant<Type = F>,
--- a/src/lib.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/lib.rs Sun Jul 19 07:34:39 2026 +0200 @@ -14,7 +14,6 @@ use alg_tools::parallelism::{set_max_threads, set_num_threads}; use clap::Parser; use serde::{Deserialize, Serialize}; -use serde_json; use serde_with::skip_serializing_none; use std::num::NonZeroUsize; @@ -47,10 +46,9 @@ pub use measures::*; } -use run::{AlgorithmConfig, DefaultAlgorithm, Named, PlotLevel, RunnableExperiment}; +use run::{DefaultAlgorithm, PlotLevel, RunnableExperiment}; use subproblem::InnerMethod; use types::{ClapFloat, Float}; -use DefaultAlgorithm::*; /// Trait for customising the experiments available from the command line pub trait ExperimentSetup: @@ -84,8 +82,7 @@ /// /// Not all algorithms are available for all the experiments. /// In particular, only PDPS is available for the experiments with L¹ data term. - #[arg(value_enum, value_name = "ALGORITHM", long, short = 'a', - default_values_t = [FB, PDPS, SlidingFB, FW, RadonFB])] + #[arg(value_enum, value_name = "ALGORITHM", long, short = 'a')] algorithm: Vec<DefaultAlgorithm>, /// Saved algorithm configration(s) to use on the experiments @@ -144,48 +141,66 @@ /// /// The first parameter is the number of bootstrap insertion iterations, and the second /// the maximum number of iterations on each of them. - bootstrap_insertions: Option<Vec<usize>>, + pub bootstrap_insertions: Option<Vec<usize>>, #[arg(long, requires = "algorithm")] /// Primal step length parameter override for --algorithm. /// /// Only use if running just a single algorithm, as different algorithms have different /// regularisation parameters. Does not affect the algorithms fw and fwrelax. - tau0: Option<F>, + pub tau0: Option<F>, #[arg(long, requires = "algorithm")] /// Second primal step length parameter override for SlidingPDPS. /// /// Only use if running just a single algorithm, as different algorithms have different /// regularisation parameters. - sigmap0: Option<F>, + pub sigmap0: Option<F>, #[arg(long, requires = "algorithm")] /// Dual step length parameter override for --algorithm. /// /// Only use if running just a single algorithm, as different algorithms have different /// regularisation parameters. Only affects PDPS. - sigma0: Option<F>, + pub sigma0: Option<F>, #[arg(long)] /// Normalised transport step length for sliding methods. - theta0: Option<F>, + pub theta0: Option<F>, + + #[arg(long)] + /// Unnormalised transport step length for sliding methods, multiplied by tau. + pub tautheta: Option<F>, #[arg(long)] - /// A posteriori transport tolerance multiplier (C_pos) - transport_tolerance_pos: Option<F>, + /// A posteriori transport tolerance multiplier + pub transport_tolerance: Option<F>, + + #[arg(long)] + /// Multiplier for rough estimate of ℓ_{∇v}. Should be ≥ 1. + /// If explicit τθ is given, ℓ_{∇v} is one divided by this number times τθ. + /// Otherwise, if τθ relative to an estimate of ℓ_{∇v}, this number multiplies that estimate. + pub gradv_lipest_mult: Option<F>, #[arg(long)] /// Transport adaptation factor. Must be in (0, 1). - transport_adaptation: Option<F>, + pub transport_adaptation: Option<F>, + + #[arg(long)] + /// Whether partially transported spikes are allowed. + pub allow_partial_transport: Option<bool>, + + #[arg(long)] + /// Use an alternative remainder control rule. + pub alt_remainder_control: Option<bool>, #[arg(long)] /// Minimal step length parameter for sliding methods. - tau0_min: Option<F>, + pub tau0_min: Option<F>, #[arg(value_enum, long)] /// PDPS acceleration, when available. - acceleration: Option<pdps::Acceleration>, + pub acceleration: Option<pdps::Acceleration>, // #[arg(long)] // /// Perform postprocess weight optimisation for saved iterations @@ -196,58 +211,57 @@ /// Merging frequency, if merging enabled (every n iterations) /// /// Only affects FB, FISTA, and PDPS. - merge_every: Option<usize>, + pub merge_every: Option<usize>, #[arg(long)] /// Enable merging (default: determined by algorithm) - merge: Option<bool>, + pub merge: Option<bool>, #[arg(long)] /// Merging radius (default: determined by experiment) - merge_radius: Option<F>, + pub merge_radius: Option<F>, #[arg(long)] /// Interpolate when merging (default : determined by algorithm) - merge_interp: Option<bool>, + pub merge_interp: Option<bool>, #[arg(long)] /// Enable final merging (default: determined by algorithm) - final_merging: Option<bool>, + pub final_merging: Option<bool>, #[arg(long)] /// Enable fitness-based merging for relevant FB-type methods. /// This has worse convergence guarantees that merging based on optimality conditions. - fitness_merging: Option<bool>, + pub fitness_merging: Option<bool>, #[arg(long, value_names = &["ε", "θ", "p"])] /// Set the tolerance to ε_k = ε/(1+θk)^p - tolerance: Option<Vec<F>>, + pub tolerance: Option<Vec<F>>, #[arg(long)] /// Method for solving inner optimisation problems - inner_method: Option<InnerMethod>, + pub inner_method: Option<InnerMethod>, #[arg(long)] /// Step length parameter for inner problem - inner_τ0: Option<F>, + pub inner_τ0: Option<F>, #[arg(long, value_names = &["τ0", "σ0"])] /// Dual step length parameter for inner problem - inner_pdps_τσ0: Option<Vec<F>>, + pub inner_pdps_τσ0: Option<Vec<F>>, #[arg(long, value_names = &["τ", "growth"])] /// Inner proximal point method step length and its growth - inner_pp_τ: Option<Vec<F>>, + pub inner_pp_τ: Option<Vec<F>>, #[arg(long)] /// Inner tolerance multiplier - inner_tol: Option<F>, + pub inner_tol: Option<F>, } /// A generic entry point for binaries based on this library pub fn common_main<E: ExperimentSetup>() -> DynResult<()> { - let full_cli = FusedCommandLineArgs::<E>::parse(); - let cli = &full_cli.general; + let cli = FusedCommandLineArgs::<E>::parse(); #[cfg(debug_assertions)] { @@ -267,32 +281,16 @@ ); } - if let Some(n_threads) = cli.num_threads { + if let Some(n_threads) = cli.general.num_threads { let n = NonZeroUsize::new(n_threads).expect("Invalid thread count"); set_num_threads(n); } else { - let m = NonZeroUsize::new(cli.max_threads).expect("Invalid maximum thread count"); + let m = NonZeroUsize::new(cli.general.max_threads).expect("Invalid maximum thread count"); set_max_threads(m); } - for experiment in full_cli.experiment_setup.runnables()? { - let mut algs: Vec<Named<AlgorithmConfig<E::FloatType>>> = cli - .algorithm - .iter() - .map(|alg| { - let cfg = alg - .default_config() - .cli_override(&experiment.algorithm_overrides(*alg)) - .cli_override(&full_cli.algorithm_overrides); - alg.to_named(cfg) - }) - .collect(); - for filename in cli.saved_algorithm.iter() { - let f = std::fs::File::open(filename)?; - let alg = serde_json::from_reader(f)?; - algs.push(alg); - } - experiment.runall(&cli, (!algs.is_empty()).then_some(algs))?; + for experiment in cli.experiment_setup.runnables()? { + experiment.runall(&cli.general, &cli.algorithm_overrides)?; } Ok(())
--- a/src/prox_penalty/radon_squared.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/prox_penalty/radon_squared.rs Sun Jul 19 07:34:39 2026 +0200 @@ -56,7 +56,11 @@ where I: AlgIterator, { - let violation = reg.find_tolerance_violation(τv, τ, ε, true, config); + // If no merging heuristic is used, let's be more conservative about spike insertion, + // and skip it after first round. If merging is done, being more greedy about spike + // insertion also seems to improve performance. + let skip_by_rough_check = !config.merging.enabled; + let violation = reg.find_tolerance_violation(τv, τ, ε, skip_by_rough_check, config); reg.solve_oc_radonsq(μ, τv, τ, ε, violation, config, stats); Ok((None, true))
--- a/src/regularisation.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/regularisation.rs Sun Jul 19 07:34:39 2026 +0200 @@ -234,8 +234,11 @@ where M: MinMaxMapping<Domain, F>; - /// Convert bound on the regulariser to a bond on the Radon norm + /// Convert bound on the regulariser to a bound on the Radon norm fn radon_norm_bound(&self, b: F) -> F; + + /// Returns true if $v$ is within the pointwise range of the subdifferential + fn subdiff_range(&self) -> Bounds<F>; } #[replace_float_literals(F::cast_from(literal))] @@ -386,8 +389,7 @@ // Solve finite-dimensional subproblem. let inner_tolerance = ε * config.inner.tolerance_mult; let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); - stats.inner_iters += - l1squared_nonneg(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it); + stats.inner_iters += l1squared_nonneg(&y, &g_na, τα, &mut x, &config.inner, inner_it); // Update masses of μ based on solution of finite-dimensional subproblem. μ.set_masses_dvector(&x); @@ -416,12 +418,12 @@ μ.both_matching(radon_μ).all(|(α, rα, x)| { let v = -d.apply(x); // TODO: observe ad hoc negation here, after minus_τv // switch to τv. - let (l1, u1) = match α.partial_cmp(&0.0).unwrap_or(Equal) { + let (l1, u1) = match α.total_cmp(&0.0) { Greater => (τα, τα), _ => (F::NEG_INFINITY, τα), // Less should not happen; treated as Equal }; - let (l2, u2) = match rα.partial_cmp(&0.0).unwrap_or(Equal) { + let (l2, u2) = match rα.total_cmp(&0.0) { Greater => (slack, slack), Equal => (-slack, slack), Less => (-slack, -slack), @@ -465,6 +467,10 @@ fn radon_norm_bound(&self, b: F) -> F { b / self.α() } + + fn subdiff_range(&self) -> Bounds<F> { + Bounds(F::NEG_INFINITY, self.α()) + } } #[replace_float_literals(F::cast_from(literal))] @@ -644,7 +650,7 @@ let inner_tolerance = ε * config.inner.tolerance_mult; let inner_it = config.inner.iterator_options.stop_target(inner_tolerance); stats.inner_iters += - l1squared_unconstrained(&y, &g_na, τα, 1.0, &mut x, &config.inner, inner_it); + l1squared_unconstrained(&y, &g_na, τα, &mut x, &config.inner, inner_it); // Update masses of μ based on solution of finite-dimensional subproblem. μ.set_masses_dvector(&x); @@ -732,4 +738,9 @@ fn radon_norm_bound(&self, b: F) -> F { b / self.α() } + + fn subdiff_range(&self) -> Bounds<F> { + let α = self.α(); + Bounds(-α, α) + } }
--- a/src/run.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/run.rs Sun Jul 19 07:34:39 2026 +0200 @@ -23,7 +23,9 @@ }; use crate::regularisation::{NonnegRadonRegTerm, RadonRegTerm, Regularisation, SlidingRegTerm}; use crate::seminorms::*; -use crate::sliding_fb::{pointsource_sliding_fb_reg, SlidingFBConfig, TransportConfig}; +use crate::sliding_fb::{ + pointsource_sliding_fb_reg, SlidingFBConfig, TransportConfig, TransportProxPenalty, +}; use crate::sliding_pdps::{ pointsource_sliding_fb_pair, pointsource_sliding_pdps_pair, SlidingPDPSConfig, }; @@ -129,8 +131,14 @@ }; let override_transport = |g: TransportConfig<F>| TransportConfig { θ0: cli.theta0.unwrap_or(g.θ0), - tolerance_mult_con: cli.transport_tolerance_pos.unwrap_or(g.tolerance_mult_con), + τθ: cli.tautheta.or(g.τθ), + tolerance_mult: cli.transport_tolerance.unwrap_or(g.tolerance_mult), + ℓ_gradv_mult: cli.gradv_lipest_mult.unwrap_or(g.ℓ_gradv_mult), adaptation: cli.transport_adaptation.unwrap_or(g.adaptation), + allow_partial_transport: cli + .allow_partial_transport + .unwrap_or(g.allow_partial_transport), + alt_remainder_control: cli.alt_remainder_control.unwrap_or(g.alt_remainder_control), ..g }; @@ -270,12 +278,20 @@ impl DefaultAlgorithm { /// Returns the algorithm configuration corresponding to the algorithm shorthand - pub fn default_config<F: Float>(&self) -> AlgorithmConfig<F> { + pub fn default_config<F: Float>( + &self, + regularisation_hint: Option<&Regularisation<F>>, + ) -> AlgorithmConfig<F> { use DefaultAlgorithm::*; let radon_insertion = InsertionConfig { - merging: SpikeMergingMethod { interp: false, ..Default::default() }, + merging: SpikeMergingMethod { enabled: true, interp: false, ..Default::default() }, + fitness_merging: true, inner: InnerSettings { - method: InnerMethod::PDPS, // SSN not implemented + method: if let Some(Regularisation::NonnegRadon(_)) = regularisation_hint { + InnerMethod::Exact + } else { + InnerMethod::PDPS // SSN not implemented + }, ..Default::default() }, ..Default::default() @@ -320,11 +336,6 @@ } } - /// Returns the [`Named`] algorithm corresponding to the algorithm shorthand - pub fn get_named<F: Float>(&self) -> Named<AlgorithmConfig<F>> { - self.to_named(self.default_config()) - } - pub fn to_named<F: Float>(self, alg: AlgorithmConfig<F>) -> Named<AlgorithmConfig<F>> { Named { name: self.name(), data: alg } } @@ -485,12 +496,9 @@ fn runall( &self, cli: &CommandLineArgs, - algs: Option<Vec<Named<AlgorithmConfig<F>>>>, + cli_algorithm_overrides: &AlgorithmOverrides<F>, ) -> DynError; - /// Return algorithm default config - fn algorithm_overrides(&self, alg: DefaultAlgorithm) -> AlgorithmOverrides<F>; - /// Experiment name fn name(&self) -> &str; } @@ -511,6 +519,7 @@ TimingIteratorFactory<BasicAlgIteratorFactory<IterInfo<F>>>, >; +/// Trait for things used by many a [`RunnableExperiment`], but that are not dyn-compatible. pub trait RunnableExperimentExtras<F: ClapFloat>: RunnableExperiment<F> + Serialize + Sized { @@ -545,6 +554,42 @@ Ok(prefix) } + /// Collect algorithms to run based on command line args. + fn collect_algs( + &self, + cli: &CommandLineArgs, + algorithm_overrides: &AlgorithmOverrides<F>, + algorithm_overrides_fn: impl Fn(DefaultAlgorithm) -> Option<AlgorithmOverrides<F>>, + regularisation_hint: Option<&Regularisation<F>>, + dflt: &[DefaultAlgorithm], + ) -> DynResult<Vec<Named<AlgorithmConfig<F>>>> + where + F: for<'a> Deserialize<'a>, + { + let proc_alg = |alg: &DefaultAlgorithm| { + let default_cfg = alg.default_config(regularisation_hint); + let cfg = if let Some(over) = &algorithm_overrides_fn(*alg) { + default_cfg.cli_override(over) + } else { + default_cfg + } + .cli_override(algorithm_overrides); + alg.to_named(cfg) + }; + + let mut algs: Vec<Named<AlgorithmConfig<F>>> = cli.algorithm.iter().map(proc_alg).collect(); + for filename in cli.saved_algorithm.iter() { + let f = std::fs::File::open(filename)?; + let alg = serde_json::from_reader(f)?; + algs.push(alg); + } + if algs.is_empty() { + Ok(dflt.iter().map(proc_alg).collect()) + } else { + Ok(algs) + } + } + /// Helper function to run all algorithms on an experiment. fn do_runall<P, Z, Plot, const N: usize>( &self, @@ -707,22 +752,10 @@ self.name.as_ref() } - fn algorithm_overrides(&self, alg: DefaultAlgorithm) -> AlgorithmOverrides<F> { - AlgorithmOverrides { - merge_radius: Some(self.data.default_merge_radius), - ..self - .data - .algorithm_overrides - .get(&alg) - .cloned() - .unwrap_or(Default::default()) - } - } - fn runall( &self, cli: &CommandLineArgs, - algs: Option<Vec<Named<AlgorithmConfig<F>>>>, + cli_algorithm_overrides: &AlgorithmOverrides<F>, ) -> DynError { // Get experiment configuration let &ExperimentV2 { @@ -741,11 +774,34 @@ } = &self.data; // Set up algorithms - let algorithms = match (algs, dataterm) { - (Some(algs), _) => algs, - (None, DataTermType::L222) => vec![DefaultAlgorithm::FB.get_named()], - (None, DataTermType::L1) => vec![DefaultAlgorithm::PDPS.get_named()], - }; + let algorithms = self.collect_algs( + cli, + cli_algorithm_overrides, + |alg| { + Some(AlgorithmOverrides { + merge_radius: Some(self.data.default_merge_radius), + ..self + .data + .algorithm_overrides + .get(&alg) + .cloned() + .unwrap_or_else(Default::default) + }) + }, + Some(®ularisation), + match dataterm { + DataTermType::L222 => &[ + DefaultAlgorithm::FB, + DefaultAlgorithm::RadonFB, + DefaultAlgorithm::SlidingFB, + DefaultAlgorithm::RadonSlidingFB, + DefaultAlgorithm::FW, + DefaultAlgorithm::FWRelax, + DefaultAlgorithm::PDPS, + ], + DataTermType::L1 => &[DefaultAlgorithm::PDPS], + }, + )?; // Set up operators let depth = DynamicDepth(8); @@ -944,10 +1000,11 @@ F: Float + ToNalgebraRealField, I: AlgIteratorFactory<IterInfo<F>>, Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>, - Dat::DerivativeDomain: DifferentiableRealMapping<N, F> + ClosedMul<F>, + Dat::DerivativeDomain: + DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>, RNDM<N, F>: SpikeMerging<F>, Reg: SlidingRegTerm<Loc<N, F>, F>, - P: ProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>, + P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>, Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>, { let pt = P::prox_type(); @@ -1067,23 +1124,10 @@ self.name.as_ref() } - fn algorithm_overrides(&self, alg: DefaultAlgorithm) -> AlgorithmOverrides<F> { - AlgorithmOverrides { - merge_radius: Some(self.data.base.default_merge_radius), - ..self - .data - .base - .algorithm_overrides - .get(&alg) - .cloned() - .unwrap_or(Default::default()) - } - } - fn runall( &self, cli: &CommandLineArgs, - algs: Option<Vec<Named<AlgorithmConfig<F>>>>, + cli_algorithm_overrides: &AlgorithmOverrides<F>, ) -> DynError { // Get experiment configuration let &ExperimentBiased { @@ -1107,10 +1151,29 @@ } = &self.data; // Set up algorithms - let algorithms = match (algs, dataterm) { - (Some(algs), _) => algs, - _ => vec![DefaultAlgorithm::SlidingPDPS.get_named()], - }; + let algorithms = self.collect_algs( + cli, + cli_algorithm_overrides, + |alg| { + Some(AlgorithmOverrides { + merge_radius: Some(self.data.base.default_merge_radius), + ..self + .data + .base + .algorithm_overrides + .get(&alg) + .cloned() + .unwrap_or_else(Default::default) + }) + }, + Some(®ularisation), + &[ + DefaultAlgorithm::SlidingPDPS, + DefaultAlgorithm::ForwardPDPS, + DefaultAlgorithm::RadonSlidingPDPS, + DefaultAlgorithm::RadonForwardPDPS, + ], + )?; // Set up operators let depth = DynamicDepth(8); @@ -1259,19 +1322,19 @@ I: AlgIteratorFactory<IterInfo<F>>, Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>> + BoundedCurvature<F>, - S: DifferentiableRealMapping<N, F> + ClosedMul<F>, + S: DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>, for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, //Pair<S, Z>: ClosedMul<F>, RNDM<N, F>: SpikeMerging<F>, Reg: SlidingRegTerm<Loc<N, F>, F>, - P: ProxPenalty<Loc<N, F>, S, Reg, F>, + P: TransportProxPenalty<Loc<N, F>, S, Reg, F>, KOpZ: BoundedLinear<Z, L2, L2, F, Codomain = Y> + GEMV<F, Z> + SimplyAdjointable<Z, Y, AdjointCodomain = Z>, KOpZ::SimpleAdjoint: GEMV<F, Y>, Y: ClosedEuclidean<F> + Clone, for<'b> &'b Y: Instance<Y>, - Z: ClosedEuclidean<F> + Clone + ClosedMul<F>, + Z: ClosedEuclidean<F> + Clone + ClosedMul<F> + AXPY, for<'b> &'b Z: Instance<Z>, R: Prox<Z, Codomain = F>, H: Conjugable<Y, F, Codomain = F>, @@ -1342,10 +1405,10 @@ I: AlgIteratorFactory<IterInfo<F>>, Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>> + BoundedCurvature<F>, - S: DifferentiableRealMapping<N, F> + ClosedMul<F>, + S: DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>, RNDM<N, F>: SpikeMerging<F>, Reg: SlidingRegTerm<Loc<N, F>, F>, - P: ProxPenalty<Loc<N, F>, S, Reg, F>, + P: TransportProxPenalty<Loc<N, F>, S, Reg, F>, for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, Z: ClosedEuclidean<F> + AXPY + Clone, for<'b> &'b Z: Instance<Z>,
--- a/src/sliding_fb.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/sliding_fb.rs Sun Jul 19 07:34:39 2026 +0200 @@ -16,13 +16,16 @@ use crate::measures::merging::SpikeMerging; use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, RNDM}; use crate::plot::Plotter; -use crate::prox_penalty::{ProxPenalty, StepLengthBound}; +use crate::prox_penalty::{ProxPenalty, RadonSquared, StepLengthBound}; use crate::regularisation::SlidingRegTerm; +use crate::seminorms::DiscreteMeasureOp; use crate::types::*; +use alg_tools::bounds::{Bounds, MinMaxMapping}; use alg_tools::error::DynResult; use alg_tools::euclidean::Euclidean; +use alg_tools::instance::Space; use alg_tools::iterate::AlgIteratorFactory; -use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping}; +use alg_tools::mapping::{DifferentiableMapping, DifferentiableRealMapping, Mapping, RealMapping}; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::norms::Norm; use anyhow::ensure; @@ -34,14 +37,24 @@ pub struct TransportConfig<F: Float> { /// Transport step length $θ$ normalised to $(0, 1)$. pub θ0: F, + /// Unnormalised transport step length $θ$. Overrides θ0. + pub τθ: Option<F>, /// Factor in $(0, 1)$ for decreasing transport to adapt to tolerance. pub adaptation: F, - /// A posteriori transport tolerance multiplier (C_pos) - pub tolerance_mult_con: F, + /// A posteriori transport tolerance multiplier + pub tolerance_mult: F, + /// Multiplier for rough estimate of ℓ_{∇v}. Should be ≥ 1. + /// If explicit τθ is given, ℓ_{∇v} is one divided by this number times τθ. + /// Otherwise, if τθ relative to an estimate of ℓ_{∇v}, this number multiplies that estimate. + pub ℓ_gradv_mult: F, /// maximum number of adaptation iterations, until cancelling transport. pub max_attempts: usize, /// Maximum number of failed transportations for a single source point pub max_fail: usize, + /// Allow points to be transported partially. + pub allow_partial_transport: bool, + /// Use an alternative remainder control rule. + pub alt_remainder_control: bool, } #[replace_float_literals(F::cast_from(literal))] @@ -50,7 +63,7 @@ pub fn check(&self) -> DynResult<()> { ensure!(self.θ0 > 0.0); ensure!(0.0 < self.adaptation && self.adaptation < 1.0); - ensure!(self.tolerance_mult_con > 0.0); + ensure!(self.tolerance_mult > 0.0); Ok(()) } } @@ -59,9 +72,13 @@ impl<F: Float> Default for TransportConfig<F> { fn default() -> Self { TransportConfig { - θ0: 0.9, + θ0: 0.99, + τθ: None, adaptation: 0.9, - tolerance_mult_con: 100.0, + allow_partial_transport: true, + alt_remainder_control: false, + tolerance_mult: 1e1, + ℓ_gradv_mult: 3.0, max_attempts: 2, max_fail: usize::MAX, } @@ -98,63 +115,514 @@ } /// Internal type of adaptive transport step length calculation -pub(crate) enum TransportStepLength<F: Float, G: Fn(F, F) -> F> { +#[derive(Clone, Debug, Serialize, Deserialize)] +pub enum TransportStepLength<F: Float> { /// Fixed, known step length #[allow(dead_code)] - Fixed(F), + Fixed { τθ: F, ℓ_gradv: F }, + /// Simple step lengths that do not depend on maximum transport + Simple { ℓ_gradv: F, τθ0: F, ℓ_base: F }, /// Adaptive step length, only wrt. maximum transport. - /// Content of `l` depends on use case, while `g` calculates the step length from `l`. - AdaptiveMax { l: F, max_transport: F, g: G }, + AdaptiveMax { + ℓ_gradv: F, + adaptive_max_transport: F, + τθ0: F, + ℓ_base: F, + ℓ_base_max_transport: F, + }, /// Adaptive step length. - /// Content of `l` depends on use case, while `g` calculates the step length from `l`. - FullyAdaptive { l: F, max_transport: F, g: G }, + FullyAdaptive { + adaptive_ℓ_gradv: F, + adaptive_max_transport: F, + τθ0: F, + ℓ_base: F, + ℓ_base_max_transport: F, + }, +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F: Float> TransportStepLength<F> { + fn get_ℓ_gradv(&self) -> F { + use TransportStepLength::*; + match *self { + Fixed { ℓ_gradv, .. } => ℓ_gradv, + Simple { ℓ_gradv, .. } => ℓ_gradv, + AdaptiveMax { ℓ_gradv, .. } => ℓ_gradv, + FullyAdaptive { adaptive_ℓ_gradv, .. } => adaptive_ℓ_gradv, + } + } + + fn new( + maybe_ℓ_gradv_est: DynResult<F>, + tconfig: &TransportConfig<F>, + ℓ_base: F, + ℓ_base_max_transport: F, + ) -> Self { + if let Some(τθ) = tconfig.τθ { + TransportStepLength::Fixed { τθ, ℓ_gradv: 1.0 / (τθ * tconfig.ℓ_gradv_mult) } + } else { + match maybe_ℓ_gradv_est { + Ok(ℓ_gradv_est) => { + if ℓ_base_max_transport == 0.0 { + TransportStepLength::Simple { + ℓ_gradv: tconfig.ℓ_gradv_mult * ℓ_gradv_est, + τθ0: tconfig.θ0, + ℓ_base, + } + } else { + TransportStepLength::AdaptiveMax { + ℓ_gradv: tconfig.ℓ_gradv_mult * ℓ_gradv_est, + adaptive_max_transport: 0.0, + τθ0: tconfig.θ0, + ℓ_base, + ℓ_base_max_transport, + } + } + } + Err(_) => TransportStepLength::FullyAdaptive { + adaptive_ℓ_gradv: 10.0 * F::EPSILON, // Start with something very small to estimate differentials + adaptive_max_transport: 0.0, + τθ0: tconfig.θ0, + ℓ_base, + ℓ_base_max_transport, + }, + } + } + } } #[derive(Clone, Debug, Serialize)] -pub struct SingleTransport<const N: usize, F: Float> { +pub struct SingleTransport<Domain, F: Float> { /// Source point - x: Loc<N, F>, + x: Domain, /// Target point - y: Loc<N, F>, + y: Domain, /// Original mass α_μ_orig: F, /// Transported mass α_γ: F, /// Helper for pruning - prune: bool, + retain: bool, /// Fail count fail_count: usize, + /// Contribution to remainder (temporary variable) + excess: F, } #[derive(Clone, Debug, Serialize)] -pub struct Transport<const N: usize, F: Float> { - vec: Vec<SingleTransport<N, F>>, +pub struct Transport<Domain, F: Float> { + vec: Vec<SingleTransport<Domain, F>>, } -/// Whether partiall transported points are allowed. +/// Whether partially transported points are allowed. /// /// Partial transport can cause spike count explosion, so full or zero /// transport is generally preferred. If this is set to `true`, different /// transport adaptation heuristics will be used. -const ALLOW_PARTIAL_TRANSPORT: bool = true; const MINIMAL_PARTIAL_TRANSPORT: bool = true; +const NEW_APPROACH: bool = true; -impl<const N: usize, F: Float> Transport<N, F> { +pub trait TransportProxPenalty<Domain, PreadjointCodomain, Reg, F = f64>: + ProxPenalty<Domain, PreadjointCodomain, Reg, F> +where + F: Float + ToNalgebraRealField, + Reg: SlidingRegTerm<Domain, F>, + Domain: Space + Clone, +{ + type TransportStepLength; + + /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)` + /// with step lengh τ and transport step length `θ_or_adaptive`. + fn initial_transport( + &self, + γ: &mut Transport<Domain, F>, + μ: &DiscreteMeasure<Domain, F>, + ε: F, + τ: F, + τθ_or_adaptive: &mut Self::TransportStepLength, + v: &PreadjointCodomain, + tconfig: &TransportConfig<F>, + ); + + /// A posteriori transport adaptation. + fn aposteriori_transport( + &self, + γ: &mut Transport<Domain, F>, + μ: &DiscreteMeasure<Domain, F>, + μ̆: &DiscreteMeasure<Domain, F>, + τv̆: &mut PreadjointCodomain, + v: &mut PreadjointCodomain, + extra: Option<F>, + ε: F, + τ: F, + τθ_or_adaptive: &Self::TransportStepLength, + reg: &Reg, + tconfig: &TransportConfig<F>, + rconfig: &RefinementSettings<F>, + attempts: &mut usize, + ) -> bool; + + fn get_transport_steplength( + &self, + lips: (DynResult<F>, DynResult<F>, DynResult<F>), + tconfig: &TransportConfig<F>, + ℓ_base: F, + ℓ_base_max_transport: F, + ) -> Self::TransportStepLength; +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F, M, Reg, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for RadonSquared +where + RadonSquared: ProxPenalty<Loc<N, F>, M, Reg, F>, + F: Float + ToNalgebraRealField, + M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>, + Reg: SlidingRegTerm<Loc<N, F>, F>, + RNDM<N, F>: SpikeMerging<F>, +{ + type TransportStepLength = TransportStepLength<F>; + + fn get_transport_steplength( + &self, + (ℓ_F, maybe_ℓ_gradv, _maybe_transport_lip): (DynResult<F>, DynResult<F>, DynResult<F>), + tconfig: &TransportConfig<F>, + ℓ: F, + ℓ_base_max_transport: F, + ) -> Self::TransportStepLength { + TransportStepLength::new( + maybe_ℓ_gradv, + tconfig, + ℓ + ℓ_F.unwrap_or(0.0), + ℓ_base_max_transport, + ) + } + + fn initial_transport( + &self, + γ: &mut Transport<Loc<N, F>, F>, + μ: &RNDM<N, F>, + _ε: F, + _τ: F, + τθ_or_adaptive: &mut TransportStepLength<F>, + v: &M, + tconfig: &TransportConfig<F>, + ) { + γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig); + } + + fn aposteriori_transport( + &self, + γ: &mut Transport<Loc<N, F>, F>, + μ: &RNDM<N, F>, + μ̆: &RNDM<N, F>, + τv̆: &mut M, + v: &mut M, + _extra: Option<F>, + ε: F, + τ: F, + τθ_or_adaptive: &TransportStepLength<F>, + reg: &Reg, + tconfig: &TransportConfig<F>, + _rconfig: &RefinementSettings<F>, + attempts: &mut usize, + ) -> bool { + *attempts += 1; + + if *attempts > tconfig.max_attempts { + // Previous round has set transport to zero + return true; + } + + let nΔ = μ.dist_matching(&μ̆); + let all_ok = γ.do_new_aposteriori_transport( + μ, + τv̆, + v, + ε, + τ, + τθ_or_adaptive, + reg, + tconfig, + |_, mass_zero, negated| { + use std::cmp::Ordering::*; + match (mass_zero, negated) { + (Less, _) => -nΔ, + (Greater, _) => nΔ, + (Equal, false) => nΔ, // pessimistic estimated for ω + (Equal, true) => -nΔ, // pessimistic estimated for -ω + } + }, + ); + + if !all_ok { + if *attempts >= tconfig.max_attempts { + for ρ in γ.iter_mut() { + ρ.α_γ = 0.0; + } + } + } else { + for ρ in γ.iter_mut() { + if ρ.α_γ == 0.0 { + ρ.fail_count += 1; + } else { + ρ.fail_count = 0; + } + } + } + all_ok + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl<F, M, Reg, 𝒟, O, const N: usize> TransportProxPenalty<Loc<N, F>, M, Reg, F> for 𝒟 +where + F: Float + ToNalgebraRealField, + 𝒟: DiscreteMeasureOp<Loc<N, F>, F>, + 𝒟::Codomain: RealMapping<N, F>, + M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>, + for<'a> &'a M: std::ops::Add<𝒟::PreCodomain, Output = O>, + O: MinMaxMapping<Loc<N, F>, F>, + Reg: SlidingRegTerm<Loc<N, F>, F>, + Self: ProxPenalty<Loc<N, F>, M, Reg, F>, + RNDM<N, F>: SpikeMerging<F>, +{ + type TransportStepLength = TransportStepLength<F>; + + fn get_transport_steplength( + &self, + (ℓ_F, maybe_ℓ_gradv, _maybe_transport_lip): (DynResult<F>, DynResult<F>, DynResult<F>), + tconfig: &TransportConfig<F>, + ℓ_base: F, + ℓ_base_max_transport: F, + ) -> Self::TransportStepLength { + TransportStepLength::new( + maybe_ℓ_gradv, + tconfig, + ℓ_base + ℓ_F.unwrap_or(0.0), + ℓ_base_max_transport, + ) + } + + fn initial_transport( + &self, + γ: &mut Transport<Loc<N, F>, F>, + μ: &RNDM<N, F>, + _ε: F, + _τ: F, + τθ_or_adaptive: &mut TransportStepLength<F>, + v: &M, + tconfig: &TransportConfig<F>, + ) { + γ.do_init_transport(v, μ, τθ_or_adaptive, tconfig); + } + + fn aposteriori_transport( + &self, + γ: &mut Transport<Loc<N, F>, F>, + μ: &RNDM<N, F>, + μ̆: &RNDM<N, F>, + τv̆: &mut M, + v: &mut M, + extra: Option<F>, + ε: F, + τ: F, + τθ_or_adaptive: &TransportStepLength<F>, + reg: &Reg, + tconfig: &TransportConfig<F>, + _rconfig: &RefinementSettings<F>, + attempts: &mut usize, + ) -> bool { + *attempts += 1; + + if *attempts > tconfig.max_attempts { + // Previous round has set transport to zero + return true; + } + + let all_ok = if NEW_APPROACH { + let ω = self.apply(μ.sub_matching(&μ̆)); + γ.do_new_aposteriori_transport( + μ, + τv̆, + v, + ε, + τ, + τθ_or_adaptive, + reg, + tconfig, + |x, _, _| ω.apply(x), + ) + } else { + let mut all_ok0 = true; + + // 1. If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not, + // then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1 + // at that point to zero, and retry. + for (δ, ρ) in izip!(μ.iter_spikes(), γ.iter_mut()) { + if δ.α == 0.0 && ρ.α_γ != 0.0 { + all_ok0 = false; + ρ.α_γ = 0.0; + } + // TODO: sign + } + + // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z). + // through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ̆^k + // which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient. + + let nγ = γ.norm(Radon); + let nΔ = μ.dist_matching(&μ̆) + extra.unwrap_or(0.0); + let t = ε * tconfig.tolerance_mult; + if nγ * nΔ > t && *attempts >= tconfig.max_attempts { + all_ok0 = false; + } else if nγ * nΔ > t { + // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1, + // this will guarantee that eventually ‖γ‖ decreases sufficiently that we + // will not enter here. + //*γ *= tconfig.adaptation * t / (nγ * nΔ); + + // We want a consistent behaviour that has the potential to set many weights to zero. + // Therefore, we find the smallest uniform reduction `chg_one`, subtracted + // from all weights, that achieves total `adapt` adaptation. + let adapt_to = tconfig.adaptation * t / nΔ; + let reduction_target = nγ - adapt_to; + assert!(reduction_target > 0.0); + if tconfig.allow_partial_transport { + if MINIMAL_PARTIAL_TRANSPORT { + // This reduces weights of transport, starting from … until `adapt` is + // exhausted. It will, therefore, only ever cause one extrap point insertion + // at the sources, unlike “full” partial transport. + //let refs = γ.vec.iter_mut().collect::<Vec<_>>(); + //refs.sort_by(|ρ1, ρ2| ρ1.α_γ.abs().partial_cmp(&ρ2.α_γ.abs()).unwrap()); + // let mut it = refs.into_iter(); + // + // Maybe sort by differential norm + // let mut refs = γ + // .vec + // .iter_mut() + // .map(|ρ| { + // let val = v.differential(&ρ.x).norm2_squared(); + // (ρ, val) + // }) + // .collect::<Vec<_>>(); + // refs.sort_by(|(_, v1), (_, v2)| v2.partial_cmp(&v1).unwrap()); + // let mut it = refs.into_iter().map(|(ρ, _)| ρ); + let mut it = γ.vec.iter_mut().rev(); + let _unused = it.try_fold(reduction_target, |left, ρ| { + let w = ρ.α_γ.abs(); + if left <= w { + ρ.α_γ = ρ.α_γ.signum() * (w - left); + ControlFlow::Break(()) + } else { + ρ.α_γ = 0.0; + ControlFlow::Continue(left - w) + } + }); + } else { + // This version equally reduces all weights. It causes partial transport, which + // has the problem that that we need to then adapt weights in both start and + // end points, in insert_and_reweigh, somtimes causing the number of spikes μ + // to explode. + let mut abs_weights = γ + .vec + .iter() + .map(|ρ| ρ.α_γ.abs()) + .filter(|t| *t > F::EPSILON) + .collect::<Vec<F>>(); + abs_weights.sort_by(|a, b| a.total_cmp(b)); + let n = abs_weights.len(); + // Cannot have partial transport; can cause spike count explosion + let chg = abs_weights.into_iter().zip((1..=n).rev()).try_fold( + 0.0, + |smaller_total, (w, m)| { + let mf = F::cast_from(m); + let reduction = w * mf + smaller_total; + if reduction >= reduction_target { + ControlFlow::Break((reduction_target - smaller_total) / mf) + } else { + ControlFlow::Continue(smaller_total + w) + } + }, + ); + match chg { + ControlFlow::Continue(_) => γ.vec.iter_mut().for_each(|δ| δ.α_γ = 0.0), + ControlFlow::Break(chg_one) => γ.vec.iter_mut().for_each(|ρ| { + let t = ρ.α_γ.abs(); + if t > 0.0 { + if tconfig.allow_partial_transport { + let new = (t - chg_one).max(0.0); + ρ.α_γ = ρ.α_γ.signum() * new; + } + } + }), + } + } + } else { + // This version zeroes smallest weights, avoiding partial transport. + let mut abs_weights_idx = γ + .vec + .iter() + .map(|ρ| ρ.α_γ.abs()) + .zip(0..) + .filter(|(w, _)| *w >= 0.0) + .collect::<Vec<(F, usize)>>(); + abs_weights_idx.sort_by(|(a, _), (b, _)| a.total_cmp(b)); + + let mut left = reduction_target; + + for (w, i) in abs_weights_idx { + left -= w; + let ρ = &mut γ.vec[i]; + ρ.α_γ = 0.0; + if left < 0.0 { + break; + } + } + } + + all_ok0 = false + } + all_ok0 + }; + + if !all_ok { + if *attempts >= tconfig.max_attempts { + for ρ in γ.iter_mut() { + ρ.α_γ = 0.0; + } + } + } else { + for ρ in γ.iter_mut() { + if ρ.α_γ == 0.0 { + ρ.fail_count += 1; + } else { + ρ.fail_count = 0; + } + } + } + all_ok + } +} + +#[replace_float_literals(F::cast_from(literal))] +impl<const N: usize, F: Float> Transport<Loc<N, F>, F> { pub(crate) fn new() -> Self { Transport { vec: Vec::new() } } - pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<N, F>> { + pub(crate) fn iter(&self) -> impl Iterator<Item = &'_ SingleTransport<Loc<N, F>, F>> { self.vec.iter() } - pub(crate) fn iter_mut(&mut self) -> impl Iterator<Item = &'_ mut SingleTransport<N, F>> { + pub(crate) fn iter_mut( + &mut self, + ) -> impl Iterator<Item = &'_ mut SingleTransport<Loc<N, F>, F>> { self.vec.iter_mut() } pub(crate) fn extend<I>(&mut self, it: I) where - I: IntoIterator<Item = SingleTransport<N, F>>, + I: IntoIterator<Item = SingleTransport<Loc<N, F>, F>>, { self.vec.extend(it) } @@ -171,7 +639,6 @@ // } /// Construct `μ̆`, replacing the contents of `μ`. - #[replace_float_literals(F::cast_from(literal))] pub(crate) fn μ̆_into(&self, μ: &mut RNDM<N, F>) { assert!(self.len() <= μ.len()); @@ -190,287 +657,28 @@ // Then source points with partial transport let mut i = self.len(); - if ALLOW_PARTIAL_TRANSPORT { - // This can cause the number of points to explode, so cannot have partial transport. - for ρ in self.iter() { - let α = ρ.α_μ_orig - ρ.α_γ; - if ρ.α_γ.abs() > F::EPSILON && α != 0.0 { - let δ = DeltaMeasure { α, x: ρ.x }; - if i < μ.len() { - μ[i] = δ; - } else { - μ.push(δ) - } - i += 1; + // This can cause the number of points to explode, so cannot have partial transport. + for ρ in self.iter() { + let α = ρ.α_μ_orig - ρ.α_γ; + if ρ.α_γ.abs() > F::EPSILON && α.abs() > F::EPSILON { + let δ = DeltaMeasure { α, x: ρ.x }; + if i < μ.len() { + μ[i] = δ; + } else { + μ.push(δ) } + i += 1; } } μ.truncate(i); } - /// Constrution of initial transport `γ1` from initial measure `μ` and `v=F'(μ)` - /// with step lengh τ and transport step length `θ_or_adaptive`. - #[replace_float_literals(F::cast_from(literal))] - pub(crate) fn initial_transport<G, D>( - &mut self, - μ: &RNDM<N, F>, - _τ: F, - τθ_or_adaptive: &mut TransportStepLength<F, G>, - v: D, - tconfig: &TransportConfig<F>, - ) where - G: Fn(F, F) -> F, - D: DifferentiableRealMapping<N, F>, - { - use TransportStepLength::*; - - // Initialise transport structure weights - for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) { - ρ.α_μ_orig = δ.α; - ρ.x = δ.x; - if ρ.fail_count > tconfig.max_fail { - ρ.α_γ = 0.0 - } else { - // If old transport has opposing sign, the new transport will be none. - ρ.α_γ = if (ρ.α_γ > 0.0 && δ.α < 0.0) || (ρ.α_γ < 0.0 && δ.α > 0.0) { - 0.0 - } else { - δ.α - } - } - } - - let γ_prev_len = self.len(); - assert!(μ.len() >= γ_prev_len); - self.extend(μ[γ_prev_len..].iter().map(|δ| SingleTransport { - x: δ.x, - y: δ.x, // Just something, will be filled properly in the next phase - α_μ_orig: δ.α, - α_γ: δ.α, - prune: false, - fail_count: 0, - })); - - // Calculate transport rays. - match *τθ_or_adaptive { - Fixed(θ) => { - for ρ in self.iter_mut() { - if ρ.fail_count <= tconfig.max_fail { - ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θ); - } - } - } - AdaptiveMax { l: ℓ_F, ref mut max_transport, g: ref calculate_θτ } => { - *max_transport = max_transport.max(self.norm(Radon)); - let θτ = calculate_θτ(ℓ_F, *max_transport); - for ρ in self.iter_mut() { - if ρ.fail_count <= tconfig.max_fail { - ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * θτ); - } - } - } - FullyAdaptive { - l: ref mut adaptive_ℓ_F, - ref mut max_transport, - g: ref calculate_θτ, - } => { - *max_transport = max_transport.max(self.norm(Radon)); - let mut θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport); - // Do two runs through the spikes to update θ, breaking if first run did not cause - // a change. - for _i in 0..=1 { - let mut changes = false; - for ρ in self.iter_mut() { - if ρ.fail_count < tconfig.max_fail { - let dv_x = v.differential(&ρ.x); - let g = &dv_x * (ρ.α_γ.signum() * θτ); - ρ.y = ρ.x - g; - let n = g.norm2(); - if n >= F::EPSILON { - // Estimate Lipschitz factor of ∇v - let this_ℓ_F = (dv_x - v.differential(&ρ.y)).norm2() / n; - *adaptive_ℓ_F = adaptive_ℓ_F.max(this_ℓ_F); - θτ = calculate_θτ(*adaptive_ℓ_F, *max_transport); - changes = true - } - } - } - if !changes { - break; - } - } - } - } - } - - /// A posteriori transport adaptation. - #[replace_float_literals(F::cast_from(literal))] - pub(crate) fn aposteriori_transport<D>( - &mut self, - μ: &RNDM<N, F>, - μ̆: &RNDM<N, F>, - _v: &mut D, - extra: Option<F>, - ε: F, - tconfig: &TransportConfig<F>, - attempts: &mut usize, - ) -> bool - where - D: DifferentiableRealMapping<N, F>, - { - *attempts += 1; - - // 1. If π_♯^1γ^{k+1} = γ1 has non-zero mass at some point y, but μ = μ^{k+1} does not, - // then the ansatz ∇w̃_x(y) = w^{k+1}(y) may not be satisfied. So set the mass of γ1 - // at that point to zero, and retry. - let mut all_ok = true; - for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) { - if δ.α == 0.0 && ρ.α_γ != 0.0 { - all_ok = false; - ρ.α_γ = 0.0; - } - } - - // 2. Through bounding ∫ B_ω(y, z) dλ(x, y, z). - // through the estimate ≤ C ‖Δ‖‖γ^{k+1}‖ for Δ := μ^{k+1}-μ̆^k - // which holds for some some C if the convolution kernel in 𝒟 has Lipschitz gradient. - let nγ = self.norm(Radon); - let nΔ = μ.dist_matching(&μ̆) + extra.unwrap_or(0.0); - let t = ε * tconfig.tolerance_mult_con; - if nγ * nΔ > t && *attempts >= tconfig.max_attempts { - all_ok = false; - } else if nγ * nΔ > t { - // Since t/(nγ*nΔ)<1, and the constant tconfig.adaptation < 1, - // this will guarantee that eventually ‖γ‖ decreases sufficiently that we - // will not enter here. - //*self *= tconfig.adaptation * t / (nγ * nΔ); - - // We want a consistent behaviour that has the potential to set many weights to zero. - // Therefore, we find the smallest uniform reduction `chg_one`, subtracted - // from all weights, that achieves total `adapt` adaptation. - let adapt_to = tconfig.adaptation * t / nΔ; - let reduction_target = nγ - adapt_to; - assert!(reduction_target > 0.0); - if ALLOW_PARTIAL_TRANSPORT { - if MINIMAL_PARTIAL_TRANSPORT { - // This reduces weights of transport, starting from … until `adapt` is - // exhausted. It will, therefore, only ever cause one extrap point insertion - // at the sources, unlike “full” partial transport. - //let refs = self.vec.iter_mut().collect::<Vec<_>>(); - //refs.sort_by(|ρ1, ρ2| ρ1.α_γ.abs().partial_cmp(&ρ2.α_γ.abs()).unwrap()); - // let mut it = refs.into_iter(); - // - // Maybe sort by differential norm - // let mut refs = self - // .vec - // .iter_mut() - // .map(|ρ| { - // let val = v.differential(&ρ.x).norm2_squared(); - // (ρ, val) - // }) - // .collect::<Vec<_>>(); - // refs.sort_by(|(_, v1), (_, v2)| v2.partial_cmp(&v1).unwrap()); - // let mut it = refs.into_iter().map(|(ρ, _)| ρ); - let mut it = self.vec.iter_mut().rev(); - let _unused = it.try_fold(reduction_target, |left, ρ| { - let w = ρ.α_γ.abs(); - if left <= w { - ρ.α_γ = ρ.α_γ.signum() * (w - left); - ControlFlow::Break(()) - } else { - ρ.α_γ = 0.0; - ControlFlow::Continue(left - w) - } - }); - } else { - // This version equally reduces all weights. It causes partial transport, which - // has the problem that that we need to then adapt weights in both start and - // end points, in insert_and_reweigh, somtimes causing the number of spikes μ - // to explode. - let mut abs_weights = self - .vec - .iter() - .map(|ρ| ρ.α_γ.abs()) - .filter(|t| *t > F::EPSILON) - .collect::<Vec<F>>(); - abs_weights.sort_by(|a, b| a.partial_cmp(b).unwrap()); - let n = abs_weights.len(); - // Cannot have partial transport; can cause spike count explosion - let chg = abs_weights.into_iter().zip((1..=n).rev()).try_fold( - 0.0, - |smaller_total, (w, m)| { - let mf = F::cast_from(m); - let reduction = w * mf + smaller_total; - if reduction >= reduction_target { - ControlFlow::Break((reduction_target - smaller_total) / mf) - } else { - ControlFlow::Continue(smaller_total + w) - } - }, - ); - match chg { - ControlFlow::Continue(_) => self.vec.iter_mut().for_each(|δ| δ.α_γ = 0.0), - ControlFlow::Break(chg_one) => self.vec.iter_mut().for_each(|ρ| { - let t = ρ.α_γ.abs(); - if t > 0.0 { - if ALLOW_PARTIAL_TRANSPORT { - let new = (t - chg_one).max(0.0); - ρ.α_γ = ρ.α_γ.signum() * new; - } - } - }), - } - } - } else { - // This version zeroes smallest weights, avoiding partial transport. - let mut abs_weights_idx = self - .vec - .iter() - .map(|ρ| ρ.α_γ.abs()) - .zip(0..) - .filter(|(w, _)| *w >= 0.0) - .collect::<Vec<(F, usize)>>(); - abs_weights_idx.sort_by(|(a, _), (b, _)| a.partial_cmp(b).unwrap()); - - let mut left = reduction_target; - - for (w, i) in abs_weights_idx { - left -= w; - let ρ = &mut self.vec[i]; - ρ.α_γ = 0.0; - if left < 0.0 { - break; - } - } - } - - all_ok = false - } - - if !all_ok && *attempts >= tconfig.max_attempts { - for ρ in self.iter_mut() { - ρ.α_γ = 0.0; - } - } - - for ρ in self.iter_mut() { - if ρ.α_γ == 0.0 { - ρ.fail_count += 1; - } else if all_ok { - ρ.fail_count = 0; - } - } - - all_ok - } - /// Returns $‖μ\^k - π\_♯\^0γ\^{k+1}‖$ pub(crate) fn μ0_minus_γ0_radon(&self) -> F { self.vec.iter().map(|ρ| (ρ.α_μ_orig - ρ.α_γ).abs()).sum() } /// Returns $∫ c_2 d|γ|$ - #[replace_float_literals(F::cast_from(literal))] pub(crate) fn c2integral(&self) -> F { self.vec .iter() @@ -478,7 +686,6 @@ .sum() } - #[replace_float_literals(F::cast_from(literal))] pub(crate) fn get_transport_stats(&self, stats: &mut IterInfo<F>, μ: &RNDM<N, F>) { // TODO: This doesn't take into account μ[i].α becoming zero in the latest tranport // attempt, for i < self.len(), when a corresponding source term also exists with index @@ -522,22 +729,229 @@ assert!(self.vec.len() <= μ.len()); let old_len = μ.len(); for (ρ, δ) in self.vec.iter_mut().zip(μ.iter_spikes()) { - ρ.prune = !(δ.α.abs() > F::EPSILON); + ρ.retain = δ.α.abs() > F::EPSILON; } μ.prune_by(|δ| δ.α.abs() > F::EPSILON); stats.pruned += old_len - μ.len(); - self.vec.retain(|ρ| !ρ.prune); + self.vec.retain(|ρ| ρ.retain); assert!(self.vec.len() <= μ.len()); } + + /// Helper for initial transport. Called from [`TransportProxPenalty::initial_transport`]. + pub(super) fn do_init_transport<D>( + &mut self, + v: &D, + μ: &RNDM<N, F>, + τθ_or_adaptive: &mut TransportStepLength<F>, + tconfig: &TransportConfig<F>, + ) where + D: DifferentiableRealMapping<N, F>, + { + use TransportStepLength::*; + + // Initialise transport structure weights + for (δ, ρ) in izip!(μ.iter_spikes(), self.iter_mut()) { + ρ.α_μ_orig = δ.α; + ρ.x = δ.x; + ρ.y = δ.x; // Later updated if no fails. + ρ.α_γ = if ρ.fail_count > tconfig.max_fail { + 0.0 + } else { + // If old transport has opposing sign, the new transport will be none. + if (ρ.α_γ > 0.0 && δ.α < 0.0) || (ρ.α_γ < 0.0 && δ.α > 0.0) { + 0.0 + } else { + δ.α + } + }; + } + + let γ_prev_len = self.len(); + assert!(μ.len() >= γ_prev_len); + self.extend(μ[γ_prev_len..].iter().map(|δ| SingleTransport { + x: δ.x, + y: δ.x, // Just something, will be filled properly in the next phase + α_μ_orig: δ.α, + α_γ: δ.α, + retain: true, + fail_count: 0, + excess: 0.0, + })); + + // Calculate transport rays. + let simple_τθ = match *τθ_or_adaptive { + Fixed { τθ, .. } => Some(τθ), + Simple { ℓ_gradv, τθ0, ℓ_base } => Some(τθ0 / (ℓ_gradv + ℓ_base)), + AdaptiveMax { + ℓ_gradv, + ref mut adaptive_max_transport, + τθ0, + ℓ_base, + ℓ_base_max_transport, + } => { + *adaptive_max_transport = adaptive_max_transport.max(self.norm(Radon)); + Some(τθ0 / (ℓ_gradv + ℓ_base + ℓ_base_max_transport * *adaptive_max_transport)) + } + FullyAdaptive { + ref mut adaptive_ℓ_gradv, + ref mut adaptive_max_transport, + τθ0, + ℓ_base, + ℓ_base_max_transport, + } => { + *adaptive_max_transport = adaptive_max_transport.max(self.norm(Radon)); + let mut τθ = τθ0 + / (*adaptive_ℓ_gradv + ℓ_base + ℓ_base_max_transport * *adaptive_max_transport); + // Do two runs through the spikes to update θ, breaking if first run did not cause + // a change. + for _i in 0..=1 { + let mut changes = false; + for ρ in self.iter_mut() { + if ρ.fail_count < tconfig.max_fail { + let dv_x = v.differential(&ρ.x); + let g = &dv_x * (ρ.α_γ.signum() * τθ); + ρ.y = ρ.x - g; + let n = g.norm2(); + if n >= F::EPSILON { + // Estimate Lipschitz factor of ∇v + let this_ℓ_gradv = (dv_x - v.differential(&ρ.y)).norm2() / n; + *adaptive_ℓ_gradv = adaptive_ℓ_gradv.max(this_ℓ_gradv); + τθ = τθ0 + / (*adaptive_ℓ_gradv + + ℓ_base + + ℓ_base_max_transport * *adaptive_max_transport); + changes = true + } + } + } + if !changes { + break; + } + } + None + } + }; + + if let Some(τθ) = simple_τθ { + for ρ in self.iter_mut() { + if ρ.fail_count <= tconfig.max_fail { + ρ.y = ρ.x - v.differential(&ρ.x) * (ρ.α_γ.signum() * τθ); + } + } + } + } + + /// Helper for a posteriori transport error control. + /// Called from [`TransportProxPenalty::aposteriori_transport`]. + fn do_new_aposteriori_transport<Reg, M>( + &mut self, + μ: &RNDM<N, F>, + τv̆: &mut M, + v: &mut M, + ε: F, + τ: F, + τθ_or_adaptive: &TransportStepLength<F>, + reg: &Reg, + tconfig: &TransportConfig<F>, + ω: impl Fn(&Loc<N, F>, std::cmp::Ordering, bool) -> F, + ) -> bool + where + Reg: SlidingRegTerm<Loc<N, F>, F>, + F: ToNalgebraRealField, + M: MinMaxMapping<Loc<N, F>, F> + DifferentiableRealMapping<N, F>, + { + let τℓ_gradv2 = τ * τθ_or_adaptive.get_ℓ_gradv() / 2.0; + let Bounds(α_lower, α_upper) = reg.subdiff_range(); + + let (all_ok, m) = izip!(self.vec.iter_mut(), μ.iter_spikes()).fold( + (true, 0.0), + |(all_ok_so_far, total_excess), (ρ, δ)| { + use std::cmp::Ordering::*; + // NOTE: The tolerances ε are commented out, because they will in any case + // be consumed by `t` below by suitably large choice of `tolerance_mult`. + // Hence, we simply implicitly adapt the `tolerance_mult` be commenting out + // the `ε` here. That way, `d` is in its entirely multiplied by τ, as is `t`, + // making all factors independent of τ. + let maybe_excess = match ρ.α_γ.total_cmp(&0.0) { + Greater => (δ.α >= 0.0).then(|| { + let gvx = v.differential(&ρ.x); + let d = if tconfig.alt_remainder_control { + let τv̆y = τv̆.apply(&ρ.y); + (/*ε +*/τ * α_upper + ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), false)) + + (τv̆y + τ * gvx.dot(&ρ.x - &ρ.y)) + - τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x) + } else { + (/*2.0 * ε*/-ω(&ρ.y, δ.α.total_cmp(&ρ.α_γ), true) + + ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), false)) + + τ * gvx.dot(&ρ.x - &ρ.y) + - τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x) + }; + d * ρ.α_γ + }), + Less => (δ.α <= 0.0).then(|| { + let gvx = v.differential(&ρ.x); + let d = if tconfig.alt_remainder_control { + let τv̆y = τv̆.apply(&ρ.y); + (/*ε*/-τ * α_lower - ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), true)) + - (τv̆y + τ * gvx.dot(&ρ.x - &ρ.y)) + - τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x) + } else { + (/*2.0 * ε +*/ω(&ρ.y, δ.α.total_cmp(&ρ.α_γ), false) + - ω(&ρ.x, ρ.α_γ.total_cmp(&ρ.α_μ_orig), true)) + - τ * gvx.dot(&ρ.x - &ρ.y) + - τℓ_gradv2 * ρ.y.dist2_squared(&ρ.x) + }; + d * (-ρ.α_γ) + }), + Equal => Some(0.0), + }; + match maybe_excess { + None => { + ρ.α_γ = 0.0; + ρ.excess = 0.0; + (false, total_excess) + } + Some(e) => { + ρ.excess = e; + (all_ok_so_far, total_excess + e) + } + } + }, + ); + + let t = τ * ε * tconfig.tolerance_mult; + + if m > t { + let mut it = self.vec.iter_mut().rev().filter(|ρ| ρ.excess > 0.0); + let reduction_target = m - tconfig.adaptation * t; + let _unused = it.try_fold(reduction_target, |left, ρ| { + let d = ρ.excess; + if d >= left { + if tconfig.allow_partial_transport { + ρ.α_γ *= (d - left) / d; + } else { + ρ.α_γ = 0.0; + } + ControlFlow::Break(()) + } else { + ρ.α_γ = 0.0; + ControlFlow::Continue(left - d) + } + }); + false + } else { + all_ok + } + } } -impl<const N: usize, F: Float> Norm<Radon, F> for Transport<N, F> { +impl<const N: usize, F: Float> Norm<Radon, F> for Transport<Loc<N, F>, F> { fn norm(&self, _: Radon) -> F { self.iter().map(|ρ| ρ.α_γ.abs()).sum() } } -impl<const N: usize, F: Float> MulAssign<F> for Transport<N, F> { +impl<const N: usize, F: Float> MulAssign<F> for Transport<Loc<N, F>, F> { fn mul_assign(&mut self, factor: F) { for ρ in self.iter_mut() { ρ.α_γ *= factor; @@ -564,11 +978,12 @@ F: Float + ToNalgebraRealField, I: AlgIteratorFactory<IterInfo<F>>, Dat: DifferentiableMapping<RNDM<N, F>, Codomain = F> + BoundedCurvature<F>, - Dat::DerivativeDomain: DifferentiableRealMapping<N, F> + ClosedMul<F>, + Dat::DerivativeDomain: + DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>, //for<'a> Dat::Differential<'a>: Lipschitz<&'a P, FloatType = F>, RNDM<N, F>: SpikeMerging<F>, Reg: SlidingRegTerm<Loc<N, F>, F>, - P: ProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>, + P: TransportProxPenalty<Loc<N, F>, Dat::DerivativeDomain, Reg, F> + StepLengthBound<F, Dat>, Plot: Plotter<P::ReturnMapping, Dat::DerivativeDomain, RNDM<N, F>>, { // Check parameters @@ -584,30 +999,14 @@ //let max_transport = config.max_transport.scale // * reg.radon_norm_bound(b.norm2_squared() / 2.0); //let ℓ = opA.transport.lipschitz_factor(L2Squared) * max_transport; - let ℓ = 0.0; let τ = config.τ0 / prox_penalty.step_length_bound(&f)?; - let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) { - (_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0), - (maybe_ℓ_F, Ok(transport_lip)) => { - let calculate_θτ = move |ℓ_F, max_transport| { - let ℓ_r = transport_lip * max_transport; - config.transport.θ0 / (ℓ + ℓ_F + ℓ_r) - }; - match maybe_ℓ_F { - Ok(ℓ_F) => TransportStepLength::AdaptiveMax { - l: ℓ_F, // TODO: could estimate computing the real reesidual - max_transport: 0.0, - g: calculate_θτ, - }, - Err(_) => TransportStepLength::FullyAdaptive { - l: 10.0 * F::EPSILON, // Start with something very small to estimate differentials - max_transport: 0.0, - g: calculate_θτ, - }, - } - } - }; + let mut τθ_or_adaptive = prox_penalty.get_transport_steplength( + f.curvature_bound_components(config.guess), + &config.transport, + 0.0, + 0.0, + ); // We multiply tolerance by τ for FB since our subproblems depending on tolerances are scaled // by τ compared to the conditional gradient approach. let tolerance = config.insertion.tolerance * τ * reg.tolerance_scaling(); @@ -625,9 +1024,18 @@ // Run the algorithm for state in iterator.iter_init(|| full_stats(&μ, ε, stats.clone())) { + let mut v = f.differential(&μ); + // Calculate initial transport - let v = f.differential(&μ); - γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v, &config.transport); + prox_penalty.initial_transport( + &mut γ, + &μ, + ε, + τ, + &mut τθ_or_adaptive, + &v, + &config.transport, + ); let mut attempts = 0; @@ -658,8 +1066,21 @@ )?; // A posteriori transport adaptation. - if γ.aposteriori_transport(&μ, &μ̆, &mut τv̆, None, ε, &config.transport, &mut attempts) - { + if prox_penalty.aposteriori_transport( + &mut γ, + &μ, + &μ̆, + &mut τv̆, + &mut v, + None, + ε, + τ, + &τθ_or_adaptive, + reg, + &config.transport, + &config.insertion.refinement, + &mut attempts, + ) { break 'adapt_transport (maybe_d, within_tolerances, τv̆, μ̆); } @@ -672,7 +1093,7 @@ // 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 γ. if config.insertion.merge_now(&state) { - stats.merged += prox_penalty.merge_spikes( + let m = prox_penalty.merge_spikes( &mut μ, &mut τv̆, &μ̆, @@ -682,6 +1103,10 @@ ®, Some(|μ̃: &RNDM<N, F>| f.apply(μ̃)), ); + if m > 0 { + stats.merged += m; + v = f.differential(&μ); + } } γ.prune_compat(&mut μ, &mut stats);
--- a/src/sliding_pdps.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/sliding_pdps.rs Sun Jul 19 07:34:39 2026 +0200 @@ -8,10 +8,11 @@ use crate::measures::merging::SpikeMerging; use crate::measures::{DiscreteMeasure, RNDM}; use crate::plot::Plotter; -use crate::prox_penalty::{ProxPenalty, StepLengthBoundPair}; +use crate::prox_penalty::StepLengthBoundPair; use crate::regularisation::SlidingRegTerm; -use crate::sliding_fb::{SlidingFBConfig, Transport, TransportConfig, TransportStepLength}; +use crate::sliding_fb::{SlidingFBConfig, Transport, TransportConfig, TransportProxPenalty}; use crate::types::*; +use alg_tools::bounds::MinMaxMapping; use alg_tools::convex::{Conjugable, Prox, Zero}; use alg_tools::direct_product::Pair; use alg_tools::error::DynResult; @@ -100,12 +101,12 @@ I: AlgIteratorFactory<IterInfo<F>>, Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>> + BoundedCurvature<F>, - S: DifferentiableRealMapping<N, F> + ClosedMul<F>, + S: DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>, for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, //Pair<S, Z>: ClosedMul<F>, RNDM<N, F>: SpikeMerging<F>, Reg: SlidingRegTerm<Loc<N, F>, F>, - P: ProxPenalty<Loc<N, F>, S, Reg, F>, + P: TransportProxPenalty<Loc<N, F>, S, Reg, F>, // KOpM : Linear<RNDM<N, F>, Codomain=Y> // + GEMV<F, RNDM<N, F>> // + Preadjointable< @@ -122,7 +123,7 @@ KOpZ::SimpleAdjoint: GEMV<F, Y>, Y: ClosedEuclidean<F>, for<'b> &'b Y: Instance<Y>, - Z: ClosedEuclidean<F>, + Z: ClosedEuclidean<F> + AXPY, for<'b> &'b Z: Instance<Z>, R: Prox<Z, Codomain = F>, H: Conjugable<Y, F, Codomain = F>, @@ -153,7 +154,6 @@ let bigM = 0.0; //opKμ.adjoint_product_bound(&op𝒟).unwrap().sqrt(); let nKz = opKz.opnorm_bound(L2, L2)?; let is_fb = nKz == 0.0; - let ℓ = 0.0; let idOpZ = IdOp::new(); let opKz_adj = opKz.adjoint(); let (l, l_z) = Pair(prox_penalty, &idOpZ).step_length_bound_pair(&f)?; @@ -183,27 +183,13 @@ // The factor two in the manuscript disappears due to the definition of 𝚹 being // for ‖x-y‖₂² instead of c_2(x, y)=‖x-y‖₂²/2. - let mut θ_or_adaptive = match f.curvature_bound_components(config.guess) { - (_, Err(_)) => TransportStepLength::Fixed(config.transport.θ0), - (maybe_ℓ_F, Ok(transport_lip)) => { - let calculate_θτ = move |ℓ_F, max_transport| { - let ℓ_r = transport_lip * max_transport; - config.transport.θ0 / ((ℓ + ℓ_F + ℓ_r) + κ * bigθ * max_transport / τ) - }; - match maybe_ℓ_F { - Ok(ℓ_F) => TransportStepLength::AdaptiveMax { - l: ℓ_F, // TODO: could estimate computing the real reesidual - max_transport: 0.0, - g: calculate_θτ, - }, - Err(_) => TransportStepLength::FullyAdaptive { - l: F::EPSILON, // Start with something very small to estimate differentials - max_transport: 0.0, - g: calculate_θτ, - }, - } - } - }; + let mut τθ_or_adaptive = prox_penalty.get_transport_steplength( + f.curvature_bound_components(config.guess), + &config.transport, + 0.0, + κ * bigθ / τ, // = 0 currently + ); + // Acceleration is not currently supported // let γ = dataterm.factor_of_strong_convexity(); let ω = 1.0; @@ -230,8 +216,8 @@ // Run the algorithm for state in iterator.iter_init(|| full_stats(&μ, &z, ε, stats.clone())) { - // Calculate initial transport - let Pair(v, _) = f.differential(Pair(&μ, &z)); + let Pair(mut v, _) = f.differential(Pair(&μ, &z)); + //opKμ.preadjoint().apply_add(&mut v, y); // We want to proceed as in Example 4.12 but with v and v̆ as in §5. // With A(ν, z) = A_μ ν + A_z z, following Example 5.1, we have @@ -242,7 +228,15 @@ //dbg!(&μ); - γ.initial_transport(&μ, τ, &mut θ_or_adaptive, v, &config.transport); + prox_penalty.initial_transport( + &mut γ, + &μ, + ε, + τ, + &mut τθ_or_adaptive, + &v, + &config.transport, + ); let mut attempts = 0; @@ -254,7 +248,8 @@ let μ̆ = μ.clone(); // Calculate τv̆ = τA_*(A[μ_transported + μ_transported_base]-b) - let Pair(mut τv̆, τz̆) = f.differential(Pair(&μ̆, &z)) * τ; + let Pair(v̆, mut z̆) = f.differential(Pair(&μ̆, &z)); + let mut τv̆ = v̆ * τ; // opKμ.preadjoint().gemv(&mut τv̆, τ, y, 1.0); // Construct μ^{k+1} by solving finite-dimensional subproblems and insert new spikes. @@ -270,18 +265,25 @@ )?; // Do z variable primal update here to able to estimate B_{v̆^k-v^{k+1}} - let mut z_new = τz̆; - opKz_adj.gemv(&mut z_new, -σ_p, &y, -σ_p / τ); - z_new = fnR.prox(σ_p, z_new + &z); + opKz_adj.apply_add(&mut z̆, &y); + z̆.axpy(1.0, &z, -σ_p); + let z_new = fnR.prox(σ_p, z̆); + //let z_new = fnR.prox(σ_p, z̆ * (-σ_p) + &z); // A posteriori transport adaptation. - if γ.aposteriori_transport( + if prox_penalty.aposteriori_transport( + &mut γ, &μ, &μ̆, &mut τv̆, - Some(z_new.dist2(&z)), + &mut v, + None, //Some(z_new.dist2(&z)), ε, + τ, + &τθ_or_adaptive, + reg, &config.transport, + &config.insertion.refinement, &mut attempts, ) { break 'adapt_transport (maybe_d, within_tolerances, τv̆, z_new, μ̆); @@ -294,7 +296,7 @@ // 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 γ. if config.insertion.merge_now(&state) { - stats.merged += prox_penalty.merge_spikes( + let m = prox_penalty.merge_spikes( &mut μ, &mut τv̆, &μ̆, @@ -302,17 +304,20 @@ ε, &config.insertion, ®, - is_fb.then_some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))), + Some(|μ̃: &RNDM<N, F>| f.apply(Pair(μ̃, &z))), ); + if m > 0 { + stats.merged += m; + } } γ.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}] - opKz.gemv(&mut y, σ_d * (1.0 + ω), &z_new, 1.0); - // opKμ.gemv(&mut y, -σ_d*ω, μ_base, 1.0);// y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b - opKz.gemv(&mut y, -σ_d * ω, z, 1.0); // y = y + σ_d K[(1+ω)(μ,z)^{k+1} - ω (μ,z)^k]-b + 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; @@ -364,10 +369,10 @@ I: AlgIteratorFactory<IterInfo<F>>, Dat: DifferentiableMapping<MeasureZ<F, Z, N>, Codomain = F, DerivativeDomain = Pair<S, Z>> + BoundedCurvature<F>, - S: DifferentiableRealMapping<N, F> + ClosedMul<F>, + S: DifferentiableRealMapping<N, F> + ClosedMul<F> + MinMaxMapping<Loc<N, F>, F>, RNDM<N, F>: SpikeMerging<F>, Reg: SlidingRegTerm<Loc<N, F>, F>, - P: ProxPenalty<Loc<N, F>, S, Reg, F>, + P: TransportProxPenalty<Loc<N, F>, S, Reg, F>, for<'a> Pair<&'a P, &'a IdOp<Z>>: StepLengthBoundPair<F, Dat>, Z: ClosedEuclidean<F> + AXPY + Clone, for<'b> &'b Z: Instance<Z>,
--- a/src/subproblem.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/subproblem.rs Sun Jul 19 07:34:39 2026 +0200 @@ -19,6 +19,8 @@ #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug, ValueEnum)] #[allow(dead_code)] pub enum InnerMethod { + /// Exact solver (when available) + Exact, /// Forward-backward FB, /// Semismooth Newton
--- a/src/subproblem/l1squared_nonneg.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/subproblem/l1squared_nonneg.rs Sun Jul 19 07:34:39 2026 +0200 @@ -3,13 +3,11 @@ */ use itertools::izip; -use nalgebra::DVector; +use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1}; use numeric_literals::replace_float_literals; -//use std::iter::zip; -use std::cmp::Ordering::*; use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState}; -use alg_tools::nalgebra_support::ToNalgebraRealField; +use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField}; use alg_tools::norms::{Dist, L1}; use super::l1squared_unconstrained::l1squared_prox; @@ -19,6 +17,7 @@ /// Return maximum of `dist` and distnce of inteval `[lb, ub]` to zero. #[replace_float_literals(F::cast_from(literal))] +#[inline] pub(super) fn max_interval_dist_to_zero<F: Float>(dist: F, lb: F, ub: F) -> F { if lb < 0.0 { if ub > 0.0 { @@ -27,7 +26,7 @@ dist.max(-ub) } } else - /* ub ≥ 0.0*/ + /* lb ≥ 0.0*/ { dist.max(lb) } @@ -37,44 +36,35 @@ /// /// `v` will be modified and cannot be trusted to contain useful values afterwards. #[replace_float_literals(F::cast_from(literal))] -fn min_subdifferential<F: Float + nalgebra::RealField>( - y: &DVector<F>, - x: &DVector<F>, - g: &DVector<F>, +fn min_subdifferential<F: Float + nalgebra::RealField, S1, S2, S3>( + y: &Vector<F, Dyn, S1>, + x: &Vector<F, Dyn, S2>, + g: &Vector<F, Dyn, S3>, λ: F, - β: F, -) -> F { +) -> F +where + S1: Storage<F, Dyn>, + S2: Storage<F, Dyn>, + S3: Storage<F, Dyn>, + ShapeConstraint: StridesOk<F, Dyn, U1, S2>, +{ let mut val = 0.0; - let tmp = β * y.dist(x, L1); + let tmp = y.dist(x, L1); for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) { - let (mut lb, mut ub) = (-g_i, -g_i); - match x_i.partial_cmp(y_i) { - Some(Greater) => { - lb += tmp; - ub += tmp - } - Some(Less) => { - lb -= tmp; - ub -= tmp - } - Some(Equal) => { - lb -= tmp; - ub += tmp - } - None => {} + let (mut lb, mut ub) = (-g_i + λ, -g_i + λ); + if num_traits::abs(x_i - *y_i) < F::EPSILON { + lb -= tmp; + ub += tmp + } else if x_i > *y_i { + lb += tmp; + ub += tmp + } else { + lb -= tmp; + ub -= tmp } - match x_i.partial_cmp(&0.0) { - Some(Greater) => { - lb += λ; - ub += λ - } - // Less should not happen - Some(Less | Equal) => { - lb = F::NEG_INFINITY; - ub += λ - } - None => {} - }; + if x_i < F::EPSILON { + lb = F::NEG_INFINITY; + } val = max_interval_dist_to_zero(val, lb, ub); } val @@ -164,11 +154,14 @@ /// (`shift` in the code), and below a value that would cause changes in the locked set /// (`max_shift` in the code). #[replace_float_literals(F::cast_from(literal))] -pub fn l1squared_nonneg_prox<F: Float + nalgebra::RealField>( - x: &mut DVector<F>, - y: &DVector<F>, +pub fn l1squared_nonneg_prox<F: Float + nalgebra::RealField, S1, S2>( + x: &mut Vector<F, Dyn, S1>, + y: &Vector<F, Dyn, S2>, β: F, -) { +) where + S2: Storage<F, Dyn>, + S1: StorageMut<F, Dyn>, +{ // nalgebra double-definition bullshit workaround let abs = alg_tools::NumTraitsFloat::abs; let min = alg_tools::NumTraitsFloat::min; @@ -218,12 +211,11 @@ /// The `λ` component of the model is handled in the proximal step instead of the gradient step /// for potential performance improvements. #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] -pub fn l1squared_nonneg_pp<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_nonneg_pp<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ_: F, - β_: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, τ_: F, θ_: F, iterator: I, @@ -231,9 +223,12 @@ where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { let λ = λ_.to_nalgebra_mixed(); - let β = β_.to_nalgebra_mixed(); let mut τ = τ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); let mut iters = 0; @@ -242,7 +237,7 @@ // Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2+δ_{≥0}+}(x^k - τ(λ𝟙^⊤-g)) x.apply(|x_i| *x_i -= τ * λ); x.axpy(τ, g, 1.0); - l1squared_nonneg_prox(x, y, τ * β); + l1squared_nonneg_prox(x, y, τ); iters += 1; // This gives O(1/N^2) rates due to monotonicity of function values. @@ -253,7 +248,7 @@ // Higher acceleration does not seem to be numerically stable. //τ + = F::cast_from(iters).to_nalgebra_mixed()*θ; - state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ))) }); iters @@ -266,12 +261,11 @@ /// for potential performance improvements. /// The parameter `θ` is used to multiply the rescale the operator (identity) of the PDPS model. #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] -pub fn l1squared_nonneg_pdps<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_nonneg_pdps<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ_: F, - β_: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, τ_: F, σ_: F, θ_: F, @@ -280,22 +274,25 @@ where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { let λ = λ_.to_nalgebra_mixed(); - let β = β_.to_nalgebra_mixed(); let τ = τ_.to_nalgebra_mixed(); let σ = σ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); let mut w = DVector::zeros(x.len()); let mut tmp = DVector::zeros(x.len()); - let mut xprev = x.clone(); + let mut xprev = x.clone_owned(); let mut iters = 0; iterator.iterate(|state| { // Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2}(x^k - τ (w^k - g)) x.axpy(-τ * θ, &w, 1.0); x.axpy(τ, g, 1.0); - l1squared_prox(&mut tmp, x, y, τ * β); + l1squared_prox(&mut tmp, x, y, τ); // Dual step: w^{k+1} = proj_{[-∞,λ]}(w^k + σ(2x^{k+1}-x^k)) w.axpy(2.0 * σ * θ, x, 1.0); @@ -305,7 +302,7 @@ iters += 1; - state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ))) }); iters @@ -323,18 +320,17 @@ /// We rewrite /// <div>$$ /// \begin{split} -/// & \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x) \\ +/// & \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x) \\ /// & = \min_{x ∈ ℝ^n} \max_{w} ⟨θ w, x⟩ - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x) -/// - \left(x ↦ \frac{β}{2θ} |x-y|_1^2 \right)^*(w). +/// - \left(x ↦ \frac{1}{2θ} |x-y|_1^2 \right)^*(w). /// \end{split} /// $$</div> #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] -pub fn l1squared_nonneg_pdps_alt<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_nonneg_pdps_alt<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ_: F, - β_: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, τ_: F, σ_: F, θ_: F, @@ -343,80 +339,282 @@ where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { let λ = λ_.to_nalgebra_mixed(); let τ = τ_.to_nalgebra_mixed(); let σ = σ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); - let β = β_.to_nalgebra_mixed(); let σθ = σ * θ; - let τθ = τ * θ; + let τλ = τ * λ; + let one_div_σθ = 1.0 / σθ; let mut w = DVector::zeros(x.len()); let mut tmp = DVector::zeros(x.len()); - let mut xprev = x.clone(); + let mut xprev = x.clone_owned(); let mut iters = 0; + let mut y_scale = y.clone_owned(); + y_scale *= σ; + iterator.iterate(|state| { // Primal step: x^{k+1} = nonnegsoft_τλ(x^k - τ(θ w^k -g)) - x.axpy(-τθ, &w, 1.0); - x.axpy(τ, g, 1.0); - x.apply(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τ * λ)); + if θ == 1.0 { + for (x_i, xprev_i, w_i, &g_i) in + izip!(x.iter_mut(), xprev.iter_mut(), w.iter_mut(), g.iter()) + { + *x_i = nonneg_soft_thresholding(*x_i - τ * (*w_i - g_i), τλ); + // Fused dual part from below + *w_i += σ * (2.0 * *x_i - *xprev_i); + *xprev_i = *w_i; + } + } else { + for (x_i, xprev_i, w_i, &g_i) in + izip!(x.iter_mut(), xprev.iter_mut(), w.iter_mut(), g.iter()) + { + *x_i = nonneg_soft_thresholding(*x_i - τ * (θ * *w_i - g_i), τλ); + // Fused dual part from below + *w_i += σ * (2.0 * *x_i - *xprev_i); + *xprev_i = *w_i; + } + } + // This is numerically unstable: + // x.axpy(-τθ, &w, 1.0); + // x.axpy(τ, g, 1.0); + // x.apply(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τ * λ)); // Dual step: with g(x) = (β/(2θ))‖x-y‖₁² and q = w^k + σ(2x^{k+1}-x^k), // we compute w^{k+1} = prox_{σg^*}(q) for // = q - σ prox_{g/σ}(q/σ) // = q - σ prox_{(β/(2θσ))‖.-y‖₁²}(q/σ) // = σ(q/σ - prox_{(β/(2θσ))‖.-y‖₁²}(q/σ)) + // ALT + // = q - prox_{(β/(2θσ))‖.-σy‖₁²}(q) // where q/σ = w^k/σ + (2x^{k+1}-x^k), - w /= σ; - w.axpy(2.0, x, 1.0); - w.axpy(-1.0, &xprev, 1.0); - xprev.copy_from(&w); // use xprev as temporary variable - l1squared_prox(&mut tmp, &mut xprev, y, β / σθ); + + // This has been fused into the loop below + // for (xprev_i, w_i, &x_i) in izip!(xprev.iter_mut(), w.iter_mut(), x.iter()) { + // *w_i += σ * (2.0 * x_i - *xprev_i); + // *xprev_i = *w_i; + // } + // xprev.axpy(2.0, x, -1.0); + // w.axpy(σ, &xprev, 1.0); + // xprev.copy_from(&w); // use xprev as temporary variable + //l1squared_prox(&mut tmp, &mut xprev, &y_scale, β_div_σθ); + l1squared_prox(&mut tmp, &mut xprev, &y_scale, one_div_σθ); w -= &xprev; - w *= σ; xprev.copy_from(x); iters += 1; - state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ))) }); iters } +/// This is an exact solver for +/// <div>$$ +/// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x). +/// $$</div> +/// i.e., +/// <div>$$ +/// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 + (λ - g)^⊤ x + δ_{≥ 0}(x), +/// $$</div> +/// which has the optimality conditions +/// <div>$$ +/// 0 ∈ |x-y|_1\sign(x-y)_i + λ - g_i + δ_{≥ 0}(x_i) +/// \quad\text{for all}\quad i. +/// $$</div> +/// +/// If $x_i > 0$ and $x_i ≠ y_i$, then this forces $|g_i -λ| = |x-y|_1$. +/// Let $i^*$ be such an index. Then if $i ≠ i^*$ does not have the same value of $|g_i-λ|$, +/// we *must* have $x_i=y_i$ or $x_i=0$. In fact, even $x_i=y_i$ is impossible for $y_i>0$ if +/// $|g_i-λ| > |g_{i^*}-λ|$, whereas, otherwise $x_i=y_i$ is possible if $y_i>=0$. +/// If $|g_i-λ| > |g_{i^*}-λ|$, we must, therefore, either have $x_i=0$ or be able to take $i^*=i$. +/// +/// If $|g_i-λ| ≤ |g_{i^*}-λ|$, we can have either $x_i=0$ or $x_i=y_i$. If $y_i<=0$, then, +/// of course $x_i=0$. If $y_i>0$, we can take $x_i=y_i$. +/// +/// We, therefore, sort the |g_i-λ|, and look for an index $i^*$ that gives a non-contradictory +/// $β|x-y|_1=|g_i -λ|$, noting that contributions to β|x-y|_1= come from, besides $i^*$, from +/// indices $i$ such that $x_i= 0 ≠ y_i$. +/// +/// Finally, if no index $i$ reaches $|g_i -λ| = |x-y|_1$, we have to $|x-y|_1$ being in an +/// intermediate interval. This can similarly done using sorting. +#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] +pub(super) fn l1squared_nonneg_solve_exact<F, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, + λ_: F, + x: &mut Vector<F::MixedType, Dyn, S3>, +) -> bool +where + F: Float + ToNalgebraRealField, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, +{ + let λ = λ_.to_nalgebra_mixed(); + + assert_eq!(y.len(), g.len()); + assert_eq!(x.len(), g.len()); + assert!(x.len() <= u32::MAX as usize); + + #[derive(Debug)] + struct Tmp<F> { + y: F, + d: F, + contrib: F, + i: u32, // Shrink size for sort + λ_le_g: bool, + } + let mut sorted = Vec::from_iter(izip!(g.iter(), y.iter(), 0..).map(|(&g_i, &y_i, i)| { + // We already precompute the comparison λ <= g_i here, since we need it for the abs, + // and would need to store g_i for that in any case. + let (d, λ_le_g) = if λ <= g_i { + (g_i - λ, true) + } else { + (λ - g_i, false) + }; + Tmp { d, λ_le_g, y: y_i, i, contrib: 0.0 } + })); + sorted + .as_mut_slice() + .sort_unstable_by(|a, b| b.d.total_cmp(&a.d)); + + // Reverse-compute contribs + sorted.iter_mut().rev().fold(0.0, |contrib, a| { + a.contrib = contrib; + if a.y < 0.0 { + contrib - a.y + } else { + contrib + } + }); + + let mut contrib0 = 0.0; + + x.fill(0.0); + + let mut it = sorted.iter(); + let mut found = false; + + // We first try to find an index m such that y_m ≠ x_m > 0.0 that satisfies ‖x-y‖₁ = |λ - g_m|. + 'search: while let Some(m) = it.next() { + let d_m = m.d; + let y_m = m.y; + let contrib = m.contrib + contrib0; + let δ = d_m - contrib; + if δ > 0.0 { + // If λ < g[m], we must have x[m]≥y[m] so (contrib + x[m]-y[m])=g[m]-λ=d + // If λ = g_m, then also d_m=0, so contrib=0, and this gtives *x_m=y_m. + // If λ > g[m], We must have x[m]≤y[m] so -(contrib + y[m]-x[m])=g[m]-λ=-d. + let x_m_prime = if m.λ_le_g { y_m + δ } else { y_m - δ }; + if x_m_prime > 0.0 { + let x_m = unsafe { x.get_unchecked_mut(m.i as usize) }; + *x_m = x_m_prime; + found = true; + break 'search; + } + } + contrib0 += num_traits::abs(y_m); + } + + if !found { + // Every x_i is either zero or equal to y_i. + // We scan intervals (bound,prev_bound) between the d-values of each component, + // for the value of ‖x-y‖₁. Then it can be decided whether x_i should be one or equal to + // y_i. Due to the sorting and pre-calculation of posterior contributions, this simplifies + // into the following: + let mut prev_bound = F::MixedType::INFINITY; + let mut contrib0 = 0.0; + it = sorted.iter(); + 'search_degenerate: while let Some(m) = it.next() { + let bound = m.d; + let contrib = m.contrib + contrib0 - m.y.min(0.0); + if prev_bound >= contrib && contrib >= bound { + if m.y > 0.0 { + let x_m = unsafe { x.get_unchecked_mut(m.i as usize) }; + *x_m = m.y; + } + found = true; + break 'search_degenerate; + } + prev_bound = bound; + if m.y > 0.0 { + contrib0 += m.y; + } + } + // Check final interval from last element to -∞, if nothing found. + if !found && !(prev_bound >= contrib0) { + panic!("l1squared_nonneg_solve_exact failure") + } + } + + // Set remaining components to their now known values + while let Some(a) = it.next() { + // Safety: size checked above. + let y_i = a.y; + if y_i >= 0.0 { + // We can leave unchanged, as ‖x-y‖₁ is guaranteed large enough if + // the current maximising index attempt succeeds. + let x_i = unsafe { x.get_unchecked_mut(a.i as usize) }; + *x_i = y_i + } + // Zero fill in `else` case done above in init already + } + + debug_assert!(min_subdifferential(y, x, g, λ) <= F::EPSILON.to_nalgebra_mixed() * 10.0); + + return true; +} + /// This function applies an iterative method for the solution of the problem /// <div>$$ -/// \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x). +/// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x). /// $$</div> /// /// This function returns the number of iterations taken. #[replace_float_literals(F::cast_from(literal))] -pub fn l1squared_nonneg<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_nonneg<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ: F, - β: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, inner: &InnerSettings<F>, iterator: I, ) -> usize where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { + if let InnerMethod::Exact = inner.method { + // Try exact solution, fall back to PDPS if it does not work. + if l1squared_nonneg_solve_exact(y, g, λ, x) { + return 1; + } + } + match inner.method { - InnerMethod::PDPS => { + InnerMethod::PDPS | InnerMethod::Exact => { let inner_θ = 1.0; - // Estimate of ‖K‖ for K=θ\Id. + //Estimate of ‖K‖ for K=θ\Id. let normest = inner_θ; let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest); - l1squared_nonneg_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator) + l1squared_nonneg_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator) } InnerMethod::PP | InnerMethod::FB => { let inner_τ = inner.pp_τ.0; let inner_θ = inner.pp_τ.1; - l1squared_nonneg_pp(y, g, λ, β, x, inner_τ, inner_θ, iterator) + l1squared_nonneg_pp(y, g, λ, x, inner_τ, inner_θ, iterator) } other => unimplemented!("${other:?} is unimplemented"), }
--- a/src/subproblem/l1squared_unconstrained.rs Fri May 15 14:40:02 2026 -0500 +++ b/src/subproblem/l1squared_unconstrained.rs Sun Jul 19 07:34:39 2026 +0200 @@ -3,13 +3,12 @@ */ use itertools::izip; -use nalgebra::DVector; +use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1}; use numeric_literals::replace_float_literals; use std::cmp::Ordering::*; use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState}; -use alg_tools::nalgebra_support::ToNalgebraRealField; -use alg_tools::nanleast::NaNLeast; +use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField}; use alg_tools::norms::{Dist, L1}; use std::iter::zip; @@ -44,31 +43,47 @@ /// Clearly, if this condition fails for $x\_i$, it will fail for all the components /// already exluced. While, if it holds, it will hold for all components not excluded. #[replace_float_literals(F::cast_from(literal))] -pub(super) fn l1squared_prox<F: Float + nalgebra::RealField>( +pub(super) fn l1squared_prox<F: Float + nalgebra::RealField, S1, S2>( sorted_abs: &mut DVector<F>, - x: &mut DVector<F>, - y: &DVector<F>, + x: &mut Vector<F, Dyn, S1>, + y: &Vector<F, Dyn, S2>, β: F, -) { +) where + S2: Storage<F, Dyn>, + S1: StorageMut<F, Dyn>, +{ + //let orig_x = x.clone(); sorted_abs.copy_from(x); sorted_abs.axpy(-1.0, y, 1.0); sorted_abs.apply(|z_i| *z_i = num_traits::abs(*z_i)); - sorted_abs - .as_mut_slice() - .sort_unstable_by(|a, b| NaNLeast(*a).cmp(&NaNLeast(*b))); + sorted_abs.as_mut_slice().sort_unstable_by(F::total_cmp); let mut n = sorted_abs.sum(); for (m, az_i) in zip((1..=x.len() as u32).rev(), sorted_abs) { - // test first - let tmp = β * n / (1.0 + β * F::cast_from(m)); - if *az_i <= tmp { + // test first. This is just *az_i <= tmp, for tmp defined below, without the division. + if *az_i <= β * (n - F::cast_from(m) * *az_i) { // Fail n -= *az_i; } else { // Success + let tmp = β * n / (1.0 + β * F::cast_from(m)); x.zip_apply(y, |x_i, y_i| { *x_i = y_i + soft_thresholding(*x_i - y_i, tmp) }); + // //Check 0 ∈ w-x + β\norm{w-y}\_1\sign (w-y). + // let n: F = izip!(x.iter(), y) + // .map(|(&w_i, &y_i)| NumTraitsFloat::abs(w_i - y_i)) + // .sum(); + // for (&mut w_i, &x_i, &y_i) in izip!(x, &orig_x, y) { + // if w_i > y_i { + // assert_lt!(NumTraitsFloat::abs(w_i - x_i + n * β), 10.0 * F::EPSILON); + // } else if w_i < y_i { + // assert_lt!(NumTraitsFloat::abs(w_i - x_i - n * β), 10.0 * F::EPSILON); + // } else { + // assert_lt!(-n * β - 10.0 * F::EPSILON, w_i - x_i); + // assert_lt!(w_i - x_i, n * β + 10.0 * F::EPSILON); + // } + // } return; } } @@ -80,15 +95,20 @@ /// /// `v` will be modified and cannot be trusted to contain useful values afterwards. #[replace_float_literals(F::cast_from(literal))] -fn min_subdifferential<F: Float + nalgebra::RealField>( - y: &DVector<F>, - x: &DVector<F>, - g: &DVector<F>, +fn min_subdifferential<F: Float + nalgebra::RealField, S1, S2, S3>( + y: &Vector<F, Dyn, S1>, + x: &Vector<F, Dyn, S2>, + g: &Vector<F, Dyn, S3>, λ: F, - β: F, -) -> F { +) -> F +where + S1: Storage<F, Dyn>, + S2: Storage<F, Dyn>, + S3: Storage<F, Dyn>, + ShapeConstraint: StridesOk<F, Dyn, U1, S2>, +{ let mut val = 0.0; - let tmp = β * y.dist(x, L1); + let tmp = y.dist(x, L1); for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) { let (mut lb, mut ub) = (-g_i, -g_i); match x_i.partial_cmp(y_i) { @@ -132,12 +152,11 @@ /// The `λ` component of the model is handled in the proximal step instead of the gradient step /// for potential performance improvements. #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] -pub fn l1squared_unconstrained_pdps<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_unconstrained_pdps<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ_: F, - β_: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, τ_: F, σ_: F, iterator: I, @@ -145,21 +164,24 @@ where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { let λ = λ_.to_nalgebra_mixed(); - let β = β_.to_nalgebra_mixed(); let τ = τ_.to_nalgebra_mixed(); let σ = σ_.to_nalgebra_mixed(); let mut w = DVector::zeros(x.len()); let mut tmp = DVector::zeros(x.len()); - let mut xprev = x.clone(); + let mut xprev = x.clone_owned(); let mut iters = 0; iterator.iterate(|state| { // Primal step: x^{k+1} = prox_{τ|.-y|_1^2}(x^k - τ (w^k - g)) x.axpy(-τ, &w, 1.0); x.axpy(τ, g, 1.0); - l1squared_prox(&mut tmp, x, y, τ * β); + l1squared_prox(&mut tmp, x, y, τ); // Dual step: w^{k+1} = proj_{[-λ,λ]}(w^k + σ(2x^{k+1}-x^k)) w.axpy(2.0 * σ, x, 1.0); @@ -169,7 +191,7 @@ iters += 1; - state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ))) }); iters @@ -193,12 +215,11 @@ /// \end{split} /// $$</div> #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] -pub fn l1squared_unconstrained_pdps_alt<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_unconstrained_pdps_alt<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ_: F, - β_: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, τ_: F, σ_: F, θ_: F, @@ -207,17 +228,20 @@ where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { let λ = λ_.to_nalgebra_mixed(); let τ = τ_.to_nalgebra_mixed(); let σ = σ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); - let β = β_.to_nalgebra_mixed(); let σθ = σ * θ; let τθ = τ * θ; let mut w = DVector::zeros(x.len()); let mut tmp = DVector::zeros(x.len()); - let mut xprev = x.clone(); + let mut xprev = x.clone_owned(); let mut iters = 0; iterator.iterate(|state| { @@ -236,14 +260,14 @@ w.axpy(2.0, x, 1.0); w.axpy(-1.0, &xprev, 1.0); xprev.copy_from(&w); // use xprev as temporary variable - l1squared_prox(&mut tmp, &mut xprev, y, β / σθ); + l1squared_prox(&mut tmp, &mut xprev, y, 1.0 / σθ); w -= &xprev; w *= σ; xprev.copy_from(x); iters += 1; - state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ))) }); iters @@ -257,18 +281,21 @@ /// /// This function returns the number of iterations taken. #[replace_float_literals(F::cast_from(literal))] -pub fn l1squared_unconstrained<F, I>( - y: &DVector<F::MixedType>, - g: &DVector<F::MixedType>, +pub fn l1squared_unconstrained<F, I, S1, S2, S3>( + y: &Vector<F::MixedType, Dyn, S1>, + g: &Vector<F::MixedType, Dyn, S2>, λ: F, - β: F, - x: &mut DVector<F::MixedType>, + x: &mut Vector<F::MixedType, Dyn, S3>, inner: &InnerSettings<F>, iterator: I, ) -> usize where F: Float + ToNalgebraRealField, I: AlgIteratorFactory<F>, + S1: Storage<F::MixedType, Dyn>, + S2: Storage<F::MixedType, Dyn>, + S3: StorageMut<F::MixedType, Dyn>, + ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>, { // Estimate of ‖K‖ for K=θ Id. let inner_θ = 1.0; @@ -278,7 +305,7 @@ match inner.method { InnerMethod::PDPS => { - l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator) + l1squared_unconstrained_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator) } other => unimplemented!("${other:?} is unimplemented"), }
