On Fri, Jan 15, 2010 at 12:24 PM, Warren Weckesser <warren.weckes...@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.p...@gmail.com wrote: >> On Fri, Jan 15, 2010 at 11:32 AM, Sebastian Walter >> <sebastian.wal...@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.p...@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 >>> >>> >> _______________________________________________ >> 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 > _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion