# HG changeset patch # User Tuomo Valkonen # Date 1745784200 18000 # Node ID 83a95f2e1f4d1bcca7f4ec0f08a5a9dd1d81cea6 # Parent 0693cc9ba9f0c97d440179bd22f808ccea6a49d2# Parent 6099ba025aac7dd50fa356f4f56d7935f63ccc16 Start work on 3.0.0 diff -r 0693cc9ba9f0 -r 83a95f2e1f4d .hgtags --- a/.hgtags Mon Feb 17 13:51:50 2025 -0500 +++ b/.hgtags Sun Apr 27 15:03:20 2025 -0500 @@ -1,2 +1,3 @@ b71edfd403aaf3045011769ce14775bfff94c4f5 v1.0.0-pre-arxiv 12cee797937ff49574f95c766dc0b56395ed9f4f v1.0.0 +b3312eee105c4af232cb90a9c5f3a4888366b203 v2.0.0-pre diff -r 0693cc9ba9f0 -r 83a95f2e1f4d Cargo.lock --- a/Cargo.lock Mon Feb 17 13:51:50 2025 -0500 +++ b/Cargo.lock Sun Apr 27 15:03:20 2025 -0500 @@ -33,7 +33,7 @@ [[package]] name = "alg_tools" -version = "0.3.0" +version = "0.3.1" dependencies = [ "anyhow", "colored", diff -r 0693cc9ba9f0 -r 83a95f2e1f4d Cargo.toml --- a/Cargo.toml Mon Feb 17 13:51:50 2025 -0500 +++ b/Cargo.toml Sun Apr 27 15:03:20 2025 -0500 @@ -1,6 +1,6 @@ [package] name = "pointsource_algs" -version = "2.0.0-pre" +version = "3.0.0-dev" edition = "2021" rust-version = "1.85" authors = ["Tuomo Valkonen "] @@ -22,7 +22,7 @@ categories = ["mathematics", "science", "computer-vision"] [dependencies.alg_tools] -version = "~0.3.0-dev" +version = "~0.3.1-dev" path = "../alg_tools" default-features = false features = ["nightly"] diff -r 0693cc9ba9f0 -r 83a95f2e1f4d README.md --- a/README.md Mon Feb 17 13:51:50 2025 -0500 +++ b/README.md Sun Apr 27 15:03:20 2025 -0500 @@ -4,7 +4,7 @@ This package contains the [Rust] codes for the numerical experiments in the articles * T. Valkonen, “_Proximal methods for point source localisation_”, Journal of Nonsmooth Analysis and Optimization 4 (2023), 10433, [doi:10.46298/jnsao-2023-10433] ([arXiv:2212.02991]) -* T. Valkonen, “_Point source localisation with unbalanced optimal transport_” (2025), submitted. +* T. Valkonen, “_Point source localisation with unbalanced optimal transport_” (2025), submitted, [arXiv:2502.12417] It concerns solution of problems of the type $$ @@ -49,6 +49,7 @@ [rust-GSL]: https://docs.rs/GSL/6.0.0/rgsl/ [Homebrew]: https://brew.sh [arXiv:2212.02991]: https://arxiv.org/abs/2212.02991 + [arXiv:2502.12417]: https://arxiv.org/abs/2502.12417 [doi:10.46298/jnsao-2023-10433]: http://doi.org/10.46298/jnsao-2023-10433 ### Building and running the experiments @@ -69,7 +70,7 @@ ```console cargo run --release -- -o results ``` -The double-dash separates the options for the Cargo build system +The double-dash separates the options for the Cargo build system and `pointsource_algs`. ### Documentation @@ -89,4 +90,3 @@ generated documentation through an ugly workaround. Unfortunately, `rustdoc`, akin to Rust largely itself, is stuck in 80's 7-bit gringo ASCII world, and does not support modern markdown features, such as mathematics. - diff -r 0693cc9ba9f0 -r 83a95f2e1f4d misc/katex-header.html --- a/misc/katex-header.html Mon Feb 17 13:51:50 2025 -0500 +++ b/misc/katex-header.html Sun Apr 27 15:03:20 2025 -0500 @@ -1,15 +1,63 @@ - - - + + + diff -r 0693cc9ba9f0 -r 83a95f2e1f4d src/seminorms.rs --- a/src/seminorms.rs Mon Feb 17 13:51:50 2025 -0500 +++ b/src/seminorms.rs Sun Apr 27 15:03:20 2025 -0500 @@ -4,31 +4,31 @@ the trait [`DiscreteMeasureOp`]. */ -use std::iter::Zip; -use std::ops::RangeFrom; -use alg_tools::types::*; +use crate::measures::{DeltaMeasure, DiscreteMeasure, Radon, SpikeIter, RNDM}; +use alg_tools::bisection_tree::*; +use alg_tools::instance::Instance; +use alg_tools::iter::{FilterMapX, Mappable}; +use alg_tools::linops::{BoundedLinear, Linear, Mapping}; use alg_tools::loc::Loc; -use alg_tools::sets::Cube; -use alg_tools::bisection_tree::*; use alg_tools::mapping::RealMapping; -use alg_tools::iter::{Mappable, FilterMapX}; -use alg_tools::linops::{Mapping, Linear, BoundedLinear}; -use alg_tools::instance::Instance; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::norms::Linfinity; -use crate::measures::{DiscreteMeasure, DeltaMeasure, SpikeIter, Radon, RNDM}; +use alg_tools::sets::Cube; +use alg_tools::types::*; +use itertools::Itertools; use nalgebra::DMatrix; +use std::iter::Zip; use std::marker::PhantomData; -use itertools::Itertools; +use std::ops::RangeFrom; /// Abstraction for operators $𝒟 ∈ 𝕃(𝒵(Ω); C_c(Ω))$. /// /// Here $𝒵(Ω) ⊂ ℳ(Ω)$ is the space of sums of delta measures, presented by [`DiscreteMeasure`]. -pub trait DiscreteMeasureOp - : BoundedLinear, Radon, Linfinity, F> +pub trait DiscreteMeasureOp: + BoundedLinear, Radon, Linfinity, F> where - F : Float + ToNalgebraRealField, - Domain : 'static + Clone + PartialEq, + F: Float + ToNalgebraRealField, + Domain: 'static + Clone + PartialEq, { /// The output type of [`Self::preapply`]. type PreCodomain; @@ -40,12 +40,13 @@ /// for a $x_1, …, x_n$ the coordinates given by the iterator `I`, and $a=(α_1,…,α_n)$. /// Here $C_* ∈ 𝕃(C_c(Ω); ℝ^n) $ stands for the preadjoint. ///

