add tolerance to is_hermitian() for sparse "double" matrices

[orlitzky]

Discovered in #33023, the implementation of is_hermitian() for sparse double matrices is the generic one that tests equality of entries, and can fail due to numerical imprecision.

We already have an is_hermitian() for DENSE double matrices, but it is not used for the sparse class, which is its sibling and not a child. A quick way to fix this is to implement the "naive" version from the sense class in the sparse class. But a better solution might be to move the "naive" algorithm for dense matrices up into the superclass, and just allow the generic method to take a tolerance. If that's stable, it would allow us to avoid reimplementing is_hermitian() three times.

When this has been fixed, the # known bug tags need to be removed from the following tests in src/sage/matrix/matrix_double_sparse.pyx:

sage: L = A.cholesky()                    # known bug, 33031        
sage: (A - L*L.T).norm(1) < 1e-10         # known bug, 33031        
True                                                                
sage: B = A.change_ring(RR)               # known bug, 33031        
sage: (B.cholesky() - L).norm(1) < 1e-10  # known bug, 33031

sage: L = A.cholesky()                    # known bug, 33031        
sage: (A - L*L.H).norm(1) < 1e-10         # known bug, 33031        
True                                                                
sage: B = A.change_ring(CC)               # known bug, 33031        
sage: (B.cholesky() - L).norm(1) < 1e-10  # known bug, 33031

Depends on #33023

CC: @collares @kliem

Component: linear algebra

Author: Michael Orlitzky

Branch/Commit: 4c14a9e

Reviewer: Dima Pasechnik

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

[orlitzky]

Description changed:

--- 
+++ 
@@ -2,3 +2,20 @@

 We already have an `is_hermitian()` for DENSE double matrices, but it is not used for the sparse class, which is its sibling and not a child. A quick way to fix this is to implement the "naive" version from the sense class in the sparse class. But a better solution might be to move the "naive" algorithm for dense matrices up into the superclass, and just allow the generic method to take a tolerance. If that's stable, it would allow us to avoid reimplementing `is_hermitian()` three times.

+When this has been fixed, the `# known bug` tags need to be removed from the following tests in `src/sage/matrix/matrix_double_sparse.pyx`:
+
+```
+sage: L = A.cholesky()                    # known bug, 33031        
+sage: (A - L*L.T).norm(1) < 1e-10         # known bug, 33031        
+True                                                                
+sage: B = A.change_ring(RR)               # known bug, 33031        
+sage: (B.cholesky() - L).norm(1) < 1e-10  # known bug, 33031
+```
+
+```
+sage: L = A.cholesky()                    # known bug, 33031        
+sage: (A - L*L.H).norm(1) < 1e-10         # known bug, 33031        
+True                                                                
+sage: B = A.change_ring(CC)               # known bug, 33031        
+sage: (B.cholesky() - L).norm(1) < 1e-10  # known bug, 33031
+```

[dimpase]

comment:2

This also needs to take into account cases with entries such as RBF and CBF, RIF amd CIF.

[orlitzky]

comment:3

Replying to @dimpase:

This also needs to take into account cases with entries such as RBF and CBF, RIF amd CIF.

is_hermitian() already fails on RBF and RIF matrices for want of a conjugate() method...

[orlitzky]

comment:4

Here's a proof-of-concept that attempts to fix sparse RDF/CDF without breaking anything else.

---

New commits:

295508e	Trac #33023: disable flaky cholesky() tests.
a2359e7	Trac #33031: add tolerance to _is_hermitian().
7ff6363	Trac #33031: use superclass for naive RBF/CDF is_hermitian().
dd5b12d	Trac #33031: add tolerance for sparse RDF/CDF is_hermitian().
d202b67	Trac #33031: speed up is_hermitian() for sparse matrices.
1dea31b	Trac #33031: fix is_hermitian() over finite fields.

[orlitzky]

Branch: u/mjo/ticket/33031

[orlitzky]

Author: Michael Orlitzky

[orlitzky]

Commit: 1dea31b

[orlitzky]

Dependencies: #33023

[orlitzky]

Work Issues: actually fix is_hermitian() and not just cholesky()

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

Changed commit from 1dea31b to 0f04c34

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

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

0752485	Trac #33031: add tolerance for sparse RDF/CDF is_hermitian().
95714d1	Trac #33031: speed up is_hermitian() for sparse matrices.
9339252	Trac #33031: fix is_hermitian() over finite fields.
0f04c34	Trac #33031: add tolerance for sparse RDF/CDF is_skew_hermitian().

[orlitzky]

Changed work issues from actually fix is_hermitian() and not just cholesky() to none

[orlitzky]

comment:8

I've also fixed is_skew_hermitian(), since doing so is trivial. First things first: can someone confirm that my tolerance of 1e-12 is large enough to avoid the issue on aarch64 (or wherever the problem was observed)?

