pointsource_algs: comparison src/subproblem/l1squared

-:b6bdb6cb4d44
+:6a7365d73e4c
 });
 iters
 }
+/// Alternative PDPS implementation of [`l1squared_unconstrained`].
+/// For detailed documentation of the inputs and outputs, refer to there.
+///
+/// By not dualising the 1-norm, this should produce more sparse solutions than
+/// [`l1squared_unconstrained_pdps`].
+///
+/// The `λ` component of the model is handled in the proximal step instead of the gradient step
+/// for potential performance improvements.
+/// The parameter `θ` is used to multiply the rescale the operator (identity) of the PDPS model.
+/// We rewrite
+/// <div>$$
+///     \begin{split}
+///     & \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ \\
+///     & = \min_{x ∈ ℝ^n} \max_{w} ⟨θ w, x⟩ - g^⊤ x + λ\|x\|₁
+///      - \left(x ↦ \frac{β}{2θ} |x-y|_1^2 \right)^*(w).
+///     \end{split}
+/// $$</div>
+#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
+pub fn l1squared_unconstrained_pdps_alt<F, I>(
+y : &DVector<F::MixedType>,
+g : &DVector<F::MixedType>,
+λ_ : F,
+β_ : F,
+x : &mut DVector<F::MixedType>,
+τ_ : F,
+σ_ : F,
+θ_ : F,
+iterator : I
+) -> usize
+where F : Float + ToNalgebraRealField,
+I : AlgIteratorFactory<F>
+{
+let λ = λ_.to_nalgebra_mixed();
+let τ = τ_.to_nalgebra_mixed();
+let σ = σ_.to_nalgebra_mixed();
+let θ = θ_.to_nalgebra_mixed();
+let β = β_.to_nalgebra_mixed();
+let σθ = σ*θ;
+let τθ = τ*θ;
+let mut w = DVector::zeros(x.len());
+let mut tmp = DVector::zeros(x.len());
+let mut xprev = x.clone();
+let mut iters = 0;
+iterator.iterate(|state| {
+// Primal step: x^{k+1} = soft_τλ(x^k - τ(θ w^k -g))
+x.axpy(-τθ, &w, 1.0);
+x.axpy(τ, g, 1.0);
+x.apply(|x_i| *x_i = soft_thresholding(*x_i, τ*λ));
+// Dual step: with g(x) = (β/(2θ))‖x-y‖₁² and q = w^k + σ(2x^{k+1}-x^k),
+// we compute w^{k+1} = prox_{σg^*}(q) for
+//                    = q - σ prox_{g/σ}(q/σ)
+//                    = q - σ prox_{(β/(2θσ))‖.-y‖₁²}(q/σ)
+//                    = σ(q/σ - prox_{(β/(2θσ))‖.-y‖₁²}(q/σ))
+// where q/σ = w^k/σ + (2x^{k+1}-x^k),
+w /= σ;
+w.axpy(2.0, x, 1.0);
+w.axpy(-1.0, &xprev, 1.0);
+xprev.copy_from(&w); // use xprev as temporary variable
+l1squared_prox(&mut tmp, &mut xprev, y, β/σθ);
+w -= &xprev;
+w *= σ;
+xprev.copy_from(x);
+iters += 1;
+state.if_verbose(|| {
+F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))
+})
+});
+iters
+}
 /// This function applies an iterative method for the solution of the problem
 /// <div>$$
 ///     \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁.
 /// $$</div>
 iterator : I
 ) -> usize
 where F : Float + ToNalgebraRealField,
 I : AlgIteratorFactory<F>
 {
-// Estimate of ‖K‖ for K=Id.
+// Estimate of ‖K‖ for K=θ Id.
-let normest = 1.0;
+let inner_θ = 1.0;
+let normest = inner_θ;
 let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest);
 match inner.method {
 InnerMethod::PDPS =>
-l1squared_unconstrained_pdps(y, g, λ, β, x, inner_τ, inner_σ, iterator),
+l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator),
 other => unimplemented!("${other:?} is unimplemented"),
 }
 }

comparison: src/subproblem/l1squared_unconstrained.rs

src/subproblem/l1squared_unconstrained.rs

Mercurial > repos > pointsource_algs / file comparison

comparison: src/subproblem/l1squared_unconstrained.rs

src/subproblem/l1squared_unconstrained.rs