mathru 0.16.2

Fundamental algorithms for scientific computing in Rust
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

use crate::algebra::abstr::{AbsDiffEq, Field, Scalar};
use crate::algebra::linear::matrix::{Diagonal, UnitUpperTriangular};
use crate::algebra::linear::matrix::{EigenDec, EigenDecomposition, General};
use crate::elementary::Power;

impl<T> EigenDecomposition<T> for UnitUpperTriangular<T>
where
    T: Field + Scalar + Power + AbsDiffEq<Epsilon = T>,
{
    /// Computes the eigen decomposition
    ///
    /// # Arguments
    ///
    /// # Example
    ///
    /// ```
    /// use mathru::algebra::linear::matrix::{EigenDec, EigenDecomposition, General, UnitUpperTriangular};
    /// use mathru::{matrix, assert_relative_eq};
    ///
    /// let a: UnitUpperTriangular<f64> = matrix![  1.0, 3.0, 6.0;
    ///                                             0.0, 1.0, -6.0;
    ///                                             0.0, 0.0, 1.0].into();
    ///
    /// let (values, vectors) = a.dec_eigen().unwrap().pair();
    ///
    /// assert_relative_eq!(
    ///    &Into::<General<f64>>::into(a) * &vectors,
    ///    &vectors * &Into::<General<f64>>::into(values),
    ///    epsilon = 1.0e-5
    /// );
    /// ```
    fn dec_eigen(&self) -> Result<EigenDec<T>, String> {
        let (m, _): (usize, usize) = self.dim();
        let eigen_values = vec![T::one(); m];
        let values: Diagonal<T> = Diagonal::new(&eigen_values);

        let vectors: General<T> = self.calc_eigenvector(&eigen_values);
        Ok(EigenDec::new(values, vectors))
    }
}

impl<T> UnitUpperTriangular<T>
where
    T: Field + Scalar + Power + AbsDiffEq<Epsilon = T>,
{
    pub fn calc_eigenvector(&self, eigen_values: &[T]) -> General<T> {
        let mut x: General<T> = General::zero(eigen_values.len(), eigen_values.len());

        for (lambda_idx, lambda) in eigen_values.iter().rev().enumerate() {
            for row in (0..eigen_values.len()).rev() {
                let sum = (row..self.matrix.n).fold(T::zero(), |s, i| {
                    s + self.matrix[[row, i]] * x[[i, lambda_idx]]
                });

                let divisor = *lambda - self.matrix[[row, row]];
                if divisor != T::zero() {
                    x[[row, lambda_idx]] = sum / divisor;
                } else if row == 0 {
                    x[[row, lambda_idx]] = T::one();
                } else {
                    x[[row, lambda_idx]] = T::zero();
                }
            }
        }

        UnitUpperTriangular::norm(x)
    }

    fn norm(mut m: General<T>) -> General<T> {
        let (rows, cols) = m.dim();
        for c in 0..cols {
            let norm = m.get_column(c).p_norm(&T::from_f32(2.0));

            for r in 0..rows {
                m[[r, c]] /= norm;
            }
        }

        m
    }
}