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.