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

[josef.pktd at gmail.com]

On Fri, Jan 15, 2010 at 12:24 PM, Warren Weckesser
<warren.weckesser at enthought.com> wrote:
> 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
>

Thanks, I thought I had tried similar versions, but I guess not with
the matrix without multiplicity of eigenvalues

>>> np.max(np.abs(np.dot(evecl, (ev * evecl).T)-omega))
3.6415315207705135e-014
>>> np.max(np.abs(np.dot(evecr, (ev * evecr).T)-omega))
2.6256774532384952e-014
>>> np.max(np.abs(np.dot(evecl, np.dot(np.diag(ev), evecl.T))-omega))
3.6415315207705135e-014

So, the my confusion was just because eig doesn't treat multiple
eigenvalues in the same way as eigh.

Josef


>
>
> 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
>>
>
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
>