Don't subtype AbstractQ <: AbstractMatrix (#46196)

Author: dkarrasch

NEWS.md

@@ -51,6 +51,15 @@ Standard library changes
 
 #### LinearAlgebra
 
+* `AbstractQ` no longer subtypes to `AbstractMatrix`. Moreover, `adjoint(Q::AbstractQ)`
+  no longer wraps `Q` in an `Adjoint` type, but instead in an `AdjointQ`, that itself
+  subtypes `AbstractQ`. This change accounts for the fact that typically `AbstractQ`
+  instances behave like function-based, matrix-backed linear operators, and hence don't
+  allow for efficient indexing. Also, many `AbstractQ` types can act on vectors/matrices
+  of different size, acting like a matrix with context-dependent size. With this change,
+  `AbstractQ` has a well-defined API that is described in detail in the
+  [Julia documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#man-linalg-abstractq)
+  ([#46196]).
 * Adjoints and transposes of `Factorization` objects are no longer wrapped in `Adjoint`
   and `Transpose` wrappers, respectively. Instead, they are wrapped in
   `AdjointFactorization` and `TranposeFactorization` types, which themselves subtype

stdlib/LinearAlgebra/docs/src/index.md

@@ -150,7 +150,10 @@ julia> sB\x
  -1.1086956521739126
  -1.4565217391304346
 ```
-The `\` operation here performs the linear solution. The left-division operator is pretty powerful and it's easy to write compact, readable code that is flexible enough to solve all sorts of systems of linear equations.
+
+The `\` operation here performs the linear solution. The left-division operator is pretty
+powerful and it's easy to write compact, readable code that is flexible enough to solve all
+sorts of systems of linear equations.
 
 ## Special matrices
 
@@ -309,6 +312,94 @@ Adjoints and transposes of [`Factorization`](@ref) objects are lazily wrapped in
 `AdjointFactorization` and `TransposeFactorization` objects, respectively. Generically,
 transpose of real `Factorization`s are wrapped as `AdjointFactorization`.
 
+## [Orthogonal matrices (`AbstractQ`)](@id man-linalg-abstractq)
+
+Some matrix factorizations generate orthogonal/unitary "matrix" factors. These
+factorizations include QR-related factorizations obtained from calls to [`qr`](@ref), i.e.,
+`QR`, `QRCompactWY` and `QRPivoted`, the Hessenberg factorization obtained from calls to
+[`hessenberg`](@ref), and the LQ factorization obtained from [`lq`](@ref). While these
+orthogonal/unitary factors admit a matrix representation, their internal representation
+is, for performance and memory reasons, different. Hence, they should be rather viewed as
+matrix-backed, function-based linear operators. In particular, reading, for instance, a
+column of its matrix representation requires running "matrix"-vector multiplication code,
+rather than simply reading out data from memory (possibly filling parts of the vector with
+structural zeros). Another clear distinction from other, non-triangular matrix types is
+that the underlying multiplication code allows for in-place modification during multiplication.
+Furthermore, objects of specific `AbstractQ` subtypes as those created via [`qr`](@ref),
+[`hessenberg`](@ref) and [`lq`](@ref) can behave like a square or a rectangular matrix
+depending on context:
+
+```julia
+julia> using LinearAlgebra
+
+julia> Q = qr(rand(3,2)).Q
+3×3 LinearAlgebra.QRCompactWYQ{Float64, Matrix{Float64}, Matrix{Float64}}
+
+julia> Matrix(Q)
+3×2 Matrix{Float64}:
+ -0.320597   0.865734
+ -0.765834  -0.475694
+ -0.557419   0.155628
+
+julia> Q*I
+3×3 Matrix{Float64}:
+ -0.320597   0.865734  -0.384346
+ -0.765834  -0.475694  -0.432683
+ -0.557419   0.155628   0.815514
+
+julia> Q*ones(2)
+3-element Vector{Float64}:
+  0.5451367118802273
+ -1.241527373086654
+ -0.40179067589600226
+
+julia> Q*ones(3)
+3-element Vector{Float64}:
+  0.16079054743832022
+ -1.674209978965636
+  0.41372375588835797
+
+julia> ones(1,2) * Q'
+1×3 Matrix{Float64}:
+ 0.545137  -1.24153  -0.401791
+
+julia> ones(1,3) * Q'
+1×3 Matrix{Float64}:
+ 0.160791  -1.67421  0.413724
+```
+
+Due to this distinction from dense or structured matrices, the abstract `AbstractQ` type
+does not subtype `AbstractMatrix`, but instead has its own type hierarchy. Custom types
+that subtype `AbstractQ` can rely on generic fallbacks if the following interface is satisfied.
+For example, for
+
+```julia
+struct MyQ{T} <: LinearAlgebra.AbstractQ{T}
+    # required fields
+end
+```
+
+provide overloads for
+
+```julia
+Base.size(Q::MyQ) # size of corresponding square matrix representation
+Base.convert(::Type{AbstractQ{T}}, Q::MyQ) # eltype promotion [optional]
+LinearAlgebra.lmul!(Q::MyQ, x::AbstractVecOrMat) # left-multiplication
+LinearAlgebra.rmul!(A::AbstractMatrix, Q::MyQ) # right-multiplication
+```
+
+If `eltype` promotion is not of interest, the `convert` method is unnecessary, since by
+default `convert(::Type{AbstractQ{T}}, Q::AbstractQ{T})` returns `Q` itself.
+Adjoints of `AbstractQ`-typed objects are lazily wrapped in an `AdjointQ` wrapper type,
+which requires its own `LinearAlgebra.lmul!` and `LinearAlgebra.rmul!` methods. Given this
+set of methods, any `Q::MyQ` can be used like a matrix, preferably in a multiplicative
+context: multiplication via `*` with scalars, vectors and matrices from left and right,
+obtaining a matrix representation of `Q` via `Matrix(Q)` (or `Q*I`) and indexing into the
+matrix representation all work. In contrast, addition and subtraction as well as more
+generally broadcasting over elements in the matrix representation fail because that would
+be highly inefficient. For such use cases, consider computing the matrix representation
+up front and cache it for future reuse.
+
 ## Standard functions
 
 Linear algebra functions in Julia are largely implemented by calling functions from [LAPACK](http://www.netlib.org/lapack/).
@@ -505,38 +596,42 @@ four methods defined, for [`Float32`](@ref), [`Float64`](@ref), [`ComplexF32`](@
 and [`ComplexF64`](@ref Complex) arrays.
 
 ### [BLAS character arguments](@id stdlib-blas-chars)
+
 Many BLAS functions accept arguments that determine whether to transpose an argument (`trans`),
 which triangle of a matrix to reference (`uplo` or `ul`),
 whether the diagonal of a triangular matrix can be assumed to
 be all ones (`dA`) or which side of a matrix multiplication
 the input argument belongs on (`side`). The possibilities are:
 
 #### [Multiplication order](@id stdlib-blas-side)
+
 | `side` | Meaning                                                             |
 |:-------|:--------------------------------------------------------------------|
 | `'L'`  | The argument goes on the *left* side of a matrix-matrix operation.  |
 | `'R'`  | The argument goes on the *right* side of a matrix-matrix operation. |
 
 #### [Triangle referencing](@id stdlib-blas-uplo)
+
 | `uplo`/`ul` | Meaning                                               |
 |:------------|:------------------------------------------------------|
 | `'U'`       | Only the *upper* triangle of the matrix will be used. |
 | `'L'`       | Only the *lower* triangle of the matrix will be used. |
 
 #### [Transposition operation](@id stdlib-blas-trans)
+
 | `trans`/`tX` | Meaning                                                 |
 |:-------------|:--------------------------------------------------------|
 | `'N'`        | The input matrix `X` is not transposed or conjugated.   |
 | `'T'`        | The input matrix `X` will be transposed.                |
 | `'C'`        | The input matrix `X` will be conjugated and transposed. |
 
 #### [Unit diagonal](@id stdlib-blas-diag)
+
 | `diag`/`dX` | Meaning                                                   |
 |:------------|:----------------------------------------------------------|
 | `'N'`       | The diagonal values of the matrix `X` will be read.       |
 | `'U'`       | The diagonal of the matrix `X` is assumed to be all ones. |
 
-
 ```@docs
 LinearAlgebra.BLAS
 LinearAlgebra.BLAS.set_num_threads
@@ -582,6 +677,7 @@ and define matrix-vector operations.
 [Dongarra-1988]: https://dl.acm.org/doi/10.1145/42288.42291
 
 **return a vector**
+
 ```@docs
 LinearAlgebra.BLAS.gemv!
 LinearAlgebra.BLAS.gemv(::Any, ::Any, ::Any, ::Any)
@@ -611,6 +707,7 @@ LinearAlgebra.BLAS.trsv
 ```
 
 **return a matrix**
+
 ```@docs
 LinearAlgebra.BLAS.ger!
 # xGERU

stdlib/LinearAlgebra/src/LinearAlgebra.jl

@@ -15,7 +15,7 @@ import Base: USE_BLAS64, abs, acos, acosh, acot, acoth, acsc, acsch, adjoint, as
     length, log, map, ndims, one, oneunit, parent, permutedims, power_by_squaring,
     print_matrix, promote_rule, real, round, sec, sech, setindex!, show, similar, sin,
     sincos, sinh, size, sqrt, strides, stride, tan, tanh, transpose, trunc, typed_hcat,
-    vec, zero
+    vec, view, zero
 using Base: IndexLinear, promote_eltype, promote_op, promote_typeof,
     @propagate_inbounds, reduce, typed_hvcat, typed_vcat, require_one_based_indexing,
     splat
@@ -431,8 +431,6 @@ include("tridiag.jl")
 include("triangular.jl")
 
 include("factorization.jl")
-include("qr.jl")
-include("lq.jl")
 include("eigen.jl")
 include("svd.jl")
 include("symmetric.jl")
@@ -443,7 +441,10 @@ include("diagonal.jl")
 include("symmetriceigen.jl")
 include("bidiag.jl")
 include("uniformscaling.jl")
+include("qr.jl")
+include("lq.jl")
 include("hessenberg.jl")
+include("abstractq.jl")
 include("givens.jl")
 include("special.jl")
 include("bitarray.jl")