sim_lib_numbers_optimize/
minimize.rs

1//! Projected multivariate objective minimization.
2
3use super::*;
4
5/// Projected BFGS/trust-region objective minimization with acceptance ratio.
6pub fn minimize<F, G>(
7    mut f: F,
8    gradient: Option<G>,
9    mut x: Vec<f64>,
10    plan: &ObjectivePlan,
11) -> Result<OptimizeResult, Error>
12where
13    F: FnMut(&[f64]) -> f64,
14    G: Fn(&[f64], &mut [f64]),
15{
16    let n = x.len();
17    if plan.bounds.lower.len() != n {
18        return Err(Error::Dimension("bounds and point differ"));
19    }
20    validate_scale(&plan.scale, n)?;
21    if plan.policy != StepPolicy::ProjectedBfgs || plan.initial_radius <= 0.0 {
22        return Err(Error::InvalidPlan(
23            "multivariate objective requires projected BFGS and positive radius",
24        ));
25    }
26    plan.bounds.project(&mut x);
27    let memory = n * n * 8 + n * 40;
28    if memory > plan.limits.memory_bytes {
29        return Ok(opt_result(
30            x,
31            f64::NAN,
32            vec![],
33            Termination::WorkLimit,
34            0,
35            0,
36            memory,
37        ));
38    }
39    let mut h = vec![vec![0.0; n]; n];
40    for (i, row) in h.iter_mut().enumerate() {
41        row[i] = 1.0
42    }
43    let mut fx = f(&x);
44    let mut evals = 1;
45    let mut radius = plan.initial_radius;
46    if !fx.is_finite() {
47        return Ok(opt_result(
48            x,
49            fx,
50            vec![],
51            Termination::NonFinite,
52            evals,
53            0,
54            memory,
55        ));
56    }
57    for iter in 0..plan.limits.iterations {
58        let mut g = vec![0.0; n];
59        if let Some(ref grad) = gradient {
60            grad(&x, &mut g)
61        } else if let Some(v) = numerical_gradient(
62            &mut f,
63            &x,
64            fx,
65            &plan.scale,
66            &mut evals,
67            plan.limits.evaluations,
68        ) {
69            g = v
70        } else {
71            return Ok(opt_result(
72                x,
73                fx,
74                g,
75                Termination::NonFinite,
76                evals,
77                iter,
78                memory,
79            ));
80        };
81        if !finite(&g) {
82            return Ok(opt_result(
83                x,
84                fx,
85                g,
86                Termination::NonFinite,
87                evals,
88                iter,
89                memory,
90            ));
91        }
92        let mut pg = g.clone();
93        for i in 0..n {
94            if (x[i] <= plan.bounds.lower[i] && g[i] > 0.0)
95                || (x[i] >= plan.bounds.upper[i] && g[i] < 0.0)
96            {
97                pg[i] = 0.0
98            }
99        }
100        if norm(&pg) <= plan.tolerances.gradient {
101            let boundary = pg != g;
102            return Ok(opt_result(
103                x,
104                fx,
105                g,
106                if boundary {
107                    Termination::BoundaryConverged
108                } else {
109                    Termination::Converged
110                },
111                evals,
112                iter,
113                memory,
114            ));
115        }
116        let mut p = vec![0.0; n];
117        for i in 0..n {
118            p[i] = -h[i].iter().zip(&pg).map(|(a, b)| a * b).sum::<f64>() / plan.scale[i]
119        }
120        let pn = norm(&p);
121        if pn > radius {
122            for z in &mut p {
123                *z *= radius / pn
124            }
125        }
126        let mut y = x.iter().zip(&p).map(|(a, b)| a + b).collect::<Vec<_>>();
127        plan.bounds.project(&mut y);
128        let step = y.iter().zip(&x).map(|(a, b)| a - b).collect::<Vec<_>>();
129        if norm(&step) <= plan.tolerances.argument {
130            return Ok(opt_result(
131                x,
132                fx,
133                g,
134                Termination::NoProgress,
135                evals,
136                iter,
137                memory,
138            ));
139        }
140        if evals >= plan.limits.evaluations {
141            return Ok(opt_result(
142                x,
143                fx,
144                g,
145                Termination::WorkLimit,
146                evals,
147                iter,
148                memory,
149            ));
150        }
151        let fy = f(&y);
152        evals += 1;
153        if !fy.is_finite() {
154            radius *= 0.25;
155            continue;
156        }
157        let predicted = (-pg.iter().zip(&step).map(|(a, b)| a * b).sum::<f64>()).max(f64::EPSILON);
158        let ratio = (fx - fy) / predicted;
159        if ratio > 0.1 {
160            let old = x;
161            x = y;
162            let old_fx = fx;
163            fx = fy;
164            let mut ng = vec![0.0; n];
165            if let Some(ref grad) = gradient {
166                grad(&x, &mut ng)
167            } else if let Some(v) = numerical_gradient(
168                &mut f,
169                &x,
170                fx,
171                &plan.scale,
172                &mut evals,
173                plan.limits.evaluations,
174            ) {
175                ng = v
176            } else {
177                return Ok(opt_result(
178                    x,
179                    fx,
180                    g,
181                    Termination::WorkLimit,
182                    evals,
183                    iter,
184                    memory,
185                ));
186            };
187            let s = x.iter().zip(&old).map(|(a, b)| a - b).collect::<Vec<_>>();
188            let q = ng.iter().zip(&g).map(|(a, b)| a - b).collect::<Vec<_>>();
189            let sq = s.iter().zip(&q).map(|(a, b)| a * b).sum::<f64>();
190            if sq > 1e-14 {
191                let rho = 1.0 / sq;
192                let hq = h
193                    .iter()
194                    .map(|r| r.iter().zip(&q).map(|(a, b)| a * b).sum::<f64>())
195                    .collect::<Vec<_>>();
196                let qhq = q.iter().zip(&hq).map(|(a, b)| a * b).sum::<f64>();
197                for i in 0..n {
198                    for j in 0..n {
199                        h[i][j] += (1.0 + qhq * rho) * rho * s[i] * s[j]
200                            - rho * (s[i] * hq[j] + hq[i] * s[j]);
201                    }
202                }
203            }
204            if (old_fx - fx).abs() <= plan.tolerances.objective {
205                return Ok(opt_result(
206                    x,
207                    fx,
208                    ng,
209                    Termination::Converged,
210                    evals,
211                    iter + 1,
212                    memory,
213                ));
214            }
215            if ratio > 0.75 {
216                radius *= 2.0
217            }
218        } else {
219            radius *= 0.25
220        }
221    }
222    Ok(opt_result(
223        x,
224        fx,
225        vec![],
226        Termination::WorkLimit,
227        evals,
228        plan.limits.iterations,
229        memory,
230    ))
231}
232pub(crate) fn opt_result(
233    x: Vec<f64>,
234    value: f64,
235    g: Vec<f64>,
236    termination: Termination,
237    evaluations: usize,
238    iterations: usize,
239    memory_bytes: usize,
240) -> OptimizeResult {
241    OptimizeResult {
242        active: vec![],
243        gradient_norm: norm(&g),
244        point: x,
245        value,
246        termination,
247        work: Work {
248            evaluations,
249            iterations,
250            memory_bytes,
251        },
252    }
253}
254
255pub(crate) fn solve(mut a: Vec<Vec<f64>>, mut b: Vec<f64>, tol: f64) -> (Vec<f64>, usize) {
256    let n = b.len();
257    let mut rank = 0;
258    for k in 0..n {
259        let mut p = k;
260        for i in k + 1..n {
261            if a[i][k].abs() > a[p][k].abs() {
262                p = i
263            }
264        }
265        if a[p][k].abs() <= tol {
266            continue;
267        }
268        a.swap(k, p);
269        b.swap(k, p);
270        rank += 1;
271        let pivot_row = a[k].clone();
272        for i in k + 1..n {
273            let q = a[i][k] / a[k][k];
274            for (value, pivot) in a[i][k..].iter_mut().zip(&pivot_row[k..]) {
275                *value -= q * pivot
276            }
277            b[i] -= q * b[k]
278        }
279    }
280    let mut x = vec![0.0; n];
281    for i in (0..n).rev() {
282        if a[i][i].abs() > tol {
283            x[i] = (b[i]
284                - a[i]
285                    .iter()
286                    .enumerate()
287                    .skip(i + 1)
288                    .map(|(j, v)| v * x[j])
289                    .sum::<f64>())
290                / a[i][i]
291        }
292    }
293    (x, rank)
294}
295pub(crate) fn normal(j: &[Vec<f64>], r: &[f64], damping: f64) -> (Vec<Vec<f64>>, Vec<f64>) {
296    let n = j.first().map_or(0, Vec::len);
297    let mut a = vec![vec![0.0; n]; n];
298    let mut b = vec![0.0; n];
299    for (row, ri) in j.iter().zip(r) {
300        for p in 0..n {
301            b[p] -= row[p] * ri;
302            for q in 0..n {
303                a[p][q] += row[p] * row[q]
304            }
305        }
306    }
307    for (i, row) in a.iter_mut().enumerate() {
308        row[i] += damping
309    }
310    (a, b)
311}
312
313pub(crate) fn reflective_step(
314    j: &[Vec<f64>],
315    r: &[f64],
316    x: &[f64],
317    bounds: &Bounds,
318    damping: f64,
319    radius: f64,
320    tolerance: f64,
321) -> Vec<f64> {
322    let (mut a, b) = normal(j, r, 0.0);
323    let distance = b
324        .iter()
325        .enumerate()
326        .map(|(i, descent)| {
327            if *descent >= 0.0 {
328                bounds.upper[i] - x[i]
329            } else {
330                x[i] - bounds.lower[i]
331            }
332        })
333        .map(|v| v.max(tolerance).sqrt())
334        .collect::<Vec<_>>();
335    for i in 0..x.len() {
336        for k in 0..x.len() {
337            a[i][k] *= distance[i] * distance[k];
338        }
339        a[i][i] += damping;
340    }
341    let scaled_rhs = b
342        .iter()
343        .zip(&distance)
344        .map(|(v, d)| v * d)
345        .collect::<Vec<_>>();
346    let (mut scaled, _) = solve(a, scaled_rhs, tolerance);
347    let scaled_norm = norm(&scaled);
348    if scaled_norm > radius {
349        for value in &mut scaled {
350            *value *= radius / scaled_norm;
351        }
352    }
353    scaled
354        .iter()
355        .zip(distance)
356        .map(|(value, d)| value * d)
357        .collect()
358}
sim_lib_numbers_optimize/minimize.rs

sim_lib_numbers_optimize/
minimize.rs
