bayes_estimate/linalg/
rcond.rs

1//! Numerical approximation of matrix reciprocal condition numbers 'rcond'.
2//!
3//! Matrices are well-conditioned if the reciprocal condition number is near 1 and ill-conditioned if it is near zero.
4//!
5//! For positive (semi-)definite matrices a simple definition is the minimum diagonal element / maximum diagonal element.
6//! This is best applied to Cholesky factorisation by using the square of the factors diagonal or for a UdU' by using 'd'
7//! directly.
8//!
9//! A rcond = 0 implies the matrix is semi definite, < 0 implies the matrix is negative.
10//!
11//! Required for all linear algebra in models and filters.
12
13use nalgebra::{allocator::Allocator, DefaultAllocator, Dim, OMatrix, RealField};
14
15/// Determine the reciprocal condition number of a symmetric matrix.
16///
17/// The 'rcond' is simply the minimum diagonal element / maximum diagonal element. A NaN result a negative result.
18pub fn rcond_symmetric<N: Copy + RealField, R: Dim, C: Dim>(sm: &OMatrix<N, R, C>) -> N
19where
20    DefaultAllocator: Allocator<R, C>,
21{
22    // Special case an empty matrix
23    let n = sm.nrows();
24    if n == 0 {
25        N::zero()
26    } else {
27        let mut mind = sm[(0, 0)];
28        let mut maxd = mind;
29
30        for i in 0..n {
31            let d = sm[(i, i)];
32            if d != d {
33                // NaN
34                mind = N::one().neg();
35                break;
36            }
37            if d < mind {
38                mind = d;
39            }
40            if d > maxd {
41                maxd = d;
42            }
43        }
44
45        rcond_min_max(mind, maxd)
46    }
47}
48
49/// Determine the 'rcond' from minimum and maximum diagonal elements.
50fn rcond_min_max<N: RealField>(mind: N, maxd: N) -> N {
51    if mind < N::zero() {
52        // matrix is negative
53        mind // mind < 0 but does not represent a rcond
54    } else {
55        assert!(mind <= maxd); // check sanity
56
57        // NOTE mind -0 will be handled here and propagate into rcond
58        let rcond = mind / maxd; // rcond from min/max norm
59        if !rcond.is_finite() {
60            // NaN, singular due to (mind == maxd) == (zero or infinity)
61            N::zero()
62        } else {
63            assert!(rcond <= N::one());
64            rcond
65        }
66    }
67}
68
69#[cfg(test)]
70mod test {
71    use super::*;
72
73    #[test]
74    fn test_simple_rconds() {
75        assert_eq!(rcond_min_max(1., 1.), 1.);
76        assert_eq!(rcond_min_max(1., 2.), 0.5);
77    }
78
79    #[test]
80    fn test_special_rconds() {
81        assert_eq!(rcond_min_max(-1., 55.), -1.);
82        assert_eq!(rcond_min_max(-2.,55.), -2.);
83        assert_eq!(rcond_min_max(-0.,55.), -0.);
84        assert_eq!(rcond_min_max(-0., 0.), 0.);
85        assert_eq!(rcond_min_max(-0., -0.), 0.);
86    }
87
88    #[should_panic(expected = "mind <= maxd")]
89    #[test]
90    fn test_greater_one_panics() {
91        rcond_min_max(2., 1.);
92    }
93
94}
95
96/// Checks the reciprocal condition number is > 0 .
97///
98/// IEC 559 NaN values are never true
99pub fn check_positive<N: RealField>(rcond: N, message: &str) -> Result<(), &str> {
100    if rcond > N::zero() {
101        Ok(())
102    } else {
103        Err(message)
104    }
105}
106
107/// Checks the reciprocal condition number is >= 0 .
108///
109/// IEC 559 NaN values are never true
110pub fn check_non_negative<N: RealField>(rcond: N, message: &str) -> Result<(), &str> {
111    if rcond >= N::zero() {
112        Ok(())
113    } else {
114        Err(message)
115    }
116}
117
118
bayes_estimate/linalg/rcond.rs

bayes_estimate/linalg/
rcond.rs
