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