flow_linalg/
condition.rs

1//! Matrix condition number (κ₂) and complexity index.
2
3use faer::MatRef;
4
5/// Condition number (κ₂ = σ_max / σ_min) and complexity index (log₁₀ κ).
6#[derive(Debug, Clone, Copy, PartialEq)]
7pub struct ConditionMetrics {
8    pub condition_number: f64,
9    pub complexity_index: f64,
10    pub singular: bool,
11}
12
13/// Compute κ₂ and complexity for a (possibly rectangular) `f64` matrix via SVD.
14///
15/// For an m×n mixing matrix, κ₂ = σ_max / σ_min over the nonzero singular values
16/// (same 2-norm condition used for square spillover matrices).
17pub fn condition_metrics(matrix: MatRef<'_, f64>) -> Result<ConditionMetrics, String> {
18    let m = matrix.nrows();
19    let n = matrix.ncols();
20    if m == 0 || n == 0 {
21        return Err("condition metrics require a non-empty matrix".into());
22    }
23    let sigma = matrix
24        .singular_values()
25        .map_err(|e| format!("SVD failed while assessing matrix condition: {e:?}"))?;
26    if sigma.is_empty() {
27        return Err("SVD returned no singular values".into());
28    }
29    let sigma_max = sigma[0];
30    let sigma_min = *sigma.last().unwrap_or(&0.0);
31    let singular = !(sigma_min.is_finite() && sigma_min > f64::EPSILON);
32    let condition_number = if singular {
33        f64::INFINITY
34    } else {
35        sigma_max / sigma_min
36    };
37    let complexity_index = if !condition_number.is_finite() {
38        f64::INFINITY
39    } else {
40        condition_number.max(1.0).log10()
41    };
42    Ok(ConditionMetrics {
43        condition_number,
44        complexity_index,
45        singular,
46    })
47}
48
49/// Compute κ₂ and complexity for a (possibly rectangular) `f32` matrix (promotes to f64 for SVD).
50pub fn condition_metrics_f32(matrix: MatRef<'_, f32>) -> Result<ConditionMetrics, String> {
51    let m = matrix.nrows();
52    let n = matrix.ncols();
53    if m == 0 || n == 0 {
54        return Err("condition metrics require a non-empty matrix".into());
55    }
56    let owned = faer::Mat::<f64>::from_fn(m, n, |i, j| f64::from(matrix[(i, j)]));
57    condition_metrics(owned.as_ref())
58}
59
60/// 2-norm condition number only (∞ when singular / empty SVD).
61pub fn condition_number_2(matrix: MatRef<'_, f64>) -> f64 {
62    condition_metrics(matrix)
63        .map(|m| m.condition_number)
64        .unwrap_or(f64::INFINITY)
65}
66
67#[cfg(test)]
68mod tests {
69    use super::*;
70    use faer::Mat;
71
72    #[test]
73    fn identity_has_kappa_one() {
74        let m = Mat::<f64>::from_fn(2, 2, |i, j| if i == j { 1.0 } else { 0.0 });
75        let metrics = condition_metrics(m.as_ref()).expect("metrics");
76        assert!(!metrics.singular);
77        assert!((metrics.condition_number - 1.0).abs() < 1e-9);
78        assert!(metrics.complexity_index.abs() < 1e-9);
79    }
80
81    #[test]
82    fn f32_identity_matches() {
83        let m = Mat::<f32>::from_fn(2, 2, |i, j| if i == j { 1.0 } else { 0.0 });
84        let metrics = condition_metrics_f32(m.as_ref()).expect("metrics");
85        assert!((metrics.condition_number - 1.0).abs() < 1e-5);
86    }
87
88    #[test]
89    fn tall_rectangular_is_finite() {
90        // 3 detectors × 2 endmembers (normalized-ish columns).
91        let m = Mat::<f64>::from_fn(3, 2, |i, j| if i == j { 1.0 } else { 0.1 });
92        let metrics = condition_metrics(m.as_ref()).expect("metrics");
93        assert!(!metrics.singular);
94        assert!(metrics.condition_number.is_finite());
95        assert!(metrics.condition_number >= 1.0);
96    }
97}
flow_linalg/condition.rs

flow_linalg/
condition.rs
