Multiplying `AbstractQ` with structured matrices can cause large allocations

[mikmoore]

Came up on https://discourse.julialang.org/t/best-way-for-linear-regression-problem-on-product-features/127786.

julia> using LinearAlgebra

julia> N = 1_000_000; X = ones(N, 1); Y = Diagonal(ones(N)); Xqr = qr(X);

julia> Xqr \ Y
ERROR: OutOfMemoryError()

julia> Xqr.Q' * Y # tries to densely materialize Y # EDIT: No, see below where we discuss that this operation fails because the result is NxN.
ERROR: OutOfMemoryError()

julia> Matrix(Xqr.Q)' * Y # the materialized rectangular Q works fine
1×1000000 Matrix{Float64}:
 -0.001  -0.001  -0.001  -0.001  -0.001  …  -0.001  -0.001  -0.001  -0.001

julia> Xqr.R \ (Matrix(Xqr.Q)' * Y) # can solve the system by manual application
1×1000000 Matrix{Float64}:
 1.0e-6  1.0e-6  1.0e-6  1.0e-6  1.0e-6  …  1.0e-6  1.0e-6  1.0e-6  1.0e-6

This happens with both QRCompactWYQ (this case) and QRPackedQ (pivoted case), at least.

I haven't wrapped my head fully around what the full solution should be (AbstractQ representations make my head spin). What's efficient may vary with context, so it's possible there are some real design decisions involved around this issue. Broadly speaking, it seems that applying a "dimension-reducing" AbstractQ to certain structured matrices (Diagonal and AbstractSparseMatrix, for example) might benefit from materializing the AbstractQ rather than the structured matrix.

[mikmoore]

@dkarrasch in case you have thoughts on this.

[araujoms]

Xqr.Q' * Y does not densely materialize Y, it keeps it as Diagonal. Instead, the problem is that Xqr.Q' * Y is an N x N matrix, and Julia becomes out of memory when you try to allocate that. On the other hand, Matrix(Xqr.Q)' * Y is an 1 x N matrix, which doesn't pose any difficulty.

The underlying issue is that the handling of non-square AbstractQ is a mess. Sometimes it is interpreted as a square matrix, sometimes as a non-square matrix. I had problems with it before, see #1172.

[mikmoore]

Indeed, you're right on that point. Xqr.Q' * Y is a square matrix and this is why that particular operation fails.

But that doesn't excuse Xqr \ Y from working. Xqr.R is 1x1 in this case, so it only needs the 1xN-shaped result of Xqr.Q' * Y.