QR factorization still producing negative diagonals?

[markusm73]

If I understand the release notes of LAPACK 3.2 correctly, QR factorization should always produce non-negative diagonals in matrix R since that version:

http://www.netlib.org/lapack/lapack-3.2.html#4

But I still see negative diagonals in R after calling DGEQRF. Multiplying Q as obtained via DORGQR and the previously copied R yields the initial matrix A again so the factorization seems otherwise correct. Is there anything I may need to do differently?

[Julien Langou]

Hi @markusm73,

We have reverted this. Yes for LAPACK 3.2, GEQRF was producing nonnegative diagonal. But we discovered a (very bad) bug. This was a numerical bug, and this was not easy to fix. We are not 100% sure of our fix, so we moved QR with nonnegative diagonal to GEQRFP (have a look!), and restored GEQRF as it was pre-3.2.

So to repeat, in 3.8,
GEQRF + or - on the diagonal of R
GEQRFP, only + (or 0) on the diagonal of R

Julien.

[markusm73]

Thanks for the fast reply, Julien, I missed that reversion. I guess I'll stick with GEQRF and manually negate Q columns and R rows as required, which is probably the next best solution.

[Julien Langou]

Question: what is your application? Why do you care about the element of R being nonnegative? (We are looking for example of application of GERFP, so it would be good for us to know why you care. We have a few examples, but not many to be fair.)

Comment: Yes, if you have Q and R, please just use GEQRF and +/- as appropriate. GEQRFP is when you want to compute Q^T x without forming Q.

[markusm73]

It's a machine learning application (Gaussian Processes) that needed to calculate the log marginal likelihood. By summing the logarithms of the entries along the diagonal of R, we obtain a log determinant, which contributes to the likelihood. Obviously, negative entries along the diagonal would cause problems with logarithms.

[markusm73]

A more general justification for producing positive entries along the diagonal of R is that R in "[Q, R] = qr(A)" corresponds to the Cholesky factor "R = chol(A'*A)". If A has full rank, the entries of the two Rs have to be exactly the same modulo signs, but Cholesky always produces positive diagonals. Forming the inner product "A'*A" just to obtain the Cholesky factorization is numerically unstable, which is why I think the QR factorization should directly produce a "sane" Cholesky factor for the inner product.

[Julien Langou]

Sounds good Markus. Thanks for the comments.

If all is needed to satisfy the user's need is achieved by a simple post processing, so a +/- on R (and on Q, if Q is needed), then I do not feel too bad to leave this up to the user. (This is not that we do not want to do it, this is that we do not know how to do it, we tried, and we failed and had to roll back the code.) In LAPACK case, we do not return Q. We return the Householder vectors. Call them V. If we return R only, yes +/- to achieve nonnegative diagonal in R is easy. If we return Q and R, yes +/- to achieve nonnegative diagonal in R is easy. We return V and R and it is not easy to return a nonnegative diagonal R. This is what GEQRFP does and we had a hard time to get it right and we are not 100% of the numerical stability. We return V and R. We want V and R to be consistent. And getting R with a nonnegative diagonal has complicated implication for V. So this is where we are. I hope this makes sense. So we know of a few applications where users really need V and R, and really need nonnegative diagonal for R. For these applications GEQRFP is really needed. These applications cannot be handle by a simple post processing. The two applications in mind: (a) when doing ``triangle on top of triangle`` QR factorization, (so take R1 and R2 return R factors of R1 on top of R2) the binary operation becomes commutative if we have nonnegative diagonal R in the QR step, (B) this has to do with random distribution of orthogonal matrices. (How to get a random orthogonal matrix.) And now, that I write both, I am not even sure if GEQRFP cannot be bypass actually.

Julien.

[markusm73]

Thanks for the insightful explanation of the problem. Post-processing of R to obtain a positive diagonal is not a big deal, of course. It's definitely surprising that the intermediate representation based on Householder vectors causes problems with numerical stability if you want to enforce a positive diagonal.