Skip to content

Eigh eigenvectors are inverted #368

Open
@pdcook

Description

@pdcook

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).

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions