Positive-definiteness of a Matrix instance in sympy

[Yuxiang Wang]

Dear all,

I am currently trying to do solid mechanics (finite deformation) in SymPy. There are a lot of matrices that are positive definite, but I do not know whether there is a way to define this property in SymPy. Any help would be deeply appreciated!

Take this code snippet for example:

```

from sympy import symbols, init_printing, sqrt, simplify, Matrix

init_printing()

# Define a positive-definite real symmetric matrix

b11, b22, b33, b12, b13, b23 = symbols('b11, b22, b33, b12, b13, b23')

b = Matrix([[b11, b12, b13], [b12, b22, b23], [b13, b23, b33]])

J = sqrt(b.det())

lambda1sq, lambda2sq, lambda3sq = b.eigenvals()

lambda1, lambda2, lambda3 = sqrt(lambda1sq), sqrt(lambda2sq), sqrt(lambda3sq)

simplify(lambda1 * lambda2 * lambda3 - J)

```

1) Supposedly, the result of the simplify() should be zero. However it is not, because we are not sure whether the eigenvalues of the matrix b are positive! They are because the matrix should be positive definite, but I do not know how to define that.

2) Also, I know that I could've used real=True in where I defined the symbols of b11, b22, b33.... But I didn't do that. The reason is, if I use real=True, then I cannot get the eigenvalues any more - it will say 

TypeError: cannot determine truth value of-b11*b22*b33 + b11*b23**2 + b12**2*b33 - 2*b12*b13*b23 + b13**2*b22 + 2*(-b11 - b22 - b33)**3/27 - (-b11 - b22 - b33)*(b11*b22 + b11*b33 - b12**2 - b13**2 + b22*b33 - b23**2)/3 < 0

Any way to tighten that up as well?

Thanks,

Shawn

[Aaron Meurer]

On Mon, Feb 2, 2015 at 1:14 PM, Yuxiang Wang <wangyux...@gmail.com> wrote:

Dear all,

I am currently trying to do solid mechanics (finite deformation) in SymPy. There are a lot of matrices that are positive definite, but I do not know whether there is a way to define this property in SymPy. Any help would be deeply appreciated!

Take this code snippet for example:

```

from sympy import symbols, init_printing, sqrt, simplify, Matrix

init_printing()

# Define a positive-definite real symmetric matrix

b11, b22, b33, b12, b13, b23 = symbols('b11, b22, b33, b12, b13, b23')

b = Matrix([[b11, b12, b13], [b12, b22, b23], [b13, b23, b33]])

J = sqrt(b.det())

lambda1sq, lambda2sq, lambda3sq = b.eigenvals()

lambda1, lambda2, lambda3 = sqrt(lambda1sq), sqrt(lambda2sq), sqrt(lambda3sq)

simplify(lambda1 * lambda2 * lambda3 - J)

```

1) Supposedly, the result of the simplify() should be zero. However it is not, because we are not sure whether the eigenvalues of the matrix b are positive! They are because the matrix should be positive definite, but I do not know how to define that.

Maybe refine(lambda1, Q.positive(lambda1sq)).

 

2) Also, I know that I could've used real=True in where I defined the symbols of b11, b22, b33.... But I didn't do that. The reason is, if I use real=True, then I cannot get the eigenvalues any more - it will say 

TypeError: cannot determine truth value of-b11*b22*b33 + b11*b23**2 + b12**2*b33 - 2*b12*b13*b23 + b13**2*b22 + 2*(-b11 - b22 - b33)**3/27 - (-b11 - b22 - b33)*(b11*b22 + b11*b33 - b12**2 - b13**2 + b22*b33 - b23**2)/3 < 0

This error almost always indicates a bug.

Aaron Meurer

 

Any way to tighten that up as well?

Thanks,

Shawn

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/02b195fa-d0ef-4523-95fb-059599803090%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Yuxiang Wang]

Hi Aaron,

Thank you for your response!

1) You are right that it is a bug. If I do this:

from sympy import symbols, sqrt, simplify, Matrix, Q

from sympy.assumptions.assume import global_assumptions

b11, b22, b33, b12, b13, b23 = symbols('b11, b22, b33, b12, b13, b23', real=True)

b = Matrix([[b11, b12, b13], [b12, b22, b23], [b13, b23, b33]])

