numra_optim/
optim_sensitivity.rs

1//! Parametric sensitivity analysis for optimization problems.
2//!
3//! Given a parameterized problem where the objective and/or constraints
4//! depend on parameters `p`, compute `dx*/dp` by finite differences:
5//! for each parameter p_j, perturb p_j, re-solve, and compute
6//! (x*(p+h) - x*(p-h)) / (2h).
7//!
8//! Author: Moussa Leblouba
9//! Date: 8 February 2026
10//! Modified: 2 May 2026
11
12use numra_core::Scalar;
13
14use crate::error::OptimError;
15use crate::problem::OptimProblem;
16use crate::types::ParamSensitivity;
17
18/// Compute parametric sensitivity of an optimization problem.
19///
20/// `build_problem` is a closure that takes a parameter vector and returns
21/// an `OptimProblem` ready to solve. For each parameter, a central finite
22/// difference is used to estimate `dx*/dp_j`.
23///
24/// # Arguments
25///
26/// * `build_problem` - Factory that builds an `OptimProblem` from a parameter vector.
27/// * `params` - Nominal parameter values.
28/// * `names` - Human-readable names for each parameter.
29/// * `eps` - Relative perturbation size (default 1e-5).
30///
31/// # Returns
32///
33/// A `ParamSensitivity` matrix of shape `n_vars x n_params`.
34pub fn compute_param_sensitivity<S, F>(
35    build_problem: F,
36    params: &[S],
37    names: &[&str],
38    eps: Option<S>,
39) -> Result<ParamSensitivity<S>, OptimError>
40where
41    S: Scalar + faer::SimpleEntity + faer::Conjugate<Canonical = S> + faer::ComplexField,
42    F: Fn(&[S]) -> OptimProblem<S>,
43{
44    let eps = eps.unwrap_or_else(|| S::EPSILON.cbrt());
45    let n_params = params.len();
46
47    // Solve the nominal problem to get x* and n_vars.
48    let nominal = build_problem(params).solve()?;
49    let x_star = &nominal.x;
50    let n_vars = x_star.len();
51
52    let mut values = vec![S::ZERO; n_vars * n_params];
53
54    for j in 0..n_params {
55        let h = eps * (S::ONE + params[j].abs());
56
57        // Forward perturbation
58        let mut p_plus = params.to_vec();
59        p_plus[j] += h;
60        let result_plus = build_problem(&p_plus).solve()?;
61        let x_plus = &result_plus.x;
62
63        // Backward perturbation
64        let mut p_minus = params.to_vec();
65        p_minus[j] -= h;
66        let result_minus = build_problem(&p_minus).solve()?;
67        let x_minus = &result_minus.x;
68
69        // Central finite difference
70        for i in 0..n_vars {
71            values[i * n_params + j] = (x_plus[i] - x_minus[i]) / (S::TWO * h);
72        }
73    }
74
75    Ok(ParamSensitivity {
76        names: names.iter().map(|s| s.to_string()).collect(),
77        values,
78        n_vars,
79        n_params,
80    })
81}
82
83#[cfg(test)]
84mod tests {
85    use super::*;
86
87    #[test]
88    fn test_sensitivity_quadratic() {
89        // min (x - p)^2 => x* = p, so dx*/dp = 1
90        let params = [3.0];
91        let sens = compute_param_sensitivity(
92            |p: &[f64]| {
93                let p_val = p[0];
94                OptimProblem::new(1)
95                    .x0(&[0.0])
96                    .objective(move |x: &[f64]| (x[0] - p_val) * (x[0] - p_val))
97            },
98            &params,
99            &["p"],
100            None,
101        )
102        .unwrap();
103
104        assert_eq!(sens.n_vars, 1);
105        assert_eq!(sens.n_params, 1);
106        assert!(
107            (sens.get(0, 0) - 1.0).abs() < 1e-3,
108            "dx/dp = {}, expected 1.0",
109            sens.get(0, 0)
110        );
111    }
112
113    #[test]
114    fn test_sensitivity_two_params() {
115        // min (x - p1)^2 + (y - p2)^2 => x* = p1, y* = p2
116        // dx/dp1 = 1, dx/dp2 = 0, dy/dp1 = 0, dy/dp2 = 1
117        let params = [3.0, 7.0];
118        let sens = compute_param_sensitivity(
119            |p: &[f64]| {
120                let p1 = p[0];
121                let p2 = p[1];
122                OptimProblem::new(2)
123                    .x0(&[0.0, 0.0])
124                    .objective(move |x: &[f64]| {
125                        (x[0] - p1) * (x[0] - p1) + (x[1] - p2) * (x[1] - p2)
126                    })
127            },
128            &params,
129            &["p1", "p2"],
130            None,
131        )
132        .unwrap();
133
134        assert_eq!(sens.n_vars, 2);
135        assert_eq!(sens.n_params, 2);
136
137        // dx/dp1 = 1
138        assert!(
139            (sens.get(0, 0) - 1.0).abs() < 1e-3,
140            "dx/dp1 = {}, expected 1.0",
141            sens.get(0, 0)
142        );
143        // dx/dp2 = 0
144        assert!(
145            sens.get(0, 1).abs() < 1e-3,
146            "dx/dp2 = {}, expected 0.0",
147            sens.get(0, 1)
148        );
149        // dy/dp1 = 0
150        assert!(
151            sens.get(1, 0).abs() < 1e-3,
152            "dy/dp1 = {}, expected 0.0",
153            sens.get(1, 0)
154        );
155        // dy/dp2 = 1
156        assert!(
157            (sens.get(1, 1) - 1.0).abs() < 1e-3,
158            "dy/dp2 = {}, expected 1.0",
159            sens.get(1, 1)
160        );
161
162        // Also test column/row accessors
163        let col0 = sens.column(0);
164        assert!((col0[0] - 1.0).abs() < 1e-3);
165        assert!(col0[1].abs() < 1e-3);
166
167        let row1 = sens.row(1);
168        assert!(row1[0].abs() < 1e-3);
169        assert!((row1[1] - 1.0).abs() < 1e-3);
170    }
171
172    #[test]
173    fn test_sensitivity_bounded() {
174        // min (x - p)^2 with x in [0, 10]
175        // At p=5 (interior): x*=5, dx*/dp = 1
176        // At p=15 (bound active): x*=10, dx*/dp ~ 0
177
178        // Interior case: p = 5
179        let sens_interior = compute_param_sensitivity(
180            |p: &[f64]| {
181                let p_val = p[0];
182                OptimProblem::new(1)
183                    .x0(&[5.0])
184                    .objective(move |x: &[f64]| (x[0] - p_val) * (x[0] - p_val))
185                    .bounds(0, (0.0, 10.0))
186            },
187            &[5.0],
188            &["p"],
189            None,
190        )
191        .unwrap();
192
193        assert!(
194            (sens_interior.get(0, 0) - 1.0).abs() < 1e-3,
195            "interior dx/dp = {}, expected 1.0",
196            sens_interior.get(0, 0)
197        );
198
199        // Bound-active case: p = 15 => x* = 10 (upper bound)
200        let sens_bound = compute_param_sensitivity(
201            |p: &[f64]| {
202                let p_val = p[0];
203                OptimProblem::new(1)
204                    .x0(&[5.0])
205                    .objective(move |x: &[f64]| (x[0] - p_val) * (x[0] - p_val))
206                    .bounds(0, (0.0, 10.0))
207            },
208            &[15.0],
209            &["p"],
210            None,
211        )
212        .unwrap();
213
214        assert!(
215            sens_bound.get(0, 0).abs() < 0.1,
216            "bound-active dx/dp = {}, expected ~0",
217            sens_bound.get(0, 0)
218        );
219    }
220
221    #[test]
222    fn test_sensitivity_canonical_default_eps_at_large_param() {
223        // Regression: with `eps = None` the FD step factor must be the
224        // canonical `cbrt(EPSILON)`, additively scaled by `(1 + |p|)`.
225        // For min (x - p)^2 the analytical sensitivity is dx*/dp = 1
226        // independent of |p|. Exercising at |p| = 100 ensures the
227        // canonical factor (~6e-6) and the additive scaling cooperate
228        // to recover the analytical answer at large parameter magnitude
229        // — pinning the property the canonical formula provides, not
230        // just non-zero output.
231        let params = [100.0];
232        let sens = compute_param_sensitivity(
233            |p: &[f64]| {
234                let p_val = p[0];
235                OptimProblem::new(1)
236                    .x0(&[100.0])
237                    .objective(move |x: &[f64]| (x[0] - p_val) * (x[0] - p_val))
238            },
239            &params,
240            &["p"],
241            None,
242        )
243        .unwrap();
244
245        let s = sens.get(0, 0);
246        assert!(s.is_finite(), "sensitivity must be finite, got {s}");
247        assert!(
248            (s - 1.0).abs() < 1e-2,
249            "canonical default at |p| = 100: dx/dp = {s}, expected ~1.0"
250        );
251    }
252}
numra_optim/optim_sensitivity.rs

numra_optim/
optim_sensitivity.rs
