Opened 11 years ago
Closed 8 months ago
#10332 closed enhancement (fixed)
isPositiveSemiDefinite not accessible
Reported by:  monniaux  Owned by:  jason, was 

Priority:  major  Milestone:  sage9.3 
Component:  linear algebra  Keywords:  psd, semidefinite positive 
Cc:  mjo, ghkliem, ghmwageringel, dimpase, rbeezer  Merged in:  
Authors:  Michael Orlitzky  Reviewers:  Dima Pasechnik 
Report Upstream:  N/A  Work issues:  
Branch:  909c1e1 (Commits, GitHub, GitLab)  Commit:  909c1e19c205fc8741c5d2934f10cd8969c7a507 
Dependencies:  Stopgaps: 
Description
I would like to test rational matrices for semidefinite positiveness. There is such a function in linbox, apparently (isPositiveSemiDefinite, from ispositivesemidefinite.h), but it is not accessible from Sage.
Change History (13)
comment:1 Changed 9 years ago by
 Component changed from linbox to linear algebra
 Owner changed from cpernet to jason, was
comment:2 Changed 13 months ago by
 Cc mjo ghkliem ghmwageringel added
comment:3 Changed 13 months ago by
I'm slowly working on a pivoted LDLT factorization at
http://gitweb.michael.orlitzky.com/?p=sage.d.git;a=blob;f=mjo/ldlt.py
At some point I will have to carefully check (unless someone has a reference) how it can be used to determine positivesemidefiniteness, especially with complex matrices.
comment:4 Changed 13 months ago by
 Branch set to u/mjo/ticket/10332
 Cc rbeezer added
 Commit set to b2b925e81cdfa95475a545d4e2b37bfad337edfb
This is nearing presentability.
My branch has a numerically stable blockLDLT factorization that works for indefinite matrices and in inexact arithmetic. The cited paper shows why we can use it to determine positivesemidefiniteness: each 2x2 diagonal block corresponds to one positive and one negative eigenvalue, and Sylvester's inertia theorem handles the rest.. I've added an is_positive_semidefinite()
method for matrices based on that.
Compared to is_positive_definite()
, the new method...
 Doesn't distinguish between real/complex algorithms. We simply insist that the matrix be Hermitian (which is true of real symmetric matrices, of course), and always use
conjugate_transpose()
. This is not hugely expensive and cleans up the user interface a lot.  Doesn't do any checks on the input matrix. So long as the matrix is Hermitian, you should be okay  but you're in charge of checking that. Running
if not A.is_hermitian()...
on your own is a lot simpler than asking the method to check it for you, and then waiting to catch an exception. And when you want it to be fast, there's no need to disable the checks; it's already fast.
This could conceivably also be used to speed up solve_left
and solve_right
for Hermitian matrices.
comment:5 Changed 9 months ago by
 Status changed from new to needs_review
Since I haven't changed anything in four months I guess this is ready for review.
comment:6 Changed 9 months ago by
 Cc dimpase added
comment:7 Changed 9 months ago by
 Commit changed from b2b925e81cdfa95475a545d4e2b37bfad337edfb to 909c1e19c205fc8741c5d2934f10cd8969c7a507
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
21e0454  Trac #10332: add the Bunch/Kaufman blockLDLT reference.

b4196fd  Trac #10332: add the Higham "Accuracy and Stability..." reference.

8cd9e60  Trac #10332: add block_ldlt() method for matrices.

1143a11  Trac #10332: crossreference block_ldlt <> indefinite_factorization.

909c1e1  Trac #10332: add is_positive_semidefinite() method for matrices.

comment:8 Changed 9 months ago by
Rebased onto develop.
comment:9 Changed 9 months ago by
matrix2.pyx is long overdue for splitting, it's 670K !!! (unless we aim for Guinness Book of Records, of course).
comment:10 Changed 9 months ago by
 Reviewers set to Dima Pasechnik
 Status changed from needs_review to positive_review
otherwise, OK.
comment:11 Changed 9 months ago by
Thanks!
comment:12 Changed 9 months ago by
 Milestone set to sage9.3
comment:13 Changed 8 months ago by
 Branch changed from u/mjo/ticket/10332 to 909c1e19c205fc8741c5d2934f10cd8969c7a507
 Resolution set to fixed
 Status changed from positive_review to closed
I went down a rabbit hole with the cholesky inverse and wound up here. The inverse of a positivedefinite matrix can actually be computed a tiiiiny bit faster through the
indefinite_factorization
method, which returns a naive (suitable only in exact arithmetic) LDLT factorization.If we modify the
indefinite_factorization
to pivot (put the zeros at the bottom ofD
), we could also perform an LDLT factorization for positiveSEMIdefinite matrices, making the goal of this ticket attainable.Note also that some people have asked for square roots of PSD matrices in #23424 and #25464. Right now that's only implemented for positivedefinite matrices, due to the aforementioned limitation of the LDLT factorization, so we could potentially improve that too.