# HG changeset patch # User Tuomo Valkonen # Date 1669728972 -7200 # Node ID 7a953a87b6c16a7ac8b665c12db3e1882e6bd8f2 # Parent eb3c7813b67acf6b6b89c71acd0b24594b0fccd6 fubar diff -r eb3c7813b67a -r 7a953a87b6c1 src/forward_model.rs --- a/src/forward_model.rs Thu Dec 01 23:07:35 2022 +0200 +++ b/src/forward_model.rs Tue Nov 29 15:36:12 2022 +0200 @@ -53,14 +53,22 @@ type Observable : Euclidean + AXPY + Clone; + /// Finite-dimensional model for weight optimisation + type FindimModel : Linear, Codomain=Self::Observable>; - /// Return A_*A and A_* b + /// Return A_*A and A_* b instantiated at the locations of the spikes of `μ`. fn findim_quadratic_model( &self, μ : &DiscreteMeasure, b : &Self::Observable ) -> (DMatrix, DVector); + /// Return A instantiated at the locations of the spikes of `μ` + fn findim_model( + &self, + μ : &DiscreteMeasure, + ) -> Self::FindimModel; + /// Write an observable into a file. fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError; @@ -517,6 +525,7 @@ ShiftedSensor : LocalAnalysis, Weighted, F> : LocalAnalysis { type Observable = DVector; + type FindimModel = DMatrix; fn findim_quadratic_model( &self, @@ -524,6 +533,15 @@ b : &Self::Observable ) -> (DMatrix, DVector) { assert_eq!(b.len(), self.n_sensors()); + let mA = self.findim_model(μ); + let mAt = mA.transpose(); + (&mAt * mA, &mAt * b) + } + + fn findim_model( + &self, + μ : &DiscreteMeasure, F>, + ) -> DMatrix { let mut mA = DMatrix::zeros(self.n_sensors(), μ.len()); let grid = self.grid(); for (mut mAcol, δ) in mA.column_iter_mut().zip(μ.iter_spikes()) { @@ -532,8 +550,7 @@ mAcol[d] += sensor.apply(&δ.x); } } - let mAt = mA.transpose(); - (&mAt * mA, &mAt * b) + mA } fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError { diff -r eb3c7813b67a -r 7a953a87b6c1 src/frank_wolfe.rs --- a/src/frank_wolfe.rs Thu Dec 01 23:07:35 2022 +0200 +++ b/src/frank_wolfe.rs Tue Nov 29 15:36:12 2022 +0200 @@ -52,10 +52,10 @@ }; use crate::forward_model::ForwardModel; #[allow(unused_imports)] // Used in documentation -use crate::subproblem::{ - quadratic_nonneg, - InnerSettings, - InnerMethod, +use crate::subproblem::InnerSettings; +use crate::weight_optim::{ + prepare_optimise_weights, + optimise_weights_l2 }; use crate::tolerance::Tolerance; use crate::plot::{ @@ -104,75 +104,6 @@ } } -/// Helper struct for pre-initialising the finite-dimensional subproblems solver -/// [`prepare_optimise_weights`]. -/// -/// The pre-initialisation is done by [`prepare_optimise_weights`]. -pub struct FindimData { - opAnorm_squared : F -} - -/// Return a pre-initialisation struct for [`prepare_optimise_weights`]. -/// -/// The parameter `opA` is the forward operator $A$. -pub fn prepare_optimise_weights(opA : &A) -> FindimData -where F : Float + ToNalgebraRealField, - A : ForwardModel, F> { - FindimData{ - opAnorm_squared : opA.opnorm_bound().powi(2) - } -} - -/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity -/// point source localisation problem. -/// -/// That is, we minimise -///
$$ -/// μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ) -/// $$
-/// only with respect to the weights of $μ$. -/// -/// The parameter `μ` is the discrete measure whose weights are to be optimised. -/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the -/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while -/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to -/// save intermediate iteration states as images. The parameter `findim_data` should be -/// prepared using [`prepare_optimise_weights`]: -/// -/// Returns the number of iterations taken by the method configured in `inner`. -pub fn optimise_weights<'a, F, A, I, const N : usize>( - μ : &mut DiscreteMeasure, F>, - opA : &'a A, - b : &A::Observable, - α : F, - findim_data : &FindimData, - inner : &InnerSettings, - iterator : I -) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory, - A : ForwardModel, F> -{ - // Form and solve finite-dimensional subproblem. - let (Ã, g̃) = opA.findim_quadratic_model(&μ, b); - let mut x = μ.masses_dvector(); - - // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to - // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the - // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to - // 2-norm, we estimate - // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2 - // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2}, - // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no - // square root is needed when we scale: - let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len())); - let iters = quadratic_nonneg(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator); - // Update masses of μ based on solution of finite-dimensional subproblem. - μ.set_masses_dvector(&x); - - iters -} - /// Solve point source localisation problem using a conditional gradient method /// for the 2-norm-squared data fidelity, i.e., the problem ///
$$ @@ -279,7 +210,7 @@ } }; - inner_iters += optimise_weights(&mut μ, opA, b, α, &findim_data, &config.inner, inner_it); + inner_iters += optimise_weights_l2(&mut μ, opA, b, α, &findim_data, &config.inner, inner_it); // Merge spikes and update residual for next step and `if_verbose` below. let n_before_merge = μ.len(); diff -r eb3c7813b67a -r 7a953a87b6c1 src/main.rs --- a/src/main.rs Thu Dec 01 23:07:35 2022 +0200 +++ b/src/main.rs Tue Nov 29 15:36:12 2022 +0200 @@ -30,6 +30,7 @@ pub mod forward_model; pub mod plot; pub mod subproblem; +pub mod weight_optim; pub mod tolerance; pub mod fb; pub mod frank_wolfe; diff -r eb3c7813b67a -r 7a953a87b6c1 src/pdps.rs --- a/src/pdps.rs Thu Dec 01 23:07:35 2022 +0200 +++ b/src/pdps.rs Tue Nov 29 15:36:12 2022 +0200 @@ -352,4 +352,4 @@ opA, α, op𝒟, τ, &config.insertion, iterator, plotter, y, pdps ) -} \ No newline at end of file +} diff -r eb3c7813b67a -r 7a953a87b6c1 src/run.rs --- a/src/run.rs Thu Dec 01 23:07:35 2022 +0200 +++ b/src/run.rs Tue Nov 29 15:36:12 2022 +0200 @@ -57,8 +57,10 @@ FWConfig, FWVariant, pointsource_fw, +}; +use crate::weight_optim::{ prepare_optimise_weights, - optimise_weights, + optimise_weights_l2, }; use crate::subproblem::InnerSettings; use crate::seminorms::*; @@ -451,7 +453,7 @@ Some(mut μ) => { match dataterm { DataTerm::L2Squared => { - optimise_weights( + optimise_weights_l2( &mut μ, &opA, &b, α, &findim_data, &inner_config, inner_it ); diff -r eb3c7813b67a -r 7a953a87b6c1 src/subproblem.rs --- a/src/subproblem.rs Thu Dec 01 23:07:35 2022 +0200 +++ b/src/subproblem.rs Tue Nov 29 15:36:12 2022 +0200 @@ -15,8 +15,10 @@ Verbose, Step, }; -use alg_tools::linops::GEMV; +use alg_tools::linops::{GEMV, Adjointable, AXPY}; use alg_tools::nalgebra_support::ToNalgebraRealField; +use alg_tools::norms::{Projection, Linfinity}; +use alg_tools::euclidean::Euclidean; use crate::types::*; @@ -100,6 +102,8 @@ let mut v = DVector::zeros(x.len()); let mut iters = 0; + assert!(λ > 0.0); + iterator.iterate(|state| { // Replace `x` with $x - τ[Ax-g]= [x + τg]- τAx$ v.copy_from(g); // v = g @@ -219,6 +223,8 @@ let mut decomp = nalgebra::linalg::LU::new(DMatrix::zeros(0, 0)); let mut iters = 0; + assert!(λ > 0.0); + let res = iterator.iterate_fallible(|state| { // 1. Perform delayed SSN-update based on previously computed step on active // coordinates. The step is delayed to the beginning of the loop because @@ -371,3 +377,107 @@ }) } } + +/// PDPSimplementation of [`l1_nonneg`]. +/// For detailed documentation of the inputs and outputs, refer to there. +#[replace_float_literals(F::cast_from(literal))] +pub fn l1_nonneg_pdps( + mA : &A, + b : &Y, + λ : F, + x : &mut DVector, + τ : F, + σ : F, + iterator : I +) -> usize +where F : Float + ToNalgebraRealField, + I : AlgIteratorFactory, + Y : Euclidean + Projection + AXPY, + A : GEMV, Codomain=Y> + + Adjointable, Y>, + for<'a> A::Adjoint<'a> : GEMV> +{ + let mAt = mA.adjoint(); + let mut xprev = x.clone(); + let τλ = τ * λ; + let mut v = DVector::zeros(x.len()); + let mut y = b.similar_origin(); + let mut iters = 0; + + assert!(λ > 0.0); + + iterator.iterate(|state| { + // Replace `x` with $x - τA^*y + xprev.copy_from(x); + mAt.gemv(x, -τ, &y, 1.0); + // Calculate the proximal map + x.iter_mut().for_each(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τλ)); + // Calculate v = 2x - xprev + v.copy_from(x); + v.axpy(-1.0, &xprev, 2.0); + // Calculate y = y + σ(A v - b) + y.axpy(-σ, b, 1.0); + mA.gemv(&mut y, σ, &v, 1.0); + y.proj_ball_mut(1.0, Linfinity); + + iters +=1; + + state.if_verbose(|| { + // The subdifferential of the objective is $A^* ∂\|.\|(Ax - g) + λ + ∂δ_{≥ 0}(x)$. + // We return the minimal ∞-norm over all subderivatives. + v.copy_from(b); // d = g + mA.gemv(&mut ỹ, 1.0, x, -1.0); // d = Ax - b + let mut val = 0.0; + izip!(ỹ.iter(), x.iter()).for_each(|(&ỹ_i, &x_i)| { + use std::cmp::Ordering::*; + let tmp = match ỹ_i.partial_cmp(&0.0) { + None => F::INFINITY, + Some(c) => match (c, x_i > 0.0) { + (Greater, true) => 1.0 + λ, + (Greater, false) | (Equal, false) => 0.0, + (Less, true) => -1.0 + λ, + (Less, false) => 0.0.min(-1.0 + λ), + (Equal, true) => 0.0.max(-1.0 + λ), // since λ > 0 + } + }; + val = val.max(tmp); + }); + val + }) + }); + + iters +} + +/// Method for L1 weight optimisation +#[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)] +#[allow(dead_code)] +pub enum L1InnerMethod { + /// PDPS + PDPS, +} + +/// This function applies an iterative method for the solution of the L1 non-negativity +/// constrained problem +///
$$ +/// \min_{x ∈ ℝ^n} \|Ax - b\|_1 + λ{\vec 1}^⊤ x + δ_{≥ 0}(x). +/// $$
+/// +/// This function returns the number of iterations taken. +pub fn l1_nonneg( + method : L1InnerMethod, + mA : &DMatrix, + b : &DVector, + λ : F, + σ : F, + x : &mut DVector, + τ : F, + iterator : I +) -> usize +where F : Float + ToNalgebraRealField, + I : AlgIteratorFactory { + + match method { + L1InnerMethod::PDPS => l1_nonneg_pdps(mA, b, λ, x, τ, σ, iterator) + } +} diff -r eb3c7813b67a -r 7a953a87b6c1 src/weight_optim.rs --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/weight_optim.rs Tue Nov 29 15:36:12 2022 +0200 @@ -0,0 +1,131 @@ +//! Weight optimisation routines + +use numeric_literals::replace_float_literals; +use alg_tools::nalgebra_support::ToNalgebraRealField; +use alg_tools::loc::Loc; +use alg_tools::iterate::AlgIteratorFactory; +use crate::types::*; +use crate::subproblem::*; +use crate::measures::*; +use crate::forward_model::ForwardModel; + +/// Helper struct for pre-initialising the finite-dimensional subproblems solver +/// [`prepare_optimise_weights_l2`]. +/// +/// The pre-initialisation is done by [`prepare_optimise_weights_l2`]. +pub struct FindimData { + opAnorm_squared : F +} + +/// Return a pre-initialisation struct for [`prepare_optimise_weights_l2`]. +/// +/// The parameter `opA` is the forward operator $A$. +pub fn prepare_optimise_weights(opA : &A) -> FindimData +where F : Float + ToNalgebraRealField, + A : ForwardModel, F> { + FindimData{ + opAnorm_squared : opA.opnorm_bound().powi(2) + } +} + +/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity +/// point source localisation problem. +/// +/// That is, we minimise +///
$$ +/// μ ↦ \frac{1}{2}\|Aμ-b\|_2^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ) +/// $$
+/// only with respect to the weights of $μ$. +/// +/// The parameter `μ` is the discrete measure whose weights are to be optimised. +/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the +/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while +/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to +/// save intermediate iteration states as images. The parameter `findim_data` should be +/// prepared using [`prepare_optimise_weights`]: +/// +/// Returns the number of iterations taken by the method configured in `inner`. +pub fn optimise_weights_l2<'a, F, A, I, const N : usize>( + μ : &mut DiscreteMeasure, F>, + opA : &'a A, + b : &A::Observable, + α : F, + findim_data : &FindimData, + inner : &InnerSettings, + iterator : I +) -> usize +where F : Float + ToNalgebraRealField, + I : AlgIteratorFactory, + A : ForwardModel, F> +{ + // Form and solve finite-dimensional subproblem. + let (Ã, g̃) = opA.findim_quadratic_model(&μ, b); + let mut x = μ.masses_dvector(); + + // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to + // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the + // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to + // 2-norm, we estimate + // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2 + // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2}, + // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no + // square root is needed when we scale: + let inner_τ = inner.τ0 / (findim_data.opAnorm_squared * F::cast_from(μ.len())); + let iters = quadratic_nonneg(inner.method, &Ã, &g̃, α, &mut x, inner_τ, iterator); + // Update masses of μ based on solution of finite-dimensional subproblem. + μ.set_masses_dvector(&x); + + iters +} + +/// Solve the finite-dimensional weight optimisation problem for the 2-norm-squared data fidelity +/// point source localisation problem. +/// +/// That is, we minimise +///
$$ +/// μ ↦ \|Aμ-b\|_1 + α\|μ\|_ℳ + δ_{≥ 0}(μ) +/// $$
+/// only with respect to the weights of $μ$. +/// +/// The parameter `μ` is the discrete measure whose weights are to be optimised. +/// The `opA` parameter is the forward operator $A$, while `b`$ and `α` are as in the +/// objective above. The method parameter are set in `inner` (see [`InnerSettings`]), while +/// `iterator` is used to iterate the steps of the method, and `plotter` may be used to +/// save intermediate iteration states as images. +/// +/// Returns the number of iterations taken by the method configured in `inner`. +#[replace_float_literals(F::cast_from(literal))] +pub fn optimise_weights_l1<'a, F, A, I, const N : usize>( + μ : &mut DiscreteMeasure, F>, + opA : &'a A, + b : &A::Observable, + α : F, + findim_data : &FindimData, + inner : &InnerSettings, + iterator : I +) -> usize +where F : Float + ToNalgebraRealField, + I : AlgIteratorFactory, + A : ForwardModel, F> +{ + // Form and solve finite-dimensional subproblem. + let à = opA.findim_model(&μ); + let mut x = μ.masses_dvector(); + + // `inner_τ1` is based on an estimate of the operator norm of $A$ from ℳ(Ω) to + // ℝ^n. This estimate is a good one for the matrix norm from ℝ^m to ℝ^n when the + // former is equipped with the 1-norm. We need the 2-norm. To pass from 1-norm to + // 2-norm, we estimate + // ‖A‖_{2,2} := sup_{‖x‖_2 ≤ 1} ‖Ax‖_2 ≤ sup_{‖x‖_1 ≤ C} ‖Ax‖_2 + // = C sup_{‖x‖_1 ≤ 1} ‖Ax‖_2 = C ‖A‖_{1,2}, + // where C = √m satisfies ‖x‖_1 ≤ C ‖x‖_2. Since we are intested in ‖A_*A‖, no + // square root is needed when we scale: + let l = (findim_data.opAnorm_squared * F::cast_from(μ.len())).sqrt(); + let inner_σ = (0.99 / inner.τ0) / l; + let inner_τ = inner.τ0 / l; + let iters = l1_nonneg(L1InnerMethod::PDPS, &Ã, b, α, &mut x, inner_τ, inner_σ, iterator); + // Update masses of μ based on solution of finite-dimensional subproblem. + μ.set_masses_dvector(&x); + + iters +}