src/subproblem.rs

Tue, 06 Dec 2022 14:12:20 +0200

author
Tuomo Valkonen <tuomov@iki.fi>
date
Tue, 06 Dec 2022 14:12:20 +0200
changeset 10
b71edfd403aa
parent 0
eb3c7813b67a
child 2
7a953a87b6c1
child 4
5aa5c279e341
child 24
d29d1fcf5423
permissions
-rw-r--r--

v1.0.0-pre-arxiv (missing arXiv links)

0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
1 //! Iterative algorithms for solving finite-dimensional subproblems.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
2
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
3 use serde::{Serialize, Deserialize};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
4 use nalgebra::{DVector, DMatrix};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
5 use numeric_literals::replace_float_literals;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
6 use itertools::{izip, Itertools};
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
7 use colored::Colorize;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
8
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
9 use alg_tools::iter::Mappable;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
10 use alg_tools::error::NumericalError;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
11 use alg_tools::iterate::{
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
12 AlgIteratorFactory,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
13 AlgIteratorState,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
14 AlgIteratorOptions,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
15 Verbose,
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::*;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
22
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
23 /// Method for solving finite-dimensional subproblems
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
24 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
25 #[allow(dead_code)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
26 pub enum InnerMethod {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
27 /// Forward-backward
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
28 FB,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
29 /// Semismooth Newton
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
30 SSN,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
31 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
32
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
33 /// Settings for the solution of finite-dimensional subproblems
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
34 #[derive(Clone, Copy, Eq, PartialEq, Serialize, Deserialize, Debug)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
35 pub struct InnerSettings<F : Float> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
36 /// Method
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
37 pub method : InnerMethod,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
38 /// Proportional step length (∈ [0, 1) for `InnerMethod::FB`).
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
39 pub τ0 : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
40 /// Fraction of `tolerance` given to inner algorithm
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
41 pub tolerance_mult : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
42 /// Iterator options
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
43 #[serde(flatten)]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
44 pub iterator_options : AlgIteratorOptions,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
45 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
46
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
47 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
48 impl<F : Float> Default for InnerSettings<F> {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
49 fn default() -> Self {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
50 InnerSettings {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
51 τ0 : 0.99,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
52 iterator_options : AlgIteratorOptions {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
53 // max_iter cannot be very small, as initially FB needs many iterations, although
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
54 // on later invocations even one or two tends to be enough
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
55 max_iter : 2000,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
56 // verbose_iter affects testing of sufficient convergence, so we set it to
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
57 // a small value…
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
58 verbose_iter : Verbose::Every(1),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
59 // … but don't print out anything
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
60 quiet : true,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
61 .. Default::default()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
62 },
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
63 method : InnerMethod::FB,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
64 tolerance_mult : 0.01,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
65 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
66 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
67 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
68
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
69 /// Compute the proximal operator of $x \mapsto x + \delta\_{[0, \infty)}$, i.e.,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
70 /// the non-negativity contrained soft-thresholding operator.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
71 #[inline]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
72 #[replace_float_literals(F::cast_from(literal))]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
73 fn nonneg_soft_thresholding<F : Float>(v : F, λ : F) -> F {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
74 (v - λ).max(0.0)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
75 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
76
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
77 /// Forward-backward splitting implementation of [`quadratic_nonneg`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
78 /// For detailed documentation of the inputs and outputs, refer to there.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
79 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
80 /// 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
81 /// for potential performance improvements.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
82 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
83 pub fn quadratic_nonneg_fb<F, I>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
84 mA : &DMatrix<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
85 g : &DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
86 //c_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
87 λ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
88 x : &mut DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
89 τ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
90 iterator : I
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
91 ) -> usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
92 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
93 I : AlgIteratorFactory<F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
94 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
95 let mut xprev = x.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
96 //let c = c_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
97 let λ = λ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
98 let τ = τ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
99 let τλ = τ * λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
100 let mut v = DVector::zeros(x.len());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
101 let mut iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
102
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
103 iterator.iterate(|state| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
104 // Replace `x` with $x - τ[Ax-g]= [x + τg]- τAx$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
105 v.copy_from(g); // v = g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
106 v.axpy(1.0, x, τ); // v = x + τ*g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
107 v.sygemv(-τ, mA, x, 1.0); // v = [x + τg]- τAx
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
108 let backup = state.if_verbose(|| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
109 xprev.copy_from(x)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
110 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
111 // Calculate the proximal map
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
112 x.iter_mut().zip(v.iter()).for_each(|(x_i, &v_i)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
113 *x_i = nonneg_soft_thresholding(v_i, τλ);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
114 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
115
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
116 iters +=1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
117
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
118 backup.map(|_| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
119 // The subdifferential of the objective is $Ax - g + λ + ∂ δ_{≥ 0}(x)$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
120 // We return the minimal ∞-norm over all subderivatives.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
121 v.copy_from(g); // d = g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
122 mA.gemv(&mut v, 1.0, x, -1.0); // d = Ax - g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
123 let mut val = 0.0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
124 for (&v_i, &x_i) in izip!(v.iter(), x.iter()) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
125 let d = v_i + λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
126 if x_i > 0.0 || d < 0.0 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
127 val = val.max(d.abs());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
128 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
129 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
130 F::from_nalgebra_mixed(val)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
131 })
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
132 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
133
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
134 iters
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
135 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
136
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
137 /// Semismooth Newton implementation of [`quadratic_nonneg`].
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
138 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
139 /// For detailed documentation of the inputs, refer to there.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
140 /// 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
141 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
142 /// ## Method derivation
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
143 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
144 /// **The below may look like garbage. Sorry, but rustdoc is obsolete rubbish
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
145 /// 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
146 /// forces one into unreliable KaTeX autorender postprocessing andescape hell and that
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
147 /// it doesn't even process correctly.**
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
148 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
149 /// <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
150 /// For the objective
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
151 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
152 /// 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
153 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
154 /// we have the optimality condition
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
155 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
156 /// x - \mathop{\mathrm{prox}}_{τλ{\vec 1}^⊤ + δ_{≥ 0}}(x - τ[Ax-g^⊤]) = 0,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
157 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
158 /// which we write as
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
159 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
160 /// x - [G ∘ F](x)=0
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
161 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
162 /// for
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
163 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
164 /// G(x) = \mathop{\mathrm{prox}}_{λ{\vec 1}^⊤ + δ_{≥ 0}}
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
165 /// \quad\text{and}\quad
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
166 /// F(x) = x - τ Ax + τ g^⊤
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
167 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
168 /// We can use Newton derivative chain rule to compute
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
169 /// $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
170 /// $D_N F(x) = \mathop{\mathrm{Id}} - τ A$,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
171 /// 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
172 /// </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
173 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
174 /// <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
175 /// 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
176 /// updating $x^{k+1} = x^k + s^k$. Consequently
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
177 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
178 /// 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
179 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
180 /// For $𝒜$ the set of active coordinates and $ℐ$ the set of inactive coordinates, this
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
181 /// expands as
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
182 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
183 /// [τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
184 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
185 /// and
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
186 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
187 /// s^k_𝒜 = - x^k_𝒜 + [G ∘ F](x^k)_𝒜.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
188 /// $$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
189 /// 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
190 /// </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
191 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
192 /// <p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
193 /// 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
194 /// 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
195 /// forward-backward step.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
196 /// </p>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
197 #[replace_float_literals(F::cast_from(literal).to_nalgebra_mixed())]
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
198 pub fn quadratic_nonneg_ssn<F, I>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
199 mA : &DMatrix<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
200 g : &DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
201 //c_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
202 λ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
203 x : &mut DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
204 τ_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
205 iterator : I
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
206 ) -> Result<usize, NumericalError>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
207 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
208 I : AlgIteratorFactory<F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
209 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
210 let n = x.len();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
211 let mut xprev = x.clone();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
212 let mut v = DVector::zeros(n);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
213 //let c = c_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
214 let λ = λ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
215 let τ = τ_.to_nalgebra_mixed();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
216 let τλ = τ * λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
217 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
218 let mut s = DVector::zeros(0);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
219 let mut decomp = nalgebra::linalg::LU::new(DMatrix::zeros(0, 0));
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
220 let mut iters = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
221
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
222 let res = iterator.iterate_fallible(|state| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
223 // 1. Perform delayed SSN-update based on previously computed step on active
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
224 // coordinates. The step is delayed to the beginning of the loop because
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
225 // 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
226 // 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
227 let mut si = s.iter();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
228 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
229 if ast {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
230 *x_i = *xprev_i + *si.next().unwrap()
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
231 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
232 *xprev_i = *x_i;
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 //xprev.copy_from(x);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
236
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
237 // 2. Calculate FB step.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
238 // 2.1. Replace `x` with $x⁻ - τ[Ax⁻-g]= [x⁻ + τg]- τAx⁻$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
239 x.axpy(τ, g, 1.0); // x = x⁻ + τ*g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
240 x.sygemv(-τ, mA, &xprev, 1.0); // x = [x⁻ + τg]- τAx⁻
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
241 // 2.2. Calculate prox and set of active coordinates at the same time
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
242 let mut act_changed = false;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
243 let mut n_inact = 0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
244 for (x_i, ast) in izip!(x.iter_mut(), inact.iter_mut()) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
245 if *x_i > τλ {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
246 *x_i -= τλ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
247 if !*ast {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
248 act_changed = true;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
249 *ast = true;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
250 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
251 n_inact += 1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
252 } else {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
253 *x_i = 0.0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
254 if *ast {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
255 act_changed = true;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
256 *ast = false;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
257 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
258 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
259 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
260
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
261 // *** x now contains forward-backward step ***
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
262
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
263 // 3. Solve SSN step `s`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
264 // 3.1 Construct [τ A_{ℐ × ℐ}] if the set of inactive coordinates has changed.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
265 if act_changed {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
266 let decomp_iter = inact.iter().cartesian_product(inact.iter()).zip(mA.iter());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
267 let decomp_constr = decomp_iter.filter_map(|((&i_inact, &j_inact), &mAij)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
268 //(i_inact && j_inact).then_some(mAij * τ)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
269 (i_inact && j_inact).then_some(mAij) // 🔺 below matches removal of τ
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 mat = DMatrix::from_iterator(n_inact, n_inact, decomp_constr);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
272 decomp = nalgebra::linalg::LU::new(mat);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
273 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
274
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
275 // 3.2 Solve `s` = $s_ℐ^k$ from
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
276 // $[τ A_{ℐ × ℐ}]s^k_ℐ = - x^k_ℐ + [G ∘ F](x^k)_ℐ - [τ A_{ℐ × 𝒜}]s^k_𝒜$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
277 // 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
278 // 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
279 // The matrix $[τ A_{ℐ × ℐ}]$ we have already LU-decomposed above into `decomp`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
280 s = if n_inact > 0 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
281 // 3.2.1 Construct `rhs` = $(x-x⁻)_ℐ - [τ A_{ℐ × 𝒜}](x-x⁻)_𝒜$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
282 let inactfilt = inact.iter().copied();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
283 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
284 let rhs_constr = rhs_iter.map(|(&x_i, &xprev_i, mAi)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
285 // Calculate row i of [τ A_{ℐ × 𝒜}]s^k_𝒜 = [τ A_{ℐ × 𝒜}](x-xprev)_𝒜
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
286 let actfilt = inact.iter().copied().map(std::ops::Not::not);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
287 let actit = izip!(x.iter(), xprev.iter(), mAi.iter()).filter_zip(actfilt);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
288 let actpart = actit.map(|(&x_j, &xprev_j, &mAij)| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
289 mAij * (x_j - xprev_j)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
290 }).sum();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
291 // Subtract it from [x-prev]_i
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
292 //x_i - xprev_i - τ * actpart
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
293 (x_i - xprev_i) / τ - actpart // 🔺 change matches removal of τ above
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
294 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
295 let mut rhs = DVector::from_iterator(n_inact, rhs_constr);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
296 assert_eq!(rhs.len(), n_inact);
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
297 // Solve the system
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
298 if !decomp.solve_mut(&mut rhs) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
299 return Step::Failure(NumericalError(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
300 "Failed to solve linear system for subproblem SSN."
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
301 ))
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
302 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
303 rhs
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
304 } else {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
305 DVector::zeros(0)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
306 };
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
307
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
308 iters += 1;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
309
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
310 // 4. Report solution quality
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
311 state.if_verbose(|| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
312 // Calculate subdifferential at the FB step `x` that hasn't yet had `s` yet added.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
313 // The subdifferential of the objective is $Ax - g + λ + ∂ δ_{≥ 0}(x)$.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
314 // We return the minimal ∞-norm over all subderivatives.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
315 v.copy_from(g); // d = g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
316 mA.gemv(&mut v, 1.0, x, -1.0); // d = Ax - g
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
317 let mut val = 0.0;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
318 for (&v_i, &x_i) in izip!(v.iter(), x.iter()) {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
319 let d = v_i + λ;
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
320 if x_i > 0.0 || d < 0.0 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
321 val = val.max(d.abs());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
322 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
323 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
324 F::from_nalgebra_mixed(val)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
325 })
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
326 });
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
327
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
328 res.map(|_| iters)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
329 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
330
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
331 /// 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
332 /// constrained problem
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
333 /// <div>$$
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
334 /// \min_{x ∈ ℝ^n} \frac{1}{2} x^⊤Ax - g^⊤ x + λ{\vec 1}^⊤ x + c + δ_{≥ 0}(x).
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
335 /// $$</div>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
336 /// Semismooth Newton or forward-backward are supported based on the setting in `method`.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
337 /// 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
338 /// 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
339 /// `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
340 /// stopping. The “verbose” value output by all methods is the $ℓ\_∞$ distance of some
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
341 /// subdifferential of the objective to zero.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
342 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
343 /// Interior point methods could offer a further alternative, for example, the one in:
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
344 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
345 /// * Valkonen T. - _A method for weighted projections to the positive definite
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
346 /// cone_, <https://doi.org/10.1080/02331934.2014.929680>.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
347 ///
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
348 /// This function returns the number of iterations taken.
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
349 pub fn quadratic_nonneg<F, I>(
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
350 method : InnerMethod,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
351 mA : &DMatrix<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
352 g : &DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
353 //c_ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
354 λ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
355 x : &mut DVector<F::MixedType>,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
356 τ : F,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
357 iterator : I
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
358 ) -> usize
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
359 where F : Float + ToNalgebraRealField,
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
360 I : AlgIteratorFactory<F>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
361 {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
362
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
363 match method {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
364 InnerMethod::FB =>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
365 quadratic_nonneg_fb(mA, g, λ, x, τ, iterator),
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
366 InnerMethod::SSN =>
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
367 quadratic_nonneg_ssn(mA, g, λ, x, τ, iterator).unwrap_or_else(|e| {
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
368 println!("{}", format!("{e}. Using FB fallback.").red());
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
369 let ins = InnerSettings::<F>::default();
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
370 quadratic_nonneg_fb(mA, g, λ, x, τ, ins.iterator_options)
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
371 })
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
372 }
eb3c7813b67a Initial version
Tuomo Valkonen <tuomov@iki.fi>
parents:
diff changeset
373 }

mercurial