diff --git a/Project.toml b/Project.toml index 3f106e9..b065c6d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "UpdatableQRFactorizations" uuid = "8b110c62-ada2-4dd7-b53a-29f29fe8f7f4" authors = ["Sebastian Ament "] -version = "1.0.0" +version = "1.0.1" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/UpdatableQRFactorizations.jl b/src/UpdatableQRFactorizations.jl index 870d802..ea9c65a 100644 --- a/src/UpdatableQRFactorizations.jl +++ b/src/UpdatableQRFactorizations.jl @@ -1,12 +1,13 @@ module UpdatableQRFactorizations using LinearAlgebra -using LinearAlgebra: Givens +using LinearAlgebra: Givens, AbstractQ export AbstractQR, GivensQ, GivensQR, UpdatableGivensQR, UpdatableQR, UQR, add_column!, remove_column! const AbstractMatOrFac{T} = Union{AbstractMatrix{T}, Factorization{T}} +const AdjointQ = isdefined(LinearAlgebra, :AdjointQ) ? LinearAlgebra.AdjointQ : Adjoint # IDEA: could have efficient version for sparse matrices include("abstract.jl") diff --git a/src/abstract.jl b/src/abstract.jl index 4a35d2b..6a68665 100644 --- a/src/abstract.jl +++ b/src/abstract.jl @@ -3,15 +3,18 @@ abstract type AbstractQR{T} <: Factorization{T} end # AbstractQR assumes existence of n, m, Q, R fields Base.size(F::AbstractQR) = (F.n, F.m) Base.size(F::AbstractQR, i::Int) = i > 2 ? 1 : size(F)[i] -Base.eltype(F::AbstractQR{T}) where T = T +Base.eltype(::AbstractQR{T}) where T = T Base.:\(F::AbstractQR, x::AbstractVecOrMat) = ldiv!(F, copy(x)) +# disambiguation +Base.:\(F::AbstractQR{T}, B::VecOrMat{Complex{T}}) where {T<:LinearAlgebra.BlasReal} = + complex.(F \ real(B), F \ imag(B)) # some helpers """ -``` + number_of_rotations(n::Int, m::Int) -``` + Computes the number of Givens rotations that are necessary to compute the QR factorization of a general matrix of size n by m. """ @@ -25,9 +28,9 @@ number_of_rotations_to_append_column(n::Int, m::Int, k::Int = 1) = k*(n-m) - (k* number_of_rotations_to_remove_column(m::Int, k::Int) = m-k """ -``` + allocate_rotations(T::DataType, n::Int, m::Int) -``` + Allocates a vector of Givens rotations of a length that is necessary to compute the QR factorization of a general matrix of size n by m. """ diff --git a/src/givensQ.jl b/src/givensQ.jl index 134b7b9..145c01e 100644 --- a/src/givensQ.jl +++ b/src/givensQ.jl @@ -1,5 +1,5 @@ # Givens representation of Q matrix, memory efficient compared to dense representation -struct GivensQ{T, QT<:AbstractVector{<:Givens{T}}} <: Factorization{T} +struct GivensQ{T, QT<:AbstractVector{<:Givens{T}}} <: AbstractQ{T} rotations::QT n::Int m::Int @@ -15,67 +15,104 @@ end Base.size(F::GivensQ) = (F.n, F.n) # could be square but also (n, m)? Base.size(F::GivensQ, i::Int) = i > 2 ? 1 : size(F)[i] Base.copy(F::GivensQ) = GivensQ(copy(F.rotations), F.n, F.m) -Base.adjoint(F::GivensQ) = Adjoint(F) +if !isdefined(LinearAlgebra, :AdjointQ) # VERSION < v"1.10-" + Base.adjoint(F::GivensQ) = Adjoint(F) +end # NOTE: this also works if rotations is empty (is identity operator) function LinearAlgebra.lmul!(F::GivensQ, X::AbstractVecOrMat) - m = size(X, 1) for G in Iterators.reverse(F.rotations) # the adjoint of a GivensQ reverses the order of the rotations lmul!(G', X) # ... and takes their adjoint end return X end -function LinearAlgebra.lmul!(A::Adjoint{<:Any, <:GivensQ}, X::AbstractVecOrMat) - F = A.parent +function LinearAlgebra.lmul!(A::AdjointQ{<:Any, <:GivensQ}, X::AbstractVecOrMat) + F = parent(A) for G in F.rotations # the adjoint of a GivensQ reverses the order of the rotations lmul!(G, X) # ... and takes their adjoint end return X end -LinearAlgebra.rmul!(X::AbstractVecOrMat, F::GivensQ) = lmul!(F', X')' -LinearAlgebra.rmul!(X::AbstractVecOrMat, F::Adjoint{<:Any, <:GivensQ}) = lmul!(F', X')' - -lmul(F, X) = lmul!(F, copy(X)) -rmul(X, F) = rmul!(copy(X), F) +# disambiguation +LinearAlgebra.lmul!(F::GivensQ, X::LinearAlgebra.AbstractTriangular) = lmul!(F, full!(X)) +LinearAlgebra.lmul!(A::AdjointQ{<:Any, <:GivensQ}, X::LinearAlgebra.AbstractTriangular) = + lmul!(F, full!(X)) -Base.:*(F::GivensQ, X::AbstractVector) = lmul(F, X) -Base.:*(F::GivensQ, X::AbstractMatrix) = lmul(F, X) -Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = lmul(F, X) -Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul(F, X) - -Base.:*(X::AbstractMatrix, F::GivensQ) = rmul(X, F) -Base.:*(X::AbstractMatrix, F::Adjoint{<:Any, <:GivensQ}) = rmul(X, F) +LinearAlgebra.rmul!(X::AbstractVecOrMat, F::GivensQ) = lmul!(F', X')' +LinearAlgebra.rmul!(X::AbstractVecOrMat, F::AdjointQ{<:Any, <:GivensQ}) = lmul!(F', X')' +# disambiguation +LinearAlgebra.rmul!(X::LinearAlgebra.AbstractTriangular, F::GivensQ) = lmul!(F', X')' +LinearAlgebra.rmul!(X::LinearAlgebra.AbstractTriangular, F::AdjointQ{<:Any, <:GivensQ}) = + lmul!(F', X')' -function Base.:*(F::GivensQ, G::Adjoint{<:Any, <:GivensQ}) +function Base.:*(F::GivensQ, G::AdjointQ{<:Any, <:GivensQ}) T = promote_type(eltype(F), eltype(G)) F === G' ? (one(T) * I)(F.n) : F * Matrix(G) end -function Base.:*(F::Adjoint{<:Any, <:GivensQ}, G::GivensQ) +function Base.:*(F::AdjointQ{<:Any, <:GivensQ}, G::GivensQ) T = promote_type(eltype(F), eltype(G)) F' === G ? (one(T) * I)(G.n) : F * Matrix(G) end + function Base.:(==)(F::GivensQ, G::GivensQ) F.rotations == G.rotations && F.n == G.n && F.m == G.m end -LinearAlgebra.ldiv!(F::GivensQ, x::AbstractVector) = lmul!(F', x) -LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector) = lmul!(F', x) -LinearAlgebra.ldiv!(F::GivensQ, X::AbstractMatrix) = lmul!(F', X) -LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul!(F', X) +if VERSION < v"1.10-" + lmul(F, X) = lmul!(F, copy(X)) + rmul(X, F) = rmul!(copy(X), F) + Base.:*(F::GivensQ, X::AbstractVector) = lmul(F, X) + Base.:*(F::GivensQ, X::AbstractMatrix) = lmul(F, X) + Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = lmul(F, X) + Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul(F, X) -function LinearAlgebra.mul!(y::AbstractVector, F::GivensQ, x::AbstractVector) - @. y = x - lmul!(F, y) -end -function LinearAlgebra.mul!(Y::AbstractMatrix, F::GivensQ, X::AbstractMatrix) - @. Y = X - lmul!(F, Y) -end -function LinearAlgebra.mul!(y::AbstractVector, F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector) - @. y = x - lmul!(F, y) -end -function LinearAlgebra.mul!(Y::AbstractMatrix, F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) - @. Y = X - lmul!(F, Y) + Base.:*(X::AbstractMatrix, F::GivensQ) = rmul(X, F) + Base.:*(X::AbstractMatrix, F::AdjointQ{<:Any, <:GivensQ}) = rmul(X, F) + + Base.:*(F::GivensQ, X::StridedVector) = lmul(F, X) + Base.:*(F::GivensQ, X::StridedMatrix) = lmul(F, X) + Base.:*(F::GivensQ, X::Adjoint{<:Any, <:StridedVecOrMat}) = lmul(F, X) + Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::StridedVector) = lmul(F, X) + Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::StridedMatrix) = lmul(F, X) + Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::Adjoint{<:Any, <:StridedVecOrMat}) = lmul(F, X) + + Base.:*(X::StridedMatrix, F::GivensQ) = rmul(X, F) + Base.:*(X::StridedMatrix, F::AdjointQ{<:Any, <:GivensQ}) = rmul(X, F) + + LinearAlgebra.mul!(Y::StridedVector, F::GivensQ, X::StridedVector) = + invoke(LinearAlgebra.mul!, Tuple{AbstractVector, GivensQ, AbstractVector}, Y, F, X) + LinearAlgebra.mul!(Y::StridedMatrix, F::GivensQ, X::StridedMatrix) = + invoke(LinearAlgebra.mul!, Tuple{AbstractMatrix, GivensQ, AbstractMatrix}, Y, F, X) + LinearAlgebra.mul!(Y::StridedVector, F::AdjointQ{<:Any, <:GivensQ}, X::StridedVector) = + invoke(LinearAlgebra.mul!, Tuple{AbstractVector, AdjointQ{<:Any, <:GivensQ}, AbstractVector}, Y, F, X) + LinearAlgebra.mul!(Y::StridedMatrix, F::AdjointQ{<:Any, <:GivensQ}, X::StridedMatrix) = + invoke(LinearAlgebra.mul!, Tuple{AbstractMatrix, AdjointQ{<:Any, <:GivensQ}, AbstractMatrix}, Y, F, X) + + if VERSION < v"1.8-" + Base.:\(F::GivensQ, X::AbstractVector) = F'X + Base.:\(F::GivensQ, X::AbstractMatrix) = F'X + Base.:\(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = F'X + Base.:\(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = F'X + end + LinearAlgebra.ldiv!(F::GivensQ, x::AbstractVector) = lmul!(F', x) + LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector) = lmul!(F', x) + LinearAlgebra.ldiv!(F::GivensQ, X::AbstractMatrix) = lmul!(F', X) + LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul!(F', X) + + function LinearAlgebra.mul!(y::AbstractVector, F::GivensQ, x::AbstractVector) + @. y = x + lmul!(F, y) + end + function LinearAlgebra.mul!(Y::AbstractMatrix, F::GivensQ, X::AbstractMatrix) + @. Y = X + lmul!(F, Y) + end + function LinearAlgebra.mul!(y::AbstractVector, F::AdjointQ{<:Any, <:GivensQ}, x::AbstractVector) + @. y = x + lmul!(F, y) + end + function LinearAlgebra.mul!(Y::AbstractMatrix, F::AdjointQ{<:Any, <:GivensQ}, X::AbstractMatrix) + @. Y = X + lmul!(F, Y) + end end diff --git a/src/updatableQR.jl b/src/updatableQR.jl index 68861dc..38d8576 100644 --- a/src/updatableQR.jl +++ b/src/updatableQR.jl @@ -1,7 +1,7 @@ """ -``` + UpdatableGivensQR <: AbstractQR <: Factorization -``` + Holds storage for Givens rotations and the R matrix for a QR factorization. Allows for efficient addition and deletion of columns. See `add_column!` and `remove_column!`. @@ -102,9 +102,9 @@ function LinearAlgebra.ldiv!(F::UQR, x::AbstractVecOrMat) end """ -``` + add_column!(F::UpdatableGivensQR, x::AbstractVector, k::Int = size(F, 2) + 1) -``` + Given an existing QR factorization `F` of a matrix, computes the factorization of the same matrix after a new column `x` has been added as the `k`ᵗʰ column. Computational complexity: O(nm) where size(F) = (n, m). @@ -195,9 +195,9 @@ function ensure_space_to_append_column!(rotations::AbstractVector{<:Givens}, end """ -``` + remove_column!(F::UpdatableGivensQR, k::Int = size(F, 2)) -``` + Updates the existing QR factorization `F` of a matrix to the factorization of the same matrix after its kᵗʰ column has been deleted. Computational complexity: O(m²) where size(F) = (n, m). diff --git a/test/UpdatableQRFactorizations.jl b/test/UpdatableQRFactorizations.jl index 4d443ac..ca095a5 100644 --- a/test/UpdatableQRFactorizations.jl +++ b/test/UpdatableQRFactorizations.jl @@ -74,9 +74,11 @@ using UpdatableQRFactorizations x = randn(m) Ax = A*x @test F \ (A*x) ≈ x + @test F \ complex(A*x) ≈ x r = 4 X = randn(m, r) @test F \ (A*X) ≈ X + @test F \ complex(A*X) ≈ X end @testset "UpdatableQR" begin