Fix Matrix(QRSparseQ) constructor

[andreasnoack]

Fixes JuliaLang/LinearAlgebra.jl#505

[KlausC]

I propose to leave SuiteSparse unchanged and change the constructor of Matrix as noted in my comment.

[andreasnoack]

To add some comments to the change here. The Q in QR of A is actually square because it is a product of elementary (Householderish) reflectors. Consequencly in the tall case, we actually don't have that A = Q*R but instead A = Q*[R; 0]. It is customary to decompose Q into two column blocks Q = [Q0 Q1] such that A = Q0*R, i.e. Q0 has size(A, 2) columns, and in many cases it is useful just to work with Q0. When Q is stored as elementary reflectors it can be interpreted as both Q and Q0 and that is indeed what we currently do for *:

julia> A = randn(4,2);

julia> F = qrfact(A);

julia> F.Q*Matrix(I,2,2)
4×2 Array{Float64,2}:
 -0.481088  -0.70301
  0.404963  -0.240095
  0.397627  -0.668398
  0.668171   0.037106

julia> F.Q*Matrix(I,3,2)
ERROR: DimensionMismatch("first dimension of matrix must have size either 4 or 2")
Stacktrace:
 [1] *(::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}, ::Array{Bool,2}) at /Users/andreasnoack/julia-dev/usr/share/julia/site/v0.7/LinearAlgebra/src/qr.jl:607
 [2] top-level scope

julia> F.Q*Matrix(I,4,2)
4×2 Array{Float64,2}:
 -0.481088  -0.70301
  0.404963  -0.240095
  0.397627  -0.668398
  0.668171   0.037106

You can argue if that is a good idea but if one of the interpretations should go, it would definately be the Q0 interpretation. When converting to Matrix, we'd have to make a choice and we've decided to return Q0 because it saves a lot of memory since it is much smaller than Q. The confusing part is that size(Q, 2) != size(Matrix(Q), 2).

Regarding the implementation then we'd need a way to get the size of Q0 from the data stored in Q type. In the dense case that is fairly simple because it is the same as the number of elementary reflectors, i.e. size(F.factors, 2) but in the sparse case that need not be the case since

1. In contrast to LAPACK, it doesn't store identity reflectors
2. It estimates rank so when the trailing matrix is considered zero, no more reflectors are produced.

As a result the number of columns in QRSparseQ.factors corresponds to the minimum number of elementary reflectors needed for A ≈ Q*[R; 0] for the specified tolerance. For the test matrix that you provided

A = sparse([0.0 1 0 0; 0 0 0 0])

the number of reflectors needed is zero. When A is rank deficient that means that we can decompose R = [R00 R01; 0 0]. Notice that this means that part of Q0, as it was defined above, is actually redundant which is the reason for JuliaLang/LinearAlgebra.jl#506 because the "R" returned by SPQR is only [R00 R01]. We therefore have the choice of using either size(A,2) or rank(A) as the number of columns in Q0 and adjust R accordingly. (I went with the latter option in #26391.) In any case, neither size(A,2) nor rank(A) is readly available in the current QRSparseQ which is the reason why I had to adjust the type. Hope this explanation clarified the matter a bit.

[KlausC]

I just added some analysis to my post in discourse: about QR-decomposition

[KlausC]

You can argue if that is a good idea but if one of the interpretations should go, it would definately be the Q0 interpretation.

I don't want to argue like that - when left multiplying with F.Q*B several row counts of B are acceptable. At the end, the rows of the rhs will have to be padded with zeroes, before applying the Householder factors.

Unfortunately I found another issue:
It is not possible to evaluate R' * Q', which is the transpose of Q*R. I think right multiplication with Q' should also be defined for AbstractQ types.

The confusing part is that size(Q, 2) != size(Matrix(Q), 2)

All I want here is clarity: that should be documented close to qrfact.

We therefore have the choice of using either size(A,2) or rank(A) as the number of columns in Q0 and adjust R accordingly.

I agree completely. In my post I suggest to make a decision about those case and also to store the row count of R in Q.

In the dense case that is fairly simple because it is the same as the number of elementary reflectors, i.e. size(F.factors, 2) but in the sparse case that need not be the case since ...

I am not sure, if it is wise to rely on an implementation detail of the upstream library here. I would prefer to store the row count of R in Q also in the LAPACK and generic cases.

It estimates rank so when the trailing matrix is considered zero, no more reflectors are produced.

It would be easy, to add the rank-estimation also in the LAPACK and generic cases with pivoting: Just remove the last rows of R, if they are small. The corresponding Householder-factors would do no harm, because they operate on the complement of image(A).

[andreasnoack]

I am not sure, if it is wise to rely on an implementation detail of the upstream library here...
...It would be easy, to add the rank-estimation also in the LAPACK and generic cases with pivoting

In the dense case, the rank estimation is currently delayed until e.g. ldiv! which takes tolerance. This is how LAPACK is structured. I like the idea of keeping as much information as long as possible but it makes the sparse and dense cases behave slightly differently. The alternative is to move the rank estimation into the factorization in the dense case. This actually would be in line with pivoted Cholesky so actually LAPACK is a bit inconsistent here. I'll open separate issues for the insufficient documention that you pointed out and an issue for when the rank should be determined and then I'll merge this since it fixes a a current bug.

[andreasnoack]

Weird. I cannot reproduce the test errors locally.

[JeffBezanson]

Is this really completely non-breaking?

[andreasnoack]

Yes. It shouldn't change any correct behavior.

[StefanKarpinski]

The policy for backports to long-term support branches is to be extremely conservative. Can you elaborate on the details of what that means?

[andreasnoack]

It's a bugfix. Before this PR, we wrongly assumed that the number of reflectors was the same as the number of columns in a matrix. That isn't the case for sparse matrices so a field is added to the struct to indicate number of columns.