causal_hub/utils/
pseudo_inverse.rs

1use ndarray::prelude::*;
2use ndarray_linalg::SVD;
3use ndarray_stats::QuantileExt;
4
5use crate::types::EPSILON;
6
7/// Moore-Penrose pseudo-inverse.
8pub trait PseudoInverse {
9    /// Computes the Moore-Penrose pseudo-inverse of a matrix.
10    ///
11    /// # Panics
12    ///
13    /// * Panics if the SVD computation fails.
14    ///
15    /// # Returns
16    ///
17    /// The pseudo-inverse of the matrix.
18    ///
19    fn pinv(&self) -> Self;
20}
21
22impl PseudoInverse for Array2<f64> {
23    fn pinv(&self) -> Self {
24        // Step 0: Scale the matrix to improve numerical stability.
25        let a = *self.abs().max().unwrap_or(&1.);
26        let m = self / a;
27        // Step 1: Compute the Single Value Decomposition (SVD).
28        let (u, s, vt) = m.svd(true, true).unwrap_or_else(|e| {
29            panic!(
30                "Failed to compute the SVD \n\
31                \t for the matrix: \n\
32                \t {m:?} \n\
33                \t with error: \n\
34                \t {e:?}."
35            )
36        });
37        let u = u.expect("Failed to get U from the SVD.");
38        let vt = vt.expect("Failed to get VT from the SVD.");
39        // Step 2: Compute the pseudo-inverse of the singular values.
40        let s_max = s.max().unwrap_or(&0.);
41        let r_tol = f64::max(EPSILON, s.len() as f64 * s_max * EPSILON);
42        let s_inv = Array2::from_diag(&s.mapv(|x| if x > r_tol { 1. / x } else { 0. }));
43        // Step 3: Compute the pseudo-inverse of S_zz.
44        vt.t().dot(&s_inv).dot(&u.t()) / a
45    }
46}
causal_hub/utils/pseudo_inverse.rs

causal_hub/utils/
pseudo_inverse.rs
