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: sage-9.3
Component: linear algebra Keywords: psd, semidefinite positive
Cc: mjo, gh-kliem, gh-mwageringel, 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:

Status badges

Description

I would like to test rational matrices for semidefinite positiveness. There is such a function in linbox, apparently (isPositiveSemiDefinite, from is-positive-semidefinite.h), but it is not accessible from Sage.

Change History (13)

comment:1 Changed 9 years ago by jdemeyer

  • Component changed from linbox to linear algebra
  • Owner changed from cpernet to jason, was

comment:2 Changed 13 months ago by mjo

  • Cc mjo gh-kliem gh-mwageringel added

I went down a rabbit hole with the cholesky inverse and wound up here. The inverse of a positive-definite 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 of D), we could also perform an LDLT factorization for positive-SEMIdefinite 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 positive-definite matrices, due to the aforementioned limitation of the LDLT factorization, so we could potentially improve that too.

comment:3 Changed 13 months ago by mjo

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 positive-semidefiniteness, especially with complex matrices.

comment:4 Changed 13 months ago by mjo

  • Authors set to Michael Orlitzky
  • Branch set to u/mjo/ticket/10332
  • Cc rbeezer added
  • Commit set to b2b925e81cdfa95475a545d4e2b37bfad337edfb

This is nearing presentability.

My branch has a numerically stable block-LDLT factorization that works for indefinite matrices and in inexact arithmetic. The cited paper shows why we can use it to determine positive-semidefiniteness: 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...

  1. 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.
  2. 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 mjo

  • 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 dimpase

  • Cc dimpase added

comment:7 Changed 9 months ago by git

  • Commit changed from b2b925e81cdfa95475a545d4e2b37bfad337edfb to 909c1e19c205fc8741c5d2934f10cd8969c7a507

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

21e0454Trac #10332: add the Bunch/Kaufman block-LDLT reference.
b4196fdTrac #10332: add the Higham "Accuracy and Stability..." reference.
8cd9e60Trac #10332: add block_ldlt() method for matrices.
1143a11Trac #10332: cross-reference block_ldlt <-> indefinite_factorization.
909c1e1Trac #10332: add is_positive_semidefinite() method for matrices.

comment:8 Changed 9 months ago by mjo

Rebased onto develop.

comment:9 Changed 9 months ago by dimpase

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 dimpase

  • Reviewers set to Dima Pasechnik
  • Status changed from needs_review to positive_review

otherwise, OK.

comment:11 Changed 9 months ago by mjo

Thanks!

comment:12 Changed 9 months ago by chapoton

  • Milestone set to sage-9.3

comment:13 Changed 8 months ago by vbraun

  • Branch changed from u/mjo/ticket/10332 to 909c1e19c205fc8741c5d2934f10cd8969c7a507
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.