pointsource_algs: comparison src/subproblem.rs

-:eb3c7813b67a
+:7a953a87b6c1
 AlgIteratorState,
 AlgIteratorOptions,
 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::*;
 /// Method for solving finite-dimensional subproblems
 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
 let λ = λ_.to_nalgebra_mixed();
 let τ = τ_.to_nalgebra_mixed();
 let τλ = τ * λ;
 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
 v.axpy(1.0, x, τ);          // v = x + τ*g
 let τλ = τ * λ;
 let mut inact : Vec<bool> = Vec::from_iter(std::iter::repeat(false).take(n));
 let mut s = DVector::zeros(0);
 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
 // the SSN step may violate constraints, so we arrange `x` to contain at the
 let ins = InnerSettings::<F>::default();
 quadratic_nonneg_fb(mA, g, λ, x, τ, ins.iterator_options)
 })
 }
 }
+/// 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 ỹ, 1.0, x, -1.0);   // d =  Ax - b
+let mut val = 0.0;
+izip!(ỹ.iter(), x.iter()).for_each(|(&ỹ_i, &x_i)| {
+use std::cmp::Ordering::*;
+let tmp = match ỹ_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)
+}
+}

Mercurial > repos > pointsource_algs / file comparison

comparison: src/subproblem.rs

src/subproblem.rs