This is only a minimal refactoring to address the immediate problem. I've opened #33053 to track the larger issue.

[collares]

comment:9

The test runs a little further, but fails after the change_ring call:

sage: n = ZZ.random_element(1,5)
....: A = matrix.random(CDF, n, sparse=True)
....: I = matrix.identity(CDF, n, sparse=True)
....: A = A*A.conjugate_transpose() + I
....: L = A.cholesky()
....: assert (A - L*L.H).norm(1) < 1e-10
....: B = A.change_ring(CC)
....: assert (B.cholesky() - L).norm(1) < 1e-10
....: 
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-193612c9eade> in <module>
      6 assert (A - L*L.H).norm(Integer(1)) < RealNumber('1e-10')
      7 B = A.change_ring(CC)
----> 8 assert (B.cholesky() - L).norm(Integer(1)) < RealNumber('1e-10')

/nix/store/zqyv54wwjaal14k1skabdk9ch6rr6v6r-python3-3.9.6-env/lib/python3.9/site-packages/sage/matrix/matrix2.pyx in sage.matrix.matrix2.Matrix.cholesky (build/cythonized/sage/matrix/matrix2.c:84977)()
  12631 
  12632         if not self.is_hermitian():
> 12633             raise ValueError("matrix is not Hermitian")
  12634 
  12635         # Use classical=True to ensure that we don't get a permuted L.

ValueError: matrix is not Hermitian
sage: A
[ 3.229039802593019 + 1.6322627083148353e-17*I  -0.021330037612755842 - 1.3824367515761828*I    0.18142736612755045 + 0.5825316936751774*I]
[ -0.021330037612755842 + 1.3824367515761828*I 4.0201724676288855 - 1.6208942612402586e-17*I    -0.2251207850381418 - 1.1919564519572095*I]
[   0.18142736612755045 - 0.5825316936751775*I    -0.2251207850381418 + 1.1919564519572095*I  1.9658297004722771 - 2.185885234400264e-17*I]
sage: B
[3.22903980259302 + 1.63226270831484e-17*I  -0.0213300376127558 - 1.38243675157618*I   0.181427366127550 + 0.582531693675177*I]
[ -0.0213300376127558 + 1.38243675157618*I 4.02017246762889 - 1.62089426124026e-17*I   -0.225120785038142 - 1.19195645195721*I]
[  0.181427366127550 - 0.582531693675178*I   -0.225120785038142 + 1.19195645195721*I 1.96582970047228 - 2.18588523440026e-17*I]

This does not seem to be due to the tolerance, though, because (B - B.H).norm(1) is 1.5474000715052094e-16.

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

Branch pushed to git repo; I updated commit sha1. New commits:

4c14a9e

Trac #33031: compare sparse RDF/CDF cholesky() answers against dense ones.

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

Changed commit from 0f04c34 to 4c14a9e

[orlitzky]

comment:11

Replying to @collares:

The test runs a little further, but fails after the change_ring call:

Argh, of course, you're now hitting the same problem over RR and CC, which we use to check the answer. Thankfully those aren't crucial to the process. I just pushed another commit that checks the answer against the dense RDF/CDF factorization instead, which will use a different but supposedly stable is_hermitian(), but the same cholesky() that we had with RR/CC.

[collares]

comment:12

matrix_double_sparse.pyx passes now. Thanks!

[dimpase]

comment:13

I would have kept skew = False as the default.

-    def _is_hermitian(self, skew = False):
+    def _is_hermitian(self, skew, tolerance):

[dimpase]

comment:14

maybe a default for tolerance may be set, too.

[orlitzky]

comment:15

Replying to @dimpase:

I would have kept skew = False as the default.

-    def _is_hermitian(self, skew = False):
+    def _is_hermitian(self, skew, tolerance):

This is mainly a style choice since it's an internal function and users won't be annoyed by the extra parameter. I think it's kind of nice to see the symmetry/contrast between uses of the internal function in the user-facing ones. i.e.:

# in the superclass
def is_hermitian(self):
    return self._is_hermitian(skew=False, tolerance=0)
def is_skew_hermitian(self):
    return self._is_hermitian(skew=True, tolerance=0)
# in the subclass
def is_hermitian(self, tolerance=1e-12):
    return self._is_hermitian(skew=False, tolerance=tolerance)
def is_skew_hermitian(self, tolerance=1e-12):
    return self._is_hermitian(skew=True, tolerance=tolerance)

The user-facing method always knows the value of the skew parameter, so there's never any real need to "default" it. Doing so would only save a total of five characters. The tolerance parameter is similar... a default of zero would work now, but would be superfluous when we finish making the tolerance (in the superclass) user-facing in #33053.

[dimpase]

comment:16

ok

[dimpase]

Reviewer: Dima Pasechnik

[vbraun]

Changed branch from u/mjo/ticket/33031 to 4c14a9e