Cholesky decomposition incorrect solution

[Andrea Fresu]

Hello,
I have this symmetric matrix:

A=Matrix([[0, 0.623433638694170, 1, 5.42879559377412], [0.623433638694170, 0, 1, 6.81321272473753], [1, 1, 0, 0], [5.42879559377412, 6.81321272473753, 0, 0]])

I did

A.cholesky()

getting:

⎡ 0      0       0    0  ⎤

⎢                           ⎥

⎢zoo   zoo    0    0  ⎥

⎢                           ⎥

⎢zoo  nan    nan  0 ⎥

⎢                           ⎥  

⎣zoo  nan  nan  nan⎦

What can I do to get the correct result?

Thanks for the time,

Andrea Fresu

[Ondřej Čertík]

It looks like a bug. You can use numpy as follows:

In [17]: a = matrix2numpy(A, dtype=float)

In [18]: numpy.linalg.cholesky(a)
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
<ipython-input-18-0367e248ee93> in <module>()
----> 1 numpy.linalg.cholesky(a)

/local/certik/bld/profile/o3tqau5cbqzb/lib/python2.7/site-packages/numpy/linalg/linalg.pyc
in cholesky(a)
529 results = lapack_routine(_L, n, a, m, 0)
530 if results['info'] > 0:
--> 531 raise LinAlgError('Matrix is not positive definite - '
532 'Cholesky decomposition cannot be computed')
533 s = triu(a, k=0).transpose()

LinAlgError: Matrix is not positive definite - Cholesky decomposition
cannot be computed


So it looks like your matrix is not positive definite.

Ondrej

> Thanks for the time,
>
> Andrea Fresu
>

> --
> 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/18051c80-1ea2-4b6e-a7ad-e9710a4f4d7c%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Andrea Fresu]

Oh! It is my fault!
So I suggest to insert an error message doing A.cholesky() when A is not positive definite, like in numpy :)
However, thank you very much Ondrej.
Good job!

[Chris Smith]

It would be interesting to see how the positive definite is known. Page https://en.wikipedia.org/wiki/Positive-definite_matrix gives an example for which v.T*M*v gives 2*a**2 - 2*a*b + 2*b**2 - 2*b*c + 2*c**2 which is factored into a sum of squares (and hence nonnegative). I don't know if we have a way of factoring such expressions into squares, however.

[Aaron Meurer]

There's a useful mathematical fact here, which is that a multivariate
polynomial with rational coefficients is nonnegative over the reals if
and only if it can be written as a sum of squares of rational
functions (see https://en.wikipedia.org/wiki/Hilbert%27s_seventeenth_problem).
Rational is important here, as there are positive polynomials that are
not the sum of squares of polynomials (there's an example of such on
the Wikipedia page).

What I don't know, and what that Wikipedia page is quiet about, is if
there is known a general algorithm to either compute such a sum of
squares or prove that none exists (or even the weaker problem of just
detecting if one exists or not, which is what we really care about
here).

Another interesting related read is
http://andrescaicedo.wordpress.com/2008/11/11/275-positive-polynomials/.

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/fe87cb3a-721b-4ff6-a188-368ba56f0a9f%40googlegroups.com.