Wrong mismatch test on LinearAlgebra.lmul!

[lrsantos11]

https://github.com/JuliaLang/julia/blob/b56a9f07948255dfbe804eef25bdbada06ec2a57/stdlib/SuiteSparse/src/spqr.jl#L207 The code should be

     size(A, 1) !=size(Q.factors,2)

as in https://github.com/JuliaLang/julia/blob/80516ca20297a67b996caa08c38786332379b6a5/stdlib/LinearAlgebra/src/qr.jl#L533

[andreasnoack]

Why do think that? Both the current versions test for mismatch of the first dimension of the arrays.

[lrsantos11]

Because I could not make the simple F.Q*F.R work, where F = qr(A) and A is a $ 27\times 2$ sparse matrix due to the fact that F.Q is $27\times 27$ and F.R is $2\times 2$. Also, why the other SPQR routines to do multiplication or solving the linear system are not used?

[andreasnoack]

I see. In lmul!, Q is square and that is the same for dense and sparse.

julia> A = randn(4,2);

julia> Asp = sparse(A);

julia> F = qr(A);

julia> Fsp = qr(Asp);

julia> lmul!(F.Q, copy(F.R))
ERROR: DimensionMismatch("first dimensions of C, 2, and V, 4 must match")
Stacktrace:
 [1] gemqrt!(::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,2}, ::Array{Float64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/lapack.jl:2842
 [2] lmul!(::LinearAlgebra.QRCompactWYQ{Float64,Array{Float64,2}}, ::Array{Float64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/qr.jl:525
 [3] top-level scope at none:0

julia> lmul!(Fsp.Q, Array(Fsp.R))
ERROR: DimensionMismatch("size(Q) = (4, 4) but size(A) = (2, 2)")
Stacktrace:
 [1] lmul!(::SuiteSparse.SPQR.QRSparseQ{Float64,Int64}, ::Array{Float64,2}) at /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/SuiteSparse/src/spqr.jl:208
 [2] top-level scope at none:0

However, in the dense case * allows both a thin and wide interpretation of Q. That is why F.Q*F.R works. What happens is that R is appended with zeros before the call to lmul!. We could make something similar work for sparse matrices.

[andreasnoack]

I should mention that some might not like the current behavior where Q can act as both square and tall so in a future version of LinearAlgebra, we might want to remove the version where Q acts as tall. Being able to multiply the thin Q with R is rarely a relevant computation in real problems.

Also, why the other SPQR routines to do multiplication or solving the linear system are not used?

Please elaborate. I'm not sure what you are asking here.

[lrsantos11]

Yeah. That would be easy to implement. However, even in this case, F.Q*F.R could be different from A, because, by default, SPQR performs pivoting both on rows and columns. I discovered that because I'm trying to solve the minimum 2-norm solution for a fat matrix A, to the problem Ax=b. If you want to use QR factorization, you should compute F = qr(transpose(A)) and the solution is x = F.Q*(F.R'\b). However, because of the permutations and because of the problem with the multiplication, I could not.

BTW, F = qr(transpose(A)) also won't work because qr is not extended to adjoints. I extended in my code using

qr(M::Adjoint{Float64,SparseMatrixCSC{Float64,Int64}}) = qr(copy(M))

[lrsantos11]

Also, why the other SPQR routines to do multiplication or solving the linear system are not used?

Please elaborate. I'm not sure what you are asking here.

SPQR has all the routines implemented in C++ to perform the multiplications with Q and R and the respective adjoints (they are avaliable on the Julia 0.6.4) as well as the routines to solve the minimum norm solutions, whether the matrix is tall (problem min ||Ax - b||^2) or being fat (problem min ||x|| s.t. Ax=b). Why not use them?

For instance, call them when A\b is used.

[chriscoey]

SPQR has all the routines implemented in C++ to perform the multiplications with Q and R and the respective adjoints (they are avaliable on the Julia 0.6.4) as well as the routines to solve the minimum norm solutions, whether the matrix is tall (problem min ||Ax - b||^2) or being fat (problem min ||x|| s.t. Ax=b). Why not use them?

Yes - why were these auxiliary SPQR methods removed? I am multiplying a sparse matrix A by a QRSparseQ (A * F.Q) and it is extremely slow, hitting some generic fallback. Presumably it would be much faster if we could use SPQR's routines.

[andreasnoack]

The representation we use in Julia was completely changed, i.e. we use to hold just a pointer to the SPQR object but based on popular demand, I changed it to storing the actual data in Julia. The idea was that we could support the required methods in Julia code. If we are hitting slow fallbacks we should just fix that.

[chriscoey]

@andreasnoack Fix that by accessing the old SPQR methods? Or by implementing our own?

[andreasnoack]

We should just implement our own. It's no longer an option to use the SPQR methods since we store the data in a Julia struct. It's probably not hard to implement and it would be great if you could prepare a PR.