is_positive_definite returns wrong results

[mwageringel]

For some double dense matrices that are not positive definite, Sage incorrectly claims that they are:

```
sage: A = matrix(CDF, [[1+I]])
sage: A.is_positive_definite()  # should be False
True
sage: A.is_positive_semidefinite()  # correct
False
```

Clearly, A is not even Hermitian. The implementation relies on testing whether a Cholesky decomposition exists or not, but .cholesky() gives an incorrect result for this matrix

```
sage: A.cholesky()
[1.0]
```

while for most other matrices this would raise an error.

The SciPy documentation is not entirely clear about this, but it seems that SciPy assumes that the matrix is already Hermitian positive-definite and thus only considers half the matrix (see SciPy issues 2116 and 15012).

This ticket adds a check whether the matrix is Hermitian before computing the Cholesky decomposition, which resolves the problem. The default algorithm for is_hermitian is changed from "orthonormal" to "naive".

Depends on #33107

CC: @orlitzky

Component: linear algebra

Author: Michael Orlitzky

Branch/Commit: 1d6bae6

Reviewer: Markus Wageringel

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

[mwageringel]

comment:1

This also fails for larger matrices, such as these linear combinations of positive definite matrices (for many seeds):

sage: set_random_seed(0)
sage: random_pd_matrix = lambda n: (B := matrix.random(CDF, n)) * B.H
sage: r = 2  # 2 or larger
sage: n = 3  # 1 or larger
sage: As = [random_pd_matrix(n) for j in range(r)]
sage: C = random_pd_matrix(r)
sage: A = sum(sum(C[i,j] for i in range(r)) * As[j] for j in range(r))
sage: A.is_positive_definite()  # should usually be False
True
sage: A.is_positive_semidefinite()  # correct, but pd would imply psd
False
sage: A.is_hermitian(tol=1e-6)  # correct, but pd would imply Hermitian
False

[orlitzky]

comment:2

My first thought was that we should be able to delete cholesky() and is_positive_definite() in the subclass to obtain something possibly slower, but that actually works. There are two small problems though.

First, the generic cholesky() doesn't work on the trivial matrix, breaking one doctest (#33107). This is an easy fix.

Second, the subclass implementation of is_hermitian() used by the superclass's is_cholesky() doesn't work out-of-the box. For example,

sage: M = matrix(RDF,[[ 1,  1,    1,     1,     1],                 
....:                 [ 1,  5,   31,   121,   341],                 
....:                 [ 1, 31,  341,  1555,  4681],                 
....:                 [ 1,121, 1555,  7381, 22621],                 
....:                 [ 1,341, 4681, 22621, 69905]])
M.is_hermitian(algorithm='orthonormal') # the default
False

while

sage: M.is_hermitian(algorithm='orthonormal', tol=1e-10)
True

or

sage: M.is_hermitian(algorithm='naive')
True

We'd have to either tweak the default tolerance there, or again maybe just delete the subclass method. The "naive" algorithm now defers to a superclass implementation where the default tolerance works.

[orlitzky]

Dependencies: #33107

[orlitzky]

Commit: c70330b

[orlitzky]

comment:3

Rather than remove the subclass is_hermitian(), I've only switched its default to the "naive" algorithm. (But note that I don't think it's any more naive than "orthonormal"). Afterwards, eliminating cholesky() and is_positive_definite() from the subclass allows the superclass methods to return the correct result.

---

New commits:

89c46c6	Trac #33107: cholesky() for trivial matrices.
38c86e5	Trac #33100: switch default RDF/CDF is_hermitian algorithm to "naive".
f6e4b64	Trac #33100: drop special RDF/CDF cholesky() and is_positive_definite().
c70330b	Trac #33100: add doctests for the example on the ticket.

[orlitzky]

Author: Michael Orlitzky

[orlitzky]

Branch: u/mjo/ticket/33100

[mwageringel]

comment:4

I agree with using the "naive" algorithm here and do not see why it would be more naive/less appropriate than the "orthonormal" one. For the orthonormal algorithm, no tolerance would be large enough to work for all such matrices that are Hermitian in an exact sense, when considered entrywise. Scipy 1.8.0 will also have a function scipy.linalg.ishermitian which works like the naive algorithm, i.e. entrywise.

Replying to @orlitzky:

My first thought was that we should be able to delete cholesky() and is_positive_definite() in the subclass to obtain something possibly slower, but that actually works.

