Givens rotation based compact QR for sparse matrices

Specific Domains Numerics
1

Is there any implementation/API for the QR factorization based on Givens rotations that don’t produce the Q matrix but gives an API for applying the Givens rotations to a vector possibly after inverting/transposing them? I believe the Givens rotation variant should be more sparsity and memory friendly than the Householder reflection variant.

Also it seems like qrfact doesn’t work at all on sparse matrices. Am I missing something?

julia> a = sprand(100,100,0.05);

julia> b = qrfact(a);

julia> b[:Q];
ERROR: MethodError: no method matching getindex(::Base.SparseArrays.SPQR.Factorization{Float64}, ::Symbol)

julia> b[:R];
ERROR: MethodError: no method matching getindex(::Base.SparseArrays.SPQR.Factorization{Float64}, ::Symbol)

julia> b[:Q] * rand(100);
ERROR: MethodError: no method matching getindex(::Base.SparseArrays.SPQR.Factorization{Float64}, ::Symbol)

julia> a = rand(100,100);

julia> b = qrfact(a);

julia> b[:Q];

julia> b[:R];

julia> b[:Q] * rand(100);
1 Like
2

I think this is not a matter of Givens vs Householder but Householder reflectors vs explicit dense storage of the Q matrix. In next version of Julia, you’ll be able to use qrfact on sparse matrices in the same way as for dense matrices and F[:Q] will store the Householder reflectors which are typically sparse whereas the explicit version of Q is typically dense. You can try it out by downloading the nightly binaries. Then you’ll be able to do things like

julia> A = sprandn(100,10,0.1);

julia> F = qrfact(A);

julia> F[:Q]'ones(100)
100-element Array{Float64,1}:
  0.41592945666376946
  0.3323257932564315
 -0.08952795654142387
 -0.46495867887952524
  1.9466537337738807
 -0.5502857954699574
...
2 Likes
3

Householder would also work but then most references stress that the Givens QR is better for sparse matrices (less intermediate fill-in) and parallel implementations, while the Householder is better for dense matrices (fewer flops). So unless one of the fancy variants of Householder QR like the one described in this paper is used, it might be a better idea to go with Givens, no? Of course the main issue here is to not generate the Q matrix fully, but why not use the fastest algorithm among the ones that don’t do that? But then I suppose it depends on what’s out there to be wrapped, and some historical benchmarking-based de-facto standards which I am not entirely aware of. Any comments?

4

If your matrices are banded you can use BandedMatrices.jl

My impression is that there is confusion re- givens vs householder: a sparse givens and a sparse householder have the exact same fill-in. The only case where Givens wins that I can think of is 2x2, where there are less operations by avoiding the sign check.

At won’t point I changed BandedMatrices QR from Givens to Householder with a significant speed up.

5

Interesting, I suppose whatever is de-facto standard is probably better than other alternatives in literature. I will need to dig in deeper to make any stronger claims, which I may leave for later if necessary.