Hi all, I am confused by the results of the following test from numpy import set_printoptions from numpy.random import seed, rand from scipy.linalg import eig from numpy import dot, diag from scipy.sparse.linalg import eigs from scipy.linalg import norm set_printoptions(precision=6) seed(10) n=5 A = rand(n,n)+ 1j*rand(n,n) B = rand(n,n)+ 1j*rand(n,n) w,vl,vr = eig(A,B,left=True, right=True) print 'Eigenvalues by eig' print w k = n-2 ws,vs = eigs(A,k,B) for i in range(k): r =dot(A,vr[:,i])-w[i]*dot(B,vr[:,i]) res = norm(r) print 'Residual', res print 'Eigenvalues by eigs' k = n-2 ws,vs = eigs(A,k,B) print ws for i in range(k): r =dot(A,vs[:,i])-ws[i]*dot(B,vs[:,i]) res = norm(r) print 'Residual', res Eigenvalues by eig [ 2.093178-1.783532j -0.202628+1.195944j 1.186895+0.56595j 0.014805-0.287354j 0.772942-0.266609j] Residual 1.9173960652e-15 Residual 2.20392640899e-15 Residual 3.21653472875e-15 Eigenvalues by eigs [ 19.064654+27.897369j 15.797549-27.544265j -0.249578-30.673974j] Residual 2.39453499962e-14 Residual 7.33698764801e-14 Residual 1.76193714267e-14 The eigenvalues returned by eig and eigs are completely different. Any idea ? I am using
scipy.__version__ '0.11.0.dev-600e81f'
Nils
On Thu, 20 Oct 2011 19:34:03 +0200 Pauli Virtanen <pav@iki.fi> wrote:
20.10.2011 18:53, Nils Wagner kirjoitti: [clip]
The eigenvalues returned by eig and eigs are completely different.
The matrix B must be hermitian positive (semi-)definite for `eigs`.
-- Pauli Virtanen
It means that there is no check for B inside eigs. IMHO, a warning should be raised if B is not hermitian positive (semi-)definite. Nils
(20.10.2011 19:57), Nils Wagner wrote: [clip]
It means that there is no check for B inside eigs. IMHO, a warning should be raised if B is not hermitian positive (semi-)definite.
That could be useful. Checking PD is may be more expensive since it requires trying to do a Cholesky decomposition. Would need some benchmarks to check whether it matters. It's also possible to do the check only for dense matrices. Scipy doesn't have a sparse Cholesky at the moment, and moreover, the linear operator can be an arbitrary function with no way to obtain the transpose. -- Pauli Virtanen
On Thu, Oct 20, 2011 at 4:16 PM, Pauli Virtanen <pav@iki.fi> wrote:
(20.10.2011 19:57), Nils Wagner wrote: [clip]
It means that there is no check for B inside eigs. IMHO, a warning should be raised if B is not hermitian positive (semi-)definite.
That could be useful.
Checking PD is may be more expensive since it requires trying to do a Cholesky decomposition. Would need some benchmarks to check whether it matters.
It's also possible to do the check only for dense matrices. Scipy doesn't have a sparse Cholesky at the moment, and moreover, the linear operator can be an arbitrary function with no way to obtain the transpose.
none of the scipy numpy eigh functions do a check. Since I use them mostly for gram, covariance matrices, where I already know it's symmetric, I wouldn't like any expensive checks, in linalg, I don't use sparse so far. user responsibility to check the doc string ? Josef
-- Pauli Virtanen
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
I'm not really an authority on these matters, but I agree with Josef. In the dense case you would have to do a Cholesky decomposition, and we have to assume that it does not break down until the last step. As far as I know, the cost for this is O(n^3). On the other hand, computing a few eigenvalues should be O(n^2) if you do not want them all (in which case you would be using eig anyway). In addition to what Pauli wrote about the sparse case, even if there was a sparse (incomplete) Cholesky decomposition implemented you could not really use it. To be sure that the matrix really was positive definite you would essentially have to set all tolerances to zero and compute the full factorization. This would (possibly) result in huge amounts of fill-in. Not good. So yes, please leave it as it is. Kind regards, Tony Stillfjord On Thu, Oct 20, 2011 at 10:39 PM, <josef.pktd@gmail.com> wrote:
(20.10.2011 19:57), Nils Wagner wrote: [clip]
It means that there is no check for B inside eigs. IMHO, a warning should be raised if B is not hermitian positive (semi-)definite.
That could be useful.
Checking PD is may be more expensive since it requires trying to do a Cholesky decomposition. Would need some benchmarks to check whether it matters.
It's also possible to do the check only for dense matrices. Scipy doesn't have a sparse Cholesky at the moment, and moreover, the linear operator can be an arbitrary function with no way to obtain the
On Thu, Oct 20, 2011 at 4:16 PM, Pauli Virtanen <pav@iki.fi> wrote: transpose.
none of the scipy numpy eigh functions do a check.
Since I use them mostly for gram, covariance matrices, where I already know it's symmetric, I wouldn't like any expensive checks, in linalg, I don't use sparse so far.
user responsibility to check the doc string ?
Josef
-- Pauli Virtanen
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
_______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
On Thu, 20 Oct 2011 16:39:28 -0400 josef.pktd@gmail.com wrote:
On Thu, Oct 20, 2011 at 4:16 PM, Pauli Virtanen <pav@iki.fi> wrote:
(20.10.2011 19:57), Nils Wagner wrote: [clip]
It means that there is no check for B inside eigs. IMHO, a warning should be raised if B is not hermitian positive (semi-)definite.
That could be useful.
Checking PD is may be more expensive since it requires trying to do a Cholesky decomposition. Would need some benchmarks to check whether it matters.
It's also possible to do the check only for dense matrices. Scipy doesn't have a sparse Cholesky at the moment, and moreover, the linear operator can be an arbitrary function with no way to obtain the transpose.
none of the scipy numpy eigh functions do a check.
Since I use them mostly for gram, covariance matrices, where I already know it's symmetric, I wouldn't like any expensive checks, in linalg, I don't use sparse so far.
user responsibility to check the doc string ?
Josef
from scipy.sparse.linalg import eigs help (eigs)
scipy.__version__ '0.10.0b2'
The docstring of eigs is misleading in that context. M must represent a real symmetric matrix. For best results, M should be of the same type as A. Additionally: * If sigma==None, M is positive definite * If sigma is specified, M is positive semi-definite If sigma==None, eigs requires an operator to compute the solution of the linear equation `M * x = b`. This is done internally via a (sparse) LU decomposition for an explicit matrix M, or via an iterative solver for a general linear operator. Alternatively, the user can supply the matrix or operator Minv, which gives x = Minv * b = M^-1 * b Nils
participants (4)
-
josef.pktd＠gmail.com -
Nils Wagner -
Pauli Virtanen -
Tony Stillfjord