spectrafit_trust_region/problem.rs
1//! The [`TrustRegionProblem`] trait — the only coupling point between this
2//! numerics-only crate and a concrete model/graph.
3//!
4//! An implementor owns the parameter state and is responsible for any
5//! parameterisation-level projection (bounds reflection, expression-tied
6//! parameters) inside [`set_params`](TrustRegionProblem::set_params). The
7//! driver treats the problem as an opaque source of weighted residuals and a
8//! weighted Jacobian written into caller-owned `faer` buffers — no allocation
9//! on the hot path.
10
11use faer::{Mat, MatMut};
12use spectrafit_types::CoreError;
13
14/// A weighted nonlinear least-squares problem driven by the trust-region core.
15///
16/// Conventions:
17/// * "weighted" means residuals and Jacobian already carry the per-point `1/σ`
18/// factor, so the driver minimises `½‖r‖²` directly.
19/// * `n_params` is the number of **free** parameters (tied/fixed excluded).
20/// * Bounds reflection and tied-parameter application live in
21/// [`set_params`](TrustRegionProblem::set_params);
22/// the driver never sees them, it only proposes raw steps.
23pub trait TrustRegionProblem {
24 /// Number of residual rows `m` (total points across datasets).
25 fn n_residuals(&self) -> usize;
26
27 /// Number of free parameters `p`.
28 fn n_params(&self) -> usize;
29
30 /// Current free-parameter vector (length `p`).
31 ///
32 /// Called after [`set_params`](Self::set_params) so the driver can read back
33 /// the value the problem actually applied (e.g. after bounds reflection).
34 fn params(&self) -> Vec<f64>;
35
36 /// Install a proposed free-parameter vector `p` (length `n_params`).
37 ///
38 /// The implementor applies any reflection / tied-plan here, so a subsequent
39 /// [`params`](Self::params) may differ from `p`.
40 fn set_params(&mut self, p: &[f64]);
41
42 /// Write the weighted residual vector `r` (shape `m × 1`) for the current
43 /// parameters into `r`.
44 fn residuals_into(&mut self, r: MatMut<'_, f64>) -> Result<(), CoreError>;
45
46 /// Write the weighted Jacobian `J` (shape `m × p`) for the current
47 /// parameters into `jac`.
48 fn jacobian_into(&mut self, jac: MatMut<'_, f64>) -> Result<(), CoreError>;
49
50 /// Apply the weighted Jacobian as a linear operator: write `out = J·v`
51 /// (length `m`) for the current parameters, where `v` has length `p`.
52 ///
53 /// This is the **matrix-free** half of the problem contract, intended for
54 /// Krylov subproblem solvers (truncated-CG / Steihaug) that need only
55 /// Jacobian–vector products, never the dense `J`. The default materializes
56 /// `J` via [`jacobian_into`](Self::jacobian_into) and multiplies — correct
57 /// but `O(m·p)` storage per call; a matrix-free implementor overrides both
58 /// operators to avoid forming `J`. Implementations must stay consistent with
59 /// `jacobian_into`: `apply_jacobian(v) == J·v`.
60 ///
61 /// **Reserved API — no driver in this workspace calls it today.** The
62 /// Newton-CG method works on the dense *scaled* `J` through
63 /// `Subproblem::hvec`, not through these operators, so the only exercise
64 /// the default implementations get is the framework's own contract test.
65 /// It is kept because it is the seam a future matrix-free implementor
66 /// needs, not because something behind it is already wired up.
67 fn apply_jacobian(&mut self, v: &[f64], out: &mut [f64]) -> Result<(), CoreError> {
68 let m = self.n_residuals();
69 let p = self.n_params();
70 let mut jac = Mat::<f64>::zeros(m, p);
71 self.jacobian_into(jac.as_mut())?;
72 for (i, o) in out.iter_mut().enumerate().take(m) {
73 let mut s = 0.0;
74 for (c, &vc) in v.iter().enumerate().take(p) {
75 s += jac[(i, c)] * vc;
76 }
77 *o = s;
78 }
79 Ok(())
80 }
81
82 /// Apply the transposed weighted Jacobian: write `out = Jᵀ·u` (length `p`)
83 /// for the current parameters, where `u` has length `m`.
84 ///
85 /// The transpose-multiply companion to [`apply_jacobian`](Self::apply_jacobian)
86 /// (forms the Krylov normal-equation operator `v ↦ Jᵀ(J·v)` without `JᵀJ`).
87 /// The default materializes `J` and multiplies; matrix-free implementors
88 /// override it. Must satisfy `apply_jacobian_transpose(u) == Jᵀ·u`.
89 ///
90 /// **Reserved API** on the same terms as
91 /// [`apply_jacobian`](Self::apply_jacobian): no driver calls it today and
92 /// the default is exercised only by tests.
93 fn apply_jacobian_transpose(&mut self, u: &[f64], out: &mut [f64]) -> Result<(), CoreError> {
94 let m = self.n_residuals();
95 let p = self.n_params();
96 let mut jac = Mat::<f64>::zeros(m, p);
97 self.jacobian_into(jac.as_mut())?;
98 for (c, o) in out.iter_mut().enumerate().take(p) {
99 let mut s = 0.0;
100 for (i, &ui) in u.iter().enumerate().take(m) {
101 s += jac[(i, c)] * ui;
102 }
103 *o = s;
104 }
105 Ok(())
106 }
107
108 /// Per-parameter scale factors (length `p`) used for column scaling of the
109 /// damping diagonal. Defaults to all-ones (Levenberg damping).
110 ///
111 /// This trait method itself has no callers in the workspace today —
112 /// `Parameter.scale` is already wired, but through a different path:
113 /// `spectrafit-solver`'s `LmProblem` (`lm_problem.rs`) holds its own
114 /// `scales` field and applies it directly when assembling the scaled
115 /// working-variable Jacobian, rather than by overriding this method.
116 fn scales(&self) -> Vec<f64> {
117 vec![1.0; self.n_params()]
118 }
119
120 /// Coleman–Li trust-region scaling `v ∈ (0, 1]` per free parameter for the
121 /// current point and gradient `grad`, or `None` when the problem is
122 /// unbounded. `v_i → 0` as parameter `i` approaches the active bound; the
123 /// driver folds `D_i ← D_i / √v_i`, so the trust region shrinks near bounds
124 /// (the "reflective" part of Trust-Region-Reflective). Default: `None`.
125 fn trust_scaling(&self, _grad: &[f64]) -> Option<Vec<f64>> {
126 None
127 }
128}