src/subproblem/nonneg.rs

Mon, 05 Dec 2022 23:50:22 +0200

author
Tuomo Valkonen <tuomov@iki.fi>
date
Mon, 05 Dec 2022 23:50:22 +0200
changeset 56
3a784e6e475a
parent 39
6316d68b58af
permissions
-rw-r--r--

Zenodo packaging hacks

24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
1 /*!
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
2 Iterative algorithms for solving the finite-dimensional subproblem with non-negativity constraints.
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
3 */
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
4
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
5 use nalgebra::{DVector, DMatrix};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
6 use numeric_literals::replace_float_literals;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
7 use itertools::{izip, Itertools};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
8 use colored::Colorize;
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
9 use std::cmp::Ordering::*;
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
10
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
11 use alg_tools::iter::Mappable;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
12 use alg_tools::error::NumericalError;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
13 use alg_tools::iterate::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
14 AlgIteratorFactory,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
15 AlgIteratorState,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
16 Step,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
17 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
18 use alg_tools::linops::GEMV;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
19 use alg_tools::nalgebra_support::ToNalgebraRealField;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
20
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
21 use crate::types::*;
24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
22 use super::{
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
23 InnerMethod,
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
24 InnerSettings
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
25 };
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
26 use super::l1squared_nonneg::max_interval_dist_to_zero;
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
27
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
28 /// Compute the proximal operator of $x \mapsto x + \delta\_{[0, \infty)}$, i.e.,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
29 /// the non-negativity contrained soft-thresholding operator.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
30 #[inline]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
31 #[replace_float_literals(F::cast_from(literal))]
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
32 pub(super) fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F {
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
33 (v - λ).max(0.0)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
34 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
35
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
36 /// Returns the ∞-norm minimal subdifferential of $x ↦ x^⊤Ax/2 - g^⊤ x + λ\vec 1^⊤ x + δ_{≥ 0}(x)$
24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
37 /// at $x$.
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
38 ///
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
39 /// `v` will be modified and cannot be trusted to contain useful values afterwards.
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
40 #[replace_float_literals(F::cast_from(literal))]
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
41 fn min_subdifferential<F : Float + nalgebra::RealField>(
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
42 v : &mut DVector<F>,
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
43 mA : &DMatrix<F>,
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
44 x : &DVector<F>,
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
45 g : &DVector<F>,
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
46 λ : F
24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
47 ) -> F {
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
48 v.copy_from(g);
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
49 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: 0
diff changeset
50 let mut val = 0.0;
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
51 for (&v_i, &x_i) in izip!(v.iter(), x.iter()) {
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
52 let (mut lb, mut ub) = (v_i, v_i);
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
53 match x_i.partial_cmp(&0.0) {
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
54 Some(Greater) => { lb += λ; ub += λ },
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
55 // Less should not happen
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
56 Some(Less|Equal) => { lb = F::NEG_INFINITY; ub += λ },
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
57 None => {},
24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
58 }
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
59 val = max_interval_dist_to_zero(val, lb, ub);
24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
60 }
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
61 val
24
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
62 }
d29d1fcf5423 Support arbitrary regularisation terms; implement non-positivity-constrained regularisation.
Tuomo Valkonen <tuomov@iki.fi>
parents: 0
diff changeset
63
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
64 /// Forward-backward splitting implementation of [`quadratic_nonneg`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
65 /// For detailed documentation of the inputs and outputs, refer to there.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
66 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
67 /// The `λ` component of the model is handled in the proximal step instead of the gradient step
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
68 /// for potential performance improvements.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
69 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
70 pub fn quadratic_nonneg_fb<F, I>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
71 mA : &DMatrix<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
72 g : &DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
73 //c_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
74 λ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
75 x : &mut DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
76 τ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
77 iterator : I
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
78 ) -> usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
79 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
80 I : AlgIteratorFactory<F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
81 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
82 let mut xprev = x.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
83 //let c = c_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
84 let λ = λ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
85 let τ = τ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
86 let τλ = τ * λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
87 let mut v = DVector::zeros(x.len());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
88 let mut iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
89
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
90 iterator.iterate(|state| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
91 // Replace `x` with $x - τ[Ax-g]= [x + τg]- τAx$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
92 v.copy_from(g); // v = g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
93 v.axpy(1.0, x, τ); // v = x + τ*g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
94 v.sygemv(-τ, mA, x, 1.0); // v = [x + τg]- τAx
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
95 let backup = state.if_verbose(|| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
96 xprev.copy_from(x)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
97 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
98 // Calculate the proximal map
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
99 x.iter_mut().zip(v.iter()).for_each(|(x_i, &v_i)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
100 *x_i = nonneg_soft_thresholding(v_i, τλ);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
101 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
102
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
103 iters +=1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
104
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
105 backup.map(|_| {
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
106 F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
107 })
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
108 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
109
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
110 iters
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
111 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
112
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
113 /// Semismooth Newton implementation of [`quadratic_nonneg`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
114 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
115 /// For detailed documentation of the inputs, refer to there.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
116 /// This function returns the number of iterations taken if there was no inversion failure,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
117 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
118 /// ## Method derivation
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
119 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
120 /// **The below may look like garbage. Sorry, but rustdoc is obsolete rubbish
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
121 /// that doesn't directly support by-now standard-in-markdown LaTeX math. Instead it
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
122 /// forces one into unreliable KaTeX autorender postprocessing andescape hell and that
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
123 /// it doesn't even process correctly.**
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
124 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
125 /// <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
126 /// For the objective
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
127 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
128 /// J(x) = \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + c + δ_{≥ 0}(x),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
129 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
130 /// we have the optimality condition
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
131 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
132 /// x - \mathop{\mathrm{prox}}_{τλ{\vec 1}^⊤ + δ_{≥ 0}}(x - τ[Ax-g^⊤]) = 0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
133 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
134 /// which we write as
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
135 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
136 /// x - [G ∘ F](x)=0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
137 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
138 /// for
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
139 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
140 /// G(x) = \mathop{\mathrm{prox}}_{λ{\vec 1}^⊤ + δ_{≥ 0}}
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
141 /// \quad\text{and}\quad
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
142 /// F(x) = x - τ Ax + τ g^⊤
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
143 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
144 /// We can use Newton derivative chain rule to compute
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
145 /// $D_N[G ∘ F](x) = D_N G(F(x)) D_N F(x)$, where
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
146 /// $D_N F(x) = \mathop{\mathrm{Id}} - τ A$,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
147 /// and $[D_N G(F(x))]_i = 1$ for inactive coordinates and $=0$ for active coordinates.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
148 /// </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
149 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
150 /// <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
151 /// The method itself involves solving $D_N[Id - G ∘ F](x^k) s^k = - [Id - G ∘ F](x^k)$ and
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
152 /// updating $x^{k+1} = x^k + s^k$. Consequently
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
153 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
154 /// s^k - D_N G(F(x^k)) [s^k - τ As^k] = - x^k + [G ∘ F](x^k)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
155 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
156 /// For $𝒜$ the set of active coordinates and $ℐ$ the set of inactive coordinates, this
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
157 /// expands as
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
158 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
159 /// [τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
160 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
161 /// and
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
162 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
163 /// s^k_𝒜 = - x^k_𝒜 + [G ∘ F](x^k)_𝒜.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
164 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
165 /// Thus on $𝒜$ the update $[x^k + s^k]_𝒜 = [G ∘ F](x^k)_𝒜$ is just the forward-backward update.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
166 /// </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
167 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
168 /// <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
169 /// We need to detect stopping by a subdifferential and return $x$ satisfying $x ≥ 0$,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
170 /// which is in general not true for the SSN. We therefore use that $[G ∘ F](x^k)$ is a valid
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
171 /// forward-backward step.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
172 /// </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
173 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
174 pub fn quadratic_nonneg_ssn<F, I>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
175 mA : &DMatrix<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
176 g : &DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
177 //c_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
178 λ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
179 x : &mut DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
180 τ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
181 iterator : I
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
182 ) -> Result<usize, NumericalError>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
183 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
184 I : AlgIteratorFactory<F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
185 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
186 let n = x.len();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
187 let mut xprev = x.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
188 let mut v = DVector::zeros(n);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
189 //let c = c_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
190 let λ = λ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
191 let τ = τ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
192 let τλ = τ * λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
193 let mut inact : Vec<bool> = Vec::from_iter(std::iter::repeat(false).take(n));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
194 let mut s = DVector::zeros(0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
195 let mut decomp = nalgebra::linalg::LU::new(DMatrix::zeros(0, 0));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
196 let mut iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
197
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
198 let res = iterator.iterate_fallible(|state| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
199 // 1. Perform delayed SSN-update based on previously computed step on active
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
200 // coordinates. The step is delayed to the beginning of the loop because
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
201 // the SSN step may violate constraints, so we arrange `x` to contain at the
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
202 // end of the loop the valid FB step that forms part of the SSN step
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
203 let mut si = s.iter();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
204 for (&ast, x_i, xprev_i) in izip!(inact.iter(), x.iter_mut(), xprev.iter_mut()) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
205 if ast {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
206 *x_i = *xprev_i + *si.next().unwrap()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
207 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
208 *xprev_i = *x_i;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
209 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
210
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
211 //xprev.copy_from(x);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
212
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
213 // 2. Calculate FB step.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
214 // 2.1. Replace `x` with $x⁻ - τ[Ax⁻-g]= [x⁻ + τg]- τAx⁻$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
215 x.axpy(τ, g, 1.0); // x = x⁻ + τ*g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
216 x.sygemv(-τ, mA, &xprev, 1.0); // x = [x⁻ + τg]- τAx⁻
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
217 // 2.2. Calculate prox and set of active coordinates at the same time
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
218 let mut act_changed = false;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
219 let mut n_inact = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
220 for (x_i, ast) in izip!(x.iter_mut(), inact.iter_mut()) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
221 if *x_i > τλ {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
222 *x_i -= τλ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
223 if !*ast {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
224 act_changed = true;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
225 *ast = true;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
226 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
227 n_inact += 1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
228 } else {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
229 *x_i = 0.0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
230 if *ast {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
231 act_changed = true;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
232 *ast = false;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
233 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
234 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
235 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
236
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
237 // *** x now contains forward-backward step ***
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
238
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
239 // 3. Solve SSN step `s`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
240 // 3.1 Construct [τ A_{ℐ × ℐ}] if the set of inactive coordinates has changed.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
241 if act_changed {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
242 let decomp_iter = inact.iter().cartesian_product(inact.iter()).zip(mA.iter());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
243 let decomp_constr = decomp_iter.filter_map(|((&i_inact, &j_inact), &mAij)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
244 //(i_inact && j_inact).then_some(mAij * τ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
245 (i_inact && j_inact).then_some(mAij) // 🔺 below matches removal of τ
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
246 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
247 let mat = DMatrix::from_iterator(n_inact, n_inact, decomp_constr);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
248 decomp = nalgebra::linalg::LU::new(mat);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
249 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
250
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
251 // 3.2 Solve `s` = $s_ℐ^k$ from
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
252 // $[τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
253 // With current variable setup we have $[G ∘ F](x^k) = $`x` and $x^k = x⁻$ = `xprev`,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
254 // so the system to solve is $[τ A_{ℐ × ℐ}]s^k_ℐ = (x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
255 // The matrix $[τ A_{ℐ × ℐ}]$ we have already LU-decomposed above into `decomp`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
256 s = if n_inact > 0 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
257 // 3.2.1 Construct `rhs` = $(x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
258 let inactfilt = inact.iter().copied();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
259 let rhs_iter = izip!(x.iter(), xprev.iter(), mA.row_iter()).filter_zip(inactfilt);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
260 let rhs_constr = rhs_iter.map(|(&x_i, &xprev_i, mAi)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
261 // Calculate row i of [τ A_{ℐ × 𝒜}]s^k_𝒜 = [τ A_{ℐ × 𝒜}](x-xprev)_𝒜
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
262 let actfilt = inact.iter().copied().map(std::ops::Not::not);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
263 let actit = izip!(x.iter(), xprev.iter(), mAi.iter()).filter_zip(actfilt);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
264 let actpart = actit.map(|(&x_j, &xprev_j, &mAij)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
265 mAij * (x_j - xprev_j)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
266 }).sum();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
267 // Subtract it from [x-prev]_i
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
268 //x_i - xprev_i - τ * actpart
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
269 (x_i - xprev_i) / τ - actpart // 🔺 change matches removal of τ above
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
270 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
271 let mut rhs = DVector::from_iterator(n_inact, rhs_constr);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
272 assert_eq!(rhs.len(), n_inact);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
273 // Solve the system
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
274 if !decomp.solve_mut(&mut rhs) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
275 return Step::Failure(NumericalError(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
276 "Failed to solve linear system for subproblem SSN."
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
277 ))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
278 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
279 rhs
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
280 } else {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
281 DVector::zeros(0)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
282 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
283
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
284 iters += 1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
285
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
286 // 4. Report solution quality
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
287 state.if_verbose(|| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
288 // Calculate subdifferential at the FB step `x` that hasn't yet had `s` yet added.
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
289 F::from_nalgebra_mixed(min_subdifferential(&mut v, mA, x, g, λ))
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
290 })
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
291 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
292
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
293 res.map(|_| iters)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
294 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
295
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
296 /// This function applies an iterative method for the solution of the quadratic non-negativity
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
297 /// constrained problem
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
298 /// <div>$$
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
299 /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + δ_{≥ 0}(x).
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
300 /// $$</div>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
301 /// Semismooth Newton or forward-backward are supported based on the setting in `method`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
302 /// The parameter `mA` is matrix $A$, and `g` and `λ` are as in the mathematical formulation.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
303 /// The constant $c$ does not need to be provided. The step length parameter is `τ` while
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
304 /// `x` contains the initial iterate and on return the final one. The `iterator` controls
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
305 /// stopping. The “verbose” value output by all methods is the $ℓ\_∞$ distance of some
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
306 /// subdifferential of the objective to zero.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
307 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
308 /// Interior point methods could offer a further alternative, for example, the one in:
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
309 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
310 /// * Valkonen T. - _A method for weighted projections to the positive definite
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
311 /// cone_, <https://doi.org/10.1080/02331934.2014.929680>.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
312 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
313 /// This function returns the number of iterations taken.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
314 pub fn quadratic_nonneg<F, I>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
315 mA : &DMatrix<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
316 g : &DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
317 λ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
318 x : &mut DVector<F::MixedType>,
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
319 mA_normest : F,
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
320 inner : &InnerSettings<F>,
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
321 iterator : I
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
322 ) -> usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
323 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
324 I : AlgIteratorFactory<F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
325 {
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
326 let inner_τ = inner.fb_τ0 / mA_normest;
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
327
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
328 match inner.method {
39
6316d68b58af Merging adjustments, parameter tuning, etc.
Tuomo Valkonen <tuomov@iki.fi>
parents: 34
diff changeset
329 InnerMethod::FB =>
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
330 quadratic_nonneg_fb(mA, g, λ, x, inner_τ, iterator),
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
331 InnerMethod::SSN =>
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
332 quadratic_nonneg_ssn(mA, g, λ, x, inner_τ, iterator).unwrap_or_else(|e| {
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
333 println!("{}", format!("{e}. Using FB fallback.").red());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
334 let ins = InnerSettings::<F>::default();
34
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
335 quadratic_nonneg_fb(mA, g, λ, x, inner_τ, ins.iterator_options)
efa60bc4f743 Radon FB + sliding improvements
Tuomo Valkonen <tuomov@iki.fi>
parents: 24
diff changeset
336 }),
39
6316d68b58af Merging adjustments, parameter tuning, etc.
Tuomo Valkonen <tuomov@iki.fi>
parents: 34
diff changeset
337 other => unimplemented!("${other:?} is unimplemented"),
0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
338 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
339 }

mercurial