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