Function russell_lab::eigen_decomp_lr

pub fn eigen_decomp_lr(
    l_real: &mut [f64], 
    l_imag: &mut [f64], 
    u_real: &mut Matrix, 
    u_imag: &mut Matrix, 
    v_real: &mut Matrix, 
    v_imag: &mut Matrix, 
    a: &mut Matrix
) -> Result<(), StrError>

Performs the eigen-decomposition of a square matrix (left and right)

Computes the eigenvalues l and left eigenvectors u, such that:

ujᴴ ⋅ a = lj ⋅ ujᴴ

where lj is the component j of l and ujᴴ is the column j of uᴴ, with uᴴ being the conjugate-transpose of u.

Also, computes the right eigenvectors v, such that:

a ⋅ vj = lj ⋅ vj

where vj is the column j of v.

Output

• l_real – (m) eigenvalues; real part
• l_imag – (m) eigenvalues; imaginary part
• u_real – (m,m) left eigenvectors (as columns); real part
• u_imag – (m,m) left eigenvectors (as columns); 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

Example

use num_complex::Complex64;
use russell_chk::assert_approx_eq;
use russell_lab::{
    complex_add_matrices, complex_mat_mat_mul, complex_mat_zip,
    complex_matrix_norm, complex_vec_zip, ComplexMatrix, NormMat,
};
use russell_lab::{eigen_decomp_lr, Matrix, StrError};

fn main() -> Result<(), StrError> {
    // set matrix
    let data = [[0.0, 1.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 0.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 u_real = Matrix::new(m, m);
    let mut u_imag = Matrix::new(m, m);
    let mut v_real = Matrix::new(m, m);
    let mut v_imag = Matrix::new(m, m);

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

    // check results
    let l_real_correct = "[-0.5, -0.5, 0.9999999999999998]";
    let l_imag_correct = "[0.8660254037844389, -0.8660254037844389, 0.0]";
    assert_eq!(format!("{:?}", l_real), l_real_correct);
    assert_eq!(format!("{:?}", l_imag), l_imag_correct);

    // check the eigen-decomposition (similarity transformation)
    // ```text
    // a⋅v = v⋅λ
    // err := a⋅v - v⋅λ
    // ```
    let a = ComplexMatrix::from(&data);
    let v = complex_mat_zip(&v_real, &v_imag)?;
    let d = complex_vec_zip(&l_real, &l_imag)?;
    let lam = ComplexMatrix::diagonal(d.as_data());
    let mut a_v = ComplexMatrix::new(m, m);
    let mut v_l = ComplexMatrix::new(m, m);
    let mut err = ComplexMatrix::filled(m, m, Complex64::new(f64::MAX, f64::MAX));
    let one = Complex64::new(1.0, 0.0);
    let m_one = Complex64::new(-1.0, 0.0);
    complex_mat_mat_mul(&mut a_v, one, &a, &v)?;
    complex_mat_mat_mul(&mut v_l, one, &v, &lam)?;
    complex_add_matrices(&mut err, one, &a_v, m_one, &v_l)?;
    assert_approx_eq!(complex_matrix_norm(&err, NormMat::Max), 0.0, 1e-15);
    Ok(())
}