Converting QRCompactWY.Q to Array changes its size

General Usage
1

I’m trying to convert the Q matrix from a QR transform into an Array. But it changes size. Is this a bug or am I doing it wrong?

julia> size(qr(ones(9,5)).Q)
(9, 9)

julia> size(Array(qr(ones(9,5)).Q))
(9, 5)

Motivation: I want an orthonormal basis for the support of a matrix, and a basis for the perpendicular space. The surrounding code looks like this:

q = qr(mat, Val(true))
take = [norm(x) >= tol for x in eachrow(q.R)]
qQ = Array(q.Q)
resize!(take, size(qQ)[2])
support = qQ[:,take]
perp = qQ[:,.!take]

It works if I don’t cast q.Q to an Array, but the indexing operations are then very slow.

2

The Q matrix is not stored explicitly in a QR decomposition. It is in a form (a set of elementary reflectors, also called Householder transformations) in which it is easy to compute Q times an array or Q’ times an array. If you really want the whole Q to create a basis, as you say, you may find that just multiplying Q by the identity is the easiest approach.

1 Like
3

That works, but is very slow.

julia> @time size(Array(qr(randn(200,200)).Q))
  0.022359 seconds (16 allocations: 1.636 MiB)
(200, 200)

julia> @time size(qr(randn(200,200)).Q * I)
  2.480357 seconds (120.01 k allocations: 207.935 MiB, 1.63% gc time)
(200, 200)
4

Oh, but it’s fast if I convert I to an Array:

julia> @time size(qr(randn(200,200)).Q * Matrix(1.0*I, (200, 200)))
  0.004260 seconds (16 allocations: 1.636 MiB)
(200, 200)

It seems to me that julia could do better here by improving support for Array(qr.Q) and qr.Q * I. I think I’ll file a bug report. Thanks for the hints.

1 Like
5

I think there’s currently an open PR for this.

6

So there is. And posted just a few hours before my question! Thanks.

https://github.com/JuliaLang/julia/issues/38972

7

You may find it is better to use lmul!(Q, Matrix(1.0*I, (200, 200))) which performs an in-place multiplication.

1 Like