LinearAlgebra.jl icon indicating copy to clipboard operation
LinearAlgebra.jl copied to clipboard

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

Open andreasnoack opened this issue 7 years ago • 6 comments

Following up on the discussion in JuliaLang/julia#26392. E.g.

julia> using LinearAlgebra

julia> A = randn(4,2);

julia> F = qr(A);

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

andreasnoack avatar Mar 23 '18 11:03 andreasnoack

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

StefanKarpinski avatar Mar 23 '18 14:03 StefanKarpinski

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.

Nosferican avatar Mar 27 '18 06:03 Nosferican

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

dlfivefifty avatar Mar 11 '19 13:03 dlfivefifty

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).

andreasnoack avatar Mar 11 '19 15:03 andreasnoack

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/LinearAlgebra.jl/issues/513 and perhaps "fix" Matrix(Q) s.t. size(F.Q) == size(Matrix(F.Q))?

tkf avatar Aug 20 '19 14:08 tkf

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

cortner avatar Aug 23 '19 20:08 cortner