Sun, 11 Dec 2022 23:25:53 +0200
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
* Fixes the conditional gradient methods that were incorrectly solving
positivity constrained subproblems although the infinite-dimensional
versions do not include such constraints.
* Introduces the `ExperimentV2` struct that has `regularisation` in place
of `α`. The `Experiment` struct is now deprecated.
* The L^2-squared experiments were switch to be unconstrained, as the
Franke-Wolfe implementations do not support constraints. (This would
be easy to add for the “fully corrective” variant, but is not
immediate for the “relaxed” variant.)
24
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
1 | /*! |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
2 | Iterative algorithms for solving the finite-dimensional subproblem without constraints. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
3 | */ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
4 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
5 | use nalgebra::{DVector, DMatrix}; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
6 | use numeric_literals::replace_float_literals; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
7 | use itertools::{izip, Itertools}; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
8 | use colored::Colorize; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
9 | use std::cmp::Ordering::*; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
10 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
11 | use alg_tools::iter::Mappable; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
12 | use alg_tools::error::NumericalError; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
13 | use alg_tools::iterate::{ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
14 | AlgIteratorFactory, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
15 | AlgIteratorState, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
16 | Step, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
17 | }; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
18 | use alg_tools::linops::GEMV; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
19 | use alg_tools::nalgebra_support::ToNalgebraRealField; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
20 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
21 | use crate::types::*; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
22 | use super::{ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
23 | InnerMethod, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
24 | InnerSettings |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
25 | }; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
26 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
27 | /// Compute the proximal operator of $x \mapsto |x|$, i.e., the soft-thresholding operator. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
28 | #[inline] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
29 | #[replace_float_literals(F::cast_from(literal))] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
30 | fn soft_thresholding<F : Float>(v : F, λ : F) -> F { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
31 | if v > λ { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
32 | v - λ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
33 | } else if v < -λ { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
34 | v + λ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
35 | } else { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
36 | 0.0 |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
37 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
38 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
39 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
40 | /// Returns the ∞-norm minimal subdifferential of $x ↦ x^⊤Ax - g^⊤ x + λ\|x\|₁$ at $x$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
41 | /// |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
42 | /// `v` will be modified and cannot be trusted to contain useful values afterwards. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
43 | #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
44 | fn min_subdifferential<F : Float + ToNalgebraRealField>( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
45 | v : &mut DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
46 | mA : &DMatrix<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
47 | x : &DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
48 | g : &DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
49 | λ : F::MixedType |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
50 | ) -> F { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
51 | v.copy_from(g); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
52 | mA.gemv(v, 1.0, x, -1.0); // v = Ax - g |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
53 | let mut val = 0.0; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
54 | for (&v_i, &x_i) in izip!(v.iter(), x.iter()) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
55 | // The subdifferential at x is $Ax - g + λ ∂‖·‖₁(x)$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
56 | val = val.max(match x_i.partial_cmp(&0.0) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
57 | Some(Greater) => v_i + λ, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
58 | Some(Less) => v_i - λ, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
59 | Some(Equal) => soft_thresholding(v_i, λ), |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
60 | None => F::MixedType::nan(), |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
61 | }) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
62 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
63 | F::from_nalgebra_mixed(val) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
64 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
65 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
66 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
67 | /// Forward-backward splitting implementation of [`quadratic_unconstrained`]. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
68 | /// For detailed documentation of the inputs and outputs, refer to there. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
69 | /// |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
70 | /// The `λ` component of the model is handled in the proximal step instead of the gradient step |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
71 | /// for potential performance improvements. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
72 | #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
73 | pub fn quadratic_unconstrained_fb<F, I>( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
74 | mA : &DMatrix<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
75 | g : &DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
76 | //c_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
77 | λ_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
78 | x : &mut DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
79 | τ_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
80 | iterator : I |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
81 | ) -> usize |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
82 | where F : Float + ToNalgebraRealField, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
83 | I : AlgIteratorFactory<F> |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
84 | { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
85 | let mut xprev = x.clone(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
86 | //let c = c_.to_nalgebra_mixed(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
87 | let λ = λ_.to_nalgebra_mixed(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
88 | let τ = τ_.to_nalgebra_mixed(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
89 | let τλ = τ * λ; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
90 | let mut v = DVector::zeros(x.len()); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
91 | let mut iters = 0; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
92 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
93 | iterator.iterate(|state| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
94 | // Replace `x` with $x - τ[Ax-g]= [x + τg]- τAx$ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
95 | v.copy_from(g); // v = g |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
96 | v.axpy(1.0, x, τ); // v = x + τ*g |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
97 | v.sygemv(-τ, mA, x, 1.0); // v = [x + τg]- τAx |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
98 | let backup = state.if_verbose(|| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
99 | xprev.copy_from(x) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
100 | }); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
101 | // Calculate the proximal map |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
102 | x.iter_mut().zip(v.iter()).for_each(|(x_i, &v_i)| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
103 | *x_i = soft_thresholding(v_i, τλ); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
104 | }); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
105 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
106 | iters +=1; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
107 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
108 | backup.map(|_| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
109 | min_subdifferential(&mut v, mA, x, g, λ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
110 | }) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
111 | }); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
112 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
113 | iters |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
114 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
115 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
116 | /// Semismooth Newton implementation of [`quadratic_unconstrained`]. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
117 | /// |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
118 | /// For detailed documentation of the inputs, refer to there. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
119 | /// This function returns the number of iterations taken if there was no inversion failure, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
120 | /// |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
121 | /// For method derivarion, see the documentation for [`super::nonneg::quadratic_nonneg`]. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
122 | #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())] |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
123 | pub fn quadratic_unconstrained_ssn<F, I>( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
124 | mA : &DMatrix<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
125 | g : &DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
126 | //c_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
127 | λ_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
128 | x : &mut DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
129 | τ_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
130 | iterator : I |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
131 | ) -> Result<usize, NumericalError> |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
132 | where F : Float + ToNalgebraRealField, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
133 | I : AlgIteratorFactory<F> |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
134 | { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
135 | let n = x.len(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
136 | let mut xprev = x.clone(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
137 | let mut v = DVector::zeros(n); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
138 | //let c = c_.to_nalgebra_mixed(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
139 | let λ = λ_.to_nalgebra_mixed(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
140 | let τ = τ_.to_nalgebra_mixed(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
141 | let τλ = τ * λ; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
142 | let mut inact : Vec<bool> = Vec::from_iter(std::iter::repeat(false).take(n)); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
143 | let mut s = DVector::zeros(0); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
144 | let mut decomp = nalgebra::linalg::LU::new(DMatrix::zeros(0, 0)); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
145 | let mut iters = 0; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
146 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
147 | let res = iterator.iterate_fallible(|state| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
148 | // 1. Perform delayed SSN-update based on previously computed step on active |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
149 | // coordinates. The step is delayed to the beginning of the loop because |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
150 | // the SSN step may violate constraints, so we arrange `x` to contain at the |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
151 | // end of the loop the valid FB step that forms part of the SSN step |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
152 | let mut si = s.iter(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
153 | for (&ast, x_i, xprev_i) in izip!(inact.iter(), x.iter_mut(), xprev.iter_mut()) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
154 | if ast { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
155 | *x_i = *xprev_i + *si.next().unwrap() |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
156 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
157 | *xprev_i = *x_i; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
158 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
159 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
160 | //xprev.copy_from(x); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
161 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
162 | // 2. Calculate FB step. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
163 | // 2.1. Replace `x` with $x⁻ - τ[Ax⁻-g]= [x⁻ + τg]- τAx⁻$ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
164 | x.axpy(τ, g, 1.0); // x = x⁻ + τ*g |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
165 | x.sygemv(-τ, mA, &xprev, 1.0); // x = [x⁻ + τg]- τAx⁻ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
166 | // 2.2. Calculate prox and set of active coordinates at the same time |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
167 | let mut act_changed = false; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
168 | let mut n_inact = 0; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
169 | for (x_i, ast) in izip!(x.iter_mut(), inact.iter_mut()) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
170 | if *x_i > τλ { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
171 | *x_i -= τλ; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
172 | if !*ast { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
173 | act_changed = true; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
174 | *ast = true; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
175 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
176 | n_inact += 1; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
177 | } else if *x_i < -τλ { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
178 | *x_i += τλ; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
179 | if !*ast { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
180 | act_changed = true; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
181 | *ast = true; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
182 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
183 | n_inact += 1; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
184 | } else { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
185 | *x_i = 0.0; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
186 | if *ast { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
187 | act_changed = true; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
188 | *ast = false; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
189 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
190 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
191 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
192 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
193 | // *** x now contains forward-backward step *** |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
194 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
195 | // 3. Solve SSN step `s`. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
196 | // 3.1 Construct [τ A_{ℐ × ℐ}] if the set of inactive coordinates has changed. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
197 | if act_changed { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
198 | let decomp_iter = inact.iter().cartesian_product(inact.iter()).zip(mA.iter()); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
199 | let decomp_constr = decomp_iter.filter_map(|((&i_inact, &j_inact), &mAij)| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
200 | //(i_inact && j_inact).then_some(mAij * τ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
201 | (i_inact && j_inact).then_some(mAij) // 🔺 below matches removal of τ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
202 | }); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
203 | let mat = DMatrix::from_iterator(n_inact, n_inact, decomp_constr); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
204 | decomp = nalgebra::linalg::LU::new(mat); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
205 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
206 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
207 | // 3.2 Solve `s` = $s_ℐ^k$ from |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
208 | // $[τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜$. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
209 | // With current variable setup we have $[G ∘ F](x^k) = $`x` and $x^k = x⁻$ = `xprev`, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
210 | // so the system to solve is $[τ A_{ℐ × ℐ}]s^k_ℐ = (x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
211 | // The matrix $[τ A_{ℐ × ℐ}]$ we have already LU-decomposed above into `decomp`. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
212 | s = if n_inact > 0 { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
213 | // 3.2.1 Construct `rhs` = $(x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
214 | let inactfilt = inact.iter().copied(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
215 | let rhs_iter = izip!(x.iter(), xprev.iter(), mA.row_iter()).filter_zip(inactfilt); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
216 | let rhs_constr = rhs_iter.map(|(&x_i, &xprev_i, mAi)| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
217 | // Calculate row i of [τ A_{ℐ × 𝒜}]s^k_𝒜 = [τ A_{ℐ × 𝒜}](x-xprev)_𝒜 |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
218 | let actfilt = inact.iter().copied().map(std::ops::Not::not); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
219 | let actit = izip!(x.iter(), xprev.iter(), mAi.iter()).filter_zip(actfilt); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
220 | let actpart = actit.map(|(&x_j, &xprev_j, &mAij)| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
221 | mAij * (x_j - xprev_j) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
222 | }).sum(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
223 | // Subtract it from [x-prev]_i |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
224 | //x_i - xprev_i - τ * actpart |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
225 | (x_i - xprev_i) / τ - actpart // 🔺 change matches removal of τ above |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
226 | }); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
227 | let mut rhs = DVector::from_iterator(n_inact, rhs_constr); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
228 | assert_eq!(rhs.len(), n_inact); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
229 | // Solve the system |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
230 | if !decomp.solve_mut(&mut rhs) { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
231 | return Step::Failure(NumericalError( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
232 | "Failed to solve linear system for subproblem SSN." |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
233 | )) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
234 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
235 | rhs |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
236 | } else { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
237 | DVector::zeros(0) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
238 | }; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
239 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
240 | iters += 1; |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
241 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
242 | // 4. Report solution quality |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
243 | state.if_verbose(|| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
244 | // Calculate subdifferential at the FB step `x` that hasn't yet had `s` yet added. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
245 | min_subdifferential(&mut v, mA, x, g, λ) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
246 | }) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
247 | }); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
248 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
249 | res.map(|_| iters) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
250 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
251 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
252 | /// This function applies an iterative method for the solution of the problem |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
253 | /// <div>$$ |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
254 | /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ\|x\|₁ + c. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
255 | /// $$</div> |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
256 | /// Semismooth Newton or forward-backward are supported based on the setting in `method`. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
257 | /// The parameter `mA` is matrix $A$, and `g` and `λ` are as in the mathematical formulation. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
258 | /// The constant $c$ does not need to be provided. The step length parameter is `τ` while |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
259 | /// `x` contains the initial iterate and on return the final one. The `iterator` controls |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
260 | /// stopping. The “verbose” value output by all methods is the $ℓ\_∞$ distance of some |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
261 | /// subdifferential of the objective to zero. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
262 | /// |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
263 | /// This function returns the number of iterations taken. |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
264 | pub fn quadratic_unconstrained<F, I>( |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
265 | method : InnerMethod, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
266 | mA : &DMatrix<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
267 | g : &DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
268 | //c_ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
269 | λ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
270 | x : &mut DVector<F::MixedType>, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
271 | τ : F, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
272 | iterator : I |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
273 | ) -> usize |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
274 | where F : Float + ToNalgebraRealField, |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
275 | I : AlgIteratorFactory<F> |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
276 | { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
277 | |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
278 | match method { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
279 | InnerMethod::FB => |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
280 | quadratic_unconstrained_fb(mA, g, λ, x, τ, iterator), |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
281 | InnerMethod::SSN => |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
282 | quadratic_unconstrained_ssn(mA, g, λ, x, τ, iterator).unwrap_or_else(|e| { |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
283 | println!("{}", format!("{e}. Using FB fallback.").red()); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
284 | let ins = InnerSettings::<F>::default(); |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
285 | quadratic_unconstrained_fb(mA, g, λ, x, τ, ins.iterator_options) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
286 | }) |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
287 | } |
d29d1fcf5423
Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff
changeset
|
288 | } |