diff --git a/Project.toml b/Project.toml index 1a0035c..ec82fd8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "ToeplitzMatrices" uuid = "c751599d-da0a-543b-9d20-d0a503d91d24" -version = "0.7.1" +version = "0.8" [deps] AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c" diff --git a/README.md b/README.md index f6c412c..8915c56 100644 --- a/README.md +++ b/README.md @@ -8,7 +8,7 @@ ToeplitzMatrices.jl Fast matrix multiplication and division for Toeplitz, Hankel and circulant matrices in Julia -# Note +## Note Multiplication of large matrices and `sqrt`, `inv`, `LinearAlgebra.eigvals`, `LinearAlgebra.ldiv!`, and `LinearAlgebra.pinv` for circulant matrices @@ -26,113 +26,137 @@ initialize the required arrays and plan the FFT only once. You can precompute the FFT factorization with `LinearAlgebra.factorize` and then use the factorization for the FFT-based computations. -# Supported matrices - -## Toeplitz +## Introduction +### Toeplitz A Toeplitz matrix has constant diagonals. It can be constructed using - ```julia Toeplitz(vc,vr) +Toeplitz{T}(vc,vr) ``` - where `vc` are the entries in the first column and `vr` are the entries in the first row, where `vc[1]` must equal `vr[1]`. For example. - ```julia Toeplitz(1:3, [1.,4.,5.]) ``` - is a sparse representation of the matrix - ```julia [ 1.0 4.0 5.0 2.0 1.0 4.0 3.0 2.0 1.0 ] ``` - -## SymmetricToeplitz - -A symmetric Toeplitz matrix is a symmetric matrix with constant diagonals. It can be constructed with - +### Special toeplitz +`SymmetricToeplitz`, `Circulant`, `UpperTriangularToeplitz` and `LowerTriangularToeplitz` only store one vector. By convention, `Circulant` stores the first column rather than the first row. They are constructed using `TYPE(v)` where `TYPE`∈{`SymmetricToeplitz`, `Circulant`, `UpperTriangularToeplitz`, `LowerTriangularToeplitz`}. + +### Hankel +A Hankel matrix has constant anti-diagonals, for example, ```julia -SymmetricToeplitz(vc) +[ 1 2 3 + 2 3 4 ] ``` +There are a few ways to construct the above `Hankel{Int}`: +- `Hankel([1,2,3,4], (2,3)) # Hankel(v, (h,w))` +- `Hankel([1,2,3,4], 2,3) # Hankel(v, h, w)` +- `Hankel([1,2], [2,3,4]) # Hankel(vc, vr)` -where `vc` are the entries of the first column. For example, +Note that the width is usually useless, since ideally, `w=length(v)-h+1`. It exists for infinite Hankel matrices. Its existence also means that `v` could be longer than necessary. `Hankel(v)`, where the size is not given, returns `Hankel(v, (l+1)÷2, (l+1)÷2)` where `l=length(v)`. -```julia -SymmetricToeplitz([1.0, 2.0, 3.0]) -``` +The `reverse` can transform between Hankel and Toeplitz. It is used to achieve fast Hankel algorithms. -is a sparse representation of the matrix +## Implemented interface -```julia -[ 1.0 2.0 3.0 - 2.0 1.0 2.0 - 3.0 2.0 1.0 ] -``` +### Operations +- ✓ implemented +- ✗ error +- _ fall back to `Matrix` -## TriangularToeplitz +||Toeplitz|Symmetric~|Circulant|UpperTriangular~|LowerTriangular~|Hankel| +|:-:|:-:|:-:|:-:|:-:|:-:|:-:| +|getindex|✓|✓|✓|✓|✓|✓| +|.vc|✓|✓|✓|✓|✓|✓| +|.vr|✓|✓|✓|✓|✓|✓| +|size|✓|✓|✓|✓|✓|✓| +|copy|✓|✓|✓|✓|✓|✓| +|similar|✓|✓|✓|✓|✓|✓| +|zero|✓|✓|✓|✓|✓|✓| +|real|✓|✓|✓|✓|✓|✓| +|imag|✓|✓|✓|✓|✓|✓| +|fill!|✓|✗|✗|✗|✗|✓| +|conj|✓|✓|✓|✓|✓|✓| +|transpose|✓|✓|✓|✓|✓|✓| +|adjoint|✓|✓|✓|✓|✓|✓| +|tril!|✓|✗|✗|✓|✓|✗| +|triu!|✓|✗|✗|✓|✓|✗| +|tril|✓|✓|✓|✓|✓|✗| +|triu|✓|✓|✓|✓|✓|✗| +|+|✓|✓|✓|✓|✓|✓| +|-|✓|✓|✓|✓|✓|✓| +|scalar
mult|✓|✓|✓|✓|✓|✓| +|==|✓|✓|✓|✓|✓|✓| +|issymmetric||||||| +|istriu||||||| +|istril||||||| +|iszero|✓|✓|✓|✓|✓|| +|isone||||||| +|copyto!|✓|✓|✓|✓|✓|✓| +|reverse|✓|✓|✓|✓|✓|✓| +|broadcast||||||| +|broadcast!||||||| -A triangular Toeplitz matrix can be constructed using +Note that scalar multiplication, `conj`, `+` and `-` could be removed once `broadcast` is implemented. -```julia -TriangularToeplitz(ve,uplo) -``` +`reverse(Hankel)` returns a `Toeplitz`, while `reverse(AbstractToeplitz)` returns a `Hankel`. -where uplo is either `:L` or `:U` and `ve` are the rows or columns, respectively. For example, +### LinearAlgebra -```julia -TriangularToeplitz([1.,2.,3.],:L) -``` +### Constructors and conversions +||Toeplitz|Symmetric~|Circulant|UpperTriangular~|LowerTriangular~|Hankel| +|:-:|:-:|:-:|:-:|:-:|:-:|:-:| +|from AbstractVector|✓|✓|✓|✓|✓|✓| +|from AbstractMatrix|✓|✓|✓|✓|✓|✓| +|from AbstractToeplitz|✓|✓|✓|✓|✓|✗| +|to supertype|✓|✓|✓|✓|✓|✓| +|to Toeplitz|-|✓|✓|✓|✓|✗| +|to another eltype|✓|✓|✓|✓|✓|✓| -is a sparse representation of the matrix +When constructing `Toeplitz` from a matrix, the first row and the first column will be considered as `vr` and `vc`. Note that `vr` and `vc` are copied in construction to avoid the cases where they share memory. If you don't want copying, construct using vectors directly. + +When constructing `SymmetricToeplitz` or `Circulant` from `AbstractMatrix`, a second argument shall specify whether the first row or the first column is used. For example, for `A = [1 2; 3 4]`, +- `SymmetricToeplitz(A,:L)` gives `[1 3; 3 1]`, while +- `SymmetricToeplitz(A,:U)` gives `[1 2; 2 1]`. +For backward compatibility and consistency with `LinearAlgebra.Symmetric`, ```julia -[ 1.0 0.0 0.0 - 2.0 1.0 0.0 - 3.0 2.0 1.0 ] +SymmetricToeplitz(A) = SymmetricToeplitz(A, :U) +Circulant(A) = Circulant(A, :L) ``` +`Hankel` constructor also accepts the second argument, `:L` denoting the first column and the last row while `:U` denoting the first row and the last column. -## Hankel - -A Hankel matrix has constant anti-diagonals. It can be constructed using - +`Symmetric`, `UpperTriangular` and `LowerTriangular` from `LinearAlgebra` are also overloaded for convenience. ```julia -Hankel(vc,vr) +Symmetric(T::Toeplitz) = SymmetricToeplitz(T) +UpperTriangular(T::Toeplitz) = UpperTriangularToeplitz(T) +LowerTriangular(T::Toeplitz) = LowerTriangularToeplitz(T) ``` -where `vc` are the entries in the first column and `vr` are the entries in the last row, where `vc[end]` must equal `vr[1]`. For example. - +### TriangularToeplitz (obsolete) +`TriangularToeplitz` is reserved for backward compatibility. ```julia -Hankel([1.,2.,3.], 3:5) +TriangularToeplitz = Union{UpperTriangularToeplitz,LowerTriangularToeplitz} ``` - -is a sparse representation of the matrix - +The old interface is implemented by ```julia -[ 1.0 2.0 3.0 - 2.0 3.0 4.0 - 3.0 4.0 5.0 ] +getproperty(UpperTriangularToeplitz,:uplo) = :U +getproperty(LowerTriangularToeplitz,:uplo) = :L ``` +This type is **obsolete** and will not be updated for features. Despite that, backward compatibility should be maintained. Codes that were using `TriangularToeplitz` should still work. -## Circulant +## Unexported interface +Methods in this section are not exported. -A circulant matrix is a special case of a Toeplitz matrix with periodic end conditions. -It can be constructed using +`_vr(A::AbstractMatrix)` returns the first row as a vector. +`_vc(A::AbstractMatrix)` returns the first column as a vector. +`_vr` and `_vc` are implemented for `AbstractToeplitz` as well. They are used to merge similar codes for `AbstractMatrix` and `AbstractToeplitz`. -```julia -Circulant(vc) -``` -where `vc` is a vector with the entries for the first column. -For example: -```julia -Circulant([1.0, 2.0, 3.0]) -``` -is a sparse representation of the matrix +`_circulate(v::AbstractVector)` converts between the `vr` and `vc` of a `Circulant`. -```julia -[ 1.0 3.0 2.0 - 2.0 1.0 3.0 - 3.0 2.0 1.0 ] -``` +`isconcrete(A::Union{AbstractToeplitz,Hankel})` decides whether the stored vector(s) are concrete. It calls `Base.isconcretetype`. diff --git a/src/ToeplitzMatrices.jl b/src/ToeplitzMatrices.jl index f22af6d..687cd4f 100644 --- a/src/ToeplitzMatrices.jl +++ b/src/ToeplitzMatrices.jl @@ -1,799 +1,50 @@ module ToeplitzMatrices -using StatsBase +# import StatsBase: levinson!, levinson +import DSP: conv +import Base: adjoint, convert, transpose, size, getindex, similar, copy, getproperty, inv, sqrt, copyto!, reverse, conj, zero, fill!, checkbounds, real, imag, isfinite, DimsInteger, iszero +import Base: ==, +, -, *, \ +import LinearAlgebra: Cholesky, Factorization +import LinearAlgebra: ldiv!, factorize, lmul!, pinv, eigvals, cholesky!, cholesky, tril!, triu!, checksquare, rmul!, dot, mul!, tril, triu +import LinearAlgebra: UpperTriangular, LowerTriangular, Symmetric, Adjoint +import AbstractFFTs: Plan, plan_fft! +import StatsBase -import Base: convert, *, \, getindex, print_matrix, size, Matrix, +, -, copy, similar, - sqrt, copyto!, adjoint, transpose -import LinearAlgebra: cholesky, cholesky!, eigvals, inv, ldiv!, mul!, pinv, tril, triu - -using LinearAlgebra: LinearAlgebra, Adjoint, Factorization, factorize, Cholesky, - DimensionMismatch, dot, rmul!, sym_uplo, Symmetric - -using AbstractFFTs -using AbstractFFTs: Plan - -using DSP: conv - -export Toeplitz, SymmetricToeplitz, Circulant, TriangularToeplitz, Hankel, chan, strang - -export durbin, durbin!, levinson, levinson!, trench, trench! - +export AbstractToeplitz, Toeplitz, SymmetricToeplitz, Circulant, LowerTriangularToeplitz, UpperTriangularToeplitz, TriangularToeplitz, Hankel +export durbin, trench, levinson include("iterativeLinearSolvers.jl") # Abstract abstract type AbstractToeplitz{T<:Number} <: AbstractMatrix{T} end +size(A::AbstractToeplitz) = (length(A.vc),length(A.vr)) +@inline _vr(A::AbstractToeplitz) = A.vr +@inline _vc(A::AbstractToeplitz) = A.vc +@inline _vr(A::AbstractMatrix) = A[1,:] +@inline _vc(A::AbstractMatrix) = A[:,1] + +AbstractArray{T}(A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) +convert(::Type{AbstractToeplitz{T}}, A::AbstractToeplitz) where T = AbstractToeplitz{T}(A) + +isconcrete(A::AbstractToeplitz) = isconcretetype(typeof(A.vc)) && isconcretetype(typeof(A.vr)) +iszero(A::AbstractToeplitz) = iszero(A.vc) && iszero(A.vr) + """ ToeplitzFactorization Factorization of a Toeplitz matrix using FFT. """ -struct ToeplitzFactorization{T<:Number,A<:AbstractToeplitz{T},S<:Number,P<:Plan{S}} <: Factorization{T} +struct ToeplitzFactorization{T,A<:AbstractToeplitz{T},S<:Number,P<:Plan{S}} <: Factorization{T} vcvr_dft::Vector{S} tmp::Vector{S} dft::P end -size(A::AbstractToeplitz) = (size(A, 1), size(A, 2)) - -convert(::Type{AbstractMatrix{T}}, S::AbstractToeplitz) where {T<:Number} = convert(AbstractToeplitz{T}, S) -convert(::Type{AbstractArray{T}}, S::AbstractToeplitz) where {T<:Number} = convert(AbstractToeplitz{T}, S) - -# Fast application of a general Toeplitz matrix to a column vector via FFT -function mul!( - y::StridedVector, A::AbstractToeplitz, x::StridedVector, α::Number, β::Number -) - m, n = size(A) - if length(y) != m - throw(DimensionMismatch( - "first dimension of A, $(m), does not match length of y, $(length(y))" - )) - end - if length(x) != n - throw(DimensionMismatch( - "second dimension of A, $(n), does not match length of x, $(length(x))" - )) - end - - # Small case: don't use FFT - N = m + n - 1 - if N < 512 - # Scale/initialize y - if iszero(β) - fill!(y, 0) - else - rmul!(y, β) - end - - @inbounds for j in 1:n - tmp = α * x[j] - for i in 1:m - y[i] = muladd(tmp, A[i,j], y[i]) - end - end - else - # Large case: use FFT - mul!(y, factorize(A), x, α, β) - end - - return y -end -function mul!( - y::StridedVector, A::ToeplitzFactorization, x::StridedVector, α::Number, β::Number -) - n = length(x) - m = length(y) - vcvr_dft = A.vcvr_dft - N = length(vcvr_dft) - if m > N || n > N - throw(DimensionMismatch( - "Toeplitz factorization does not match size of input and output vector" - )) - end - - T = Base.promote_eltype(y, A, x, α, β) - tmp = A.tmp - dft = A.dft - @inbounds begin - copyto!(tmp, 1, x, 1, n) - for i in (n+1):N - tmp[i] = zero(eltype(tmp)) - end - mul!(tmp, dft, tmp) - for i in 1:N - tmp[i] *= vcvr_dft[i] - end - dft \ tmp - if iszero(β) - for i in 1:m - y[i] = α * maybereal(T, tmp[i]) - end - else - for i in 1:m - y[i] = muladd(α, maybereal(T, tmp[i]), β * y[i]) - end - end - end - - return y -end - -# Application of a general Toeplitz matrix to a general matrix -function mul!( - C::StridedMatrix, A::AbstractToeplitz, B::StridedMatrix, α::Number, β::Number -) - return mul!(C, factorize(A), B, α, β) -end -function mul!( - C::StridedMatrix, A::ToeplitzFactorization, B::StridedMatrix, α::Number, β::Number -) - l = size(B, 2) - if size(C, 2) != l - throw(DimensionMismatch("input and output matrices must have same number of columns")) - end - for j = 1:l - mul!(view(C, :, j), A, view(B, :, j), α, β) - end - return C -end - -# Left division of a general matrix B by a general Toeplitz matrix A, i.e. the solution x of Ax=B. -function ldiv!(A::AbstractToeplitz, B::StridedMatrix) - if size(A, 1) != size(A, 2) - error("Division: Rectangular case is not supported.") - end - for j = 1:size(B, 2) - ldiv!(A, view(B, :, j)) - end - return B -end - -function (\)(A::AbstractToeplitz, b::AbstractVector) - T = promote_type(eltype(A), eltype(b)) - if T != eltype(A) - throw(ArgumentError("promotion of Toeplitz matrices not handled yet")) - end - bb = similar(b, T) - copyto!(bb, b) - ldiv!(A, bb) -end - -# General Toeplitz matrix -""" - Toeplitz - -A Toeplitz matrix. -""" -struct Toeplitz{T<:Number} <: AbstractToeplitz{T} - vc::Vector{T} - vr::Vector{T} - - function Toeplitz{T}(vc::Vector{T}, vr::Vector{T}) where {T<:Number} - if first(vc) != first(vr) - error("First element of the vectors must be the same") - end - return new{T}(vc, vr) - end -end - -""" - Toeplitz(vc::AbstractVector, vr::AbstractVector) - -Create a `Toeplitz` matrix from its first column `vc` and first row `vr` where -`vc[1] == vr[1]`. -""" -function Toeplitz(vc::AbstractVector, vr::AbstractVector) - return Toeplitz{Base.promote_eltype(vc, vr)}(vc, vr) -end -function Toeplitz{T}(vc::AbstractVector, vr::AbstractVector) where {T<:Number} - return Toeplitz{T}(convert(Vector{T}, vc), convert(Vector{T}, vr)) -end - -""" - Toeplitz(A::AbstractMatrix) - -"Project" matrix `A` onto its Toeplitz part using the first row/col of `A`. -""" -Toeplitz(A::AbstractMatrix) = Toeplitz{eltype(A)}(A) -function Toeplitz{T}(A::AbstractMatrix) where {T<:Number} - return Toeplitz{T}(A[:,1], A[1,:]) -end - -function LinearAlgebra.factorize(A::Toeplitz) - T = eltype(A) - m, n = size(A) - S = promote_type(T, Complex{Float32}) - tmp = Vector{S}(undef, m + n - 1) - copyto!(tmp, A.vc) - copyto!(tmp, m + 1, Iterators.reverse(A.vr), 1, n - 1) - dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) -end - -convert(::Type{AbstractToeplitz{T}}, A::Toeplitz) where {T} = convert(Toeplitz{T}, A) -convert(::Type{Toeplitz{T}}, A::Toeplitz) where {T} = Toeplitz(convert(Vector{T}, A.vc), - convert(Vector{T}, A.vr)) - -adjoint(A::Toeplitz) = Toeplitz(conj(A.vr), conj(A.vc)) -transpose(A::Toeplitz) = Toeplitz(A.vr, A.vc) - -# Size of a general Toeplitz matrix -function size(A::Toeplitz, dim::Int) - if dim == 1 - return length(A.vc) - elseif dim == 2 - return length(A.vr) - elseif dim > 2 - return 1 - else - error("arraysize: dimension out of range") - end -end - -# Retrieve an entry -function getindex(A::Toeplitz, i::Integer, j::Integer) - m, n = size(A) - @boundscheck if i < 1 || i > m || j < 1 || j > n - throw(BoundsError(A, (i,j))) - end - d = i - j - if d >= 0 - return A.vc[d + 1] - else - return A.vr[1 - d] - end -end - -# Form a lower triangular Toeplitz matrix by annihilating all entries above the k-th diaganal -function tril(A::Toeplitz, k::Integer = 0) - if k > 0 - error("Second argument cannot be positive") - end - Al = TriangularToeplitz(copy(A.vc), :L) - if k < 0 - for i in 1:-k - Al.ve[i] = zero(eltype(A)) - end - end - return Al -end - -# Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diagonal -function triu(A::Toeplitz, k::Integer = 0) - if k < 0 - error("Second argument cannot be negative") - end - Al = TriangularToeplitz(copy(A.vr), :U) - if k > 0 - for i in 1:k - Al.ve[i] = zero(eltype(A)) - end - end - return Al -end - -function ldiv!(A::Toeplitz, b::StridedVector) - preconditioner = factorize(strang(A)) - copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, preconditioner, 1000, 100eps())[1]) -end - -# Symmetric -""" - SymmetricToeplitz - -A symmetric Toeplitz matrix. -""" -struct SymmetricToeplitz{T<:Number} <: AbstractToeplitz{T} - vc::Vector{T} -end - -SymmetricToeplitz(vc::AbstractVector) = SymmetricToeplitz{eltype(vc)}(vc) - -SymmetricToeplitz(A::AbstractMatrix) = SymmetricToeplitz{eltype(A)}(A) -SymmetricToeplitz{T}(A::AbstractMatrix) where {T<:Number} = SymmetricToeplitz{T}(A[1, :]) - -Base.:*(a::Number, T::SymmetricToeplitz) = SymmetricToeplitz(a * T.vc) -Base.:*(T::SymmetricToeplitz, a::Number) = SymmetricToeplitz(T.vc * a) -function LinearAlgebra.lmul!(a::Number, T::SymmetricToeplitz) - T.vc .= a .* T.vc - return T -end -function LinearAlgebra.rmul!(T::SymmetricToeplitz, a::Number) - T.vc .*= a - return T -end - -function LinearAlgebra.factorize(A::SymmetricToeplitz{T}) where {T<:Number} - vc = A.vc - m = length(vc) - S = promote_type(T, Complex{Float32}) - tmp = Vector{S}(undef, 2 * m) - copyto!(tmp, vc) - @inbounds tmp[m + 1] = zero(T) - copyto!(tmp, m + 2, Iterators.reverse(vc), 1, m - 1) - dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) -end - -convert(::Type{AbstractToeplitz{T}}, A::SymmetricToeplitz) where {T} = convert(SymmetricToeplitz{T},A) -convert(::Type{SymmetricToeplitz{T}}, A::SymmetricToeplitz) where {T} = SymmetricToeplitz(convert(Vector{T},A.vc)) - -adjoint(A::SymmetricToeplitz) = SymmetricToeplitz(conj(A.vc)) -transpose(A::SymmetricToeplitz) = A - -function size(A::SymmetricToeplitz, dim::Int) - if dim < 1 - error("arraysize: dimension out of range") - end - return dim < 3 ? length(A.vc) : 1 -end - -function getindex(A::SymmetricToeplitz, i::Integer, j::Integer) - m = size(A, 1) - @boundscheck if i < 1 || i > m || j < 1 || j > m - throw(BoundsError(A, (i,j))) - end - return A.vc[abs(i - j) + 1] -end - -ldiv!(A::SymmetricToeplitz, b::StridedVector) = - copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1]) - -# Circulant -""" - Circulant - -A circulant matrix. -""" -struct Circulant{T<:Number} <: AbstractToeplitz{T} - vc::Vector{T} -end - -""" - Circulant(vc::AbstractVector{<:Number}) - -Create a circulant matrix from its first column `vc`. -""" -Circulant(vc::AbstractVector) = Circulant{eltype(vc)}(vc) - -""" - Circulant(A::AbstractMatrix) - -Create a circulant matrix from the first column of matrix `A`. -""" -Circulant(A::AbstractMatrix) = Circulant{eltype(A)}(A) -Circulant{T}(A::AbstractMatrix) where {T<:Number} = Circulant{T}(A[:,1]) - -const CirculantFactorization{T<:Number} = ToeplitzFactorization{T,Circulant{T}} -function LinearAlgebra.factorize(C::Circulant) - T = eltype(C) - vc = C.vc - S = promote_type(T, Complex{Float32}) - tmp = Vector{S}(undef, length(vc)) - copyto!(tmp, vc) - dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft) -end - -convert(::Type{AbstractToeplitz{T}}, A::Circulant) where {T} = convert(Circulant{T}, A) -convert(::Type{Circulant{T}}, A::Circulant) where {T} = Circulant(convert(Vector{T}, A.vc)) - - -function size(C::Circulant, dim::Integer) - if dim < 1 - error("arraysize: dimension out of range") - end - return dim < 3 ? length(C.vc) : 1 -end - -function getindex(C::Circulant, i::Integer, j::Integer) - n = size(C, 1) - @boundscheck if i < 1 || i > n || j < 1 || j > n - throw(BoundsError(C, (i,j))) - end - d = i - j - return C.vc[d < 0 ? n+d+1 : d+1] -end - -LinearAlgebra.ldiv!(C::Circulant, b::AbstractVector) = ldiv!(factorize(C), b) -function LinearAlgebra.ldiv!(C::CirculantFactorization, b::AbstractVector) - n = length(b) - tmp = C.tmp - vcvr_dft = C.vcvr_dft - if !(length(tmp) == length(vcvr_dft) == n) - throw(DimensionMismatch( - "size of Toeplitz factorization does not match the length of the output vector" - )) - end - dft = C.dft - @inbounds begin - copyto!(tmp, b) - dft * tmp - tmp ./= vcvr_dft - dft \ tmp - T = eltype(C) - b .= maybereal.(T, tmp) - end - return b -end - -function Base.inv(C::Circulant) - F = factorize(C) - vdft = map(inv, F.vcvr_dft) - vc = F.dft \ vdft - return Circulant(maybereal(eltype(C), vc)) -end -function Base.inv(C::CirculantFactorization{T,S,P}) where {T,S,P} - vdft = map(inv, C.vcvr_dft) - return CirculantFactorization{T,S,P}(vdft, similar(vdft), C.dft) -end - -function strang(A::AbstractMatrix{T}) where T - n = size(A, 1) - v = Vector{T}(undef, n) - n2 = div(n, 2) - for i = 1:n2 + 1 - v[i] = A[i,1] - end - for i in n2+2:n - v[i] = A[1, n - i + 2] - end - return Circulant(v) -end -function chan(A::AbstractMatrix{T}) where T - n = size(A, 1) - v = Vector{T}(undef, n) - for i = 1:n - v[i] = ((n - i + 1) * A[i, 1] + (i - 1) * A[1, min(n - i + 2, n)]) / n - end - return Circulant(v) -end - -function LinearAlgebra.pinv(C::Circulant, tol::Real=eps(real(float(one(eltype(C)))))) - F = factorize(C) - vdft = map(F.vcvr_dft) do x - z = inv(x) - return abs(x) < tol ? zero(z) : z - end - vc = F.dft \ vdft - return Circulant(maybereal(eltype(C), vc)) -end - -LinearAlgebra.eigvals(C::Circulant) = eigvals(factorize(C)) -LinearAlgebra.eigvals(C::CirculantFactorization) = copy(C.vcvr_dft) - -sqrt(C::Circulant) = sqrt(factorize(C)) -function Base.sqrt(C::CirculantFactorization) - vc = C.dft \ sqrt.(C.vcvr_dft) - return Circulant(maybereal(eltype(C), vc)) -end - -copy(C::Circulant) = Circulant(copy(C.vc)) -similar(C::Circulant) = Circulant(similar(C.vc)) -function copyto!(dest::Circulant, src::Circulant) - copyto!(dest.vc, src.vc) - return dest -end - -(+)(C1::Circulant, C2::Circulant) = Circulant(C1.vc + C2.vc) -(-)(C1::Circulant, C2::Circulant) = Circulant(C1.vc - C2.vc) -(-)(C::Circulant) = Circulant(-C.vc) - -Base.:*(A::Circulant, B::Circulant) = factorize(A) * factorize(B) -Base.:*(A::CirculantFactorization, B::Circulant) = A * factorize(B) -Base.:*(A::Circulant, B::CirculantFactorization) = factorize(A) * B -function Base.:*(A::CirculantFactorization, B::CirculantFactorization) - A_vcvr_dft = A.vcvr_dft - B_vcvr_dft = B.vcvr_dft - m = length(A_vcvr_dft) - n = length(B_vcvr_dft) - if m != n - throw(DimensionMismatch( - "size of matrix A, $(m)x$(m), does not match size of matrix B, $(n)x$(n)" - )) - end - - vc = A.dft \ (A_vcvr_dft .* B_vcvr_dft) - - return Circulant(maybereal(Base.promote_eltype(A, B), vc)) -end - -# Make an eager adjoint, similar to adjoints of Diagonal in LinearAlgebra -adjoint(C::CirculantFactorization{T,S,P}) where {T,S,P} = - CirculantFactorization{T,S,P}(conj.(C.vcvr_dft), C.tmp, C.dft) -Base.:*(A::Adjoint{<:Any,<:Circulant}, B::Circulant) = factorize(parent(A))' * factorize(B) -Base.:*(A::Adjoint{<:Any,<:Circulant}, B::CirculantFactorization) = - factorize(parent(A))' * B -Base.:*(A::Adjoint{<:Any,<:Circulant}, B::Adjoint{<:Any,<:Circulant}) = - factorize(parent(A))' * factorize(parent(B))' -Base.:*(A::Circulant, B::Adjoint{<:Any,<:Circulant}) = factorize(A) * factorize(parent(B))' -Base.:*(A::CirculantFactorization, B::Adjoint{<:Any,<:Circulant}) = A * factorize(parent(B))' - -(*)(scalar::Number, C::Circulant) = Circulant(scalar * C.vc) -(*)(C::Circulant,scalar::Number) = Circulant(C.vc * scalar) - - -# Triangular -struct TriangularToeplitz{T<:Number} <: AbstractToeplitz{T} - ve::Vector{T} - uplo::Char -end - -function TriangularToeplitz(ve::AbstractVector, uplo::Symbol) - return TriangularToeplitz{eltype(ve)}(ve, uplo) -end -function TriangularToeplitz{T}(ve::AbstractVector, uplo::Symbol) where {T<:Number} - UL = LinearAlgebra.char_uplo(uplo) - return TriangularToeplitz{T}(convert(Vector{T}, ve), UL) -end - -function TriangularToeplitz(A::AbstractMatrix, uplo::Symbol) - return TriangularToeplitz{eltype(A)}(A, uplo) -end -function TriangularToeplitz{T}(A::AbstractMatrix, uplo::Symbol) where {T<:Number} - ve = uplo === :U ? A[1, :] : A[:, 1] - return TriangularToeplitz{T}(ve, uplo) -end - -function LinearAlgebra.factorize(A::TriangularToeplitz) - T = eltype(A) - ve = A.ve - n = length(ve) - S = promote_type(T, Complex{Float32}) - tmp = zeros(S, 2 * n - 1) - if A.uplo === 'L' - copyto!(tmp, ve) - else - tmp[1] = ve[1] - copyto!(tmp, n + 1, Iterators.reverse(ve), 1, n - 1) - end - dft = plan_fft!(tmp) - return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) -end - -function convert(::Type{Toeplitz}, A::TriangularToeplitz) - if A.uplo == 'L' - Toeplitz(A.ve, [A.ve[1]; zeros(length(A.ve) - 1)]) - else - @assert A.uplo == 'U' - Toeplitz([A.ve[1]; zeros(length(A.ve) - 1)], A.ve) - end -end - -convert(::Type{AbstractToeplitz{T}}, A::TriangularToeplitz) where {T} = convert(TriangularToeplitz{T},A) -convert(::Type{TriangularToeplitz{T}}, A::TriangularToeplitz) where {T} = - TriangularToeplitz(convert(Vector{T}, A.ve), A.uplo) - - -function size(A::TriangularToeplitz, dim::Int) - if dim < 1 - error("arraysize: dimension out of range") - end - return dim < 3 ? length(A.ve) : 1 -end - -function getindex(A::TriangularToeplitz{T}, i::Integer, j::Integer) where T - n = size(A, 1) - @boundscheck if i < 1 || i > n || j < 1 || j > n - throw(BoundsError(A, (i,j))) - end - if A.uplo == 'L' - return i >= j ? A.ve[i - j + 1] : zero(T) - else - return i <= j ? A.ve[j - i + 1] : zero(T) - end -end - -function (*)(A::TriangularToeplitz, B::TriangularToeplitz) - n = size(A, 2) - if n != size(B, 1) - throw(DimensionMismatch("")) - end - if A.uplo == B.uplo - return TriangularToeplitz(conv(A.ve, B.ve)[1:n], A.uplo) - end - return A * Matrix(B) -end - -adjoint(A::TriangularToeplitz) = TriangularToeplitz(conj(A.ve), A.uplo == 'U' ? (:L) : (:U)) -transpose(A::TriangularToeplitz) = TriangularToeplitz(A.ve, A.uplo == 'U' ? (:L) : (:U)) - -# NB! only valid for lower triangular -function smallinv(A::TriangularToeplitz{T}) where T - n = size(A, 1) - b = zeros(T, n) - b[1] = 1 ./ A.ve[1] - for k = 2:n - tmp = zero(T) - for i = 1:k-1 - tmp += A.uplo == 'L' ? A.ve[k - i + 1]*b[i] : A.ve[i + 1] * b[k - i] - end - b[k] = -tmp/A.ve[1] - end - return TriangularToeplitz(b, A.uplo) -end - -function inv(A::TriangularToeplitz{T}) where T - n = size(A, 1) - if n <= 64 - return smallinv(A) - end - np2 = nextpow(2, n) - if n != np2 - return TriangularToeplitz( - inv(TriangularToeplitz([A.ve; zeros(T, np2 - n)], sym_uplo(A.uplo))).ve[1:n], - sym_uplo(A.uplo) - ) - end - nd2 = div(n, 2) - a1 = inv(TriangularToeplitz(A.ve[1:nd2], sym_uplo(A.uplo))).ve - return TriangularToeplitz( - [a1; -(TriangularToeplitz(a1, sym_uplo(A.uplo)) * (Toeplitz(A.ve[nd2 + 1:end], A.ve[nd2 + 1:-1:2]) * a1))], - sym_uplo(A.uplo) - ) -end - -# ldiv!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b -function ldiv!(A::TriangularToeplitz, b::StridedVector) - preconditioner = factorize(chan(A)) - copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, preconditioner, 1000, 100eps())[1]) -end - -# extend levinson -StatsBase.levinson!(x::StridedVector, A::SymmetricToeplitz, b::StridedVector) = - StatsBase.levinson!(A.vc, b, x) -function StatsBase.levinson!(C::StridedMatrix, A::SymmetricToeplitz, B::StridedMatrix) - n = size(B, 2) - if n != size(C, 2) - throw(DimensionMismatch("input and output matrices must have same number of columns")) - end - for j = 1:n - StatsBase.levinson!(view(C, :, j), A, view(B, :, j)) - end - C -end -StatsBase.levinson(A::AbstractToeplitz, B::StridedVecOrMat) = - StatsBase.levinson!(zeros(size(B)), A, copy(B)) - -# BlockTriangular -# type BlockTriangularToeplitz{T<:BlasReal} <: AbstractMatrix{T} -# Mc::Array{T,3} -# uplo::Char -# Mc_dft::Array{Complex{T},3} -# tmp::Vector{Complex{T}} -# dft::Plan -# end -# function BlockTriangularToeplitz{T<:BlasReal}(Mc::Array{T,3}, uplo::Symbol) -# n, p, _ = size(Mc) -# tmp = zeros(Complex{T}, 2n) -# dft = plan_fft!(tmp) -# Mc_dft = Array{Complex{T}}(2n, p, p) -# for j = 1:p -# for i = 1:p -# Mc_dft[1,i,j] = complex(Mc[1,i,j]) -# for t = 2:n -# Mc_dft[t,i,j] = uplo == :L ? complex(Mc[t,i,j]) : zero(Complex{T}) -# end -# Mc_dft[n+1,i,j] = zero(Complex{T}) -# for t = n+2:2n -# Mc_dft[t,i,j] = uplo == :L ? zero(Complex{T}) : complex(Mc[2n-t+2,i,j]) -# end -# dft(view(Mc_dft, 1:2n, 1:p, 1:p)) -# end -# end -# return BlockTriangularToeplitz(Mc, string(uplo)[1], Mc_dft, tmp, dft, idft) -# end - -#= Hankel Matrix - A Hankel matrix is a matrix that is constant across the anti-diagonals: - - [a_0 a_1 a_2 a_3 a_4 - a_1 a_2 a_3 a_4 a_5 - a_2 a_3 a_4 a_5 a_6] - - This is precisely a Toeplitz matrix with the columns reversed: - [0 0 0 0 1 - [a_4 a_3 a_2 a_1 a_0 0 0 0 1 0 - a_5 a_4 a_3 a_2 a_1 * 0 0 1 0 0 - a_6 a_5 a_4 a_3 a_2] 0 1 0 0 0 - 1 0 0 0 0] - We represent the Hankel matrix by wrapping the corresponding Toeplitz matrix.=# - - -# Hankel Matrix, use _Hankel as Hankel(::Toeplitz) should project to Hankel -function _Hankel end - -# Hankel Matrix -struct Hankel{TT<:Number} <: AbstractMatrix{TT} - T::Toeplitz{TT} - global _Hankel(T::Toeplitz{TT}) where {TT<:Number} = new{TT}(T) -end - -# Ctor: vc is the leftmost column and vr is the bottom row. -function Hankel{T}(vc::AbstractVector, vr::AbstractVector) where T - if vc[end] != vr[1] - error("First element of rows must equal last element of columns") - end - n = length(vr) - p = [vc; vr[2:end]] - _Hankel(Toeplitz{T}(p[n:end],p[n:-1:1])) -end - -Hankel(vc::AbstractVector, vr::AbstractVector) = - Hankel{promote_type(eltype(vc), eltype(vr))}(vc, vr) - -Hankel{T}(A::AbstractMatrix) where T = Hankel{T}(A[:,1], A[end,:]) -Hankel(A::AbstractMatrix) = Hankel(A[:,1], A[end,:]) - -convert(::Type{AbstractArray{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) -convert(::Type{AbstractMatrix{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) -convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = _Hankel(convert(Toeplitz{T}, A.T)) - - - -# Size -size(H::Hankel, k...) = size(H.T, k...) - - -# Retrieve an entry by two indices -function getindex(A::Hankel, i::Integer, j::Integer) - m, n = size(A) - @boundscheck if i < 1 || i > m || j < 1 || j > n - throw(BoundsError(A, (i,j))) - end - return A.T[i,n-j+1] -end - -# Fast application of a general Hankel matrix to a strided vector -*(A::Hankel, b::StridedVector) = A.T * reverse(b) -mul!(y::StridedVector, A::Hankel, x::StridedVector, α::Number, β::Number) = - mul!(y, A.T, view(x, reverse(axes(x, 1))), α, β) -# Fast application of a general Hankel matrix to a strided matrix -*(A::Hankel, B::StridedMatrix) = A.T * reverse(B, dims=1) -mul!(Y::StridedMatrix, A::Hankel, X::StridedMatrix, α::Number, β::Number) = - mul!(Y, A.T, view(X, reverse(axes(X, 1)), :), α, β) - -## BigFloat support - -(*)(A::Toeplitz{T}, b::AbstractVector) where {T<:BigFloat} = irfft( - rfft([ - A.vc; - reverse(A.vr[2:end])] - ) .* rfft([ - b; - zeros(length(b) - 1) - ]), - 2 * length(b) - 1 -)[1:length(b)] - -function cholesky!(L::AbstractMatrix, T::SymmetricToeplitz) - - L[:, 1] .= T.vc ./ sqrt(T.vc[1]) - v = copy(L[:, 1]) - N = size(T, 1) - - @inbounds for n in 1:N-1 - sinθn = v[n + 1] / L[n, n] - cosθn = sqrt(1 - sinθn^2) - - for n′ in n+1:N - v[n′] = (v[n′] - sinθn * L[n′ - 1, n]) / cosθn - L[n′, n + 1] = -sinθn * v[n′] + cosθn * L[n′ - 1, n] - end - end - return Cholesky(L, 'L', 0) -end - -""" - cholesky(T::SymmetricToeplitz) - -Implementation of the Bareiss Algorhithm, adapted from "On the stability of the Bareiss and -related Toeplitz factorization algorithms", Bojanczyk et al, 1993. -""" -function cholesky(T::SymmetricToeplitz) - return cholesky!(Matrix{eltype(T)}(undef, size(T, 1), size(T, 1)), T) -end +include("toeplitz.jl") +include("special.jl") +include("hankel.jl") +include("linearalgebra.jl") """ maybereal(::Type{T}, x) diff --git a/src/directLinearSolvers.jl b/src/directLinearSolvers.jl index c16f049..a496dc1 100644 --- a/src/directLinearSolvers.jl +++ b/src/directLinearSolvers.jl @@ -69,7 +69,7 @@ end # uses B to store the inverse # NOTE: only populates entries on upper triangle since the inverse is symmetric function trench!(B::AbstractMatrix, r::AbstractVector) - n = LinearAlgebra.checksquare(B) + n = checksquare(B) n == length(r) + 1 || throw(DimensionMismatch()) y = durbin(r) γ = inv(1 + dot(r, y)) diff --git a/src/hankel.jl b/src/hankel.jl new file mode 100644 index 0000000..43fefd5 --- /dev/null +++ b/src/hankel.jl @@ -0,0 +1,132 @@ +# Hankel +struct Hankel{T, V<:AbstractVector{T}, S<:DimsInteger} <: AbstractMatrix{T} + v::V + size::S # size + + function Hankel{T,V,S}(v::V, (m,n)::DimsInteger) where {T, V<:AbstractVector{T}, S<:DimsInteger} + (m < 0 || n < 0) && throw(ArgumentError("negative size: $s")) + length(v) ≥ m+n-1 || throw(ArgumentError("inconsistency between size and number of anti-diagonals")) + new{T,V,S}(v, (m,n)) + end +end +Hankel{T}(v::V, s::S) where {T, V<:AbstractVector{T}, S<:DimsInteger} = Hankel{T,V,S}(v,s) + +Hankel{T}(v::AbstractVector, s::DimsInteger) where T = Hankel{T}(convert(AbstractVector{T},v),s) +Hankel{T}(v::AbstractVector, h::Integer, w::Integer) where T = Hankel{T}(v,(h,w)) +Hankel(v::AbstractVector, s::DimsInteger) = Hankel{eltype(v)}(v,s) +Hankel(v::AbstractVector, h::Integer, w::Integer) = Hankel{eltype(v)}(v,h,w) +Hankel(v::AbstractVector) = (l=length(v);Hankel(v,((l+1)÷2,(l+1)÷2))) # square by default +function Hankel(vc::AbstractVector, vr::AbstractVector) + if vc[end] != vr[1] + throw(ArgumentError("vc[end] != vr[1]")) + end + Hankel(vcat(vc,vr[2:end]), (length(vc), length(vr))) +end + +function getproperty(A::Hankel, s::Symbol) + m,_ = getfield(A, :size) + if s == :vc + A.v[1:m] + elseif s == :vr + A.v[m:end] + else + getfield(A,s) + end +end + +# convert from general matrix to Hankel matrix +Hankel{T}(A::AbstractMatrix, uplo::Symbol = :L) where T<:Number = convert(Hankel{T}, Hankel(A,uplo)) +# using the first column and last row +function Hankel(A::AbstractMatrix, uplo::Symbol = :L) + m,n = size(A) + if uplo == :L + if isfinite(m) # InfiniteArrays.jl supports infinite + Hankel(A[:,1], A[end,:]) + else + Hankel(A[:,1], (m,n)) + end + elseif uplo == :U + if isfinite(n) + Hankel(A[1,:], A[:,end]) + else + Hankel(A[1,:], (m,n)) + end + else + throw(ArgumentError("expected :L or :U. got $uplo.")) + end +end + +convert(::Type{AbstractArray{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) +convert(::Type{AbstractMatrix{T}}, A::Hankel) where {T<:Number} = convert(Hankel{T}, A) +convert(::Type{Hankel{T}}, A::Hankel) where {T<:Number} = Hankel{T}(convert(AbstractVector{T}, A.v), size(A)) + +# Size +size(H::Hankel) = H.size + +# Retrieve an entry by two indices +Base.@propagate_inbounds function getindex(A::Hankel, i::Integer, j::Integer) + @boundscheck checkbounds(A, i, j) + return A.v[i+j-1] +end +similar(A::Hankel, T::Type, dims::DimsInteger{2}) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) +similar(A::Hankel, T::Type, dims::Tuple{Int64,Int64}) = Hankel{T}(similar(A.v, T, dims[1]+dims[2]-true), dims) # for ambiguity with `similar(a::AbstractArray, ::Type{T}, dims::Tuple{Vararg{Int64, N}}) where {T, N}` in Base +for fun in (:zero, :conj, :copy, :-, :similar, :real, :imag) + @eval $fun(A::Hankel) = Hankel($fun(A.v), size(A)) +end +for op in (:+, :-) + @eval function $op(A::Hankel, B::Hankel) + promote_shape(A,B) + Hankel($op(A.v,B.v), size(A)) + end +end +function copyto!(A::Hankel, B::Hankel) + promote_shape(A,B) + copyto!(A.v,B.v) +end +for fun in (:lmul!,) + @eval function $fun(x::Number, A::Hankel) + $fun(x, A.v) + A + end +end +for fun in (:fill!, :rmul!) + @eval function $fun(A::Hankel, x::Number) + $fun(A.v, x) + A + end +end + +transpose(A::Hankel) = Hankel(A.v, reverse(size(A))) +adjoint(A::Hankel) = transpose(conj(A)) +(==)(A::Hankel, B::Hankel) = A.v == B.v && size(A) == size(B) +(*)(scalar::Number, C::Hankel) = Hankel(scalar * C.v, size(C)) +(*)(C::Hankel,scalar::Number) = Hankel(C.v * scalar, size(C)) + +isconcrete(A::Hankel) = isconcretetype(A.v) + +function reverse(A::Hankel; dims=1) + _,n = size(A) + if dims==1 + Toeplitz(reverse(A.vc), A.vr) + elseif dims==2 + Toeplitz(A.v[n:end], A.v[n:-1:1]) + else + throw(ArgumentError("invalid dimension $dims in reverse")) + end +end +function reverse(A::AbstractToeplitz; dims=1) + if dims==1 + Hankel(reverse(A.vc), A.vr) + elseif dims==2 + Hankel(vcat(reverse(A.vr), A.vc), size(A)) + else + throw(ArgumentError("invalid dimension $dims in reverse")) + end +end + +# Fast application of a general Hankel matrix to a strided vector +*(A::Hankel, b::StridedVector) = reverse(A, dims=2) * reverse(b) +mul!(y::StridedVector, A::Hankel, x::StridedVector, α::Number, β::Number) = mul!(y, reverse(A,dims=2), view(x, reverse(axes(x, 1))), α, β) +# Fast application of a general Hankel matrix to a strided matrix +*(A::Hankel, B::StridedMatrix) = reverse(A,dims=2) * reverse(B, dims=1) +mul!(Y::StridedMatrix, A::Hankel, X::StridedMatrix, α::Number, β::Number) = mul!(Y, reverse(A,dims=2), view(X, reverse(axes(X, 1)), :), α, β) \ No newline at end of file diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl new file mode 100644 index 0000000..4f02d9d --- /dev/null +++ b/src/linearalgebra.jl @@ -0,0 +1,413 @@ +# general Toeplitz + +# Fast application of a general Toeplitz matrix to a column vector via FFT +function mul!( + y::StridedVector, A::AbstractToeplitz, x::StridedVector, α::Number, β::Number +) + m, n = size(A) + if length(y) != m + throw(DimensionMismatch( + "first dimension of A, $(m), does not match length of y, $(length(y))" + )) + end + if length(x) != n + throw(DimensionMismatch( + "second dimension of A, $(n), does not match length of x, $(length(x))" + )) + end + + # Small case: don't use FFT + N = m + n - 1 + if N < 512 + # Scale/initialize y + if iszero(β) + fill!(y, 0) + else + rmul!(y, β) + end + + @inbounds for j in 1:n + tmp = α * x[j] + for i in 1:m + y[i] = muladd(tmp, A[i,j], y[i]) + end + end + else + # Large case: use FFT + mul!(y, factorize(A), x, α, β) + end + + return y +end +function mul!( + y::StridedVector, A::ToeplitzFactorization, x::StridedVector, α::Number, β::Number +) + n = length(x) + m = length(y) + vcvr_dft = A.vcvr_dft + N = length(vcvr_dft) + if m > N || n > N + throw(DimensionMismatch( + "Toeplitz factorization does not match size of input and output vector" + )) + end + + T = Base.promote_eltype(y, A, x, α, β) + tmp = A.tmp + dft = A.dft + @inbounds begin + copyto!(tmp, 1, x, 1, n) + for i in (n+1):N + tmp[i] = zero(eltype(tmp)) + end + mul!(tmp, dft, tmp) + for i in 1:N + tmp[i] *= vcvr_dft[i] + end + dft \ tmp + if iszero(β) + for i in 1:m + y[i] = α * maybereal(T, tmp[i]) + end + else + for i in 1:m + y[i] = muladd(α, maybereal(T, tmp[i]), β * y[i]) + end + end + end + + return y +end + +# Application of a general Toeplitz matrix to a general matrix +function mul!( + C::StridedMatrix, A::AbstractToeplitz, B::StridedMatrix, α::Number, β::Number +) + return mul!(C, factorize(A), B, α, β) +end +function mul!( + C::StridedMatrix, A::ToeplitzFactorization, B::StridedMatrix, α::Number, β::Number +) + l = size(B, 2) + if size(C, 2) != l + throw(DimensionMismatch("input and output matrices must have same number of columns")) + end + for j = 1:l + mul!(view(C, :, j), A, view(B, :, j), α, β) + end + return C +end + +# Left division of a general matrix B by a general Toeplitz matrix A, i.e. the solution x of Ax=B. +function ldiv!(A::AbstractToeplitz, B::StridedMatrix) + if size(A, 1) != size(A, 2) + error("Division: Rectangular case is not supported.") + end + for j = 1:size(B, 2) + ldiv!(A, view(B, :, j)) + end + return B +end + +function (\)(A::AbstractToeplitz, b::AbstractVector) + T = promote_type(eltype(A), eltype(b)) + if T != eltype(A) + throw(ArgumentError("promotion of Toeplitz matrices not handled yet")) + end + bb = similar(b, T) + copyto!(bb, b) + ldiv!(A, bb) +end + +function factorize(A::Toeplitz) + T = eltype(A) + m, n = size(A) + S = promote_type(T, Complex{Float32}) + tmp = Vector{S}(undef, m + n - 1) + copyto!(tmp, A.vc) + copyto!(tmp, m + 1, Iterators.reverse(A.vr), 1, n - 1) + dft = plan_fft!(tmp) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) +end + +function ldiv!(A::Toeplitz, b::StridedVector) + preconditioner = factorize(strang(A)) + copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, preconditioner, 1000, 100eps())[1]) +end + +StatsBase.levinson(A::AbstractToeplitz, B::AbstractVecOrMat) = StatsBase.levinson!(zeros(size(B)), A, copy(B)) + +# SymmetricToeplitz + +function factorize(A::SymmetricToeplitz{T}) where {T<:Number} + vc = A.vc + m = length(vc) + S = promote_type(T, Complex{Float32}) + tmp = Vector{S}(undef, 2 * m) + copyto!(tmp, vc) + @inbounds tmp[m + 1] = zero(T) + copyto!(tmp, m + 2, Iterators.reverse(vc), 1, m - 1) + dft = plan_fft!(tmp) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) +end + +ldiv!(A::SymmetricToeplitz, b::StridedVector) = copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1]) + +function cholesky!(L::AbstractMatrix, T::SymmetricToeplitz) + + L[:, 1] .= T.vc ./ sqrt(T.vc[1]) + v = copy(L[:, 1]) + N = size(T, 1) + + @inbounds for n in 1:N-1 + sinθn = v[n + 1] / L[n, n] + cosθn = sqrt(1 - sinθn^2) + + for n′ in n+1:N + v[n′] = (v[n′] - sinθn * L[n′ - 1, n]) / cosθn + L[n′, n + 1] = -sinθn * v[n′] + cosθn * L[n′ - 1, n] + end + end + return Cholesky(L, 'L', 0) +end + +""" + cholesky(T::SymmetricToeplitz) + +Implementation of the Bareiss Algorhithm, adapted from "On the stability of the Bareiss and +related Toeplitz factorization algorithms", Bojanczyk et al, 1993. +""" +function cholesky(T::SymmetricToeplitz) + return cholesky!(Matrix{eltype(T)}(undef, size(T, 1), size(T, 1)), T) +end + +# extend levinson +StatsBase.levinson!(x::StridedVector, A::SymmetricToeplitz, b::StridedVector) = StatsBase.levinson!(A.vc, b, x) +function StatsBase.levinson!(C::StridedMatrix, A::SymmetricToeplitz, B::StridedMatrix) + n = size(B, 2) + if n != size(C, 2) + throw(DimensionMismatch("input and output matrices must have same number of columns")) + end + for j = 1:n + StatsBase.levinson!(view(C, :, j), A, view(B, :, j)) + end + C +end + +# circulant +const CirculantFactorization{T, V<:AbstractVector{T}} = ToeplitzFactorization{T,Circulant{T,V}} +function factorize(C::Circulant) + T = eltype(C) + vc = C.vc + S = promote_type(T, Complex{Float32}) + tmp = Vector{S}(undef, length(vc)) + copyto!(tmp, vc) + dft = plan_fft!(tmp) + return ToeplitzFactorization{T,typeof(C),S,typeof(dft)}(dft * tmp, similar(tmp), dft) +end + +Base.:*(A::Circulant, B::Circulant) = factorize(A) * factorize(B) +Base.:*(A::CirculantFactorization, B::Circulant) = A * factorize(B) +Base.:*(A::Circulant, B::CirculantFactorization) = factorize(A) * B +function Base.:*(A::CirculantFactorization, B::CirculantFactorization) + A_vcvr_dft = A.vcvr_dft + B_vcvr_dft = B.vcvr_dft + m = length(A_vcvr_dft) + n = length(B_vcvr_dft) + if m != n + throw(DimensionMismatch( + "size of matrix A, $(m)x$(m), does not match size of matrix B, $(n)x$(n)" + )) + end + + vc = A.dft \ (A_vcvr_dft .* B_vcvr_dft) + + return Circulant(maybereal(Base.promote_eltype(A, B), vc)) +end + +# Make an eager adjoint, similar to adjoints of Diagonal in LinearAlgebra +adjoint(C::CirculantFactorization{T,V,S,P}) where {T,V,S,P} = + CirculantFactorization{T,V,S,P}(conj.(C.vcvr_dft), C.tmp, C.dft) +Base.:*(A::Adjoint{<:Any,<:Circulant}, B::Circulant) = factorize(parent(A))' * factorize(B) +Base.:*(A::Adjoint{<:Any,<:Circulant}, B::CirculantFactorization) = + factorize(parent(A))' * B +Base.:*(A::Adjoint{<:Any,<:Circulant}, B::Adjoint{<:Any,<:Circulant}) = + factorize(parent(A))' * factorize(parent(B))' +Base.:*(A::Circulant, B::Adjoint{<:Any,<:Circulant}) = factorize(A) * factorize(parent(B))' +Base.:*(A::CirculantFactorization, B::Adjoint{<:Any,<:Circulant}) = A * factorize(parent(B))' + +ldiv!(C::Circulant, b::AbstractVector) = ldiv!(factorize(C), b) +function ldiv!(C::CirculantFactorization, b::AbstractVector) + n = length(b) + tmp = C.tmp + vcvr_dft = C.vcvr_dft + if !(length(tmp) == length(vcvr_dft) == n) + throw(DimensionMismatch( + "size of Toeplitz factorization does not match the length of the output vector" + )) + end + dft = C.dft + @inbounds begin + copyto!(tmp, b) + dft * tmp + tmp ./= vcvr_dft + dft \ tmp + T = eltype(C) + b .= maybereal.(T, tmp) + end + return b +end + +function inv(C::Circulant) + F = factorize(C) + vdft = map(inv, F.vcvr_dft) + vc = F.dft \ vdft + return Circulant(maybereal(eltype(C), vc)) +end +function inv(C::CirculantFactorization{T,V,S,P}) where {T,V,S,P} + vdft = map(inv, C.vcvr_dft) + return CirculantFactorization{T,V,S,P}(vdft, similar(vdft), C.dft) +end + +function strang(A::AbstractMatrix{T}) where T + n = size(A, 1) + v = Vector{T}(undef, n) + n2 = div(n, 2) + for i = 1:n2 + 1 + v[i] = A[i,1] + end + for i in n2+2:n + v[i] = A[1, n - i + 2] + end + return Circulant(v) +end +function chan(A::AbstractMatrix{T}) where T + n = size(A, 1) + v = Vector{T}(undef, n) + for i = 1:n + v[i] = ((n - i + 1) * A[i, 1] + (i - 1) * A[1, min(n - i + 2, n)]) / n + end + return Circulant(v) +end + +function pinv(C::Circulant, tol::Real=eps(real(float(one(eltype(C)))))) + F = factorize(C) + vdft = map(F.vcvr_dft) do x + z = inv(x) + return abs(x) < tol ? zero(z) : z + end + vc = F.dft \ vdft + return Circulant(maybereal(eltype(C), vc)) +end + +eigvals(C::Circulant) = eigvals(factorize(C)) +eigvals(C::CirculantFactorization) = copy(C.vcvr_dft) + +sqrt(C::Circulant) = sqrt(factorize(C)) +function sqrt(C::CirculantFactorization) + vc = C.dft \ sqrt.(C.vcvr_dft) + return Circulant(maybereal(eltype(C), vc)) +end + +# Triangular + +# NB! only valid for lower triangular +function smallinv(A::LowerTriangularToeplitz{T}) where T + n = size(A, 1) + b = zeros(T, n) + b[1] = 1 ./ A.v[1] + for k = 2:n + tmp = zero(T) + for i = 1:k-1 + tmp += A.v[k - i + 1]*b[i] + end + b[k] = -tmp/A.v[1] + end + return LowerTriangularToeplitz(b) +end +function inv(A::LowerTriangularToeplitz{T}) where T + n = size(A, 1) + if n <= 64 + return smallinv(A) + end + np2 = nextpow(2, n) + if n != np2 + return LowerTriangularToeplitz( + inv(LowerTriangularToeplitz([A.v; zeros(T, np2 - n)])).v[1:n]) + end + nd2 = div(n, 2) + a1 = inv(LowerTriangularToeplitz(A.v[1:nd2])).v + return LowerTriangularToeplitz( + [a1; -(LowerTriangularToeplitz(a1) * (Toeplitz(A.v[nd2 + 1:end], A.v[nd2 + 1:-1:2]) * a1))] + ) +end + +function smallinv(A::UpperTriangularToeplitz{T}) where T + n = size(A, 1) + b = zeros(T, n) + b[1] = 1 ./ A.v[1] + for k = 2:n + tmp = zero(T) + for i = 1:k-1 + tmp += A.v[i + 1] * b[k - i] + end + b[k] = -tmp/A.v[1] + end + return UpperTriangularToeplitz(b) +end +function inv(A::UpperTriangularToeplitz{T}) where T + n = size(A, 1) + if n <= 64 + return smallinv(A) + end + np2 = nextpow(2, n) + if n != np2 + return UpperTriangularToeplitz( + inv(UpperTriangularToeplitz([A.v; zeros(T, np2 - n)])).v[1:n] + ) + end + nd2 = div(n, 2) + a1 = inv(UpperTriangularToeplitz(A.v[1:nd2])).v + return UpperTriangularToeplitz( + [a1; -(UpperTriangularToeplitz(a1) * (Toeplitz(A.v[nd2 + 1:end], A.v[nd2 + 1:-1:2]) * a1))] + ) +end + +# ldiv!(A::TriangularToeplitz,b::StridedVector) = inv(A)*b +function ldiv!(A::TriangularToeplitz, b::StridedVector) + preconditioner = factorize(chan(A)) + copyto!(b, IterativeLinearSolvers.cgs(A, zeros(eltype(b), length(b)), b, preconditioner, 1000, 100eps())[1]) +end + +function (*)(A::TriangularToeplitz, B::TriangularToeplitz) + n = size(A, 2) + if n != size(B, 1) + throw(DimensionMismatch("")) + end + if A.uplo == B.uplo + return TriangularToeplitz(conv(A.v, B.v)[1:n], A.uplo) + end + return A * Matrix(B) +end + +function factorize(A::LowerTriangularToeplitz) + T = eltype(A) + v = A.v + n = length(v) + S = promote_type(T, Complex{Float32}) + tmp = zeros(S, 2 * n - 1) + copyto!(tmp, v) + dft = plan_fft!(tmp) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) +end +function factorize(A::UpperTriangularToeplitz) + T = eltype(A) + v = A.v + n = length(v) + S = promote_type(T, Complex{Float32}) + tmp = zeros(S, 2 * n - 1) + tmp[1] = v[1] + copyto!(tmp, n + 1, Iterators.reverse(v), 1, n - 1) + dft = plan_fft!(tmp) + return ToeplitzFactorization{T,typeof(A),S,typeof(dft)}(dft * tmp, similar(tmp), dft) +end \ No newline at end of file diff --git a/src/special.jl b/src/special.jl new file mode 100644 index 0000000..6e271c8 --- /dev/null +++ b/src/special.jl @@ -0,0 +1,228 @@ +# special Toeplitz types that can be represented by a single vector +# Symmetric, Circulant, LowerTriangular, UpperTriangular +for TYPE in (:SymmetricToeplitz, :Circulant, :LowerTriangularToeplitz, :UpperTriangularToeplitz) + @eval begin + struct $TYPE{T, V<:AbstractVector{T}} <: AbstractToeplitz{T} + v::V + end + $TYPE{T}(v::V) where {T,V<:AbstractVector{T}} = $TYPE{T,V}(v) + $TYPE{T}(v::AbstractVector) where T = $TYPE{T}(convert(AbstractVector{T},v)) + + AbstractToeplitz{T}(A::$TYPE) where T = $TYPE{T}(A) + $TYPE{T}(A::$TYPE) where T = $TYPE{T}(convert(AbstractVector{T},A.v)) + convert(::Type{$TYPE{T}}, A::$TYPE) where {T} = $TYPE{T}(A) + + size(A::$TYPE) = (length(A.v),length(A.v)) + + adjoint(A::$TYPE) = transpose(conj(A)) + (*)(scalar::Number, C::$TYPE) = $TYPE(scalar * C.v) + (*)(C::$TYPE, scalar::Number) = $TYPE(C.v * scalar) + (==)(A::$TYPE,B::$TYPE) = A.v==B.v + function zero!(A::$TYPE) + if isconcrete(A) + fill!(A.v,zero(eltype(A))) + else + A.v=zero(A.v) + end + end + + function copyto!(A::$TYPE,B::$TYPE) + copyto!(A.v,B.v) + A + end + similar(A::$TYPE, T::Type) = $TYPE{T}(similar(A.v, T)) + function lmul!(x::Number, A::$TYPE) + lmul!(x,A.v) + A + end + function rmul!(A::$TYPE, x::Number) + rmul!(A.v,x) + A + end + end + for fun in (:zero, :conj, :copy, :-, :real, :imag) + @eval $fun(A::$TYPE) = $TYPE($fun(A.v)) + end + for fun in (:iszero,) + @eval $fun(A::$TYPE) = $fun(A.v) + end + for op in (:+, :-) + @eval $op(A::$TYPE,B::$TYPE) = $TYPE($op(A.v,B.v)) + end + for TY in (:AbstractMatrix, :AbstractToeplitz) + @eval $TYPE(v::$TY) = $TYPE{eltype(v)}(v) + end +end +TriangularToeplitz{T}=Union{UpperTriangularToeplitz{T},LowerTriangularToeplitz{T}} + +# vc and vr +function getproperty(A::SymmetricToeplitz, s::Symbol) + if s == :vc || s == :vr + getfield(A,:v) + else + getfield(A,s) + end +end +function getproperty(A::Circulant, s::Symbol) + if s == :vc + getfield(A,:v) + elseif s == :vr + _circulate(getfield(A,:v)) + else + getfield(A,s) + end +end +function getproperty(A::LowerTriangularToeplitz, s::Symbol) + if s == :vc + getfield(A,:v) + elseif s == :vr + _firstnonzero(getfield(A,:v)) + elseif s == :uplo + :L + else + getfield(A,s) + end +end +function getproperty(A::UpperTriangularToeplitz, s::Symbol) + if s == :vr + getfield(A,:v) + elseif s == :vc + _firstnonzero(getfield(A,:v)) + elseif s == :uplo + :U + else + getfield(A,s) + end +end + +_circulate(v::AbstractVector) = vcat(v[1],v[end:-1:2]) +_firstnonzero(v::AbstractVector) = vcat(v[1],zero(view(v,2:lastindex(v)))) + +# transpose +transpose(A::SymmetricToeplitz) = A +transpose(A::Circulant) = Circulant(A.vr) +transpose(A::LowerTriangularToeplitz) = UpperTriangularToeplitz(A.v) +transpose(A::UpperTriangularToeplitz) = LowerTriangularToeplitz(A.v) + +# getindex +Base.@propagate_inbounds function getindex(A::SymmetricToeplitz, i::Integer, j::Integer) + @boundscheck checkbounds(A,i,j) + return A.v[abs(i - j) + 1] +end +Base.@propagate_inbounds function getindex(C::Circulant, i::Integer, j::Integer) + @boundscheck checkbounds(C,i,j) + d = i - j + return C.v[d < 0 ? size(C,1)+d+1 : d+1] +end +Base.@propagate_inbounds function getindex(A::LowerTriangularToeplitz{T}, i::Integer, j::Integer) where T + @boundscheck checkbounds(A,i,j) + return i >= j ? A.v[i - j + 1] : zero(T) +end +Base.@propagate_inbounds function getindex(A::UpperTriangularToeplitz{T}, i::Integer, j::Integer) where T + @boundscheck checkbounds(A,i,j) + return i <= j ? A.v[j - i + 1] : zero(T) +end + +# constructors +function Circulant{T}(v::AbstractVector, uplo::Symbol) where T + if uplo == :L + Circulant{T}(v) + elseif uplo == :U + Circulant{T}(_circulate(v)) + else + throw(ArgumentError("expected :L or :U; got $uplo.")) + end +end +Circulant(v::AbstractVector, uplo::Symbol) = Circulant{eltype(v)}(v,uplo) +# from AbstractMatrix, AbstractToeplitz +function Circulant{T}(A::AbstractMatrix, uplo::Symbol = :L) where T + checksquare(A) + if uplo == :L + Circulant{T}(_vc(A), uplo) + elseif uplo == :U + Circulant{T}(_vr(A), uplo) + else + throw(ArgumentError("expected :L or :U; get $uplo.")) + end +end +Circulant(A::AbstractMatrix, uplo::Symbol) = Circulant{eltype(A)}(A,uplo) +function SymmetricToeplitz{T}(A::AbstractMatrix, uplo::Symbol = :U) where T + checksquare(A) + if uplo == :L + SymmetricToeplitz{T}(_vc(A)) + elseif uplo == :U + SymmetricToeplitz{T}(_vr(A)) + else + throw(ArgumentError("expected :L or :U; got $uplo.")) + end +end +SymmetricToeplitz(A::AbstractMatrix, uplo::Symbol) = SymmetricToeplitz{eltype(A)}(A,uplo) +Symmetric(A::AbstractToeplitz, uplo::Symbol = :U) = SymmetricToeplitz(A,uplo) + +function UpperTriangularToeplitz{T}(A::AbstractMatrix) where T + checksquare(A) + UpperTriangularToeplitz{T}(_vr(A)) +end +function LowerTriangularToeplitz{T}(A::AbstractMatrix) where T + checksquare(A) + LowerTriangularToeplitz{T}(_vc(A)) +end +_toeplitztype(s::Symbol) = Symbol(s,"Toeplitz") +for TYPE in (:UpperTriangular, :LowerTriangular) + @eval begin + $TYPE{T}(A::AbstractToeplitz) where T = $(_toeplitztype(TYPE)){T}(A) + $TYPE(A::AbstractToeplitz) = $TYPE{eltype(A)}(A) + convert(::Type{TriangularToeplitz{T}},A::$(_toeplitztype(TYPE))) where T<:Number = convert($(_toeplitztype(TYPE)){T},A) + end +end + +# Triangular +for TYPE in (:AbstractMatrix, :AbstractVector) + @eval begin + TriangularToeplitz(A::$TYPE, uplo::Symbol) = TriangularToeplitz{eltype(A)}(A, uplo) + function TriangularToeplitz{T}(A::$TYPE, uplo::Symbol) where T + if uplo == :L + LowerTriangularToeplitz{T}(A) + elseif uplo == :U + UpperTriangularToeplitz{T}(A) + else + throw(ArgumentError("expected :L or :U. got $uplo.")) + end + end + end +end + +# tril and triu +tril(A::Union{SymmetricToeplitz,Circulant}, k::Integer=0) = tril!(Toeplitz(A),k) +triu(A::Union{SymmetricToeplitz,Circulant}, k::Integer=0) = triu!(Toeplitz(A),k) +function _tridiff!(A::TriangularToeplitz, k::Integer) + if k >= 0 + if isconcretetype(typeof(A.v)) + for i in k+2:lastindex(A.v) + A.v[i] = zero(eltype(A)) + end + else + A.v = vcat(A.v[1:k+1], zero(A.v[k+2:end])) + end + else + zero!(A) + end + A +end +tril!(A::UpperTriangularToeplitz, k::Integer) = _tridiff!(A,k) +triu!(A::LowerTriangularToeplitz, k::Integer) = _tridiff!(A,-k) + +function _trisame!(A::TriangularToeplitz, k::Integer) + if k < 0 + if isconcretetype(typeof(A.v)) + for i in 1:-k + A.v[i] = zero(eltype(A)) + end + else + A.v=vcat(A.v[1:-k+1], zero(A.v[-k+2:end])) + end + end + A +end +tril!(A::LowerTriangularToeplitz, k::Integer) = _trisame!(A,k) +triu!(A::UpperTriangularToeplitz, k::Integer) = _trisame!(A,-k) diff --git a/src/toeplitz.jl b/src/toeplitz.jl new file mode 100644 index 0000000..5bb97fd --- /dev/null +++ b/src/toeplitz.jl @@ -0,0 +1,146 @@ +# General Toeplitz matrix +""" + Toeplitz + +A Toeplitz matrix. +""" +struct Toeplitz{T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} <: AbstractToeplitz{T} + vc::VC + vr::VR + + function Toeplitz{T, VC, VR}(vc::VC, vr::VR) where {T, VC<:AbstractVector{T}, VR<:AbstractVector{T}} + if first(vc) != first(vr) + error("First element of the vectors must be the same") + end + return new{T,VC,VR}(vc, vr) + end +end +Toeplitz{T}(vc::AbstractVector{T}, vr::AbstractVector{T}) where T = Toeplitz{T,typeof(vc),typeof(vr)}(vc,vr) + +""" + Toeplitz(vc::AbstractVector, vr::AbstractVector) + +Create a `Toeplitz` matrix from its first column `vc` and first row `vr` where +`vc[1] == vr[1]`. +""" +function Toeplitz(vc::AbstractVector, vr::AbstractVector) + return Toeplitz{Base.promote_eltype(vc, vr)}(vc, vr) +end +function Toeplitz{T}(vc::AbstractVector, vr::AbstractVector) where T + return Toeplitz{T}(convert(AbstractVector{T}, vc), convert(AbstractVector{T}, vr)) +end + +""" + Toeplitz(A::AbstractMatrix) + +"Project" matrix `A` onto its Toeplitz part using the first row/col of `A`. +""" +Toeplitz(A::AbstractMatrix) = Toeplitz{eltype(A)}(A) +Toeplitz{T}(A::AbstractMatrix) where {T} = Toeplitz{T}(copy(_vc(A)), copy(_vr(A))) + +AbstractToeplitz{T}(A::Toeplitz) where T = Toeplitz{T}(A) +convert(::Type{Toeplitz{T}}, A::AbstractToeplitz) where {T} = Toeplitz{T}(A) +convert(::Type{Toeplitz}, A::AbstractToeplitz) = Toeplitz(A) + +# Retrieve an entry +Base.@propagate_inbounds function getindex(A::AbstractToeplitz, i::Integer, j::Integer) + @boundscheck checkbounds(A,i,j) + d = i - j + if d >= 0 + return A.vc[d + 1] + else + return A.vr[1 - d] + end +end + +checknonaliased(A::Toeplitz) = Base.mightalias(A.vc, A.vr) && throw(ArgumentError("Cannot modify Toeplitz matrices in place with aliased data")) + +function tril!(A::Toeplitz, k::Integer=0) + checknonaliased(A) + + if k >= 0 + if isconcretetype(typeof(A.vr)) + for i in k+2:lastindex(A.vr) + A.vr[i] = zero(eltype(A)) + end + else + A.vr=vcat(A.vr[1:k+1], zero(A.vr[k+2:end])) + end + else + fill!(A.vr, zero(eltype(A))) + if isconcretetype(typeof(A.vc)) + for i in 1:-k + A.vc[i]=zero(eltype(A)) + end + else + A.vc=vcat(zero(A.vc[1:-k]), A.vc[-k+1:end]) + end + end + A +end +function triu!(A::Toeplitz, k::Integer=0) + checknonaliased(A) + + if k <= 0 + if isconcretetype(typeof(A.vc)) + for i in -k+2:lastindex(A.vc) + A.vc[i] = zero(eltype(A)) + end + else + A.vc=vcat(A.vc[1:-k+1], zero(A.vc[-k+2:end])) + end + else + fill!(A.vc, zero(eltype(A))) + if isconcretetype(typeof(A.vr)) + for i in 1:k + A.vr[i]=zero(eltype(A)) + end + else + A.vr=vcat(zero(A.vr[1:k]), A.vr[k+1:end]) + end + end + A +end + +adjoint(A::AbstractToeplitz) = transpose(conj(A)) +transpose(A::AbstractToeplitz) = Toeplitz(A.vr, A.vc) +function similar(A::AbstractToeplitz, T::Type, dims::Dims{2}) + vc=similar(A.vc, T, dims[1]) + vr=similar(A.vr, T, dims[2]) + vr[1]=vc[1] + Toeplitz{T}(vc,vr) +end +for fun in (:zero, :conj, :copy, :-, :real, :imag) + @eval $fun(A::AbstractToeplitz)=Toeplitz($fun(A.vc),$fun(A.vr)) +end +for op in (:+, :-) + @eval $op(A::AbstractToeplitz,B::AbstractToeplitz)=Toeplitz($op(A.vc,B.vc),$op(A.vr,B.vr)) +end +function copyto!(A::Toeplitz, B::AbstractToeplitz) + checknonaliased(A) + copyto!(A.vc,B.vc) + copyto!(A.vr,B.vr) + A +end +(==)(A::AbstractToeplitz,B::AbstractToeplitz) = A.vr==B.vr && A.vc==B.vc + +function fill!(A::Toeplitz, x::Number) + fill!(A.vc,x) + fill!(A.vr,x) + A +end +(*)(scalar::Number, C::AbstractToeplitz) = Toeplitz(scalar * C.vc, scalar * C.vr) +(*)(C::AbstractToeplitz,scalar::Number) = Toeplitz(C.vc * scalar, C.vr * scalar) + +function lmul!(x::Number, A::Toeplitz) + checknonaliased(A) + lmul!(x, A.vc) + lmul!(x, A.vr) + A +end +function rmul!(A::Toeplitz, x::Number) + checknonaliased(A) + rmul!(A.vc, x) + rmul!(A.vr, x) + A +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 79ce02d..4e0580f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,8 @@ using Pkg Pkg.instantiate() end -using ToeplitzMatrices, StatsBase, Test, LinearAlgebra, Aqua +using ToeplitzMatrices, Test, LinearAlgebra, Aqua +import StatsBase using FFTW: fft @@ -39,17 +40,17 @@ cases = [ (Circulant(complex(0.9.^(0:ns - 1))), Circulant(complex(0.9.^(0:nl - 1))), "Complex circulant"), - (TriangularToeplitz(0.9.^(0:ns - 1), :U), - TriangularToeplitz(0.9.^(0:nl - 1), :U), + (UpperTriangularToeplitz(0.9.^(0:ns - 1)), + UpperTriangularToeplitz(0.9.^(0:nl - 1)), "Real upper triangular"), - (TriangularToeplitz(complex(0.9.^(0:ns - 1)), :U), - TriangularToeplitz(complex(0.9.^(0:nl - 1)), :U), + (UpperTriangularToeplitz(complex(0.9.^(0:ns - 1))), + UpperTriangularToeplitz(complex(0.9.^(0:nl - 1))), "Complex upper triangular"), - (TriangularToeplitz(0.9.^(0:ns - 1), :L), - TriangularToeplitz(0.9.^(0:nl - 1), :L), + (LowerTriangularToeplitz(0.9.^(0:ns - 1)), + LowerTriangularToeplitz(0.9.^(0:nl - 1)), "Real lower triangular"), - (TriangularToeplitz(complex(0.9.^(0:ns - 1)), :L), - TriangularToeplitz(complex(0.9.^(0:nl - 1)), :L), + (LowerTriangularToeplitz(complex(0.9.^(0:ns - 1))), + LowerTriangularToeplitz(complex(0.9.^(0:nl - 1))), "Complex lower triangular"), ] @@ -74,7 +75,7 @@ end @test Circulant(Float32.(1:3)) * ones(Float64, 3) == fill(6, 3) @test Matrix(Toeplitz(vc, vr)) == Matrix(Toeplitz(vv, vr)) @test Matrix(Circulant(vc)) == Matrix(Circulant(vv)) - @test Matrix(TriangularToeplitz(vc,:U)) == Matrix(TriangularToeplitz(vv,:U)) + @test Matrix(UpperTriangularToeplitz(vc)) == Matrix(UpperTriangularToeplitz(vv)) end @testset "Real general rectangular" begin @@ -214,7 +215,14 @@ end @test Hankel(A) == [1 3; 3 4] T = Toeplitz([1.0,2,3,4,5],[1.0,6,7,8,0]) @test Hankel(T) == Hankel([1.0,2,3,4,5],[5.0,4,3,2,1]) - @test Hankel(T) ≠ ToeplitzMatrices._Hankel(T) + @test isa(reverse(T),Hankel) + @test isa(reverse(T,dims=1),Hankel) + @test isa(reverse(T,dims=2),Hankel) + end + + @testset "v too small" begin + @test_throws ArgumentError Hankel(Int[], (3,4)) + @test_throws ArgumentError Hankel(1:5, (3,4)) end end @@ -275,25 +283,32 @@ end @testset "Constructors" begin A = ones(10, 10) - @test Matrix(Toeplitz(A)) == Matrix(Toeplitz{Float64}(A)) == A - @test Matrix(SymmetricToeplitz(A)) == Matrix(SymmetricToeplitz{Float64}(A)) == A - @test Matrix(Circulant(A)) == Matrix(Circulant{Float64}(A)) == A - @test Matrix(Hankel(A)) == Matrix(Hankel{Float64}(A)) == A + @test Toeplitz(A) == Toeplitz{Float64}(A) == A + @test SymmetricToeplitz(A) == SymmetricToeplitz{Float64}(A) == A + @test Circulant(A) == Circulant{Float64}(A) == A + @test Hankel(A) == Hankel{Float64}(A) == A A = [1.0 2.0; 3.0 4.0] - @test Toeplitz(A) == Toeplitz([1.,3.], [1.,2.]) - @test Toeplitz{Float64}(A) == Toeplitz([1.,3.], [1.,2.]) - @test Matrix(SymmetricToeplitz(A)) == Matrix(SymmetricToeplitz{Float64}(A)) == - Matrix(Toeplitz(Symmetric(A))) == Matrix(Symmetric(Toeplitz(A))) == [1. 2.; 2. 1.] - @test Matrix(Circulant(A)) == [1 3; 3 1] + @test Toeplitz(A) == Toeplitz([1.,3.], [1.,2.]) == Toeplitz{Float64}(A) + @test SymmetricToeplitz(A) == SymmetricToeplitz{Float64}(A) == + Toeplitz(Symmetric(A)) == Symmetric(Toeplitz(A)) == [1. 2.; 2. 1.] + @test Circulant(A,:L) == [1 3; 3 1] == Circulant(A) == SymmetricToeplitz(A,:L) + @test Circulant(A,:U) == [1 2; 2 1] == SymmetricToeplitz(A,:U) - @test TriangularToeplitz(A, :U) == TriangularToeplitz{Float64}(A, :U) == Toeplitz(UpperTriangular(A)) == UpperTriangular(Toeplitz(A)) - @test TriangularToeplitz(A, :L) == TriangularToeplitz{Float64}(A, :L) == Toeplitz(LowerTriangular(A)) == LowerTriangular(Toeplitz(A)) + @test TriangularToeplitz(A, :U) == TriangularToeplitz{Float64}(A, :U) == Toeplitz(UpperTriangular(A)) == UpperTriangular(Toeplitz(A)) == UpperTriangularToeplitz(A) == UpperTriangularToeplitz{Float64}(A) + @test TriangularToeplitz(A, :L) == TriangularToeplitz{Float64}(A, :L) == Toeplitz(LowerTriangular(A)) == LowerTriangular(Toeplitz(A)) == LowerTriangularToeplitz(A) == LowerTriangularToeplitz{Float64}(A) - @test Matrix(Hankel(A)) == Matrix(Hankel{Float64}(A)) == [1.0 3; 3 4] + @test Hankel(A) == Hankel{Float64}(A) == [1.0 3; 3 4] == Hankel([1.0,3],[3,4]) == Hankel([1.0,3,4],(2,2)) == Hankel([1.0,3,4],2,2) == Hankel{Float64}([1,3,4],(2,2)) == Hankel{Float64}([1.0,3,4],2,2) == Hankel([1.0,3,4]) + @test Hankel(A,:U) == [1.0 2;2 4] + + @test_throws ArgumentError Hankel(A,:🤣) + @test_throws ArgumentError SymmetricToeplitz(A,:🤣) + @test_throws ArgumentError Circulant(A,:🤣) + @test_throws ArgumentError Circulant(1:5,:🤣) + @test_throws ArgumentError TriangularToeplitz(A,:🤣) # Constructors should be projections @test Toeplitz(Toeplitz(A)) == Toeplitz(A) @@ -302,6 +317,62 @@ end @test TriangularToeplitz(TriangularToeplitz(A, :U), :U) == TriangularToeplitz(A, :U) @test TriangularToeplitz(TriangularToeplitz(A, :L), :L) == TriangularToeplitz(A, :L) @test Hankel(Hankel(A)) == Hankel(A) + + @test_throws ArgumentError Hankel(1:2,1:2) + @test_throws ErrorException Toeplitz(1:2,2:1) +end + +@testset "General Interface" begin + for Toep in (:Toeplitz, :Circulant, :SymmetricToeplitz, :UpperTriangularToeplitz, :LowerTriangularToeplitz, :Hankel) + @eval (A = [1.0 3.0; 3.0 4.0]; TA=$Toep(A); A = Matrix(TA)) + @eval (B = [2 1 ; 1 5 ]; TB=$Toep(B); B = Matrix(TB)) + + for fun in (:zero, :conj, :copy, :-, :real, :imag, :adjoint, :transpose, :iszero, :size) + @eval @test $fun(TA) == $fun(A) + end + + @test 2*TA == 2*A == lmul!(2,copy(TA)) + @test TA*2 == A*2 == rmul!(copy(TA),2) + @test TA+TB == A+B + @test TA-TB == A-B + + @test_throws ArgumentError reverse(TA,dims=3) + if isa(TA,AbstractToeplitz) + @test isa(reverse(TA),Hankel) + @test isa(reverse(TA,dims=1),Hankel) + @test isa(reverse(TA,dims=2),Hankel) + @test isa(tril(TA),AbstractToeplitz) && tril(TA)==tril(A) + @test isa(triu(TA),AbstractToeplitz) && triu(TA)==triu(A) + @test isa(tril(TA,1),AbstractToeplitz) && tril(TA,1)==tril(A,1) + @test isa(triu(TA,1),AbstractToeplitz) && triu(TA,1)==triu(A,1) + @test isa(tril(TA,-1),AbstractToeplitz) && tril(TA,-1)==tril(A,-1) + @test isa(triu(TA,-1),AbstractToeplitz) && triu(TA,-1)==triu(A,-1) + else + @test isa(reverse(TA),Toeplitz) + @test isa(reverse(TA,dims=1),Toeplitz) + @test isa(reverse(TA,dims=2),Toeplitz) + end + + T=copy(TA) + copyto!(T,TB) + @test T == B + + T=copy(TA) + end + @test fill!(Toeplitz(zeros(2,2)),1) == ones(2,2) + + @testset "aliasing" begin + v = [1,2,3] + T = Toeplitz(v, v) + @test_throws ArgumentError triu!(T) + @test_throws ArgumentError tril!(T) + @test_throws ArgumentError copyto!(T, Toeplitz(1:3, 1:3)) + @test_throws ArgumentError lmul!(2, T) + @test_throws ArgumentError rmul!(T, 2) + + @test triu(T) == triu(Matrix(T)) + @test tril(T) == tril(Matrix(T)) + end end @testset "Circulant mathematics" begin @@ -399,30 +470,30 @@ end @testset "TriangularToeplitz" begin A = [1.0 2.0; 3.0 4.0] - TU = TriangularToeplitz(A, :U) - TL = TriangularToeplitz(A, :L) - @test (TU * TU)::TriangularToeplitz ≈ Matrix(TU)*Matrix(TU) - @test (TL * TL)::TriangularToeplitz ≈ Matrix(TL)*Matrix(TL) + TU = UpperTriangularToeplitz(A) + TL = LowerTriangularToeplitz(A) + @test (TU * TU)::UpperTriangularToeplitz ≈ Matrix(TU)*Matrix(TU) + @test (TL * TL)::LowerTriangularToeplitz ≈ Matrix(TL)*Matrix(TL) @test (TU * TL) ≈ Matrix(TU)*Matrix(TL) for T in (TU, TL) @test inv(T)::TriangularToeplitz ≈ inv(Matrix(T)) end A = randn(ComplexF64, 3, 3) T = Toeplitz(A) - TU = triu(T) + TU = UpperTriangular(T) @test TU isa TriangularToeplitz @test istriu(TU) - @test TU == Toeplitz(triu(A)) + @test TU == Toeplitz(triu(A)) == triu(T) @test TU'ones(3) == Matrix(TU)'ones(3) @test transpose(TU)*ones(3) == transpose(Matrix(TU))*ones(3) - @test triu(T, 1)::TriangularToeplitz == triu(Matrix(T), 1) - TL = tril(T) + @test triu(TU, 1)::TriangularToeplitz == triu(Matrix(T), 1) == triu(T,1) + TL = LowerTriangular(T) @test TL isa TriangularToeplitz @test istril(TL) - @test TL == Toeplitz(tril(A)) + @test TL == Toeplitz(tril(A)) == tril(T) @test TL'ones(3) == Matrix(TL)'ones(3) @test transpose(TL)*ones(3) == transpose(Matrix(TL))*ones(3) - @test tril(T, -1)::TriangularToeplitz == tril(Matrix(T), -1) + @test tril(TL, -1)::TriangularToeplitz == tril(Matrix(T), -1) == tril(T,-1) for n in (65, 128) A = randn(n, n) TU = TriangularToeplitz(A, :U)