--- a/src/subproblem/l1squared_nonneg.rs Fri May 15 14:40:02 2026 -0500
+++ b/src/subproblem/l1squared_nonneg.rs Sun Jul 19 07:34:39 2026 +0200
@@ -3,13 +3,11 @@
*/
use itertools::izip;
-use nalgebra::DVector;
+use nalgebra::{constraint::ShapeConstraint, DVector, Dyn, Storage, StorageMut, Vector, U1};
use numeric_literals::replace_float_literals;
-//use std::iter::zip;
-use std::cmp::Ordering::*;
use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState};
-use alg_tools::nalgebra_support::ToNalgebraRealField;
+use alg_tools::nalgebra_support::{StridesOk, ToNalgebraRealField};
use alg_tools::norms::{Dist, L1};
use super::l1squared_unconstrained::l1squared_prox;
@@ -19,6 +17,7 @@
/// Return maximum of `dist` and distnce of inteval `[lb, ub]` to zero.
#[replace_float_literals(F::cast_from(literal))]
+#[inline]
pub(super) fn max_interval_dist_to_zero<F: Float>(dist: F, lb: F, ub: F) -> F {
if lb < 0.0 {
if ub > 0.0 {
@@ -27,7 +26,7 @@
dist.max(-ub)
}
} else
- /* ub ≥ 0.0*/
+ /* lb ≥ 0.0*/
{
dist.max(lb)
}
@@ -37,44 +36,35 @@
///
/// `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) {
- Some(Greater) => {
- lb += tmp;
- ub += tmp
- }
- Some(Less) => {
- lb -= tmp;
- ub -= tmp
- }
- Some(Equal) => {
- lb -= tmp;
- ub += tmp
- }
- None => {}
+ let (mut lb, mut ub) = (-g_i + λ, -g_i + λ);
+ if num_traits::abs(x_i - *y_i) < F::EPSILON {
+ lb -= tmp;
+ ub += tmp
+ } else if x_i > *y_i {
+ lb += tmp;
+ ub += tmp
+ } else {
+ lb -= tmp;
+ ub -= tmp
}
- match x_i.partial_cmp(&0.0) {
- Some(Greater) => {
- lb += λ;
- ub += λ
- }
- // Less should not happen
- Some(Less | Equal) => {
- lb = F::NEG_INFINITY;
- ub += λ
- }
- None => {}
- };
+ if x_i < F::EPSILON {
+ lb = F::NEG_INFINITY;
+ }
val = max_interval_dist_to_zero(val, lb, ub);
}
val
@@ -164,11 +154,14 @@
/// (`shift` in the code), and below a value that would cause changes in the locked set
/// (`max_shift` in the code).
#[replace_float_literals(F::cast_from(literal))]
-pub fn l1squared_nonneg_prox<F: Float + nalgebra::RealField>(
- x: &mut DVector<F>,
- y: &DVector<F>,
+pub fn l1squared_nonneg_prox<F: Float + nalgebra::RealField, S1, S2>(
+ x: &mut Vector<F, Dyn, S1>,
+ y: &Vector<F, Dyn, S2>,
β: F,
-) {
+) where
+ S2: Storage<F, Dyn>,
+ S1: StorageMut<F, Dyn>,
+{
// nalgebra double-definition bullshit workaround
let abs = alg_tools::NumTraitsFloat::abs;
let min = alg_tools::NumTraitsFloat::min;
@@ -218,12 +211,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_nonneg_pp<F, I>(
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
+pub fn l1squared_nonneg_pp<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,
@@ -231,9 +223,12 @@
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 mut τ = τ_.to_nalgebra_mixed();
let θ = θ_.to_nalgebra_mixed();
let mut iters = 0;
@@ -242,7 +237,7 @@
// Primal step: x^{k+1} = prox_{(τβ/2)|.-y|_1^2+δ_{≥0}+}(x^k - τ(λ𝟙^⊤-g))
x.apply(|x_i| *x_i -= τ * λ);
x.axpy(τ, g, 1.0);
- l1squared_nonneg_prox(x, y, τ * β);
+ l1squared_nonneg_prox(x, y, τ);
iters += 1;
// This gives O(1/N^2) rates due to monotonicity of function values.
@@ -253,7 +248,7 @@
// Higher acceleration does not seem to be numerically stable.
//τ + = F::cast_from(iters).to_nalgebra_mixed()*θ;
- 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
@@ -266,12 +261,11 @@
/// for potential performance improvements.
/// The parameter `θ` is used to multiply the rescale the operator (identity) of the PDPS model.
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn l1squared_nonneg_pdps<F, I>(
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
+pub fn l1squared_nonneg_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,
θ_: F,
@@ -280,22 +274,25 @@
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 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_{(τβ/2)|.-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);
@@ -305,7 +302,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
@@ -323,18 +320,17 @@
/// We rewrite
/// <div>$$
/// \begin{split}
-/// & \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x) \\
+/// & \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x) \\
/// & = \min_{x ∈ ℝ^n} \max_{w} ⟨θ w, x⟩ - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x)
-/// - \left(x ↦ \frac{β}{2θ} |x-y|_1^2 \right)^*(w).
+/// - \left(x ↦ \frac{1}{2θ} |x-y|_1^2 \right)^*(w).
/// \end{split}
/// $$</div>
#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
-pub fn l1squared_nonneg_pdps_alt<F, I>(
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
+pub fn l1squared_nonneg_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,
@@ -343,80 +339,282 @@
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 τλ = τ * λ;
+ let one_div_σθ = 1.0 / σθ;
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;
+ let mut y_scale = y.clone_owned();
+ y_scale *= σ;
+
iterator.iterate(|state| {
// Primal step: x^{k+1} = nonnegsoft_τλ(x^k - τ(θ w^k -g))
- x.axpy(-τθ, &w, 1.0);
- x.axpy(τ, g, 1.0);
- x.apply(|x_i| *x_i = nonneg_soft_thresholding(*x_i, τ * λ));
+ if θ == 1.0 {
+ for (x_i, xprev_i, w_i, &g_i) in
+ izip!(x.iter_mut(), xprev.iter_mut(), w.iter_mut(), g.iter())
+ {
+ *x_i = nonneg_soft_thresholding(*x_i - τ * (*w_i - g_i), τλ);
+ // Fused dual part from below
+ *w_i += σ * (2.0 * *x_i - *xprev_i);
+ *xprev_i = *w_i;
+ }
+ } else {
+ for (x_i, xprev_i, w_i, &g_i) in
+ izip!(x.iter_mut(), xprev.iter_mut(), w.iter_mut(), g.iter())
+ {
+ *x_i = nonneg_soft_thresholding(*x_i - τ * (θ * *w_i - g_i), τλ);
+ // Fused dual part from below
+ *w_i += σ * (2.0 * *x_i - *xprev_i);
+ *xprev_i = *w_i;
+ }
+ }
+ // This is numerically unstable:
+ // x.axpy(-τθ, &w, 1.0);
+ // x.axpy(τ, g, 1.0);
+ // x.apply(|x_i| *x_i = nonneg_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/σ))
+ // ALT
+ // = 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, β / σθ);
+
+ // This has been fused into the loop below
+ // for (xprev_i, w_i, &x_i) in izip!(xprev.iter_mut(), w.iter_mut(), x.iter()) {
+ // *w_i += σ * (2.0 * x_i - *xprev_i);
+ // *xprev_i = *w_i;
+ // }
+ // xprev.axpy(2.0, x, -1.0);
+ // w.axpy(σ, &xprev, 1.0);
+ // xprev.copy_from(&w); // use xprev as temporary variable
+ //l1squared_prox(&mut tmp, &mut xprev, &y_scale, β_div_σθ);
+ l1squared_prox(&mut tmp, &mut xprev, &y_scale, one_div_σθ);
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
}
+/// This is an exact solver for
+/// <div>$$
+/// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x).
+/// $$</div>
+/// i.e.,
+/// <div>$$
+/// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 + (λ - g)^⊤ x + δ_{≥ 0}(x),
+/// $$</div>
+/// which has the optimality conditions
+/// <div>$$
+/// 0 ∈ |x-y|_1\sign(x-y)_i + λ - g_i + δ_{≥ 0}(x_i)
+/// \quad\text{for all}\quad i.
+/// $$</div>
+///
+/// If $x_i > 0$ and $x_i ≠ y_i$, then this forces $|g_i -λ| = |x-y|_1$.
+/// Let $i^*$ be such an index. Then if $i ≠ i^*$ does not have the same value of $|g_i-λ|$,
+/// we *must* have $x_i=y_i$ or $x_i=0$. In fact, even $x_i=y_i$ is impossible for $y_i>0$ if
+/// $|g_i-λ| > |g_{i^*}-λ|$, whereas, otherwise $x_i=y_i$ is possible if $y_i>=0$.
+/// If $|g_i-λ| > |g_{i^*}-λ|$, we must, therefore, either have $x_i=0$ or be able to take $i^*=i$.
+///
+/// If $|g_i-λ| ≤ |g_{i^*}-λ|$, we can have either $x_i=0$ or $x_i=y_i$. If $y_i<=0$, then,
+/// of course $x_i=0$. If $y_i>0$, we can take $x_i=y_i$.
+///
+/// We, therefore, sort the |g_i-λ|, and look for an index $i^*$ that gives a non-contradictory
+/// $β|x-y|_1=|g_i -λ|$, noting that contributions to β|x-y|_1= come from, besides $i^*$, from
+/// indices $i$ such that $x_i= 0 ≠ y_i$.
+///
+/// Finally, if no index $i$ reaches $|g_i -λ| = |x-y|_1$, we have to $|x-y|_1$ being in an
+/// intermediate interval. This can similarly done using sorting.
+#[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
+pub(super) fn l1squared_nonneg_solve_exact<F, S1, S2, S3>(
+ y: &Vector<F::MixedType, Dyn, S1>,
+ g: &Vector<F::MixedType, Dyn, S2>,
+ λ_: F,
+ x: &mut Vector<F::MixedType, Dyn, S3>,
+) -> bool
+where
+ F: Float + ToNalgebraRealField,
+ 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();
+
+ assert_eq!(y.len(), g.len());
+ assert_eq!(x.len(), g.len());
+ assert!(x.len() <= u32::MAX as usize);
+
+ #[derive(Debug)]
+ struct Tmp<F> {
+ y: F,
+ d: F,
+ contrib: F,
+ i: u32, // Shrink size for sort
+ λ_le_g: bool,
+ }
+ let mut sorted = Vec::from_iter(izip!(g.iter(), y.iter(), 0..).map(|(&g_i, &y_i, i)| {
+ // We already precompute the comparison λ <= g_i here, since we need it for the abs,
+ // and would need to store g_i for that in any case.
+ let (d, λ_le_g) = if λ <= g_i {
+ (g_i - λ, true)
+ } else {
+ (λ - g_i, false)
+ };
+ Tmp { d, λ_le_g, y: y_i, i, contrib: 0.0 }
+ }));
+ sorted
+ .as_mut_slice()
+ .sort_unstable_by(|a, b| b.d.total_cmp(&a.d));
+
+ // Reverse-compute contribs
+ sorted.iter_mut().rev().fold(0.0, |contrib, a| {
+ a.contrib = contrib;
+ if a.y < 0.0 {
+ contrib - a.y
+ } else {
+ contrib
+ }
+ });
+
+ let mut contrib0 = 0.0;
+
+ x.fill(0.0);
+
+ let mut it = sorted.iter();
+ let mut found = false;
+
+ // We first try to find an index m such that y_m ≠ x_m > 0.0 that satisfies ‖x-y‖₁ = |λ - g_m|.
+ 'search: while let Some(m) = it.next() {
+ let d_m = m.d;
+ let y_m = m.y;
+ let contrib = m.contrib + contrib0;
+ let δ = d_m - contrib;
+ if δ > 0.0 {
+ // If λ < g[m], we must have x[m]≥y[m] so (contrib + x[m]-y[m])=g[m]-λ=d
+ // If λ = g_m, then also d_m=0, so contrib=0, and this gtives *x_m=y_m.
+ // If λ > g[m], We must have x[m]≤y[m] so -(contrib + y[m]-x[m])=g[m]-λ=-d.
+ let x_m_prime = if m.λ_le_g { y_m + δ } else { y_m - δ };
+ if x_m_prime > 0.0 {
+ let x_m = unsafe { x.get_unchecked_mut(m.i as usize) };
+ *x_m = x_m_prime;
+ found = true;
+ break 'search;
+ }
+ }
+ contrib0 += num_traits::abs(y_m);
+ }
+
+ if !found {
+ // Every x_i is either zero or equal to y_i.
+ // We scan intervals (bound,prev_bound) between the d-values of each component,
+ // for the value of ‖x-y‖₁. Then it can be decided whether x_i should be one or equal to
+ // y_i. Due to the sorting and pre-calculation of posterior contributions, this simplifies
+ // into the following:
+ let mut prev_bound = F::MixedType::INFINITY;
+ let mut contrib0 = 0.0;
+ it = sorted.iter();
+ 'search_degenerate: while let Some(m) = it.next() {
+ let bound = m.d;
+ let contrib = m.contrib + contrib0 - m.y.min(0.0);
+ if prev_bound >= contrib && contrib >= bound {
+ if m.y > 0.0 {
+ let x_m = unsafe { x.get_unchecked_mut(m.i as usize) };
+ *x_m = m.y;
+ }
+ found = true;
+ break 'search_degenerate;
+ }
+ prev_bound = bound;
+ if m.y > 0.0 {
+ contrib0 += m.y;
+ }
+ }
+ // Check final interval from last element to -∞, if nothing found.
+ if !found && !(prev_bound >= contrib0) {
+ panic!("l1squared_nonneg_solve_exact failure")
+ }
+ }
+
+ // Set remaining components to their now known values
+ while let Some(a) = it.next() {
+ // Safety: size checked above.
+ let y_i = a.y;
+ if y_i >= 0.0 {
+ // We can leave unchanged, as ‖x-y‖₁ is guaranteed large enough if
+ // the current maximising index attempt succeeds.
+ let x_i = unsafe { x.get_unchecked_mut(a.i as usize) };
+ *x_i = y_i
+ }
+ // Zero fill in `else` case done above in init already
+ }
+
+ debug_assert!(min_subdifferential(y, x, g, λ) <= F::EPSILON.to_nalgebra_mixed() * 10.0);
+
+ return true;
+}
+
/// 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\|₁ + δ_{≥ 0}(x).
+/// \min_{x ∈ ℝ^n} \frac{1}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁ + δ_{≥ 0}(x).
/// $$</div>
///
/// This function returns the number of iterations taken.
#[replace_float_literals(F::cast_from(literal))]
-pub fn l1squared_nonneg<F, I>(
- y: &DVector<F::MixedType>,
- g: &DVector<F::MixedType>,
+pub fn l1squared_nonneg<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>,
{
+ if let InnerMethod::Exact = inner.method {
+ // Try exact solution, fall back to PDPS if it does not work.
+ if l1squared_nonneg_solve_exact(y, g, λ, x) {
+ return 1;
+ }
+ }
+
match inner.method {
- InnerMethod::PDPS => {
+ InnerMethod::PDPS | InnerMethod::Exact => {
let inner_θ = 1.0;
- // Estimate of ‖K‖ for K=θ\Id.
+ //Estimate of ‖K‖ for K=θ\Id.
let normest = inner_θ;
let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest);
- l1squared_nonneg_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator)
+ l1squared_nonneg_pdps_alt(y, g, λ, x, inner_τ, inner_σ, inner_θ, iterator)
}
InnerMethod::PP | InnerMethod::FB => {
let inner_τ = inner.pp_τ.0;
let inner_θ = inner.pp_τ.1;
- l1squared_nonneg_pp(y, g, λ, β, x, inner_τ, inner_θ, iterator)
+ l1squared_nonneg_pp(y, g, λ, x, inner_τ, inner_θ, iterator)
}
other => unimplemented!("${other:?} is unimplemented"),
}