shifted Cholesky QR algorithm

[artpelling]

Adds the Cholesky based QR algorithm from this paper with an oblique inner product. Shifts are applied if necessary (see Algorithm 4.1 in the paper).

[artpelling]

I thought a bit about how to add an offset option to this. I think it is not as straight forward as we thought.
Consider the Gramian matrix $X$ of $A$, where the first $r$ columns are orthonormal, then

$$ X=A^TA=\begin{bmatrix}I_r&B^T\\B&Y\end{bmatrix}=R_X^TR_X, $$

where the (upper) Cholesky factor $R_X$ of $X$ can be written as

$$ R_X=\begin{bmatrix}I_r&B^T\\0&R_Z\end{bmatrix}, $$

where $R_Z$ is the (upper) Cholesky factor of the Schur complement $Z=Y-BB^T$.
@sdrave @drittelhacker @pmli Does this make sense to you? Should I move forward and implement this? I don't think it would become a problem but there are no proven error bounds for this case, yet.

[pmli]

@sdrave @drittelhacker @pmli Does this make sense to you? Should I move forward and implement this? I don't think it would become a problem but there are no proven error bounds for this case, yet.

I guess the important part in terms of the (RR)QR interface is that the first offset columns of $A$ do not change. So, if the implementation only using A.inner(A[offset:], product) does not work for some reason (BTW, your math looks right to me), an implementation using A.gramian(product) is also fine as long as it doesn't change A[:offset].

[artpelling]

I've pushed a sketch. This is the best I could come up with so far. Not sure, if the copying and appending of VectorArrays can be circumvented. Also, I don't see a way how we can do the algorithm in place with axpy or scal...

[artpelling]

I am not sure, if the logger should produce a warning, if scipy.linalg.cholesky breaks down - maybe it is too bewildering?

[pmli]

I am not sure, if the logger should produce a warning, if scipy.linalg.cholesky breaks down - maybe it is too bewildering?

Adding more logs might be nice, I guess. It shouldn't have too much, but it will probably have less than gram_schmidt :)

[sdrave]

To make this easier to read, I'd probably:

• merge _shifted_choslesky into the main function
• replace _CholeskyShifter by _compute_shift(A, X, product, product_norm, check_finite)
• replace _ProductNorm by product_norm = None at the beginning of the main loop and

  if product_norm is None:
      product_norm = 1 if product is None else np.sqrt(np.abs(eigs(self.product, k=1)[0][0]))

[artpelling]

To make this easier to read, I'd probably:

* merge `_shifted_choslesky` into the main function

* replace `_CholeskyShifter` by `_compute_shift(A, X, product, product_norm, check_finite)`

* replace `_ProductNorm` by `product_norm = None` at the beginning of the main loop and
  ```python
  if product_norm is None:
      product_norm = 1 if product is None else np.sqrt(np.abs(eigs(self.product, k=1)[0][0]))
  ```

I went ahead and put everything into the main function.

[pmli]

Just minor comments for now. Would be good to have tests (similar to pymortests/algorithms/gram_schmidt).

[sdrave]

I like the code much better now!

Tests are still missing. These should also include testing the offset argument.

Docs build is failing with

/builds/pymor/pymor/docs/source/autoapi/pymor/algorithms/cholesky_qr/index.rst:19: WARNING: Bullet list ends without a blank line; unexpected unindent.

Further, the algorithm does not seem to be tested for complex arrays:

.../pymor/src/pymor/algorithms/cholesky_qr.py:112: ComplexWarning: Casting complex values to real discards the imaginary part
  Rinv = spla.lapack.dtrtri(Rx)[0].T
00:47 shifted_chol_qr: Iteration 2
.../src/pymor/algorithms/cholesky_qr.py:124: ComplexWarning: Casting complex values to real discards the imaginary part
  spla.blas.dtrmm(1, Rx, Ri, overwrite_b=True)
00:47 shifted_chol_qr: Iteration 3
.../pymor/src/pymor/algorithms/cholesky_qr.py:140: ComplexWarning: Casting complex values to real discards the imaginary part
  R[:offset, offset:] = Bi
.../pymor/src/pymor/algorithms/cholesky_qr.py:141: ComplexWarning: Casting complex values to real discards the imaginary part
  R[offset:, offset:] = Ri

[sdrave]

@artpelling, any ideas why the your tests seem to be hanging in CI. Have you tried locally?

[sdrave]

@artpelling, @pmli, tests are really hanging because shift is computed to be 0 as the norm of XX is. This could be mitigated by setting the norm to 1 in this case. However, if the input data is zero, this will still lead to an 'orthonormal' basis of zero vectors. I'd say the only practical approach would be to fail in this case?

[pmli]

However, if the input data is zero, this will still lead to an 'orthonormal' basis of zero vectors. I'd say the only practical approach would be to fail in this case?

