fetish_lib/
sherman_morrison.rs

1extern crate ndarray;
2extern crate ndarray_linalg;
3
4use ndarray::*;
5use crate::linalg_utils::*;
6use std::ops::AddAssign;
7use std::ops::SubAssign;
8
9///Compute A + w * uu^T and its inverse in one stroke
10///for the case where A is already symmetric, updating both `A` and `A_inv` in place.
11///This uses the sherman-morrison matrix inverse formula.
12pub fn sherman_morrison_update(A : &mut Array2<f32>, A_inv : &mut Array2<f32>,
13                        w : f32, u : ArrayView1<f32>) {
14    let w_u_outer_u = w * outer(u, u);
15    A.add_assign(&w_u_outer_u);
16
17    let A_inv_u = A_inv.dot(&u);
18    let numerator = w * outer(A_inv_u.view(), A_inv_u.view());
19
20    let w_u_A_inv_u = w * u.dot(&A_inv_u);
21    let denominator = 1.0f32 + w_u_A_inv_u;
22
23    let scale = 1.0f32 / denominator;
24    let scaled = scale * &numerator;
25
26    A_inv.sub_assign(&scaled);
27}
28
29#[cfg(test)]
30mod tests {
31    use super::*;
32    use crate::pseudoinverse::*;
33    use crate::test_utils::*;
34
35    #[test]
36    fn test_sherman_morrison() {
37        let dim = 5;
38        let A = random_psd_matrix(dim);
39        let A_inv = pseudoinverse_h(&A);
40        let w = random_scalar();
41        let u = random_vector(dim); 
42
43        let mut A_updated = A.clone(); 
44        let mut A_inv_updated = A_inv.clone();
45        sherman_morrison_update(&mut A_updated, &mut A_inv_updated, w, u.view());
46
47        let expected_A_inv = pseudoinverse_h(&A_updated);
48
49        assert_equal_matrices_to_within(A_inv_updated.view(), expected_A_inv.view(), 0.001f32);
50    }
51}
fetish_lib/sherman_morrison.rs

fetish_lib/
sherman_morrison.rs
