Function russell_lab::eigen_decomp

pub fn eigen_decomp(
    l_real: &mut [f64], 
    l_imag: &mut [f64], 
    v_real: &mut Matrix, 
    v_imag: &mut Matrix, 
    a: &mut Matrix
) -> Result<(), &'static str>

Performs the eigen-decomposition of a square matrix

Computes the eigenvalues l and right eigenvectors v, such that:

a ⋅ vj = lj ⋅ vj

where lj is the component j of l and vj is the column j of v.

Output

• l_real – (m) eigenvalues; real part
• l_imag – (m) eigenvalues; imaginary part
• v_real – (m,m) right eigenvectors (as columns); real part
• v_imag – (m,m) right eigenvectors (as columns); imaginary part

Input

• a – (m,m) general matrix [will be modified]

Note

• The matrix a will be modified

Similarity transformation

The eigen-decomposition leads to a similarity transformation like so:

a = v⋅λ⋅v⁻¹

where v is a matrix whose columns are the m linearly independent eigenvectors of a, and λ is a matrix whose diagonal are the eigenvalues of a. Thus, the following is valid:

a⋅v = v⋅λ

Let us define the error err as follows:

err := a⋅v - v⋅λ

Example

// import
use russell_lab::*;
use russell_chk::*;

// set matrix
let data = [
    [2.0, 0.0, 0.0],
    [0.0, 3.0, 4.0],
    [0.0, 4.0, 9.0],
];
let mut a = Matrix::from(&data);

// allocate output arrays
let m = a.nrow();
let mut l_real = vec![0.0; m];
let mut l_imag = vec![0.0; m];
let mut v_real = Matrix::new(m, m);
let mut v_imag = Matrix::new(m, m);

// perform the eigen-decomposition
eigen_decomp(
    &mut l_real,
    &mut l_imag,
    &mut v_real,
    &mut v_imag,
    &mut a,
)?;

// check results
let l_real_correct = "[11.0, 1.0, 2.0]";
let l_imag_correct = "[0.0, 0.0, 0.0]";
let v_real_correct = "┌                      ┐\n\
                      │  0.000  0.000  1.000 │\n\
                      │  0.447  0.894  0.000 │\n\
                      │  0.894 -0.447  0.000 │\n\
                      └                      ┘";
let v_imag_correct = "┌       ┐\n\
                      │ 0 0 0 │\n\
                      │ 0 0 0 │\n\
                      │ 0 0 0 │\n\
                      └       ┘";
assert_eq!(format!("{:?}", l_real), l_real_correct);
assert_eq!(format!("{:?}", l_imag), l_imag_correct);
assert_eq!(format!("{:.3}", v_real), v_real_correct);
assert_eq!(format!("{}", v_imag), v_imag_correct);

// check eigen-decomposition (similarity transformation) of a
// symmetric matrix with real-only eigenvalues and eigenvectors
let a_copy = Matrix::from(&data);
let lam = Matrix::diagonal(&l_real);
let mut a_v = Matrix::new(m, m);
let mut v_l = Matrix::new(m, m);
let mut err = Matrix::filled(m, m, f64::MAX);
mat_mat_mul(&mut a_v, 1.0, &a_copy, &v_real)?;
mat_mat_mul(&mut v_l, 1.0, &v_real, &lam)?;
add_matrices(&mut err, 1.0, &a_v, -1.0, &v_l)?;
assert_approx_eq!(err.norm(EnumMatrixNorm::Max), 0.0, 1e-15);