spectrafit_trust_region/
step.rs

1//! The trust-region **subproblem** contract shared by the Δ-radius methods
2//! (dogleg, Newton-CG/Steihaug).
3//!
4//! At each outer iteration the [`driver`](crate::driver) hands a subproblem
5//! solver the local quadratic model in **scaled coordinates** `δ̃ = D·δ`:
6//! minimise
7//! ```text
8//!   m(δ̃) = g̃ᵀδ̃ + ½ δ̃ᵀ J̃ᵀJ̃ δ̃     subject to   ‖δ̃‖ ≤ Δ
9//! ```
10//! where `J̃ = J·diag(1/D)` is the column-scaled Jacobian and `g̃ = diag(1/D)·g`
11//! the scaled gradient (`g = Jᵀr`). Working in scaled coordinates turns the
12//! trust region into a plain Euclidean ball, so each method solves the standard
13//! subproblem and the driver maps the result back with `δ = diag(1/D)·δ̃`.
14//!
15//! The Gauss–Newton Hessian `J̃ᵀJ̃` is never formed: [`Subproblem`] exposes only
16//! the matrix-free products `J̃·v` and `J̃ᵀ·u`, which is exactly what a Krylov
17//! method (Steihaug-CG) needs and what keeps the large-residual case cheap.
18
19use faer::{Mat, MatRef};
20
21/// The trust-region subproblem at one outer iteration, in scaled coordinates.
22///
23/// Holds the column-scaled Jacobian `J̃` and scaled gradient `g̃`; a solver
24/// returns a scaled step `δ̃` with `‖δ̃‖ ≤ Δ`.
25pub struct Subproblem<'a> {
26    /// Column-scaled Jacobian `J̃ = J·diag(1/D)` (`m×p`).
27    jacobian: MatRef<'a, f64>,
28    /// Scaled gradient `g̃ = diag(1/D)·Jᵀr` (`p×1`).
29    gradient: MatRef<'a, f64>,
30}
31
32impl<'a> Subproblem<'a> {
33    /// Build a subproblem view over the scaled Jacobian and scaled gradient.
34    pub fn new(jacobian: MatRef<'a, f64>, gradient: MatRef<'a, f64>) -> Self {
35        Self { jacobian, gradient }
36    }
37
38    /// Number of free parameters `p`.
39    pub fn n_params(&self) -> usize {
40        self.jacobian.ncols()
41    }
42
43    /// The scaled gradient `g̃` (`p×1`).
44    pub fn gradient(&self) -> MatRef<'_, f64> {
45        self.gradient
46    }
47
48    /// The column-scaled Jacobian `J̃` (`m×p`). Dense methods (dogleg) form the
49    /// Gauss–Newton system from it; Krylov methods should prefer [`hvec`](Self::hvec).
50    pub fn jacobian(&self) -> MatRef<'_, f64> {
51        self.jacobian
52    }
53
54    /// Matrix-free product `J̃·v` (`m×1`) — never forms `J̃ᵀJ̃`.
55    pub fn jvec(&self, v: MatRef<'_, f64>) -> Mat<f64> {
56        self.jacobian * v
57    }
58
59    /// Matrix-free product `J̃ᵀ·u` (`p×1`).
60    pub fn jtvec(&self, u: MatRef<'_, f64>) -> Mat<f64> {
61        self.jacobian.transpose() * u
62    }
63
64    /// Hessian–vector product `H·v = J̃ᵀ(J̃·v)` (`p×1`) without forming `H`.
65    pub fn hvec(&self, v: MatRef<'_, f64>) -> Mat<f64> {
66        let jv = self.jacobian * v;
67        self.jacobian.transpose() * jv.as_ref()
68    }
69
70    /// Predicted reduction of `½‖r‖²` for a scaled step `δ̃`:
71    /// `−g̃ᵀδ̃ − ½‖J̃δ̃‖²` (positive for a descent step).
72    pub fn predicted_reduction(&self, s: MatRef<'_, f64>) -> f64 {
73        let gs = self.gradient.transpose() * s; // 1×1
74        let js = self.jacobian * s;
75        -gs[(0, 0)] - 0.5 * js.as_ref().squared_norm_l2()
76    }
77}
78
79/// A scaled step produced by a [`SubproblemStep`].
80pub struct StepResult {
81    /// The scaled step `δ̃` (`p×1`); the driver maps it back via `δ = δ̃/D`.
82    pub step: Mat<f64>,
83    /// Predicted reduction `−g̃ᵀδ̃ − ½‖J̃δ̃‖²` (≥ 0 for a useful step).
84    pub predicted_reduction: f64,
85    /// Whether the step reached the trust-region boundary `‖δ̃‖ = Δ` (the driver
86    /// expands Δ on a very good boundary step).
87    pub hit_boundary: bool,
88}
89
90/// A trust-region subproblem solver: minimise the local quadratic model within
91/// the scaled radius `Δ`. One implementation per method (dogleg, Steihaug-CG).
92pub trait SubproblemStep {
93    /// Solve the subproblem for trust radius `radius` (in scaled units).
94    fn solve(&self, sub: &Subproblem<'_>, radius: f64) -> StepResult;
95}
spectrafit_trust_region/step.rs

spectrafit_trust_region/
step.rs
