isPositiveSemiDefinite not accessible

[sagetrac-monniaux]

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.

CC: @orlitzky @kliem @mwageringel @dimpase @rbeezer

Component: linear algebra

Keywords: psd, semidefinite positive

Author: Michael Orlitzky

Branch/Commit: 909c1e1

Reviewer: Dima Pasechnik

Issue created by migration from https://trac.sagemath.org/ticket/10332

[orlitzky]

comment:2

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.

[orlitzky]

comment:3

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.

[orlitzky]

Branch: u/mjo/ticket/10332

[orlitzky]

Author: Michael Orlitzky

[orlitzky]

Commit: b2b925e

[orlitzky]

comment:4

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.

[orlitzky]

comment:5

Since I haven't changed anything in four months I guess this is ready for review.

[sagetrac-git]

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

21e0454	Trac #10332: add the Bunch/Kaufman block-LDLT reference.
b4196fd	Trac #10332: add the Higham "Accuracy and Stability..." reference.
8cd9e60	Trac #10332: add block_ldlt() method for matrices.
1143a11	Trac #10332: cross-reference block_ldlt <-> indefinite_factorization.
909c1e1	Trac #10332: add is_positive_semidefinite() method for matrices.

[sagetrac-git]

Changed commit from b2b925e to 909c1e1

[orlitzky]

comment:8

Rebased onto develop.

[dimpase]

comment:9

matrix2.pyx is long overdue for splitting, it's 670K !!! (unless we aim for Guinness Book of Records, of course).

[dimpase]

comment:10

otherwise, OK.

[dimpase]

Reviewer: Dima Pasechnik

[orlitzky]

comment:11

Thanks!

[vbraun]

Changed branch from u/mjo/ticket/10332 to 909c1e1