Julia: Make size(F.Q) == size(Matrix(F.Q))?

Created on 23 Mar 2018  Â·  6Comments  Â·  Source: JuliaLang/julia

Following up on the discussion in #26392. E.g.

julia> using LinearAlgebra

julia> A = randn(4,2);

julia> F = qrfact(A);

julia> size(F.Q), size(Matrix(F.Q))
((4, 4), (4, 2))
doc linear algebra

Most helpful comment

What about exposing the truncated Matrix(F.Q) as a different accessor instead? Maybe as F.Q0?

All 6 comments

What about exposing the truncated Matrix(F.Q) as a different accessor instead? Maybe as F.Q0?

It has been bugging me for a while that F.Q behaves like a Matrix in a is not really Q way.

using LinearAlgebra, Random
srand(0)
X = rand(100, 3)
B = rand(100, 4)
F = qrfact(X)
F.Q * B
Matrix(F.Q) * B

I prefer exposing Q more often than the compact Householder elementary reflectors.

I'm a bit confused what the use of the square Q is. For example, I was surprised this fails:

julia> Q,R = qr(randn(5,3));

julia> b = randn(5);

julia> R \ (Q'b)
ERROR: DimensionMismatch("B has first dimension 5 but needs 3")
Stacktrace:
 [1] trtrs!(::Char, ::Char, ::Char, ::Array{Float64,2}, ::Array{Float64,1}) at /Users/solver/Projects/julia-1.1/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/lapack.jl:3340
 [2] ldiv! at /Users/solver/Projects/julia-1.1/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:583 [inlined]
 [3] \ at /Users/solver/Projects/julia-1.1/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:1854 [inlined]
 [4] \(::Array{Float64,2}, ::Array{Float64,1}) at /Users/solver/Projects/julia-1.1/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/generic.jl:903
 [5] top-level scope at none:0

I'm a bit confused what the use of the square Q is

The way I see it is that a compact representation of Q naturally is square. That is also how LAPACK applies a compact QR. An application where the square behavior can be useful is the statistical application of least squares where you'd also like sum of the squared residuals which can be computed without computing the residual since

norm(A*x - b)^2 = norm(Q*[R; 0]*x - b)^2 = norm([R*x; 0] - Q'*b)^2 =
  norm(R*x - Q0'*b)^2 + norm(Q1'*b)^2 = norm(Q1'*b)^2

when x is the least squares solution R\(Q0'*b).

The issue in the OP seems to be already addressed by https://github.com/JuliaLang/julia/pull/28446. But is it still an option to add a different accessor as suggested by @StefanKarpinski https://github.com/JuliaLang/julia/issues/26591#issuecomment-375677651 and perhaps "fix" Matrix(Q) s.t. size(F.Q) == size(Matrix(F.Q))?

+1 for enabling ‘R \ (Q'b)’ or something similar.

Was this page helpful?
0 / 5 - 0 ratings

Related issues

arshpreetsingh picture arshpreetsingh  Â·  3Comments

m-j-w picture m-j-w  Â·  3Comments

StefanKarpinski picture StefanKarpinski  Â·  3Comments

TotalVerb picture TotalVerb  Â·  3Comments

dpsanders picture dpsanders  Â·  3Comments