scipy: generalized eigenvalues are sometimes wrong (on some hardware)

On some machines, in the following example of a singular real 4×4 matrix pencil, SciPy fails to find the generalized eigenvalues 4.0 and 8.0:

>>> import numpy, scipy.linalg
>>> A = numpy.array([[12.0, 28.0, 76.0, 220.0], [16.0, 32.0, 80.0, 224.0], [24.0, 40.0, 88.0, 232.0], [40.0, 56.0, 104.0, 248.0]], dtype='float64')
>>> B = numpy.array([[2.0, 4.0, 10.0, 28.0], [3.0, 5.0, 11.0, 29.0], [5.0, 7.0, 13.0, 31.0], [9.0, 11.0, 17.0, 35.0]], dtype='float64')
>>> D, V = scipy.linalg.eig(A, B); D  # wrong, D should contain 4.0 and 8.0
array([ 4.57142857+0.j,         inf+0.j, 10.98276305+0.j, 10.24256785+0.j])

The computed eigenvectors are linearly dependent and the homogeneous coordinates suggest that there are numerical problems:

>>> scipy.linalg.svdvals(V)  # V is only of rank 1
array([2.00000000e+00, 1.54059899e-15, 3.60419288e-16, 1.01691334e-16])
>>> scipy.linalg.eigvals(A, B, homogeneous_eigvals=True)
array([[5.68434189e-14+0.j, 4.86292581e+00+0.j, 7.31801688e-15+0.j,
        1.18745064e-13+0.j],
       [1.24344979e-14+0.j, 0.00000000e+00+0.j, 6.66318380e-16+0.j,
        1.15932904e-14+0.j]])

I understand that the problem may be ill-conditioned. However, the same kind of computation usually works, so it would be nice to understand why this example sometimes fails.

If one swaps A and B, the eigenvalues are computed correctly:

>>> scipy.linalg.eigvals(B, A)  # correctly finds 1/8 and 1/4
array([ 0.125     +0.j,  0.25      +0.j, -0.10038221+0.j,  0.10494098+0.j])

Over the complex numbers, the computation is correct as well:

>>> A = numpy.array(A, dtype='complex128')
>>> B = numpy.array(B, dtype='complex128')
>>> scipy.linalg.eigvals(A, B)  # correctly finds 4 and 8
array([ 4.        +0.j,  8.        +0.j, 12.16957991+0.j, -4.46153846+0.j])

I have observed this issue on two different machines:

macOS 10.13.6, MacBookPro8,1 Intel® Core™ i5-2415M CPU @ 2.30GHz:
- scipy 1.4.1, numpy 1.18.1, openblas 0.3.7 (all from Homebrew)
- scipy 1.2.0, numpy 1.16.1, openblas 0.3.6 (all from Sage 9.0)
CentOS 7.7, Intel® Xeon® CPU E5-2660 0 @ 2.20GHz:
- scipy 1.2.0, numpy 1.16.1, openblas 0.3.6 (all from Sage 9.0)

On a third and newer machine, the computation returns the correct eigenvalues:

CentOS 7.7, Intel® Xeon® Gold 6152 CPU @ 2.10GHz:
- scipy 1.2.0, numpy 1.16.1, openblas 0.3.6 (all from Sage 9.0)

>>> import numpy, scipy.linalg
>>> A = numpy.array([[12.0, 28.0, 76.0, 220.0], [16.0, 32.0, 80.0, 224.0], [24.0, 40.0, 88.0, 232.0], [40.0, 56.0, 104.0, 248.0]], dtype='float64')
>>> B = numpy.array([[2.0, 4.0, 10.0, 28.0], [3.0, 5.0, 11.0, 29.0], [5.0, 7.0, 13.0, 31.0], [9.0, 11.0, 17.0, 35.0]], dtype='float64')
>>> scipy.linalg.eigvals(A, B)  # correctly finds 4 and 8
array([ 4.        +0.j,  8.        +0.j, 12.16957991+0.j, -4.46153846+0.j])

Moreover, I have tried Octave, Matlab and Mathematica on machine 2, and all of them found the correct eigenvalues.

Any ideas where this problem is coming from? Thank you.

About this issue

Original URL
State: closed
Created 4 years ago
Comments: 23 (13 by maintainers)

Commits related to this issue

TST: linalg: add a regression test for a gen eig problem At the moment, nobody managed to reproduce it, so add a test to see if it shows up again. closes gh-11577 — committed to ev-br/scipy by ev-br 4 months ago
TST: linalg: add a regression test for a gen eig problem At the moment, nobody managed to reproduce it, so add a test to make sure we catch it if it shows up again. closes gh-11577 — committed to ev-br/scipy by ev-br 4 months ago
TST: linalg: add a regression test for a gen eig problem At the moment, nobody managed to reproduce it, so add a test to make sure we catch it if it shows up again. closes gh-11577 — committed to ev-br/scipy by ev-br 4 months ago

Most upvoted comments

To close the loop, here’s a minimal C version of the Fortran reproducer, https://github.com/scipy/scipy/issues/11577#issuecomment-1980332874

$ cat ggev_repro_gh_11577.c

#include<stdio.h>
#include"lapacke.h"

#define n 4

int main()
{
    int lda=n, ldb=n, ldvr=n, ldvl=n, info;
    char jobvl='V', jobvr='V';
    double alphar[n], alphai[n], beta[n];

    double vl[n*n], vr[n*n];
    int lwork = 156;
    double work[156];    /* cheat: 156 is the optimal lwork from the actual lapack call*/

    double a[n*n] = {12.0, 28.0, 76.0, 220.0,
                     16.0, 32.0, 80.0, 224.0,
                     24.0, 40.0, 88.0, 232.0,
                     40.0, 56.0, 104.0, 248.0};

    double b[n*n] = {2.0, 4.0, 10.0, 28.0,
                     3.0, 5.0, 11.0, 29.0,
                     5.0, 7.0, 13.0, 31.0,
                     9.0, 11.0, 17.0, 35.0};

    info = LAPACKE_dggev(LAPACK_ROW_MAJOR, jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr) ; //, work, lwork, info);

    printf("info = %d\n", info);

    printf("Re(eigv) = ");
    for(int i=0; i < n; i++){
        printf("%f , ", alphar[i] / beta[i] );
    }
    printf("\nIm(eigv = ");
    for(int i=0; i < n; i++){
        printf("%f , ", alphai[i] / beta[i] );
    }
    printf("\n");
}

To build it, need to provide the include path for OpenBLAS headers. The location may vary; on my system OpenBLAS comes from mambaforge so it’s <mamba root path>/envs/<conda env name>/include:

$ gcc ggev_repro_gh_11577.c -I/home/br/mambaforge/envs/scipy-dev/include/  -lblas -llapack
$ ./a.out
info = 0
Re(eigv) = 8.000000 , 4.000000 , 7.966917 , -10.146065 , 
Im(eigv = 0.000000 , 0.000000 , 0.000000 , 0.000000 ,

ev-br on Mar 15, 2024