spectrafit_dogleg/
step.rs

1//! Powell's dogleg subproblem solver.
2//!
3//! Given the scaled subproblem (`J̃`, `g̃`, radius `Δ`), the dogleg method picks a
4//! step on the piecewise-linear path Cauchy → Gauss–Newton:
5//!
6//! 1. If the Gauss–Newton step `δ_GN = −(J̃ᵀJ̃)⁻¹ g̃` lies inside the trust region
7//!    (`‖δ_GN‖ ≤ Δ`), take it — full Newton convergence near the optimum.
8//! 2. Otherwise compute the Cauchy point `δ_C = −(‖g̃‖²/‖J̃g̃‖²) g̃` (the
9//!    unconstrained steepest-descent minimiser). If it is already outside `Δ`,
10//!    take the steepest-descent direction clipped to the boundary.
11//! 3. Otherwise interpolate along `δ_C → δ_GN` to the point where `‖δ‖ = Δ`.
12//!
13//! A rank-deficient `J̃ᵀJ̃` (no Gauss–Newton step) degrades gracefully to the
14//! Cauchy / steepest-descent step.
15
16use faer::prelude::*;
17use faer::{Mat, Side};
18use spectrafit_trust_region::{StepResult, Subproblem, SubproblemStep};
19
20/// Powell's dogleg trust-region subproblem solver (stateless).
21pub struct DoglegStep;
22
23#[inline]
24fn norm2(v: &Mat<f64>) -> f64 {
25    v.as_ref().squared_norm_l2()
26}
27
28#[inline]
29fn dot(a: &Mat<f64>, b: &Mat<f64>) -> f64 {
30    (a.as_ref().transpose() * b.as_ref())[(0, 0)]
31}
32
33impl SubproblemStep for DoglegStep {
34    fn solve(&self, sub: &Subproblem<'_>, radius: f64) -> StepResult {
35        let p = sub.n_params();
36        let g = sub.gradient(); // g̃ (p×1)
37        let jt = sub.jacobian(); // J̃ (m×p)
38
39        let finish = |step: Mat<f64>, hit_boundary: bool| {
40            let predicted_reduction = sub.predicted_reduction(step.as_ref());
41            StepResult {
42                step,
43                predicted_reduction,
44                hit_boundary,
45            }
46        };
47
48        // Gauss–Newton step: (J̃ᵀJ̃) δ = −g̃ (None if not positive-definite).
49        let h = jt.transpose() * jt;
50        let neg_g = Mat::from_fn(p, 1, |i, _| -g[(i, 0)]);
51        let gn = h
52            .as_ref()
53            .llt(Side::Lower)
54            .ok()
55            .map(|llt| llt.solve(neg_g.as_ref()));
56
57        // 1. GN inside the trust region ⇒ take it.
58        if let Some(gn) = &gn {
59            if norm2(gn).sqrt() <= radius {
60                return finish(gn.clone(), false);
61            }
62        }
63
64        // Cauchy point along −g̃.
65        let g_norm2 = g.squared_norm_l2();
66        let g_norm = g_norm2.sqrt();
67        if g_norm == 0.0 {
68            return finish(Mat::zeros(p, 1), false); // stationary: no descent direction
69        }
70        let jg = sub.jvec(g); // J̃ g̃
71        let jg_norm2 = jg.as_ref().squared_norm_l2();
72
73        // Steepest-descent direction clipped to the boundary (used when the
74        // curvature along −g̃ is ~0 or the Cauchy point is already outside Δ).
75        let boundary_sd = || Mat::from_fn(p, 1, |i, _| -(radius / g_norm) * g[(i, 0)]);
76
77        if jg_norm2 <= 0.0 {
78            return finish(boundary_sd(), true);
79        }
80
81        let t_star = g_norm2 / jg_norm2;
82        let cauchy = Mat::from_fn(p, 1, |i, _| -t_star * g[(i, 0)]);
83        if norm2(&cauchy).sqrt() >= radius {
84            return finish(boundary_sd(), true);
85        }
86
87        // 3. Dogleg interpolation Cauchy → GN to the boundary ‖δ‖ = Δ.
88        let Some(gn) = gn else {
89            // Rank-deficient: no GN leg — accept the interior Cauchy step.
90            return finish(cauchy, false);
91        };
92        let d = Mat::from_fn(p, 1, |i, _| gn[(i, 0)] - cauchy[(i, 0)]);
93        let a = norm2(&d);
94        if a <= 0.0 {
95            return finish(cauchy, false); // GN ≈ Cauchy
96        }
97        // Solve a·τ² + b·τ + c = 0 for τ ∈ [0,1], c = ‖cauchy‖² − Δ² < 0.
98        let b = 2.0 * dot(&cauchy, &d);
99        let c = norm2(&cauchy) - radius * radius;
100        let disc = (b * b - 4.0 * a * c).max(0.0);
101        let tau = ((-b + disc.sqrt()) / (2.0 * a)).clamp(0.0, 1.0);
102        let step = Mat::from_fn(p, 1, |i, _| cauchy[(i, 0)] + tau * d[(i, 0)]);
103        finish(step, true)
104    }
105}
spectrafit_dogleg/step.rs

spectrafit_dogleg/
step.rs