I am a bit hesitant to replace a Scipy implementation by a pure Sage implementation. Upstream has confirmed here that scipy.linalg.cholesky essentially first constructs a Hermitian matrix by only considering the given triangular information of the matrix, ignoring the other half. So if we first check whether the matrix is Hermitian in .cholesky(), we will obtain a valid Cholesky decomposition.

This would make is_positive_definite work as intended by the implementation. Though, it would in most cases produce the same result as is_positive_semidefinite, as a Cholesky decomposition can exist for semidefinite matrices, but for matrices that are actually singular this can easily fail. Is there an intended difference between these methods? Perhaps debatable, but from a usage point of view I would expect something like

• is_positive_definite = all eigenvalues strictly positive (or > tolerance)
• is_positive_semidefinite = all eigenvalues > -tolerance (implemented in a numerically sane way).

[orlitzky]

comment:5

Replying to @mwageringel:

I agree with using the "naive" algorithm here and do not see why it would be more naive/less appropriate than the "orthonormal" one. For the orthonormal algorithm, no tolerance would be large enough to work for all such matrices that are Hermitian in an exact sense, when considered entrywise. Scipy 1.8.0 will also have a function scipy.linalg.ishermitian which works like the naive algorithm, i.e. entrywise.

Replying to @orlitzky:

My first thought was that we should be able to delete cholesky() and is_positive_definite() in the subclass to obtain something possibly slower, but that actually works.

I am a bit hesitant to replace a Scipy implementation by a pure Sage implementation. Upstream has confirmed here that scipy.linalg.cholesky essentially first constructs a Hermitian matrix by only considering the given triangular information of the matrix, ignoring the other half. So if we first check whether the matrix is Hermitian in .cholesky(), we will obtain a valid Cholesky decomposition.

Yeah, the subclass should be doing that anyway, since the superclass does it. (Users shouldn't see different qualitative behaviors from a method when you change the underlying representation of the matrix's entries.) I haven't checked but scipy is almost certainly faster even though I did my best to make it fast in sage. On the other hand, it's nice not to have copy & pasted the entire method's documentation.

This would make is_positive_definite work as intended by the implementation. Though, it would in most cases produce the same result as is_positive_semidefinite, as a Cholesky decomposition can exist for semidefinite matrices, but for matrices that are actually singular this can easily fail. Is there an intended difference between these methods? Perhaps debatable, but from a usage point of view I would expect something like

• is_positive_definite = all eigenvalues strictly positive (or > tolerance)
• is_positive_semidefinite = all eigenvalues > -tolerance (implemented in a numerically sane way).

