[Numpy-discussion] linalg.eig getting the original matrix back ?
Warren Weckesser
warren.weckesser at enthought.com
Fri Jan 15 12:24:10 EST 2010
For the case where all the eigenvalues are simple, this works for me:
In [1]: import numpy as np
In [2]: a = np.array([[1.0, 2.0, 3.0],[2.0, 3.0, 0.0], [3.0, 0.0, 4.0]])
In [3]: eval, evec = np.linalg.eig(a)
In [4]: eval
Out[4]: array([-1.51690942, 6.24391817, 3.27299125])
In [5]: a2 = np.dot(evec, eval[:,np.newaxis] * evec.T)
In [6]: np.allclose(a, a2)
Out[6]: True
In [7]:
Warren
josef.pktd at gmail.com wrote:
> On Fri, Jan 15, 2010 at 11:32 AM, Sebastian Walter
> <sebastian.walter at gmail.com> wrote:
>
>> numpy.linalg.eig guarantees to return right eigenvectors.
>> evec is not necessarily an orthonormal matrix when there are
>> eigenvalues with multiplicity >1.
>> For symmetrical matrices you'll have mutually orthogonal eigenspaces
>> but each eigenspace might be spanned by
>> vectors that are not orthogonal to each other.
>>
>> Your omega has eigenvalue 1 with multiplicity 3.
>>
>
> Yes, I thought about the multiplicity. However, even for random
> symmetric matrices, I don't get the result
> I change the example matrix to
> omega0 = np.random.randn(20,8)
> omega = np.dot(omega0.T, omega0)
> print np.max(np.abs(omega == omega.T))
>
> I have been playing with left and right eigenvectors, but I cannot
> figure out how I could compose my original matrix with them either.
>
> I checked with wikipedia, to make sure I remember my (basic) linear algebra
> http://en.wikipedia.org/wiki/Eigendecomposition_(matrix)#Symmetric_matrices
>
> The left and right eigenvectors are almost orthogonal
> ev, evecl, evecr = sp.linalg.eig(omega, left=1, right=1)
>
>>>> np.abs(np.dot(evecl.T, evecl) - np.eye(8))>1e-10
>>>> np.abs(np.dot(evecr.T, evecr) - np.eye(8))>1e-10
>>>>
>
> shows three non-orthogonal pairs
>
>
>>>> ev
>>>>
> array([ 6.27688862, 8.45055356, 15.03789945, 19.55477818,
> 20.33315408, 24.58589363, 28.71796764, 42.88603728])
>
>
> I always thought eigenvectors are always orthogonal, at least in the
> case without multiple roots
>
> I had assumed that eig will treat symmetric matrices in the same way as eigh.
> Since I'm mostly or always working with symmetric matrices, I will
> stick to eigh which does what I expect.
>
> Still, I'm currently not able to reproduce any of the composition
> result on the wikipedia page with linalg.eig which is puzzling.
>
> Josef
>
>
>>
>>
>> On Fri, Jan 15, 2010 at 4:31 PM, <josef.pktd at gmail.com> wrote:
>>
>>> I had a problem because linal.eig doesn't rebuild the original matrix,
>>> linalg.eigh does, see script below
>>>
>>> Whats the trick with linalg.eig to get the original (or the inverse)
>>> back ? None of my variations on the formulas worked.
>>>
>>> Thanks,
>>> Josef
>>>
>>>
>>> import numpy as np
>>> import scipy as sp
>>> import scipy.linalg
>>>
>>> omega = np.array([[ 6., 2., 2., 0., 0., 3., 0., 0.],
>>> [ 2., 6., 2., 3., 0., 0., 3., 0.],
>>> [ 2., 2., 6., 0., 3., 0., 0., 3.],
>>> [ 0., 3., 0., 6., 2., 0., 3., 0.],
>>> [ 0., 0., 3., 2., 6., 0., 0., 3.],
>>> [ 3., 0., 0., 0., 0., 6., 2., 2.],
>>> [ 0., 3., 0., 3., 0., 2., 6., 2.],
>>> [ 0., 0., 3., 0., 3., 2., 2., 6.]])
>>>
>>> for fun in [np.linalg.eig, np.linalg.eigh, sp.linalg.eig, sp.linalg.eigh]:
>>> print fun.__module__, fun
>>> ev, evec = fun(omega)
>>> omegainv = np.dot(evec, (1/ev * evec).T)
>>> omegainv2 = np.linalg.inv(omega)
>>> omegacomp = np.dot(evec, (ev * evec).T)
>>> print 'composition',
>>> print np.max(np.abs(omegacomp - omega))
>>> print 'inverse',
>>> print np.max(np.abs(omegainv - omegainv2))
>>>
>>> this prints:
>>>
>>> numpy.linalg.linalg <function eig at 0x017EDDF0>
>>> composition 0.405241032278
>>> inverse 0.405241032278
>>>
>>> numpy.linalg.linalg <function eigh at 0x017EDE30>
>>> composition 3.5527136788e-015
>>> inverse 7.21644966006e-016
>>>
>>> scipy.linalg.decomp <function eig at 0x01DB14F0>
>>> composition 0.238386662463
>>> inverse 0.238386662463
>>>
>>> scipy.linalg.decomp <function eigh at 0x01DB1530>
>>> composition 3.99680288865e-015
>>> inverse 4.99600361081e-016
>>> _______________________________________________
>>> NumPy-Discussion mailing list
>>> NumPy-Discussion at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>>
>>>
>> _______________________________________________
>> NumPy-Discussion mailing list
>> NumPy-Discussion at scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>
More information about the NumPy-Discussion
mailing list