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<(), StrError>

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::{add_matrices, eigen_decomp, mat_mat_mul, matrix_norm, Matrix, NormMat, StrError};
use russell_chk::assert_approx_eq;

fn main() -> Result<(), StrError> {
    // 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]";
    assert_eq!(format!("{:?}", l_real), l_real_correct);
    assert_eq!(format!("{:?}", l_imag), l_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!(matrix_norm(&err, NormMat::Max), 0.0, 1e-15);
    Ok(())
}