Re: [Numpy-discussion] linalg.eig getting the original matrix back ?
On Fri, Jan 15, 2010 at 12:24 PM, Warren Weckesser 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 >> 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, 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 composition 0.405241032278 inverse 0.405241032278 numpy.linalg.linalg composition 3.5527136788e-015 inverse 7.21644966006e-016 scipy.linalg.decomp composition 0.238386662463 inverse 0.238386662463 scipy.linalg.decomp 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-D
Re: [Numpy-discussion] linalg.eig getting the original matrix back ?
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.p...@gmail.com wrote: > On Fri, Jan 15, 2010 at 11:32 AM, Sebastian Walter > 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, 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 >>> composition 0.405241032278 >>> inverse 0.405241032278 >>> >>> numpy.linalg.linalg >>> composition 3.5527136788e-015 >>> inverse 7.21644966006e-016 >>> >>> scipy.linalg.decomp >>> composition 0.238386662463 >>> inverse 0.238386662463 >>> >>> scipy.linalg.decomp >>> 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
Re: [Numpy-discussion] linalg.eig getting the original matrix back ?
On Fri, Jan 15, 2010 at 12:07 PM, wrote: > On Fri, Jan 15, 2010 at 11:32 AM, Sebastian Walter > 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, 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 >>> composition 0.405241032278 >>> inverse 0.405241032278 >>> >>> numpy.linalg.linalg >>> composition 3.5527136788e-015 >>> inverse 7.21644966006e-016 >>> >>> scipy.linalg.decomp >>> composition 0.238386662463 >>> inverse 0.238386662463 >>> >>> scipy.linalg.decomp >>> 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
Re: [Numpy-discussion] linalg.eig getting the original matrix back ?
On Fri, Jan 15, 2010 at 11:32 AM, Sebastian Walter 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, 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 >> composition 0.405241032278 >> inverse 0.405241032278 >> >> numpy.linalg.linalg >> composition 3.5527136788e-015 >> inverse 7.21644966006e-016 >> >> scipy.linalg.decomp >> composition 0.238386662463 >> inverse 0.238386662463 >> >> scipy.linalg.decomp >> 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
Re: [Numpy-discussion] linalg.eig getting the original matrix back ?
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. On Fri, Jan 15, 2010 at 4:31 PM, 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 > composition 0.405241032278 > inverse 0.405241032278 > > numpy.linalg.linalg > composition 3.5527136788e-015 > inverse 7.21644966006e-016 > > scipy.linalg.decomp > composition 0.238386662463 > inverse 0.238386662463 > > scipy.linalg.decomp > 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
