Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Don't subtype AbstractQ <: AbstractMatrix #46196

Merged
merged 27 commits into from
Feb 16, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
27 commits
Select commit Hold shift + click to select a range
f93e32c
Don't subtype `AbstractQ <: AbstractMatrix`
dkarrasch Jul 27, 2022
81acd44
tiny fixes
dkarrasch Nov 30, 2022
a60cf20
Merge branch 'master' into dk/q
dkarrasch Nov 30, 2022
53f5c4f
Update special.jl
dkarrasch Nov 30, 2022
8701761
Merge branch 'master' into dk/q
dkarrasch Dec 3, 2022
baf6035
get multiple columns, add tests
dkarrasch Dec 3, 2022
0a9628c
fix merge issue
dkarrasch Dec 3, 2022
9b3a3f4
fix getindex
dkarrasch Dec 4, 2022
416d3e3
Merge branch 'master' into dk/q
dkarrasch Dec 5, 2022
13f0b13
Merge branch 'master' into dk/q
dkarrasch Dec 7, 2022
90811fc
minor additions
dkarrasch Dec 7, 2022
d4d7b8f
Merge branch 'master' into dk/q
dkarrasch Dec 7, 2022
efa47de
redirect view to getindex
dkarrasch Dec 14, 2022
183db1b
Merge branch 'master' into dk/q
dkarrasch Dec 14, 2022
58beb65
small additions/fixes
dkarrasch Dec 21, 2022
3218bcd
code layout
dkarrasch Dec 30, 2022
c6d37c0
Merge branch 'master' into dk/q
dkarrasch Jan 3, 2023
6021102
Merge branch 'master' into dk/q
dkarrasch Jan 8, 2023
dbedfea
add NEWS entry
dkarrasch Jan 14, 2023
621394b
add link to docs
dkarrasch Jan 17, 2023
c16d1c1
Merge branch 'master' into dk/q
dkarrasch Feb 7, 2023
708f459
add `logabsdet` and `det`
dkarrasch Feb 10, 2023
18ea0a8
Merge branch 'master' into dk/q
dkarrasch Feb 12, 2023
2419d63
Merge branch 'master' into dk/q
dkarrasch Feb 14, 2023
214dd04
Update hessenberg.jl
dkarrasch Feb 15, 2023
6310777
Merge branch 'master' into dk/q
dkarrasch Feb 15, 2023
78858d3
fix test
dkarrasch Feb 16, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
101 changes: 99 additions & 2 deletions stdlib/LinearAlgebra/docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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.
dkarrasch marked this conversation as resolved.
Show resolved Hide resolved
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/).
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -611,6 +707,7 @@ LinearAlgebra.BLAS.trsv
```

**return a matrix**

```@docs
LinearAlgebra.BLAS.ger!
# xGERU
Expand Down
7 changes: 4 additions & 3 deletions stdlib/LinearAlgebra/src/LinearAlgebra.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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")
Expand All @@ -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")
Expand Down
Loading