--- a/src/subproblem/l1squared_unconstrained.rs Fri May 15 14:40:02 2026 -0500
+++ b/src/subproblem/l1squared_unconstrained.rs Sun Jul 19 07:34:39 2026 +0200
@@ -3,13 +3,12 @@
*/
use itertools::izip;
-use nalgebra::DVector;
+use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1};
use numeric_literals::replace_float_literals;
use std::cmp::Ordering::*;
use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
-use alg_tools::nanleast::NaNLeast;
+use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField};
use alg_tools::norms::{Dist, L1};
use std::iter::zip;
@@ -44,31 +43,47 @@
/// Clearly, if this condition fails for $x\_i$, it will fail for all the components
/// already exluced. While, if it holds, it will hold for all components not excluded.
#[replace_float_literals(F::cast_from(literal))]
-pub(super) fn l1squared_prox<F: Float + nalgebra::RealField>(
+pub(super) fn l1squared_prox<F: Float + nalgebra::RealField, S1, S2>(
sorted_abs: &mut DVector<F>,
- x: &mut DVector<F>,
- y: &DVector<F>,
+ x: &mut Vector<F, Dyn, S1>,
+ y: &Vector<F, Dyn, S2>,
β: F,
-) {
+) where
+ S2: Storage<F, Dyn>,
+ S1: StorageMut<F, Dyn>,
+{
+ //let orig_x = x.clone();
sorted_abs.copy_from(x);
sorted_abs.axpy(-1.0, y, 1.0);
sorted_abs.apply(|z_i| *z_i = num_traits::abs(*z_i));
- sorted_abs
- .as_mut_slice()
- .sort_unstable_by(|a, b| NaNLeast(*a).cmp(&NaNLeast(*b)));
+ sorted_abs.as_mut_slice().sort_unstable_by(F::total_cmp);
let mut n = sorted_abs.sum();
for (m, az_i) in zip((1..=x.len() as u32).rev(), sorted_abs) {
- // test first
- let tmp = β * n / (1.0 + β * F::cast_from(m));
- if *az_i <= tmp {
+ // test first. This is just *az_i <= tmp, for tmp defined below, without the division.
+ if *az_i <= β * (n - F::cast_from(m) * *az_i) {
// Fail
n -= *az_i;
} else {
// Success
+ let tmp = β * n / (1.0 + β * F::cast_from(m));
x.zip_apply(y, |x_i, y_i| {
*x_i = y_i + soft_thresholding(*x_i - y_i, tmp)
});
+ // //Check 0 ∈ w-x + β\norm{w-y}\_1\sign (w-y).
+ // let n: F = izip!(x.iter(), y)
+ // .map(|(&w_i, &y_i)| NumTraitsFloat::abs(w_i - y_i))
+ // .sum();
+ // for (&mut w_i, &x_i, &y_i) in izip!(x, &orig_x, y) {
+ // if w_i > y_i {
+ // assert_lt!(NumTraitsFloat::abs(w_i - x_i + n * β), 10.0 * F::EPSILON);
+ // } else if w_i < y_i {
+ // assert_lt!(NumTraitsFloat::abs(w_i - x_i - n * β), 10.0 * F::EPSILON);
+ // } else {
+ // assert_lt!(-n * β - 10.0 * F::EPSILON, w_i - x_i);
+ // assert_lt!(w_i - x_i, n * β + 10.0 * F::EPSILON);
+ // }
+ // }
return;
}
}
@@ -80,15 +95,20 @@
///
/// `v` will be modified and cannot be trusted to contain useful values afterwards.
#[replace_float_literals(F::cast_from(literal))]
-fn min_subdifferential<F: Float + nalgebra::RealField>(
- y: &DVector<F>,
- x: &DVector<F>,
- g: &DVector<F>,
+fn min_subdifferential<F: Float + nalgebra::RealField, S1, S2, S3>(
+ y: &Vector<F, Dyn, S1>,
+ x: &Vector<F, Dyn, S2>,
+ g: &Vector<F, Dyn, S3>,
λ: F,
- β: F,
-) -> F {
+) -> F
+where
+ S1: Storage<F, Dyn>,
+ S2: Storage<F, Dyn>,
+ S3: Storage<F, Dyn>,
+ ShapeConstraint: StridesOk<F, Dyn, U1, S2>,
+{
let mut val = 0.0;
- let tmp = β * y.dist(x, L1);
+ let tmp = y.dist(x, L1);
for (&g_i, &x_i, y_i) in izip!(g.iter(), x.iter(), y.iter()) {
let (mut lb, mut ub) = (-g_i, -g_i);
match x_i.partial_cmp(y_i) {
@@ -132,12 +152,11 @@
/// The `λ` component of the model is handled in the proximal step instead of the gradient step
/// for potential performance improvements.
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn l1squared_unconstrained_pdps<F, I>(
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
+pub fn l1squared_unconstrained_pdps<F, I, S1, S2, S3>(
+ y: &Vector<F::MixedType, Dyn, S1>,
+ g: &Vector<F::MixedType, Dyn, S2>,
λ_: F,
- β_: F,
- x: &mut DVector<F::MixedType>,
+ x: &mut Vector<F::MixedType, Dyn, S3>,
τ_: F,
σ_: F,
iterator: I,
@@ -145,21 +164,24 @@
where
F: Float + ToNalgebraRealField,
I: AlgIteratorFactory<F>,
+ S1: Storage<F::MixedType, Dyn>,
+ S2: Storage<F::MixedType, Dyn>,
+ S3: StorageMut<F::MixedType, Dyn>,
+ ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
let λ = λ_.to_nalgebra_mixed();
- let β = β_.to_nalgebra_mixed();
let τ = τ_.to_nalgebra_mixed();
let σ = σ_.to_nalgebra_mixed();
let mut w = DVector::zeros(x.len());
let mut tmp = DVector::zeros(x.len());
- let mut xprev = x.clone();
+ let mut xprev = x.clone_owned();
let mut iters = 0;
iterator.iterate(|state| {
// Primal step: x^{k+1} = prox_{τ|.-y|_1^2}(x^k - τ (w^k - g))
x.axpy(-τ, &w, 1.0);
x.axpy(τ, g, 1.0);
- l1squared_prox(&mut tmp, x, y, τ * β);
+ l1squared_prox(&mut tmp, x, y, τ);
// Dual step: w^{k+1} = proj_{[-λ,λ]}(w^k + σ(2x^{k+1}-x^k))
w.axpy(2.0 * σ, x, 1.0);
@@ -169,7 +191,7 @@
iters += 1;
- state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)))
+ state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
});
iters
@@ -193,12 +215,11 @@
/// \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>,
+pub fn l1squared_unconstrained_pdps_alt<F, I, S1, S2, S3>(
+ y: &Vector<F::MixedType, Dyn, S1>,
+ g: &Vector<F::MixedType, Dyn, S2>,
λ_: F,
- β_: F,
- x: &mut DVector<F::MixedType>,
+ x: &mut Vector<F::MixedType, Dyn, S3>,
τ_: F,
σ_: F,
θ_: F,
@@ -207,17 +228,20 @@
where
F: Float + ToNalgebraRealField,
I: AlgIteratorFactory<F>,
+ S1: Storage<F::MixedType, Dyn>,
+ S2: Storage<F::MixedType, Dyn>,
+ S3: StorageMut<F::MixedType, Dyn>,
+ ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
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 xprev = x.clone_owned();
let mut iters = 0;
iterator.iterate(|state| {
@@ -236,14 +260,14 @@
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, β / σθ);
+ l1squared_prox(&mut tmp, &mut xprev, y, 1.0 / σθ);
w -= &xprev;
w *= σ;
xprev.copy_from(x);
iters += 1;
- state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)))
+ state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ)))
});
iters
@@ -257,18 +281,21 @@
///
/// This function returns the number of iterations taken.
#[replace_float_literals(F::cast_from(literal))]
-pub fn l1squared_unconstrained<F, I>(
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
+pub fn l1squared_unconstrained<F, I, S1, S2, S3>(
+ y: &Vector<F::MixedType, Dyn, S1>,
+ g: &Vector<F::MixedType, Dyn, S2>,
λ: F,
- β: F,
- x: &mut DVector<F::MixedType>,
+ x: &mut Vector<F::MixedType, Dyn, S3>,
inner: &InnerSettings<F>,
iterator: I,
) -> usize
where
F: Float + ToNalgebraRealField,
I: AlgIteratorFactory<F>,
+ S1: Storage<F::MixedType, Dyn>,
+ S2: Storage<F::MixedType, Dyn>,
+ S3: StorageMut<F::MixedType, Dyn>,
+ ShapeConstraint: StridesOk<F::MixedType, Dyn, U1, S3>,
{
// Estimate of ‖K‖ for K=θ Id.
let inner_θ = 1.0;
@@ -278,7 +305,7 @@
match inner.method {
InnerMethod::PDPS => {
- l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator)
+ l1squared_unconstrained_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator)
}
other => unimplemented!("${other:?} is unimplemented"),
}