The superclass cholesky() and is_*definite() methods support both exact and inexact rings in a numerically stable way (caveat: is_hermitian() itself is not yet reliable, that's #33053) via block_ldlt(). For inexact rings, positive-definite and positive-semidefinite are indeed basically the same up to machine precision, but the tests use > 0 and >= 0.

[mwageringel]

comment:6

Replying to @orlitzky:

I haven't checked but scipy is almost certainly faster even though I did my best to make it fast in sage. On the other hand, it's nice not to have copy & pasted the entire method's documentation.

Unless you strongly prefer otherwise, I tend towards staying with using Scipy for the Cholesky decomposition.

The superclass cholesky() and is_*definite() methods support both exact and inexact rings in a numerically stable way (caveat: is_hermitian() itself is not yet reliable, that's #33053) via block_ldlt(). For inexact rings, positive-definite and positive-semidefinite are indeed basically the same up to machine precision, but the tests use > 0 and >= 0.

Ok, dropping is_positive_definite from the subclass sounds good to me then, as then the is_*definite methods are at least based on the same implementation (as numerically the distinction between these methods is not so relevant).

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from c70330b to 2bf92ff

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

64d9028	Trac #33107: return immutable trivial Cholesky factors.
f8abd61	Trac #33100: switch default RDF/CDF is_hermitian algorithm to "naive".
e5b681a	Trac #33100: check is_hermitian() within cholesky() over RDF/CDF.
2bf92ff	Trac #33100: add doctests for the example on the ticket.

[orlitzky]

comment:8

I left is_positive_definite() alone because cholesky() and is_positive_definite() cache one another in the subclass. Using the superclass for is_positive_definite() would cache a block_ldlt() factorization, but then it would not be used for cholesky().

[mwageringel]

Reviewer: Markus Wageringel

[mwageringel]

comment:9

Ok, this looks good to me then. Thank you for implementing this.

[mwageringel]

Description changed:

--- 
+++ 
@@ -17,5 +17,6 @@

 while for most other matrices this would raise an error.

-The [SciPy documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky.html) is not entirely clear about this, but it seems that SciPy assumes that the matrix is already Hermitian positive-definite and thus only considers half the matrix (see [SciPy issues 2116](https://github.com/scipy/scipy/issues/2116#issuecomment-990487066) and [15012](https://github.com/scipy/scipy/issues/15012)). Therefore, Sage needs to check at least whether the matrix is Hermitian before calling `scipy.linalg.cholesky`, but it is not clear from the Scipy documentation whether that is enough to determine positive-definiteness.
+The [SciPy documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.cholesky.html) is not entirely clear about this, but it seems that SciPy assumes that the matrix is already Hermitian positive-definite and thus only considers half the matrix (see [SciPy issues 2116](https://github.com/scipy/scipy/issues/2116#issuecomment-990487066) and [15012](https://github.com/scipy/scipy/issues/15012)).

+This ticket adds a check whether the matrix is Hermitian before computing the Cholesky decomposition, which resolves the problem. The default algorithm for `is_hermitian` is changed from "orthonormal" to "naive".

[mkoeppe]

comment:11

Merge failure on top of:

429528ecd0 Trac #20343: Adding sage/plot/tikzpicture.py

d3f74f63be Trac #33085: adjust dochtml label so doctests pass when html doc is not built/installed

d70747d13f Trac #33094: update URLs in SPKG.rst for ECL

03f57e6ad1 Trac #32922: Change parameter name

1057236215 Trac #30999: expose method hyperbolicity in graphs

9b4a9198e2 Trac #33035: random doctest failure in src/sage/tests/book_stein_ent.py

fee8e65057 Trac #32264: Convenience syntax for quaternion ideals

c2bcc90ea8 Trac #31005: Add traceback information to exceptions during docbuild

3f697607e8 Trac #33159: wrong result solving equation system with symbolic matrix

bf9673f1a5 Trac #33158: kRegularSequence: minimization wrong

01ff471e9f Trac #33157: Issue with SR → [RC]BF conversion

379c88d769 Trac #33154: random doctest failure in sage.plot.matrix_plot

9fa34b02a7 Trac #33151: sage-conf_pypi does not build wheelhouse

ffe0b98c81 Trac #33148: fix doctest in quantum_group_gap

0bf6a86003 Trac #33107: Generic cholesky() fails on the trivial matrix

3b02b82340 Trac #33103: gitpod integration using Docker images from portability testing workflow

af5c4c48ad Trac #33098: Add more void packages to distros/void.txt

9b61187847 Trac #33096: src/tox.ini (relint): Exclude editor temp files etc.

217156c781 Trac #33064: sage_docbuild: fails when cache cannot be saved

c2097f4654 Trac #33039: Random doctest failure in continued_fraction_gosper

53781ca1e5 Trac #33015: Random failure in number_field_element.pyx

0ecfd27027 Trac #33005: Add feature for pdftocairo

961f91527a Trac #32970: fix parent of 0th Bernoulli polynomial

527fbc473d Trac #32877: Remove superfluous set_random_seed() calls

d3aec55fe5 Trac #32873: sage.features, sage_setup: Replace use of distutils.errors by setuptools

3154d1d97d Trac #32856: Get rid of "# optional - build"

d87eeb3e7f Trac #30362: Add symplectic structures

030af8db7e Trac #29865: Modularization of sagelib: Break out separate packages sagemath-objects, sagemath-categories

1fed3d635c Trac #33189: Make tests pass with arb 2.22

1ea6801ffb Trac #33131: Installation manual: Add decision tree, remove mention of Sage-mirror-hosted binary distributions

209b71330d Trac #33101: lrslib: fix doctest in game_theory/parser.py

73aa4f3e4d Trac #33077: pari-jupyter: Reinstate

180b9ee2fb Trac #30933: GH Actions: Repair upload of docker images to docker.pkg.github.com

e67e1ee810 Trac #8450: intermediate complex expression in real functions make many plot functions fail

d03870b5f2 Trac #33206: PDF documentation links in Documentation from Jupyter notebook are broken

8ea92d580a Updated SageMath version to 9.5.rc3

merge was not clean: conflicts in src/sage/matrix/matrix_double_dense.pyx

[mkoeppe]

comment:12

(brought to you by https://github.com/sagemath/git-trac-command/pull/57)

[mwageringel]

comment:13

For some reason, I did not receive an email notification for comment 11. Also, this does not seem so useful without knowing which of the three dozen tickets causes a merge conflict, so I have to wait for the new beta anyway. It would be more helpful if the CI identifies merge conflicts between positively reviewed tickets, so that one can resolve conflicts before it even gets any attention by the release manager.

[mkoeppe]

comment:14

Thanks for the input! Yes, this can certainly be refined.

[slel]

comment:15

Set milestone to sage-9.6 after Sage 9.5 release.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 2bf92ff to 137da8a

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

a1bc2c9	Trac #33107: cholesky() for trivial matrices.
ab8ba9b	Trac #33107: return immutable trivial Cholesky factors.
7756167	Trac #33100: switch default RDF/CDF is_hermitian algorithm to "naive".
b2150ef	Trac #33100: check is_hermitian() within cholesky() over RDF/CDF.
d32ee08	Trac #33100: add doctests for the example on the ticket.
30c5ce9	Trac #33100: use a more reliable abs() replacement in is_hermitian().
137da8a	Trac #33100: allow nonpositive tolerance for RDF/CDF is_hermitian().

[mwageringel]

comment:17

-                # avoid abs() which is missing for finite fields.
-                if d > tolerance or -d > tolerance:
+                # abs() support is spotty
+                if (d*d.conjugate()).sqrt() > tolerance:

I am not completely sure about theses changes. Square roots do not always exist in the same field, such as cyclotomic fields for which sqrt returns a result in SR. Granted, this probably did not work before either, as comparisons are broken in number fields (#20028).

For finite fields, the documentation advises against relying on comparisons as well. Though, technically they should work in our case, i.e. if the integer representation of a finite field element is compared to a tolerance of 0.

It might be better to explicitly check whether d is equal to zero when the tolerance is zero.

[orlitzky]

comment:18

Replying to @mwageringel:

I am not completely sure about theses changes. Square roots do not always exist in the same field, such as cyclotomic fields for which sqrt returns a result in SR. Granted, this probably did not work before either, as comparisons are broken in number fields (#20028).

This lurking doubt is what kept me from setting it back to needs_review right away. A priori we could square both sides to avoid roots, but I still have some questions. For the moment, a lot of this stuff is broken in corner cases anyway, so I don't feel too bad about making incremental but imperfect improvements so long as they don't hurt any existing use cases.

But for finite fields I think you're right, and I vaguely remember coming to the same conclusion. If the fields are exact, the tolerance should be zero unless you override it, and it's only by accident that the comparison works but it does work. If a tolerance is given, though, I'm not even sure what the user could expect to happen.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

9dc1cb2	Trac #33107: cholesky() for trivial matrices.
9cef1c9	Trac #33107: return immutable trivial Cholesky factors.
5d36b63	Trac #33100: switch default RDF/CDF is_hermitian algorithm to "naive".
6439830	Trac #33100: check is_hermitian() within cholesky() over RDF/CDF.
b553988	Trac #33100: add doctests for the example on the ticket.
c4977c4	Trac #33100: use a more reliable abs() replacement in is_hermitian().
e5e8a0f	Trac #33100: allow nonpositive tolerance for RDF/CDF is_hermitian().
3056eea	Trac #33100: more reliable tolerance tests in _is_hermitian().
7cac4e4	Trac #33100: combine sparse and dense _is_hermitian() branches.
b166f0b	Trac #33100: access sparse nonzero entries direcly in _is_hermitian().

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 137da8a to b166f0b

[orlitzky]

comment:20

Let's see what patchbot says about this one...

[mwageringel]

comment:21

There is a merge failure, but otherwise this looks good to me, as long as the tests pass.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

94ab5f1	Trac #33100: switch default RDF/CDF is_hermitian algorithm to "naive".
9e77b84	Trac #33100: check is_hermitian() within cholesky() over RDF/CDF.
1817400	Trac #33100: add doctests for the example on the ticket.
f9547e2	Trac #33100: use a more reliable abs() replacement in is_hermitian().
4b31a25	Trac #33100: allow nonpositive tolerance for RDF/CDF is_hermitian().
61fb01c	Trac #33100: more reliable tolerance tests in _is_hermitian().
32d9a91	Trac #33100: combine sparse and dense _is_hermitian() branches.
1d6bae6	Trac #33100: access sparse nonzero entries direcly in _is_hermitian().

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from b166f0b to 1d6bae6

[orlitzky]

comment:23

It rebased cleanly, I think I maybe it wanted me to drop the commits from #33107.

[mwageringel]

comment:24

Thanks.

[vbraun]

Changed branch from u/mjo/ticket/33100 to 1d6bae6