Skip to content

Commit

Permalink
Support Toeplitz(A::AbstractMatrix), SymmetricToepliz(A::AbstractMatr…
Browse files Browse the repository at this point in the history 
…ix), etc.. for projecting on to Toeplitz matrices.
  • Loading branch information
@dlfivefifty
dlfivefifty committed Nov 27, 2017
1 parent ff88029 commit e97fdf4
Show file tree
Hide file tree
Showing 2 changed files with 136 additions and 65 deletions.
106 changes: 71 additions & 35 deletions src/ToeplitzMatrices.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,13 +143,12 @@ mutable struct Toeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
end

# Ctor
function Toeplitz(vc::Vector, vr::Vector)
function Toeplitz{T}(vc::Vector, vr::Vector) where {T}
m, n = length(vc), length(vr)
if vc[1] != vr[1]
error("First element of the vectors must be the same")
end

T = promote_type(eltype(vc), eltype(vr), Float32)
vcp, vrp = Vector{T}(vc), Vector{T}(vr)

tmp = Vector{promote_type(T, Complex{Float32})}(m + n - 1)
Expand All @@ -161,9 +160,17 @@ function Toeplitz(vc::Vector, vr::Vector)
return Toeplitz(vcp, vrp, 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))
Toeplitz(vc::Vector, vr::Vector) =
Toeplitz{promote_type(eltype(vc), eltype(vr), Float32)}(vc, vr)

# Toeplitz(A::AbstractMatrix) projects onto Toeplitz part using the first row/col
Toeplitz(A::AbstractMatrix) = Toeplitz(A[:,1], A[1,:])
Toeplitz{T}(A::AbstractMatrix) where T = Toeplitz{T}(A[:,1], A[1,:])


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))

# Size of a general Toeplitz matrix
function size(A::Toeplitz, dim::Int)
Expand Down Expand Up @@ -230,15 +237,26 @@ mutable struct SymmetricToeplitz{T<:BlasReal} <: AbstractToeplitz{T}
vcvr_dft::Vector{Complex{T}}
tmp::Vector{Complex{T}}
dft::Plan

function SymmetricToeplitz{T}(vc::Vector{T}) where T<:BlasReal
tmp = convert(Array{Complex{T}}, [vc; zero(T); reverse(vc[2:end])])
dft = plan_fft!(tmp)
return new(vc, dft*tmp, similar(tmp), dft)
end
end
function SymmetricToeplitz(vc::Vector{T}) where T<:BlasReal
tmp = convert(Array{Complex{T}}, [vc; zero(T); reverse(vc[2:end])])
dft = plan_fft!(tmp)
return SymmetricToeplitz(vc, 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))
SymmetricToeplitz{T}(vc::AbstractVector) where T<:BlasReal = SymmetricToeplitz{T}(convert(Vector{T}, vc))
SymmetricToeplitz{T}(vc::AbstractVector{T}) where T = SymmetricToeplitz{promote_type(Float32, T)}(vc)
SymmetricToeplitz{T}(vc::AbstractVector) where T = SymmetricToeplitz{T}(convert(Vector{T}, vc))
SymmetricToeplitz(vc::AbstractVector{T}) where T = SymmetricToeplitz{T}(vc)

SymmetricToeplitz{T}(A::AbstractMatrix) where T = SymmetricToeplitz{T}(A[1, :])
SymmetricToeplitz(A::AbstractMatrix) = SymmetricToeplitz{eltype(A)}(A)



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))

function size(A::SymmetricToeplitz, dim::Int)
if 1 <= dim <= 2
Expand All @@ -261,19 +279,20 @@ mutable struct Circulant{T<:Number,S<:Number} <: AbstractToeplitz{T}
dft::Plan
end

function Circulant(vc::Vector{T}) where T<:Real
TT = promote_type(eltype(vc), Float32)
tmp = zeros(promote_type(TT, Complex{Float32}), length(vc))
return Circulant(real(vc), fft(vc), tmp, plan_fft!(tmp))
end
function Circulant(vc::Vector)
T = promote_type(eltype(vc), Float32)
function Circulant{T}(vc::Vector{T}) where T
tmp = zeros(promote_type(T, Complex{Float32}), length(vc))
return Circulant(vc, fft(vc), tmp, plan_fft!(tmp))
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))
Circulant{T}(vc::AbstractVector) where T = Circulant{T}(convert(Vector{T}, vc))
Circulant(vc::AbstractVector) = Circulant{promote_type(eltype(vc), Float32)}(vc)
Circulant{T}(A::AbstractMatrix) where T = Circulant{T}(A[:,1])
Circulant(A::AbstractMatrix) = Circulant(A[:,1])



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)
Expand Down Expand Up @@ -366,25 +385,35 @@ mutable struct TriangularToeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
dft::Plan
end

