[Numpy-discussion] linalg.eig getting the original matrix back ?

[josef.pktd at gmail.com]

On Fri, Jan 15, 2010 at 12:07 PM,  <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

This doesn't seem to be correct. I think, I had an old omega with
multiplicity of eigenvalues in the interpreter. Writing it as a clean
script, I get orthogonal left and right eigenvectors.

Thanks for the reply,

Josef

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