lambda1sq, lambda2sq, lambda3sq = b.eigenvals()

Then an error would pop up. Instead, if I do this:

from sympy import symbols, sqrt, simplify, Matrix, Q

from sympy.assumptions.assume import global_assumptions

b11, b22, b33, b12, b13, b23 = symbols('b11, b22, b33, b12, b13, b23')

b = Matrix([[b11, b12, b13], [b12, b22, b23], [b13, b23, b33]])

for elem in b:

    global_assumptions.add(Q.real(elem))

lambda1sq, lambda2sq, lambda3sq = b.eigenvals()

It would calculate alright.

2) I still have a problem, as can be viewed more clearly on here:

http://nbviewer.ipython.org/github/yw5aj/ipynb/blob/master/SymbolicOgden.ipynb

A brief description is - the sqrt(matrix determinant) is not equal to the product of the eigenvalues of matrix**(1/2).

Could you please help a little bit in there...? Would what I am doing even be mathematically correct?

Thanks!

Shawn

[Aaron Meurer]

You've discovered that SymPy is really bad at dealing with algebraic expressions.  I would recommend just sticking to what you did the first time, instead of dealing with sqrt(lambda), just use lambda. 

By the way, I don't think simplify() uses the new assumptions. You have to literally use refine() if you want to simplify things with them.

Aaron Meurer

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/182edef1-c430-46f2-b1a1-d92ac343f781%40googlegroups.com.

[Chris Smith]

A brief description is - the sqrt(matrix determinant) is not equal to the product of the eigenvalues of matrix**(1/2).

Could you please help a little bit in there...? Would what I am doing even be mathematically correct?

>>> var('a:d')

(a, b, c, d)

>>> m=Matrix(2,2,[a,b,c,d])

>>> m.det()

a*d - b*c

>>> m.eigenvals()

{a/2 + d/2 - sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1, a/2 + d/2 + sqrt(a**2 - 2*a*d + 4*b*c + d**2)/2: 1}

>>> Mul(*_).simplify()

a*d - b*c

>>> m3=Matrix(3,3,var('x:9'))

>>> d=m3.det()

>>> e=Mul(*m3.eigenvals()).simplify()

>>> d==e

True 

Does someone have the identity wrong? above it appears that det(M) = Mul(*M.eigenvals())

[Yuxiang Wang]

Chris,

Thank you for looking into this! I really appreciate it.

And you are right that det(M) == Mul(*M.eigenvals())

However, if M is positive definite, then sqrt(det(M)) should also equal to sqrt(lambda1)*sqrt(lambda2)*sqrt(lambda3), where lambda1 - 3 are the eigenvalues of the matrix. (Because positive-definiteness would mean det(M)>0, lambda1-3 > 0)

And this is currently where I got stuck...

Shawn

[Yuxiang Wang]

Aaron,

Thanks for your response! I did find a detour by using juse lambda (not taking sqrt). Appreciate your help!

Shawn

[Aaron Meurer]

That's a true fact. However, if the matrix is positive definite, then (by definition) all the eigenvalues are positive, so you should be able to take the square root on both sides. But sqrt(det(M)) != sqrt(eigen1)*sqrt(eigen2)*... in general, due to the fact that sqrt(x*y) != sqrt(x)*sqrt(y) in general.

Aaron Meurer 

--
You received this message because you are subscribed to the Google Groups "sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
To post to this group, send email to sy...@googlegroups.com.
Visit this group at http://groups.google.com/group/sympy.

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/f91aea76-e95a-4dfc-b861-06e4b2ac367f%40googlegroups.com.

[Yuxiang Wang]

Aaron,

Thanks for your response! And sorry for the late reply.

You mentioned: sqrt(det(M)) != sqrt(eigen1) * sqrt(eigen2) * sqrt(eigen3) is True if M is positive definite, so that det(M) is positive, eigen1-3 are positive. But I did use 

global_assumptions.add(Q.positive(J2))
global_assumptions.add(Q.positive(lambda1sq)) 

global_assumptions.add(Q.positive(lambda2sq)) 

global_assumptions.add(Q.positive(lambda3sq))

in the ipynb to declare that both sides of the equations are positive. Does that mean, somewhere along the code path, we didn't check the assumptions (where we should)?

Shawn