function TriangularToeplitz(ve::Vector, uplo::Symbol)
n = length(ve)

T = promote_type(eltype(ve), Float32)
vep = Vector{T}(ve)
function TriangularToeplitz{T}(vep::Vector{T}, uplo::Symbol) where T
n = length(vep)

tmp = zeros(promote_type(T, Complex{Float32}), 2n - 1)
if uplo == :L
copy!(tmp, vep)
else
tmp[1] = vep[1]
for i = 1:n - 1
tmp[n + i] = ve[n - i + 1]
tmp[n + i] = vep[n - i + 1]
end
end
dft = plan_fft!(tmp)
return TriangularToeplitz(vep, string(uplo)[1], dft * tmp, similar(tmp), dft)
end

TriangularToeplitz{T}(ve::AbstractVector, uplo::Symbol) where T =
TriangularToeplitz(convert(Vector{T}, ve), uplo)

TriangularToeplitz(ve::AbstractVector, uplo::Symbol) =
TriangularToeplitz{promote_type(eltype(ve), Float32)}(ve, uplo)

TriangularToeplitz{T}(A::AbstractMatrix, uplo::Symbol) where T =
TriangularToeplitz{T}(uplo == :U ? A[1,:] : A[:,1], uplo)

TriangularToeplitz(A::AbstractMatrix, uplo::Symbol) =
TriangularToeplitz(uplo == :U ? A[1,:] : A[:,1], uplo)


function convert(::Type{Toeplitz}, A::TriangularToeplitz)
if A.uplo == 'L'
Toeplitz(A.ve, [A.ve[1]; zeros(length(A.ve) - 1)])
Expand All @@ -394,8 +423,8 @@ function convert(::Type{Toeplitz}, A::TriangularToeplitz)
end
end

convert(::Type{AbstractToeplitz{T}},A::TriangularToeplitz) where {T} = convert(TriangularToeplitz{T},A)
convert(::Type{TriangularToeplitz{T}},A::TriangularToeplitz) where {T} =
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=='U' ? (:U) : (:L))


Expand Down Expand Up @@ -527,22 +556,29 @@ StatsBase.levinson(A::AbstractToeplitz, B::StridedVecOrMat) =
We represent the Hankel matrix by wrapping the corresponding Toeplitz matrix.=#

# Hankel Matrix
mutable struct Hankel{T<:Number} <: AbstractMatrix{T}
T::Toeplitz{T}
mutable struct Hankel{TT<:Number} <: AbstractMatrix{TT}
T::Toeplitz{TT}
Hankel{TT}(T::Toeplitz{TT}) where TT<:Number = new{TT}(T)
end

# Ctor: vc is the leftmost column and vr is the bottom row.
function Hankel(vc,vr)
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(p[n:end],p[n:-1:1]))
Hankel{T}(Toeplitz{T}(p[n:end],p[n:-1:1]))
end

convert(::Type{Array},A::Hankel) = convert(Matrix,A)
convert(::Type{Matrix},A::Hankel) = full(A)
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{Array}, A::Hankel) = convert(Matrix, A)
convert(::Type{Matrix}, A::Hankel) = full(A)

convert(::Type{AbstractMatrix{T}},A::Hankel) where {T} = convert(Hankel{T},A)
convert(::Type{Hankel{T}},A::Hankel) where {T} = Hankel(convert(Toeplitz{T},A.T))
Expand Down
95 changes: 65 additions & 30 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,36 +8,37 @@ xs = randn(ns, 5)
xl = randn(nl, 5)

@testset("Toeplitz: $st",
for (As, Al, st) in ((Toeplitz(0.9.^(0:ns-1), 0.4.^(0:ns-1)),
Toeplitz(0.9.^(0:nl-1), 0.4.^(0:nl-1)),
"Real general square"),
(Toeplitz(complex(0.9.^(0:ns-1)), complex(0.4.^(0:ns-1))),
Toeplitz(complex(0.9.^(0:nl-1)), complex(0.4.^(0:nl-1))),
"Complex general square"),
(Circulant(0.9.^(0:ns - 1)),
Circulant(0.9.^(0:nl - 1)),
"Real circulant"),
(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),
"Real upper triangular"),
(TriangularToeplitz(complex(0.9.^(0:ns - 1)), :U),
TriangularToeplitz(complex(0.9.^(0:nl - 1)), :U),
"Complex upper triangular"),
(TriangularToeplitz(0.9.^(0:ns - 1), :L),
TriangularToeplitz(0.9.^(0:nl - 1), :L),
"Real lower triangular"),
(TriangularToeplitz(complex(0.9.^(0:ns - 1)), :L),
TriangularToeplitz(complex(0.9.^(0:nl - 1)), :L),
"Complex lower triangular"))

