spectrafit_newton_cg/
step.rs

1//! Steihaug–Toint truncated conjugate-gradient subproblem solver.
2//!
3//! Approximately minimises the scaled quadratic model `g̃ᵀδ̃ + ½ δ̃ᵀHδ̃` (with
4//! `H = J̃ᵀJ̃`) inside the trust region `‖δ̃‖ ≤ Δ` using conjugate gradients, and
5//! truncates early on two events (Nocedal & Wright, *Numerical Optimization*,
6//! Algorithm 7.2):
7//!
8//! * **negative curvature** (`dᵀHd ≤ 0`) — follow the search direction to the
9//!   trust-region boundary;
10//! * **trust-region exit** (`‖z + αd‖ ≥ Δ`) — stop at the boundary.
11//!
12//! Otherwise CG runs until the model gradient is small (a forcing-sequence
13//! tolerance giving superlinear convergence). The Gauss–Newton Hessian is **never
14//! formed**: every iteration uses one matrix-free product `H·v = J̃ᵀ(J̃·v)` from
15//! the [`Subproblem`], so the cost scales with the *residual* count, not `p²` —
16//! the lever for large nD fits.
17
18use faer::Mat;
19use spectrafit_trust_region::{StepResult, Subproblem, SubproblemStep};
20
21/// Matrix-free Newton-CG (Steihaug–Toint) trust-region subproblem solver (stateless).
22pub struct SteihaugStep;
23
24#[inline]
25fn dot(a: &Mat<f64>, b: &Mat<f64>) -> f64 {
26    (a.as_ref().transpose() * b.as_ref())[(0, 0)]
27}
28
29#[inline]
30fn norm(v: &Mat<f64>) -> f64 {
31    v.as_ref().squared_norm_l2().sqrt()
32}
33
34/// Largest `τ ≥ 0` with `‖z + τ·d‖ = Δ` (positive root; `z` is interior so it exists).
35fn boundary_tau(z: &Mat<f64>, d: &Mat<f64>, radius: f64) -> f64 {
36    let a = dot(d, d);
37    if a <= 0.0 {
38        return 0.0;
39    }
40    let b = 2.0 * dot(z, d);
41    let c = dot(z, z) - radius * radius;
42    let disc = (b * b - 4.0 * a * c).max(0.0);
43    ((-b + disc.sqrt()) / (2.0 * a)).max(0.0)
44}
45
46impl SubproblemStep for SteihaugStep {
47    fn solve(&self, sub: &Subproblem<'_>, radius: f64) -> StepResult {
48        let p = sub.n_params();
49        let g = sub.gradient(); // g̃
50
51        let finish = |step: Mat<f64>, hit_boundary: bool| {
52            let predicted_reduction = sub.predicted_reduction(step.as_ref());
53            StepResult {
54                step,
55                predicted_reduction,
56                hit_boundary,
57            }
58        };
59
60        // CG on H δ = −g̃, starting from z₀ = 0 ⇒ residual r₀ = g̃, direction d₀ = −g̃.
61        let mut z = Mat::<f64>::zeros(p, 1);
62        let mut r = Mat::from_fn(p, 1, |i, _| g[(i, 0)]);
63        let r0_norm = norm(&r);
64        if r0_norm == 0.0 {
65            return finish(z, false); // already stationary
66        }
67        // Forcing-sequence tolerance ‖r‖ ≤ ‖r₀‖·min(0.5, √‖r₀‖) (superlinear).
68        let tol = r0_norm * 0.5_f64.min(r0_norm.sqrt());
69        let mut d = Mat::from_fn(p, 1, |i, _| -r[(i, 0)]);
70        let mut rr = dot(&r, &r);
71
72        let max_iter = 2 * p + 10;
73        for _ in 0..max_iter {
74            let hd = sub.hvec(d.as_ref());
75            let dhd = dot(&d, &hd);
76            if dhd <= 0.0 {
77                // Negative curvature: ride d to the boundary.
78                let tau = boundary_tau(&z, &d, radius);
79                let step = Mat::from_fn(p, 1, |i, _| z[(i, 0)] + tau * d[(i, 0)]);
80                return finish(step, true);
81            }
82            let alpha = rr / dhd;
83            let z_next = Mat::from_fn(p, 1, |i, _| z[(i, 0)] + alpha * d[(i, 0)]);
84            if norm(&z_next) >= radius {
85                // Exited the trust region: stop at the boundary along d.
86                let tau = boundary_tau(&z, &d, radius);
87                let step = Mat::from_fn(p, 1, |i, _| z[(i, 0)] + tau * d[(i, 0)]);
88                return finish(step, true);
89            }
90            z = z_next;
91            let r_next = Mat::from_fn(p, 1, |i, _| r[(i, 0)] + alpha * hd[(i, 0)]);
92            if norm(&r_next) <= tol {
93                return finish(z, false); // interior CG convergence
94            }
95            let rr_next = dot(&r_next, &r_next);
96            let beta = rr_next / rr;
97            d = Mat::from_fn(p, 1, |i, _| -r_next[(i, 0)] + beta * d[(i, 0)]);
98            r = r_next;
99            rr = rr_next;
100        }
101        finish(z, false) // budget exhausted — return the best interior iterate
102    }
103}
spectrafit_newton_cg/step.rs

spectrafit_newton_cg/
step.rs
