[SciPy-User] Eigenvalue decomposition bug

[Anne Archibald]

On 14 May 2010 14:05, Xavier Saint-Mleux <saintmlx at apstat.com> wrote:
> Sebastian Walter wrote:
>> Hello Pauli,
>> On what kind of matrix did you observe such unstable behavior?
>> Were there repeated eigenvalues?
>>
>
> It happens to me a lot with complex covariance matrices (Hermitian).
> Here's a simple example that returns non-real eigenvalues for an
> Hermitian matrix:

Uh, not to be difficult, but these values are not actually complex.
The complex component is within a floating-point epsilon of zero. The
only way to do better than this is to explicitly notice that the
matrix is Hermitian and branch to special-case code. And if the matrix
is only numerically Hermitian, i.e. values that should be equal differ
by a floating-point epsilon, even this won't help.

Anne

>
>>>> np.random.seed(0)
>>>> x = np.random.random((3,3)) + np.random.random((3,3)) * 1j
>>>> x = (x+x.T.conj())/2 # make it Hermitian
>>>> x == x.T.conj() # ensure it is Hermitian
> array([[ True,  True,  True],
>       [ True,  True,  True],
>       [ True,  True,  True]], dtype=bool)
>>>> np.linalg.eigvals(x) # returns complex values
> array([ 1.99062044 -4.98523579e-17j,  0.18062978 -9.36952928e-19j,
>       -0.23511915 -2.19606549e-17j])
>>>> np.linalg.eigvalsh(x) # imag always zero
> array([-0.23511915+0.j,  0.18062978+0.j,  1.99062044+0.j])
>>>>
>
>
>
> Xavier
>
>
>
>> Sebastian
>>
>>
>>
>>
>>
>> On Tue, May 11, 2010 at 10:39 PM, Pauli Virtanen <pav at iki.fi> wrote:
>>
>>> Tue, 11 May 2010 16:04:00 -0400, Ian Goodfellow wrote:
>>>
>>>> I've find that (scipy/numpy).linalg.eig have a problem where given a
>>>> symmetric matrix they return complex eigenvalues. I can use scipy.io to
>>>> save this matrix in matlab format, load it in matlab, and use matlab's
>>>> eig function to succesfully decompose it with real eigenvalues, so the
>>>> problem seems to be with scipy/numpy or their dependencies, not with my
>>>> matrix. Is this a known issue? And is there a good workaround?
>>>>
>>> Use the eigh function if you know your matrix is symmetric.
>>>
>>> Matlab IIRC checks first if the matrix is symmetric, and if yes, uses a
>>> symmetric-specific eigensolver. Numpy and Scipy don't do this automatic
>>> check.
>>>
>>> A nonsymmetric eigensolver cannot know that your matrix is supposed to
>>> have real eigenvalues, so it's possible some of them explode to complex
>>> pairs because of minuscule numerical error. The imaginary part, however,
>>> is typically small.
>>>
>>> --
>>> Pauli Virtanen
>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>