Or return Q of zero length and R with zero rows.

[sdrave]

However, if the input data is zero, this will still lead to an 'orthonormal' basis of zero vectors. I'd say the only practical approach would be to fail in this case?

Or return Q of zero length and R with zero rows.

You mean if $A = 0 \in \mathbb{R}^{n \times m}$, then $Q \in \mathbb{R}^{n \times 0}$, $R \in \mathbb{R}^{0 \times m}$?

[pmli]

However, if the input data is zero, this will still lead to an 'orthonormal' basis of zero vectors. I'd say the only practical approach would be to fail in this case?

Or return Q of zero length and R with zero rows.

You mean if $A = 0 \in \mathbb{R}^{n \times m}$, then $Q \in \mathbb{R}^{n \times 0}$, $R \in \mathbb{R}^{0 \times m}$?

Yes, exactly.

[artpelling]

The fixes from @maxbindhak are now merged. @pymor/main I think this is ready for another (final) round of review?

[pmli]

It seems that the tests are failing because of an infinite loop when a zero vector array is given. I guess shifted Cholesky can cheaply check if the input is zero and exit early.

[artpelling]

It seems that the tests are failing because of an infinite loop when a zero vector array is given. I guess shifted Cholesky can cheaply check if the input is zero and exit early.

Ah thats why. I pushed a fix thank you!

[pmli]

It looks like tests fail now because hypothesis is producing rank-1 matrices, e.g.,

from pymor.algorithms.chol_qr import shifted_chol_qr
from pymor.vectorarrays.numpy import NumpyVectorSpace

A = NumpyVectorSpace(5).ones(2)
print(A)

Q, R = shifted_chol_qr(A)
print(Q)
print(R)

produces

[[1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1.]]
00:00 shifted_chol_qr: Iteration 1
00:00 |   |WARNING|shifted_chol_qr: Cholesky factorization broke down! Matrix is ill-conditioned.
00:00 |   shifted_chol_qr: Applying shift: 3.907985046680553e-12
00:00 shifted_chol_qr: Iteration 2
00:00 |   |WARNING|shifted_chol_qr: Cholesky factorization broke down! Matrix is ill-conditioned.
00:00 |   shifted_chol_qr: Applying shift: 3.907985046677496e-14
00:00 shifted_chol_qr: Iteration 3
00:00 |   |WARNING|shifted_chol_qr: Cholesky factorization broke down! Matrix is ill-conditioned.
00:00 |   shifted_chol_qr: Applying shift: 3.9079850466802494e-14
[[4.47213595e-01 4.47213595e-01 4.47213595e-01 4.47213595e-01
  4.47213595e-01]
 [1.08446117e-20 1.08446117e-20 1.08446117e-20 1.08446117e-20
  1.08446117e-20]]
[[2.23606798e+00 2.23606798e+00]
 [0.00000000e+00 1.09243343e-19]]

Fixing this doesn't seem simple to me.

Also, hypothesis found that A = NumpyVectorSpace(5).ones(2) * 1.9e-7j causes an infinite loop.

[artpelling]

It looks like tests fail now because hypothesis is producing rank-1 matrices, e.g.,

Fixing this doesn't seem simple to me.

Indeed, I looked into it but am not sure how to deal with this. The question is what we want the algorithm to do in this case. I might need some help here to develop a strategy.

Also, hypothesis found that A = NumpyVectorSpace(5).ones(2) * 1.9e-7j causes an infinite loop.

This one could be fixed by prescribing a minimum shift.

[pmli]

Fixing this doesn't seem simple to me.

Indeed, I looked into it but am not sure how to deal with this. The question is what we want the algorithm to do in this case. I might need some help here to develop a strategy.

Based on the discussion in the previous community meeting, and since the "correct" approach is unknown (should use unpublished pivoted Cholesky, I guess), shifted Cholesky QR should probably raise an exception if the computed $Q$ is not orthonormal. Then testing should probably be restricted to well-conditioned matrices.

[sdrave]

Regarding testing, I'd suggest to call gram_schmidt in the test with some reasonable defaults to detect if the input array is rank deficient. If it is, one can assert that cholesky_qr raises an error. Otherwise it should work as expected. (With the transition to NiAS, we will probably revamp the unit testing infrasturcture, so I wouldn't invest time into improving the parametrization of the test itself.)

[artpelling]

I hope I fixed the tests now as discussed. @maxbindhak noticed that the computation of the 2-norm is the bottleneck of the method as it stands. In a discussion with @drittelhacker, we realized that we can speed up this computation by exploiting the Hermitian structure of the gramian (using dsyevr/cheevr). I added this as well :)

[artpelling]

I just fixed another error in the shift computation when product is given. This means now that we don't have to raise the rtol in the tests anymore.

@sdrave This PR is ready now from my side.