@test As * xs ≈ full(As) * xs
@test Al * xl ≈ full(Al) * xl
@test A_ldiv_B!(As, LinAlg.copy_oftype(xs, eltype(As))) ≈ full(As) \ xs
@test A_ldiv_B!(Al, LinAlg.copy_oftype(xl, eltype(Al))) ≈ full(Al) \ xl
end)
for (As, Al, st) in ((Toeplitz(0.9.^(0:ns-1), 0.4.^(0:ns-1)),
Toeplitz(0.9.^(0:nl-1), 0.4.^(0:nl-1)),
"Real general square"),
(Toeplitz(complex(0.9.^(0:ns-1)), complex(0.4.^(0:ns-1))),
Toeplitz(complex(0.9.^(0:nl-1)), complex(0.4.^(0:nl-1))),
"Complex general square"),
(Circulant(0.9.^(0:ns - 1)),
Circulant(0.9.^(0:nl - 1)),
"Real circulant"),
(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),
"Real upper triangular"),
(TriangularToeplitz(complex(0.9.^(0:ns - 1)), :U),
TriangularToeplitz(complex(0.9.^(0:nl - 1)), :U),
"Complex upper triangular"),
(TriangularToeplitz(0.9.^(0:ns - 1), :L),
TriangularToeplitz(0.9.^(0:nl - 1), :L),
"Real lower triangular"),
(TriangularToeplitz(complex(0.9.^(0:ns - 1)), :L),
TriangularToeplitz(complex(0.9.^(0:nl - 1)), :L),
"Complex lower triangular"))

@test As * xs ≈ full(As) * xs
@test Al * xl ≈ full(Al) * xl
@test A_ldiv_B!(As, LinAlg.copy_oftype(xs, eltype(As))) ≈ full(As) \ xs
@test A_ldiv_B!(Al, LinAlg.copy_oftype(xl, eltype(Al))) ≈ full(Al) \ xl
end
)

@testset "Real general rectangular" begin
Ar1 = Toeplitz(0.9.^(0:nl-1), 0.4.^(0:ns-1))
Expand Down Expand Up @@ -184,4 +185,38 @@ end
@test isa(convert(AbstractMatrix{Complex128},T),Hankel{Complex128})
@test isa(convert(AbstractArray{Complex128},T),Hankel{Complex128})
@test isa(convert(ToeplitzMatrices.Hankel{Complex128},T),Hankel{Complex128})


@test Circulant(1:5) == Circulant(Vector(1.0:5))

end


A = ones(10, 10)
@test full(Toeplitz(A)) == full(Toeplitz{Float64}(A)) == A
@test full(SymmetricToeplitz(A)) == full(SymmetricToeplitz{Float64}(A)) == A
@test full(Circulant(A)) == full(Circulant{Float64}(A)) == A
@test full(Hankel(A)) == full(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 full(SymmetricToeplitz(A)) == full(SymmetricToeplitz{Float64}(A)) ==
full(Toeplitz(Symmetric(A))) == full(Symmetric(Toeplitz(A))) == [1. 2.; 2. 1.]
@test full(Circulant(A)) == [1 3; 3 1]

@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 full(Hankel(A)) == full(Hankel{Float64}(A)) == [1.0 3; 3 4]

# Constructors should be projections
@test Toeplitz(Toeplitz(A)) == Toeplitz(A)
@test SymmetricToeplitz(SymmetricToeplitz(A)) == SymmetricToeplitz(A)
@test Circulant(Circulant(A)) == Circulant(A)
@test TriangularToeplitz(TriangularToeplitz(A, :U), :U) == TriangularToeplitz(A, :U)
@test TriangularToeplitz(TriangularToeplitz(A, :L), :L) == TriangularToeplitz(A, :L)
@test Hankel(Hankel(A)) == Hankel(A)

0 comments on commit e97fdf4

Please sign in to comment.