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))
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.
Most helpful comment
What about exposing the truncated
Matrix(F.Q)
as a different accessor instead? Maybe asF.Q0
?