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

[josef.pktd@gmai...]

On Fri, Jan 15, 2010 at 11:32 AM, Sebastian Walter
<sebastian.walter@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@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@scipy.org
>> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>