Mercurial > repos > pointsource_algs / file comparison
comparison: src/subproblem/l1squared_nonneg.rs
src/subproblem/l1squared_nonneg.rs
- branch
- dev
- changeset 64
- d3be4f7ffd60
- parent 63
- 7a8a55fd41c0
equal
deleted
inserted
replaced
| 63:7a8a55fd41c0 | 64:d3be4f7ffd60 |
|---|---|
| 98 /// \prox_f(x) = \prox\_{f\_0}(x - y) + y, | 98 /// \prox_f(x) = \prox\_{f\_0}(x - y) + y, |
| 99 /// $$ | 99 /// $$ |
| 100 /// which is also used internally. | 100 /// which is also used internally. |
| 101 /// The third output indicates whether the output is “locked” to either $0$ (resp. $-y$ in | 101 /// The third output indicates whether the output is “locked” to either $0$ (resp. $-y$ in |
| 102 /// the reformulation), or $y$ ($0$ in the reformulation). | 102 /// the reformulation), or $y$ ($0$ in the reformulation). |
| 103 /// The fourth output indicates the next highest $λ$ where such locking would happen. | 103 /// The fourth output indicates $λ$ where such locking to y would happen. It may be negative, |
| 104 /// in which case the prox is always locked to y. | |
| 104 #[replace_float_literals(F::cast_from(literal))] | 105 #[replace_float_literals(F::cast_from(literal))] |
| 105 #[inline] | 106 #[inline] |
| 106 fn shifted_nonneg_soft_thresholding<F: Float>(x: F, y: F, λ: F) -> (F, F, bool, Option<F>) { | 107 fn shifted_nonneg_soft_thresholding<F: Float>(x: F, y: F, λ: F) -> (F, F, bool, F) { |
| 108 // The OC is | |
| 109 // x ∈ w + λsign(w-y) + N_{≥ 0}(w) | |
| 110 // x-y ∈ w' + λsign(w') + N_{≥ -y}(w') | |
| 111 // if w' > -y, then soft-thresholding. So if soft-thresholding stays above -y. | |
| 112 // If w' = -y, and w'≠0, then x ≤ λsign(-y) | |
| 113 // If w' = -y and w'=0, then x ≤ λ. | |
| 114 // w_i=max{sign(x-y)*max(|x-y| - λ, 0), -y} + y = max{sign(x-y)*max(|x-y| - λ, 0)+y, 0} | |
| 115 let (w, locked, lock) = if x >= y { | |
| 116 // w_i=max{max(x-y - λ, 0)+y, 0}=max{max(x - λ, y), 0} = max{x-λ, max(y, 0)} | |
| 117 let l = 0.0.max(y); | |
| 118 // x - l = x + min(0.0, y) = min(x, x + y) | |
| 119 if x - λ >= l { | |
| 120 (x - λ, false, x - l) | |
| 121 } else { | |
| 122 (l, true, F::NEG_INFINITY) // λ is increasing, so NEG_INFINITY is ok. So is x-l < λ | |
| 123 } | |
| 124 } else { | |
| 125 // w_i=max{-max(-(x-y) - λ, 0)+y, 0}=max{min((x-y) + λ, 0)+y, 0}=max{min(x + λ, y), 0} | |
| 126 if x + λ <= y { | |
| 127 // = max(x+λ, 0) | |
| 128 if x >= -λ { | |
| 129 (x + λ, false, y - x) | |
| 130 } else { | |
| 131 (0.0, true, 0.0.min(y) - x) | |
| 132 } | |
| 133 } else { | |
| 134 (0.0.max(y), true, F::NEG_INFINITY) | |
| 135 } | |
| 136 }; | |
| 137 return (w, w - y, locked, lock); | |
| 138 | |
| 139 /* | |
| 107 let z = x - y; // The shift to f_0 | 140 let z = x - y; // The shift to f_0 |
| 108 if -y >= 0.0 { | 141 if -y >= 0.0 { |
| 109 if x > λ { | 142 if x > λ { |
| 110 (x - λ, z - λ, false, Some(x)) | 143 (x - λ, z - λ, false, Some(x)) |
| 111 } else { | 144 } else { |
| 121 } | 154 } |
| 122 } else if z > λ { | 155 } else if z > λ { |
| 123 (x - λ, z - λ, false, Some(z)) | 156 (x - λ, z - λ, false, Some(z)) |
| 124 } else { | 157 } else { |
| 125 (y, 0.0, true, None) | 158 (y, 0.0, true, None) |
| 126 } | 159 }*/ |
| 127 } | 160 } |
| 128 | 161 |
| 129 /// Calculate $\prox\_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}\_1\^2 + δ\_{≥0}(x)$. | 162 /// Calculate $\prox\_f(x)$ for $f(x)=\frac{β}{2}\norm{x-y}\_1\^2 + δ\_{≥0}(x)$. |
| 130 /// | 163 /// |
| 131 /// To derive an algorithm for this, we can use | 164 /// To derive an algorithm for this, we can use |
| 178 let mut w_locked = 0.0; | 211 let mut w_locked = 0.0; |
| 179 let mut n_unlocked = 0; // m | 212 let mut n_unlocked = 0; // m |
| 180 let mut x_prime = 0.0; | 213 let mut x_prime = 0.0; |
| 181 let mut max_shift = F::INFINITY; | 214 let mut max_shift = F::INFINITY; |
| 182 for (&x_i, &y_i) in izip!(x.iter(), y.iter()) { | 215 for (&x_i, &y_i) in izip!(x.iter(), y.iter()) { |
| 183 let (_, a_shifted, locked, next_lock) = shifted_nonneg_soft_thresholding(x_i, y_i, λ); | 216 let (_, w_i, locked, lock) = shifted_nonneg_soft_thresholding(x_i, y_i, λ); |
| 184 | 217 |
| 185 if let Some(t) = next_lock { | 218 // If not already locked (lock <= y), consider this coordinate in next level for λ. |
| 186 max_shift = min(max_shift, t); | 219 if lock > λ { |
| 187 assert!(max_shift > λ); | 220 max_shift = min(max_shift, lock); |
| 188 } | 221 } |
| 222 | |
| 189 if locked { | 223 if locked { |
| 190 w_locked += abs(a_shifted); | 224 w_locked += abs(w_i); |
| 191 } else { | 225 } else { |
| 192 n_unlocked += 1; | 226 n_unlocked += 1; |
| 193 x_prime += abs(x_i - y_i); | 227 x_prime += abs(x_i - y_i); |
| 194 } | 228 } |
| 195 } | 229 } |
| 201 λ = max_shift; | 235 λ = max_shift; |
| 202 } else { | 236 } else { |
| 203 assert!(λ_new >= λ); | 237 assert!(λ_new >= λ); |
| 204 // success | 238 // success |
| 205 x.zip_apply(y, |x_i, y_i| { | 239 x.zip_apply(y, |x_i, y_i| { |
| 206 let (a, _, _, _) = shifted_nonneg_soft_thresholding(*x_i, y_i, λ_new); | 240 (*x_i, _, _, _) = shifted_nonneg_soft_thresholding(*x_i, y_i, λ_new); |
| 207 //*x_i = y_i + lbd_soft_thresholding(*x_i, λ_new, -y_i) | |
| 208 *x_i = a; | |
| 209 }); | 241 }); |
| 210 return; | 242 return; |
| 211 } | 243 } |
| 212 } | 244 } |
| 213 } | 245 } |
