pointsource_algs: comparison src/subproblem/nonneg.rs

-:aec67cdd6b14
+:efa60bc4f743
 use nalgebra::{DVector, DMatrix};
 use numeric_literals::replace_float_literals;
 use itertools::{izip, Itertools};
 use colored::Colorize;
+use std::cmp::Ordering::*;
 use alg_tools::iter::Mappable;
 use alg_tools::error::NumericalError;
 use alg_tools::iterate::{
 AlgIteratorFactory,
 use crate::types::*;
 use super::{
 InnerMethod,
 InnerSettings
 };
+use super::l1squared_nonneg::max_interval_dist_to_zero;
 /// Compute the proximal operator of $x \mapsto x + \delta\_{[0, \infty)}$, i.e.,
 /// the non-negativity contrained soft-thresholding operator.
 #[inline]
 #[replace_float_literals(F::cast_from(literal))]
-fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F {
+pub(super) fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F {
 (v - λ).max(0.0)
 }
-/// Returns the ∞-norm minimal subdifferential of $x ↦ x^⊤Ax - g^⊤ x + λ\vec 1^⊤ x δ_{≥ 0}(x)$
+/// Returns the ∞-norm minimal subdifferential of $x ↦ x^⊤Ax/2 - g^⊤ x + λ\vec 1^⊤ x + δ_{≥ 0}(x)$
 /// at $x$.
 ///
 /// `v` will be modified and cannot be trusted to contain useful values afterwards.
-#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
+#[replace_float_literals(F::cast_from(literal))]
-fn min_subdifferential<F : Float + ToNalgebraRealField>(
+fn min_subdifferential<F : Float + nalgebra::RealField>(
-v : &mut DVector<F::MixedType>,
+v : &mut DVector<F>,
-mA : &DMatrix<F::MixedType>,
+mA : &DMatrix<F>,
-x : &DVector<F::MixedType>,
+x : &DVector<F>,
-g : &DVector<F::MixedType>,
+g : &DVector<F>,
-λ : F::MixedType
+λ : F
 ) -> F {
 v.copy_from(g);
 mA.gemv(v, 1.0, x, -1.0);   // v =  Ax - g
 let mut val = 0.0;
 for (&v_i, &x_i) in izip!(v.iter(), x.iter()) {
-// The subdifferential of the objective is $Ax - g + λ + ∂ δ_{≥ 0}(x)$.
+let (mut lb, mut ub) = (v_i, v_i);
-let d = v_i + λ;
+match x_i.partial_cmp(&0.0) {
-if x_i > 0.0 || d < 0.0 {
+Some(Greater) => { lb += λ; ub += λ },
-val = val.max(d.abs());
+// Less should not happen
+Some(Less|Equal) => { lb = F::NEG_INFINITY; ub += λ },
+None => {},
 }
+val = max_interval_dist_to_zero(val, lb, ub);
 }
-F::from_nalgebra_mixed(val)
+val
 }
 /// Forward-backward splitting implementation of [`quadratic_nonneg`].
 /// For detailed documentation of the inputs and outputs, refer to there.
 ///
 });
 iters +=1;
 backup.map(|_| {
-min_subdifferential(&mut v, mA, x, g, λ)
+F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
 })
 });
 iters
 }
 iters += 1;
 // 4. Report solution quality
 state.if_verbose(|| {
 // Calculate subdifferential at the FB step `x` that hasn't yet had `s` yet added.
-min_subdifferential(&mut v, mA, x, g, λ)
+F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
 })
 });
 res.map(|_| iters)
 }
 /// This function applies an iterative method for the solution of the quadratic non-negativity
 /// constrained problem
 /// <div>$$
-///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + c + δ_{≥ 0}(x).
+///     \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + δ_{≥ 0}(x).
 /// $$</div>
 /// Semismooth Newton or forward-backward are supported based on the setting in `method`.
 /// The parameter `mA` is matrix $A$, and `g` and `λ` are as in the mathematical formulation.
 /// The constant $c$ does not need to be provided. The step length parameter is `τ` while
 /// `x` contains the initial iterate and on return the final one. The `iterator` controls
 ///  * Valkonen T. - _A method for weighted projections to the positive definite
 ///    cone_, <https://doi.org/10.1080/02331934.2014.929680>.
 ///
 /// This function returns the number of iterations taken.
 pub fn quadratic_nonneg<F, I>(
-method : InnerMethod,
 mA : &DMatrix<F::MixedType>,
 g : &DVector<F::MixedType>,
-//c_ : F,
 λ : F,
 x : &mut DVector<F::MixedType>,
-τ : F,
+mA_normest : F,
+inner : &InnerSettings<F>,
 iterator : I
 ) -> usize
 where F : Float + ToNalgebraRealField,
 I : AlgIteratorFactory<F>
 {
+let inner_τ = inner.fb_τ0 / mA_normest;
-match method {
-InnerMethod::FB =>
+match inner.method {
-quadratic_nonneg_fb(mA, g, λ, x, τ, iterator),
+InnerMethod::FB | InnerMethod::Default =>
+quadratic_nonneg_fb(mA, g, λ, x, inner_τ, iterator),
 InnerMethod::SSN =>
-quadratic_nonneg_ssn(mA, g, λ, x, τ, iterator).unwrap_or_else(|e| {
+quadratic_nonneg_ssn(mA, g, λ, x, inner_τ, iterator).unwrap_or_else(|e| {
 println!("{}", format!("{e}. Using FB fallback.").red());
 let ins = InnerSettings::<F>::default();
-quadratic_nonneg_fb(mA, g, λ, x, τ, ins.iterator_options)
+quadratic_nonneg_fb(mA, g, λ, x, inner_τ, ins.iterator_options)
-})
+}),
+InnerMethod::PDPS => unimplemented!(),
 }
 }

Mercurial > repos > pointsource_algs / file comparison

comparison: src/subproblem/nonneg.rs

src/subproblem/nonneg.rs