Mercurial > repos > pointsource_algs / changeset

--- 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<F, Output=Self::Observable>
                       + AXPY<F>
                       + Clone;
+    /// Finite-dimensional model for weight optimisation
+    type FindimModel : Linear<DVector<F::MixedType>, 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<Domain, F>,
         b : &Self::Observable
     ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>);
 
+    /// Return A instantiated at the locations of the spikes of `μ`
+    fn findim_model(
+        &self,
+        μ : &DiscreteMeasure<Domain, F>,
+    ) -> Self::FindimModel;
+
     /// Write an observable into a file.
     fn write_observable(&self, b : &Self::Observable, prefix : String) -> DynError;
 
@@ -517,6 +525,7 @@
       ShiftedSensor<F, S, P, N> : LocalAnalysis<F, BT::Agg, N>,
       Weighted<ShiftedSensor<F, S, P, N>, F> : LocalAnalysis<F, BT::Agg, N> {
     type Observable = DVector<F>;
+    type FindimModel = DMatrix<F>;
 
     fn findim_quadratic_model(
         &self,
@@ -524,6 +533,15 @@
         b : &Self::Observable
     ) -> (DMatrix<F::MixedType>, DVector<F::MixedType>) {
         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<Loc<F, N>, F>,
+    ) -> DMatrix<F::MixedType> {
         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 {

--- 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<F : Float> {
-    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<F, A, const N : usize>(opA : &A) -> FindimData<F>
-where F : Float + ToNalgebraRealField,
-      A : ForwardModel<Loc<F, N>, 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
-/// <div>$$
-///     μ ↦ \frac{1}{2}\|Aμ-b\|_w^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
-/// $$</div>
-/// 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<Loc<F, N>, F>,
-    opA : &'a A,
-    b : &A::Observable,
-    α : F,
-    findim_data : &FindimData<F>,
-    inner : &InnerSettings<F>,
-    iterator : I
-) -> usize
-where F : Float + ToNalgebraRealField,
-      I : AlgIteratorFactory<F>,
-      A : ForwardModel<Loc<F, N>, 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
 /// <div>$$
@@ -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();

--- 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;

--- 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
+}

--- 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
                                 );

--- 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<F, I, A, Y>(
+    mA : &A,
+    b : &Y,
+    λ : F,
+    x : &mut DVector<F>,
+    τ : F,
+    σ : F,
+    iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+      I : AlgIteratorFactory<F>,
+      Y : Euclidean<F, Output = Y> + Projection<F, Linfinity> + AXPY<F>,
+      A : GEMV<F, DVector<F>, Codomain=Y>
+          + Adjointable<DVector<F>, Y>,
+      for<'a> A::Adjoint<'a> : GEMV<F, Y, Codomain=DVector<F>>
+{
+    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
+/// <div>$$
+///     \min_{x ∈ ℝ^n} \|Ax - b\|_1 + λ{\vec 1}^⊤ x + δ_{≥ 0}(x).
+/// $$</div>
+///
+/// This function returns the number of iterations taken.
+pub fn l1_nonneg<F, I>(
+    method : L1InnerMethod,
+    mA : &DMatrix<F::MixedType>,
+    b : &DVector<F::MixedType>,
+    λ : F,
+    σ : F,
+    x : &mut DVector<F::MixedType>,
+    τ : F,
+    iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+      I : AlgIteratorFactory<F> {
+
+    match method {
+        L1InnerMethod::PDPS => l1_nonneg_pdps(mA, b, λ, x, τ, σ, iterator)
+    }
+}

--- /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<F : Float> {
+    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<F, A, const N : usize>(opA : &A) -> FindimData<F>
+where F : Float + ToNalgebraRealField,
+      A : ForwardModel<Loc<F, N>, 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
+/// <div>$$
+///     μ ↦ \frac{1}{2}\|Aμ-b\|_2^2 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
+/// $$</div>
+/// 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<Loc<F, N>, F>,
+    opA : &'a A,
+    b : &A::Observable,
+    α : F,
+    findim_data : &FindimData<F>,
+    inner : &InnerSettings<F>,
+    iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+      I : AlgIteratorFactory<F>,
+      A : ForwardModel<Loc<F, N>, 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
+/// <div>$$
+///     μ ↦ \|Aμ-b\|_1 + α\|μ\|_ℳ + δ_{≥ 0}(μ)
+/// $$</div>
+/// 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<Loc<F, N>, F>,
+    opA : &'a A,
+    b : &A::Observable,
+    α : F,
+    findim_data : &FindimData<F>,
+    inner : &InnerSettings<F>,
+    iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+      I : AlgIteratorFactory<F>,
+      A : ForwardModel<Loc<F, N>, 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
+}