rust-ndarray / ndarray-linalg

Linear algebra package for rust-ndarray using LAPACK binding
Other
376 stars 77 forks source link

Eigh eigenvectors are inverted #368

Open pdcook opened 1 year ago

pdcook commented 1 year ago

Related: #307

MWE:

Cargo.toml:

[package]
name = "EighBug"
version = "0.1.0"
edition = "2021"

# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html

[dependencies]
# ndarray
ndarray = "0.15.6"
# ndarray_linalg
ndarray-linalg = { version = "0.16.0", features = ["openblas-static"] }
# num
num = "0.4.1"
# rustfmt
rustfmt = "0.10.0"

main.rs:

use ndarray::*;
use ndarray_linalg::*;
use num::complex::Complex64;

fn main() {
    let A = arr2(&[
        [
            Complex64::new(3., 0.),
            Complex64::new(2., 2.),
            Complex64::new(2., 2.),
        ],
        [
            Complex64::new(2., -2.),
            Complex64::new(3., 0.),
            Complex64::new(2., 2.),
        ],
        [
            Complex64::new(2., -2.),
            Complex64::new(2., -2.),
            Complex64::new(3., 0.),
        ],
    ]);

    println!("A = \n{}", A);

    let (eigenvalues, eigenvectors) = A.eigh(UPLO::Lower).unwrap();

    println!(
        "eigenvalues = \n{}",
        eigenvalues.mapv(|x| format!("{:.3}", x))
    );

    println!(
        "eigenvectors = \n{}",
        eigenvectors.mapv(|x| format!("{:.3}", x))
    );

    println!(
        "A @ eigenvectors = \n{}",
        A.dot(&eigenvectors).mapv(|x| format!("{:.3}", x))
    );

    println!(
        "eigenvectors @ eigenvalues = \n{}",
        eigenvectors
            .dot(&Array2::from_diag(
                &eigenvalues.mapv(|x| Complex64::from(x))
            ))
            .mapv(|x| format!("{:.3}", x))
    );

    println!("==============================");
    println!("WITH REVERSED ROWS");

    let eigenvectors = eigenvectors.slice(s![..;-1, ..]).to_owned();
    println!(
        "eigenvectors = \n{}",
        eigenvectors.mapv(|x| format!("{:.3}", x))
    );

    println!(
        "A @ eigenvectors = \n{}",
        A.dot(&eigenvectors).mapv(|x| format!("{:.3}", x))
    );

    println!(
        "eigenvectors @ eigenvalues = \n{}",
        eigenvectors
            .dot(&Array2::from_diag(
                &eigenvalues.mapv(|x| Complex64::from(x))
            ))
            .mapv(|x| format!("{:.3}", x))
    );
}

Output

A = 
[[3+0i, 2+2i, 2+2i],
 [2-2i, 3+0i, 2+2i],
 [2-2i, 2-2i, 3+0i]]
eigenvalues = 
[-1.000, 1.536, 8.464]
eigenvectors = 
[[-0.577+0.000i, -0.577+0.000i, -0.577+0.000i],
 [-0.000+0.577i, 0.500-0.289i, -0.500-0.289i],
 [0.577-0.000i, -0.289+0.500i, -0.289-0.500i]]
A @ eigenvectors = 
[[-1.732+2.309i, -1.732+0.845i, -1.732-3.155i],
 [0.000+4.041i, -1.232+0.711i, -2.232-1.289i],
 [1.732+2.309i, -1.598+1.077i, -3.598+0.077i]]
eigenvectors @ eigenvalues = 
[[0.577+0.000i, -0.887+0.000i, -4.887+0.000i],
 [0.000-0.577i, 0.768-0.443i, -4.232-2.443i],
 [-0.577+0.000i, -0.443+0.768i, -2.443-4.232i]]
==============================
WITH REVERSED ROWS
eigenvectors = 
[[0.577-0.000i, -0.289+0.500i, -0.289-0.500i],
 [-0.000+0.577i, 0.500-0.289i, -0.500-0.289i],
 [-0.577+0.000i, -0.577+0.000i, -0.577+0.000i]]
A @ eigenvectors = 
[[-0.577+0.000i, -0.443+0.768i, -2.443-4.232i],
 [0.000-0.577i, 0.768-0.443i, -4.232-2.443i],
 [0.577+0.000i, -0.887-0.000i, -4.887+0.000i]]
eigenvectors @ eigenvalues = 
[[-0.577+0.000i, -0.443+0.768i, -2.443-4.232i],
 [0.000-0.577i, 0.768-0.443i, -4.232-2.443i],
 [0.577+0.000i, -0.887+0.000i, -4.887+0.000i]]

Expected output

A @ eigenvectors should equal eigenvectors @ eigenvalues without needing to invert the row ordering.


This could very likely be a misunderstanding on my end. But it may be related to the reversal of UPLO::Upper and UPLO::Lower. For example, for the array:

A = [[3+0i, 2+2i, 2+2i],
     [   0, 3+0i, 2+2i],
     [   0,    0, 3+0i]]

eigh(UPLO::upper) returns [3, 3, 3] for the eigenvalues, which is the expected result for eigh(UPLO::lower).