- fn findim_matrix<'a, I>(&self, points : I) -> DMatrix - where I : ExactSizeIterator + Clone; + fn findim_matrix<'a, I>(&self, points: I) -> DMatrix + where + I: ExactSizeIterator + Clone; /// [`Mapping`] that typically returns an uninitialised [`PreBTFN`] /// instead of a full [`BTFN`]. - fn preapply(&self, μ : DiscreteMeasure) -> Self::PreCodomain; + fn preapply(&self, μ: DiscreteMeasure) -> Self::PreCodomain; } // Blanket implementation of a measure as a linear functional over a predual @@ -67,28 +68,37 @@ // /// A trait alias for simple convolution kernels. -pub trait SimpleConvolutionKernel -: RealMapping + Support + Bounded + Clone + 'static {} +pub trait SimpleConvolutionKernel: + RealMapping + Support + Bounded + Clone + 'static +{ +} -impl SimpleConvolutionKernel for T -where T : RealMapping + Support + Bounded + Clone + 'static {} +impl SimpleConvolutionKernel for T where + T: RealMapping + Support + Bounded + Clone + 'static +{ +} /// [`SupportGenerator`] for [`ConvolutionOp`]. -#[derive(Clone,Debug)] -pub struct ConvolutionSupportGenerator -where K : SimpleConvolutionKernel { - kernel : K, - centres : RNDM, +#[derive(Clone, Debug)] +pub struct ConvolutionSupportGenerator +where + K: SimpleConvolutionKernel, +{ + kernel: K, + centres: RNDM, } -impl ConvolutionSupportGenerator -where K : SimpleConvolutionKernel { - +impl ConvolutionSupportGenerator +where + K: SimpleConvolutionKernel, +{ /// Construct the convolution kernel corresponding to `δ`, i.e., one centered at `δ.x` and /// weighted by `δ.α`. #[inline] - fn construct_kernel<'a>(&'a self, δ : &'a DeltaMeasure, F>) - -> Weighted, F> { + fn construct_kernel<'a>( + &'a self, + δ: &'a DeltaMeasure, F>, + ) -> Weighted, F> { self.kernel.clone().shift(δ.x).weigh(δ.α) } @@ -98,22 +108,27 @@ #[inline] fn construct_kernel_and_id_filtered<'a>( &'a self, - (id, δ) : (usize, &'a DeltaMeasure, F>) + (id, δ): (usize, &'a DeltaMeasure, F>), ) -> Option<(usize, Weighted, F>)> { (δ.α != F::ZERO).then(|| (id.into(), self.construct_kernel(δ))) } } -impl SupportGenerator -for ConvolutionSupportGenerator -where K : SimpleConvolutionKernel { +impl SupportGenerator for ConvolutionSupportGenerator +where + K: SimpleConvolutionKernel, +{ type Id = usize; type SupportType = Weighted, F>; - type AllDataIter<'a> = FilterMapX<'a, Zip, SpikeIter<'a, Loc, F>>, - Self, (Self::Id, Self::SupportType)>; + type AllDataIter<'a> = FilterMapX< + 'a, + Zip, SpikeIter<'a, Loc, F>>, + Self, + (Self::Id, Self::SupportType), + >; #[inline] - fn support_for(&self, d : Self::Id) -> Self::SupportType { + fn support_for(&self, d: Self::Id) -> Self::SupportType { self.construct_kernel(&self.centres[d]) } @@ -124,51 +139,53 @@ #[inline] fn all_data(&self) -> Self::AllDataIter<'_> { - (0..).zip(self.centres.iter_spikes()) - .filter_mapX(self, Self::construct_kernel_and_id_filtered) + (0..) + .zip(self.centres.iter_spikes()) + .filter_mapX(self, Self::construct_kernel_and_id_filtered) } } /// Representation of a convolution operator $𝒟$. -#[derive(Clone,Debug)] -pub struct ConvolutionOp -where F : Float + ToNalgebraRealField, - BT : BTImpl, - K : SimpleConvolutionKernel { +#[derive(Clone, Debug)] +pub struct ConvolutionOp +where + F: Float + ToNalgebraRealField, + BT: BTImpl, + K: SimpleConvolutionKernel, +{ /// Depth of the [`BT`] bisection tree for the outputs [`Mapping::apply`]. - depth : BT::Depth, + depth: BT::Depth, /// Domain of the [`BT`] bisection tree for the outputs [`Mapping::apply`]. - domain : Cube, + domain: Cube, /// The convolution kernel - kernel : K, - _phantoms : PhantomData<(F,BT)>, + kernel: K, + _phantoms: PhantomData<(F, BT)>, } -impl ConvolutionOp -where F : Float + ToNalgebraRealField, - BT : BTImpl, - K : SimpleConvolutionKernel { - +impl ConvolutionOp +where + F: Float + ToNalgebraRealField, + BT: BTImpl, + K: SimpleConvolutionKernel, +{ /// Creates a new convolution operator $𝒟$ with `kernel` on `domain`. /// /// The output of [`Mapping::apply`] is a [`BT`] of given `depth`. - pub fn new(depth : BT::Depth, domain : Cube, kernel : K) -> Self { + pub fn new(depth: BT::Depth, domain: Cube, kernel: K) -> Self { ConvolutionOp { - depth : depth, - domain : domain, - kernel : kernel, - _phantoms : PhantomData + depth: depth, + domain: domain, + kernel: kernel, + _phantoms: PhantomData, } } /// Returns the support generator for this convolution operator. - fn support_generator(&self, μ : RNDM) - -> ConvolutionSupportGenerator { - + fn support_generator(&self, μ: RNDM) -> ConvolutionSupportGenerator { // TODO: can we avoid cloning μ? ConvolutionSupportGenerator { - kernel : self.kernel.clone(), - centres : μ + kernel: self.kernel.clone(), + centres: μ, } } @@ -178,43 +195,43 @@ } } -impl Mapping> -for ConvolutionOp +impl Mapping> for ConvolutionOp where - F : Float + ToNalgebraRealField, - BT : BTImpl, - K : SimpleConvolutionKernel, - Weighted, F> : LocalAnalysis + F: Float + ToNalgebraRealField, + BT: BTImpl, + K: SimpleConvolutionKernel, + Weighted, F>: LocalAnalysis, { - type Codomain = BTFN, BT, N>; - fn apply(&self, μ : I) -> Self::Codomain - where I : Instance> { + fn apply(&self, μ: I) -> Self::Codomain + where + I: Instance>, + { let g = self.support_generator(μ.own()); BTFN::construct(self.domain.clone(), self.depth, g) } } /// [`ConvolutionOp`]s as linear operators over [`DiscreteMeasure`]s. -impl Linear> -for ConvolutionOp +impl Linear> for ConvolutionOp where - F : Float + ToNalgebraRealField, - BT : BTImpl, - K : SimpleConvolutionKernel, - Weighted, F> : LocalAnalysis -{ } + F: Float + ToNalgebraRealField, + BT: BTImpl, + K: SimpleConvolutionKernel, + Weighted, F>: LocalAnalysis, +{ +} -impl -BoundedLinear, Radon, Linfinity, F> -for ConvolutionOp -where F : Float + ToNalgebraRealField, - BT : BTImpl, - K : SimpleConvolutionKernel, - Weighted, F> : LocalAnalysis { - - fn opnorm_bound(&self, _ : Radon, _ : Linfinity) -> F { +impl BoundedLinear, Radon, Linfinity, F> + for ConvolutionOp +where + F: Float + ToNalgebraRealField, + BT: BTImpl, + K: SimpleConvolutionKernel, + Weighted, F>: LocalAnalysis, +{ + fn opnorm_bound(&self, _: Radon, _: Linfinity) -> F { // With μ = ∑_i α_i δ_{x_i}, we have // |𝒟μ|_∞ // = sup_z |∑_i α_i φ(z - x_i)| @@ -225,31 +242,33 @@ } } - -impl DiscreteMeasureOp, F> -for ConvolutionOp -where F : Float + ToNalgebraRealField, - BT : BTImpl, - K : SimpleConvolutionKernel, - Weighted, F> : LocalAnalysis { +impl DiscreteMeasureOp, F> for ConvolutionOp +where + F: Float + ToNalgebraRealField, + BT: BTImpl, + K: SimpleConvolutionKernel, + Weighted, F>: LocalAnalysis, +{ type PreCodomain = PreBTFN, N>; - fn findim_matrix<'a, I>(&self, points : I) -> DMatrix - where I : ExactSizeIterator> + Clone { + fn findim_matrix<'a, I>(&self, points: I) -> DMatrix + where + I: ExactSizeIterator> + Clone, + { // TODO: Preliminary implementation. It be best to use sparse matrices or // possibly explicit operators without matrices let n = points.len(); let points_clone = points.clone(); let pairs = points.cartesian_product(points_clone); let kernel = &self.kernel; - let values = pairs.map(|(x, y)| kernel.apply(y-x).to_nalgebra_mixed()); + let values = pairs.map(|(x, y)| kernel.apply(y - x).to_nalgebra_mixed()); DMatrix::from_iterator(n, n, values) } /// A version of [`Mapping::apply`] that does not instantiate the [`BTFN`] codomain with /// a bisection tree, instead returning a [`PreBTFN`]. This can improve performance when /// the output is to be added as the right-hand-side operand to a proper BTFN. - fn preapply(&self, μ : RNDM) -> Self::PreCodomain { + fn preapply(&self, μ: RNDM) -> Self::PreCodomain { BTFN::new_pre(self.support_generator(μ)) } } @@ -258,47 +277,46 @@ /// for [`ConvolutionSupportGenerator`]. macro_rules! make_convolutionsupportgenerator_scalarop_rhs { ($trait:ident, $fn:ident, $trait_assign:ident, $fn_assign:ident) => { - impl, const N : usize> - std::ops::$trait_assign - for ConvolutionSupportGenerator { - fn $fn_assign(&mut self, t : F) { + impl, const N: usize> std::ops::$trait_assign + for ConvolutionSupportGenerator + { + fn $fn_assign(&mut self, t: F) { self.centres.$fn_assign(t); } } - impl, const N : usize> - std::ops::$trait - for ConvolutionSupportGenerator { + impl, const N: usize> std::ops::$trait + for ConvolutionSupportGenerator + { type Output = ConvolutionSupportGenerator; - fn $fn(mut self, t : F) -> Self::Output { + fn $fn(mut self, t: F) -> Self::Output { std::ops::$trait_assign::$fn_assign(&mut self.centres, t); self } } - impl<'a, F : Float, K : SimpleConvolutionKernel, const N : usize> - std::ops::$trait - for &'a ConvolutionSupportGenerator { + impl<'a, F: Float, K: SimpleConvolutionKernel, const N: usize> std::ops::$trait + for &'a ConvolutionSupportGenerator + { type Output = ConvolutionSupportGenerator; - fn $fn(self, t : F) -> Self::Output { - ConvolutionSupportGenerator{ - kernel : self.kernel.clone(), - centres : (&self.centres).$fn(t), + fn $fn(self, t: F) -> Self::Output { + ConvolutionSupportGenerator { + kernel: self.kernel.clone(), + centres: (&self.centres).$fn(t), } } } - } + }; } make_convolutionsupportgenerator_scalarop_rhs!(Mul, mul, MulAssign, mul_assign); make_convolutionsupportgenerator_scalarop_rhs!(Div, div, DivAssign, div_assign); - /// Generates an unary operation (e.g. [`std::ops::Neg`]) for [`ConvolutionSupportGenerator`]. macro_rules! make_convolutionsupportgenerator_unaryop { ($trait:ident, $fn:ident) => { - impl, const N : usize> - std::ops::$trait - for ConvolutionSupportGenerator { + impl, const N: usize> std::ops::$trait + for ConvolutionSupportGenerator + { type Output = ConvolutionSupportGenerator; fn $fn(mut self) -> Self::Output { self.centres = self.centres.$fn(); @@ -306,19 +324,18 @@ } } - impl<'a, F : Float, K : SimpleConvolutionKernel, const N : usize> - std::ops::$trait - for &'a ConvolutionSupportGenerator { + impl<'a, F: Float, K: SimpleConvolutionKernel, const N: usize> std::ops::$trait + for &'a ConvolutionSupportGenerator + { type Output = ConvolutionSupportGenerator; fn $fn(self) -> Self::Output { - ConvolutionSupportGenerator{ - kernel : self.kernel.clone(), - centres : (&self.centres).$fn(), + ConvolutionSupportGenerator { + kernel: self.kernel.clone(), + centres: (&self.centres).$fn(), } } } - } + }; } make_convolutionsupportgenerator_unaryop!(Neg, neg); - diff -r 0693cc9ba9f0 -r 83a95f2e1f4d src/subproblem/l1squared_unconstrained.rs --- a/src/subproblem/l1squared_unconstrained.rs Mon Feb 17 13:51:50 2025 -0500 +++ b/src/subproblem/l1squared_unconstrained.rs Sun Apr 27 15:03:20 2025 -0500 @@ -2,76 +2,74 @@ Iterative algorithms for solving the finite-dimensional subproblem without constraints. */ +use itertools::izip; use nalgebra::DVector; use numeric_literals::replace_float_literals; -use itertools::izip; use std::cmp::Ordering::*; -use std::iter::zip; -use alg_tools::iterate::{ - AlgIteratorFactory, - AlgIteratorState, -}; +use alg_tools::iterate::{AlgIteratorFactory, AlgIteratorState}; use alg_tools::nalgebra_support::ToNalgebraRealField; use alg_tools::nanleast::NaNLeast; use alg_tools::norms::{Dist, L1}; +use std::iter::zip; +use super::l1squared_nonneg::max_interval_dist_to_zero; +use super::unconstrained::soft_thresholding; +use super::{InnerMethod, InnerSettings}; use crate::types::*; -use super::{ - InnerMethod, - InnerSettings -}; -use super::unconstrained::soft_thresholding; -use super::l1squared_nonneg::max_interval_dist_to_zero; -/// Calculate $prox_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}_1^2$. +/// Calculate $\prox_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}_1^2$. /// /// To derive an algorithm for this, we can assume that $y=0$, as -/// $prox_f(x) = prox_{f_0}(x - y) - y$ for $f_0=\frac{β}{2}\norm{x}_1^2$. -/// Now, the optimality conditions for $w = prox_f(x)$ are -/// $$\tag{*} -/// 0 ∈ w-x + β\norm{w}_1\sign w. +/// $\prox\_f(x) = \prox\_{f_0}(x - y) - y$ for $f\_0=\frac{β}{2}\norm{x}\_1^2$. +/// Now, the optimality conditions for $w = \prox\_f(x)$ are /// $$ -/// Clearly then $w = \soft_{β\norm{w}_1}(x)$. +/// 0 ∈ w-x + β\norm{w}\_1\sign w. +/// $$ +/// Clearly then $w = \soft\_{β\norm{w}\_1}(x)$. /// Thus the components of $x$ with smallest absolute value will be zeroed out. /// Denoting by $w'$ the non-zero components, and by $x'$ the corresponding components /// of $x$, and by $m$ their count, multipying the corresponding lines of (*) by $\sign x'$, /// we obtain /// $$ -/// \norm{x'}_1 = (1+βm)\norm{w'}_1. +/// \norm{x'}\_1 = (1+βm)\norm{w'}\_1. /// $$ -/// That is, $\norm{w}_1=\norm{w'}_1=\norm{x'}_1/(1+βm)$. +/// That is, $\norm{w}\_1=\norm{w'}\_1=\norm{x'}\_1/(1+βm)$. /// Thus, sorting $x$ by absolute value, and sequentially in order eliminating the smallest -/// elements, we can easily calculate what $\norm{w}_1$ should be for that choice, and -/// then easily calculate $w = \soft_{β\norm{w}_1}(x)$. We just have to verify that +/// elements, we can easily calculate what $\norm{w}\_1$ should be for that choice, and +/// then easily calculate $w = \soft_{β\norm{w}\_1}(x)$. We just have to verify that /// the resulting $w$ has the same norm. There's a shortcut to this, as we work /// sequentially: just check that the smallest assumed-nonzero component $i$ satisfies the -/// condition of soft-thresholding to remain non-zero: $|x_i|>τ\norm{x'}/(1+τm)$. -/// Clearly, if this condition fails for x_i, it will fail for all the components +/// condition of soft-thresholding to remain non-zero: $|x\_i|>τ\norm{x'}/(1+τm)$. +/// 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( - sorted_abs : &mut DVector, - x : &mut DVector, - y : &DVector, - β : F +pub(super) fn l1squared_prox( + sorted_abs: &mut DVector, + x: &mut DVector, + y: &DVector, + β: F, ) { 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(|a, b| NaNLeast(*a).cmp(&NaNLeast(*b))); 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)); + let tmp = β * n / (1.0 + β * F::cast_from(m)); if *az_i <= tmp { // Fail n -= *az_i; } else { // Success - x.zip_apply(y, |x_i, y_i| *x_i = y_i + soft_thresholding(*x_i-y_i, tmp)); - return + x.zip_apply(y, |x_i, y_i| { + *x_i = y_i + soft_thresholding(*x_i - y_i, tmp) + }); + return; } } // m = 0 should always work, but x is zero. @@ -82,35 +80,52 @@ /// /// `v` will be modified and cannot be trusted to contain useful values afterwards. #[replace_float_literals(F::cast_from(literal))] -fn min_subdifferential( - y : &DVector, - x : &DVector, - g : &DVector, - λ : F, - β : F +fn min_subdifferential( + y: &DVector, + x: &DVector, + g: &DVector, + λ: F, + β: F, ) -> F { 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 => {}, + Some(Greater) => { + lb += tmp; + ub += tmp + } + Some(Less) => { + lb -= tmp; + ub -= tmp + } + Some(Equal) => { + lb -= tmp; + ub += tmp + } + None => {} } match x_i.partial_cmp(&0.0) { - Some(Greater) => { lb += λ; ub += λ }, - Some(Less) => { lb -= λ; ub -= λ }, - Some(Equal) => { lb -= λ; ub += λ }, - None => {}, + Some(Greater) => { + lb += λ; + ub += λ + } + Some(Less) => { + lb -= λ; + ub -= λ + } + Some(Equal) => { + lb -= λ; + ub += λ + } + None => {} }; val = max_interval_dist_to_zero(val, lb, ub); } val } - /// PDPS implementation of [`l1squared_unconstrained`]. /// For detailed documentation of the inputs and outputs, refer to there. /// @@ -118,17 +133,18 @@ /// for potential performance improvements. #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] pub fn l1squared_unconstrained_pdps( - y : &DVector, - g : &DVector, - λ_ : F, - β_ : F, - x : &mut DVector, - τ_ : F, - σ_ : F, - iterator : I + y: &DVector, + g: &DVector, + λ_: F, + β_: F, + x: &mut DVector, + τ_: F, + σ_: F, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory, { let λ = λ_.to_nalgebra_mixed(); let β = β_.to_nalgebra_mixed(); @@ -143,19 +159,17 @@ // 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); + w.axpy(2.0 * σ, x, 1.0); w.axpy(-σ, &xprev, 1.0); w.apply(|w_i| *w_i = num_traits::clamp(*w_i, -λ, λ)); xprev.copy_from(x); - - iters +=1; - state.if_verbose(|| { - F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β)) - }) + iters += 1; + + state.if_verbose(|| F::from_nalgebra_mixed(min_subdifferential(y, x, g, λ, β))) }); iters @@ -180,26 +194,27 @@ /// $$ #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] pub fn l1squared_unconstrained_pdps_alt( - y : &DVector, - g : &DVector, - λ_ : F, - β_ : F, - x : &mut DVector, - τ_ : F, - σ_ : F, - θ_ : F, - iterator : I + y: &DVector, + g: &DVector, + λ_: F, + β_: F, + x: &mut DVector, + τ_: F, + σ_: F, + θ_: F, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory, { let λ = λ_.to_nalgebra_mixed(); let τ = τ_.to_nalgebra_mixed(); let σ = σ_.to_nalgebra_mixed(); let θ = θ_.to_nalgebra_mixed(); let β = β_.to_nalgebra_mixed(); - let σθ = σ*θ; - let τθ = τ*θ; + let σθ = σ * θ; + let τθ = τ * θ; let mut w = DVector::zeros(x.len()); let mut tmp = DVector::zeros(x.len()); let mut xprev = x.clone(); @@ -209,8 +224,8 @@ // 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, τ*λ)); - + 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/σ) @@ -221,22 +236,19 @@ 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, β / σθ); 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 function applies an iterative method for the solution of the problem ///
$$ /// \min_{x ∈ ℝ^n} \frac{β}{2} |x-y|_1^2 - g^⊤ x + λ\|x\|₁. @@ -246,16 +258,17 @@ /// This function returns the number of iterations taken. #[replace_float_literals(F::cast_from(literal))] pub fn l1squared_unconstrained( - y : &DVector, - g : &DVector, - λ : F, - β : F, - x : &mut DVector, - inner : &InnerSettings, - iterator : I + y: &DVector, + g: &DVector, + λ: F, + β: F, + x: &mut DVector, + inner: &InnerSettings, + iterator: I, ) -> usize -where F : Float + ToNalgebraRealField, - I : AlgIteratorFactory +where + F: Float + ToNalgebraRealField, + I: AlgIteratorFactory, { // Estimate of ‖K‖ for K=θ Id. let inner_θ = 1.0; @@ -264,8 +277,9 @@ let (inner_τ, inner_σ) = (inner.pdps_τσ0.0 / normest, inner.pdps_τσ0.1 / normest); match inner.method { - InnerMethod::PDPS => - l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator), + InnerMethod::PDPS => { + l1squared_unconstrained_pdps_alt(y, g, λ, β, x, inner_τ, inner_σ, inner_θ, iterator) + } other => unimplemented!("${other:?} is unimplemented"), } }