verify whether a matrix is positive definite or not
Hi all, I have a parameter-dependent matrix B(x) = B_0 + x B_1, 0 \le x \le 1 where B_0 and B_1 are symmetric. How can I determine critical values x* (if any) such that B(x*) is not positive definite ? from scipy import * def B(x): return array(([[11.,8.],[8.,7.]])) - x*array(([[20.,1.],[1.,26]])) X = linspace(0,1,100) for x in X: print x L=linalg.cholesky(B(x),lower=1) I mean it would be nice if cholesky could return info=1 if the matrix is not spd. The current behaviour is Traceback (most recent call last): File "test_spd.py", line 11, in ? L=linalg.cholesky(B(x),lower=1) File "/usr/lib64/python2.4/site-packages/scipy/linalg/decomp.py", line 552, in cholesky if info>0: raise LinAlgError, "matrix not positive definite" numpy.linalg.linalg.LinAlgError: matrix not positive definite Helpful suggestions would be appreciated. Nils
On 28/06/07, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
I have a parameter-dependent matrix
B(x) = B_0 + x B_1, 0 \le x \le 1
where B_0 and B_1 are symmetric. How can I determine critical values x* (if any) such that B(x*) is not positive definite ?
from scipy import *
def B(x):
return array(([[11.,8.],[8.,7.]])) - x*array(([[20.,1.],[1.,26]]))
X = linspace(0,1,100)
for x in X: print x L=linalg.cholesky(B(x),lower=1)
I mean it would be nice if cholesky could return info=1 if the matrix is not spd. The current behaviour is
Traceback (most recent call last): File "test_spd.py", line 11, in ? L=linalg.cholesky(B(x),lower=1) File "/usr/lib64/python2.4/site-packages/scipy/linalg/decomp.py", line 552, in cholesky if info>0: raise LinAlgError, "matrix not positive definite" numpy.linalg.linalg.LinAlgError: matrix not positive definite
Helpful suggestions would be appreciated.
Well, you can always use try/except to catch the LinAlgError. It's remotely possible that cholesky might fail to converge and throw a different LinAlgError which you would want to re-raise. You can also look at the eigenvalues - the matrix is positive definite if and only if they're all positive. So making a function that takes a parameter and returns the least eigenvalue should give you a relatively smooth function to do root-finding on. With symmetric matrices, eigenvalue finding ought to be fairly reliable. For this particular case, note that the positive-definite matrices form a cone, that is, the sum of two or any positive multiple of a positive-definite matrix is also positive definite. In particular this means it's convex. So if you're tracing the line between two endpoints, as you are here (the endpoints are B_0 and B_0+B_1), you can check the endpoints and know that the matrix is positive definite between them if they're both positive definite. If one is positive definite and the other isn't, then clearly there's some point in between where they stop being positive definite. If neither is positive definite, then it's possible that some positive definite matrix lies between them (good luck finding it; you could try numerically maximizing the least eigenvalue). Anne M. Archibald
Hi Nils, is this not similar to the eigenvalue problems we discussed off-list? You can numerically find the smallest x (if that's what you want) by solving minimize x s.t. Bo - B1*x >= 0 In CVXOPT you can solve it as follows: from cvxopt.base import matrix from cvxopt.solvers import sdp from cvxopt.lapack import syev c = matrix(-1.0) B0 = matrix([ [11.,8.], [8.,7.] ]) B1 = matrix([ [20.,1.] ,[1.,26.] ]) Gs = [ B1[:] ] hs = [ B0 ] sol = sdp(c, Gs=Gs, hs=hs) x = sol['x'] v = matrix([0., 0.]) syev(B0 - x*B1, v) print v - Joachim On 6/28/07, Nils Wagner <nwagner@iam.uni-stuttgart.de> wrote:
Hi all,
I have a parameter-dependent matrix
B(x) = B_0 + x B_1, 0 \le x \le 1
where B_0 and B_1 are symmetric. How can I determine critical values x* (if any) such that B(x*) is not positive definite ?
from scipy import *
def B(x):
return array(([[11.,8.],[8.,7.]])) - x*array(([[20.,1.],[1.,26]]))
X = linspace(0,1,100)
for x in X: print x L=linalg.cholesky(B(x),lower=1)
I mean it would be nice if cholesky could return info=1 if the matrix is not spd. The current behaviour is
Traceback (most recent call last): File "test_spd.py", line 11, in ? L=linalg.cholesky(B(x),lower=1) File "/usr/lib64/python2.4/site-packages/scipy/linalg/decomp.py", line 552, in cholesky if info>0: raise LinAlgError, "matrix not positive definite" numpy.linalg.linalg.LinAlgError: matrix not positive definite
Helpful suggestions would be appreciated.
Nils
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.org http://projects.scipy.org/mailman/listinfo/scipy-user
Joachim Dahl wrote:
Hi Nils,
is this not similar to the eigenvalue problems we discussed off-list?
You can numerically find the smallest x (if that's what you want) by solving minimize x s.t. Bo - B1*x >= 0
In CVXOPT you can solve it as follows:
from cvxopt.base import matrix from cvxopt.solvers import sdp from cvxopt.lapack import syev
c = matrix(-1.0) B0 = matrix([ [11.,8.], [8.,7.] ]) B1 = matrix([ [20.,1.] ,[1.,26.] ]) Gs = [ B1[:] ] hs = [ B0 ]
sol = sdp(c, Gs=Gs, hs=hs) x = sol['x']
v = matrix([0., 0.]) syev(B0 - x*B1, v) print v
- Joachim
Hi Joachim, Yes indeed. It is related to my previous problem. Thank you very much for your solution ! Lieven sent me a randomly generated matrix pair (A(x), B(x)) off-list. However, the matrix B(x) corresponds to the mass matrix in structural dynamics (my background), which is almost always positive definite. Hence I was confused by the solution of his example for maximizing the smallest eigenvalue of (A(x), B(x)) subjected to some constraints. Anyway I found another way to detect the "border" by bisection. Thanks to Anne ! I haven't tested the code on other examples, so there could be mistakes. Nils S.M. Rump. Verification of Positive Definiteness. BIT Numerical Mathematics, 46:433–452, 2006.
participants (3)
-
Anne Archibald -
Joachim Dahl -
Nils Wagner