Skip to content

Commit

Permalink
relax mul! type constraints (#50)
Browse files Browse the repository at this point in the history 
* relax mul! type constraints

hope this solves #40

* Fix realness test

hopefully..

* Bump version

* Bump StatsBase version

* mul! work

* mul!

* allreal

* Handle mixed tupes

* Test mixed types

* cut toreal

* Remove toreal; add docstrings, AbstractVectors

* Add codecov
  • Loading branch information
@JeffFessler
JeffFessler authored Aug 14, 2020
1 parent c2b06c6 commit 3c2ae6e
Show file tree
Hide file tree
Showing 5 changed files with 106 additions and 54 deletions.
1 change: 1 addition & 0 deletions .travis.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,4 +9,5 @@ julia:
notifications:
email: false
after_success:
- julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(process_folder())'
- julia -e 'cd(Pkg.dir("ToeplitzMatrices")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
4 changes: 2 additions & 2 deletions Project.toml
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
name = "ToeplitzMatrices"
uuid = "c751599d-da0a-543b-9d20-d0a503d91d24"
version = "0.6.1"
version = "0.6.2"

[deps]
AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
Expand All @@ -11,7 +11,7 @@ StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
[compat]
AbstractFFTs = "0.4, 0.5"
FFTW = "1"
StatsBase = "0.32"
StatsBase = "0.32, 0.33"
julia = "1"

[extras]
Expand Down
8 changes: 5 additions & 3 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,9 +2,11 @@ ToeplitzMatrices.jl
===========

[![Build Status](https://travis-ci.org/JuliaMatrices/ToeplitzMatrices.jl.svg?branch=master)](https://travis-ci.org/JuliaMatrices/ToeplitzMatrices.jl)
[![codecov](https://codecov.io/gh/JuliaMatrices/ToeplitzMatrices.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaMatrices/ToeplitzMatrices.jl)
[![Coverage Status](https://coveralls.io/repos/github/JuliaMatrices/ToeplitzMatrices.jl/badge.svg?branch=master&bust=1)](https://coveralls.io/github/JuliaMatrices/ToeplitzMatrices.jl?branch=master)

Fast matrix multiplication and division for Toeplitz and Hankel matrices in Julia
Fast matrix multiplication and division
for Toeplitz, Hankel and circulant matrices in Julia


## ToeplitzMatrix
Expand All @@ -18,7 +20,7 @@ Toeplitz(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.,2.,3.],[1.,4.,5.])
Toeplitz(1:3, [1.,4.,5.])
```

is a sparse representation of the matrix
Expand Down Expand Up @@ -62,7 +64,7 @@ where uplo is either `:L` or `:U` and `ve` are the rows or columns, respectively
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.

```julia
Hankel([1.,2.,3.],[3.,4.,5.])
Hankel([1.,2.,3.], 3:5)
```

is a sparse representation of the matrix
Expand Down
74 changes: 53 additions & 21 deletions src/ToeplitzMatrices.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,8 @@ end
convert(::Type{Matrix}, A::AbstractToeplitz) = Matrix(A)

# Fast application of a general Toeplitz matrix to a column vector via FFT
function mul!(y::StridedVector{T}, A::AbstractToeplitz{T}, x::StridedVector, α::T, β::T) where T
function mul!(y::StridedVector, A::AbstractToeplitz, x::StridedVector, α::Number, β::Number)
T = promote_type(eltype.([y, A, x, α, β])...)
m = size(A,1)
n = size(A,2)
N = length(A.vcvr_dft)
Expand Down Expand Up @@ -86,18 +87,14 @@ function mul!(y::StridedVector{T}, A::AbstractToeplitz{T}, x::StridedVector, α:
end
A.dft \ A.tmp
for i in 1:m
y[i] += α * (T <: Real ? real(A.tmp[i]) : A.tmp[i])
y[i] += α * ((T <: Real) ? real(A.tmp[i]) : A.tmp[i])
end
return y
end
end
# Avoid ambiguity error
mul!(y::StridedVector{T}, A::AbstractToeplitz{T}, x::StridedVector, α::T, β::T) where {T<:Number} =
invoke(mul!, Tuple{StridedVector{T},ToeplitzMatrices.AbstractToeplitz{T},StridedVector,T,T} where T,
y, A, x, α, β)

# Application of a general Toeplitz matrix to a general matrix
function mul!(C::StridedMatrix{T}, A::AbstractToeplitz{T}, B::StridedMatrix, α::T, β::T) where T
function mul!(C::StridedMatrix, A::AbstractToeplitz, 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"))
Expand All @@ -107,10 +104,6 @@ function mul!(C::StridedMatrix{T}, A::AbstractToeplitz{T}, B::StridedMatrix, α:
end
return C
end
# Avoid ambiguity error
mul!(C::StridedMatrix{T}, A::AbstractToeplitz{T}, B::StridedMatrix, α::T, β::T) where {T<:Number} =
invoke(mul!, Tuple{StridedMatrix{T},ToeplitzMatrices.AbstractToeplitz{T},StridedMatrix,T,T} where T,
C, A, B, α, β)

# Translate three to five argument mul!
mul!(y::StridedVecOrMat, A::AbstractToeplitz, x::StridedVecOrMat) =
Expand Down Expand Up @@ -138,6 +131,13 @@ function (\)(A::AbstractToeplitz, b::AbstractVector)
end

# General Toeplitz matrix
"""
Toeplitz{T<:Number, S<:Number}
Subtype of `AbstractToeplitz{T}` for representing a Toeplitz matrix
where the matrix elements have `eltype T`
and the underlying DFT has `eltype S`.
"""
mutable struct Toeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
vc::Vector{T}
vr::Vector{T}
Expand All @@ -147,7 +147,13 @@ mutable struct Toeplitz{T<:Number,S<:Number} <: AbstractToeplitz{T}
end

# Ctor
function Toeplitz{T}(vc::Vector, vr::Vector) where {T}
"""
T = Toeplitz(vc::AbstractVector, vr::AbstractVector)
Create a `Toeplitz` matrix `T` from its first column `vc` and first row `vr`
where `vc[1] == vr[1]`.
"""
function Toeplitz{T}(vc::AbstractVector, vr::AbstractVector) where {T}
m, n = length(vc), length(vr)
if vc[1] != vr[1]
error("First element of the vectors must be the same")
Expand All @@ -164,10 +170,14 @@ function Toeplitz{T}(vc::Vector, vr::Vector) where {T}
return Toeplitz(vcp, vrp, dft*tmp, similar(tmp), dft)
end

Toeplitz(vc::Vector, vr::Vector) =
Toeplitz(vc::AbstractVector, vr::AbstractVector) =
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)
"Project" matrix `A` onto its Toeplitz part using the first row/col of `A`.
"""
Toeplitz(A::AbstractMatrix) = Toeplitz(A[:,1], A[1,:])
Toeplitz{T}(A::AbstractMatrix) where T = Toeplitz{T}(A[:,1], A[1,:])

Expand Down Expand Up @@ -222,7 +232,7 @@ function tril(A::Toeplitz, k = 0)
return Al
end

# Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diaganal
# Form a lower triangular Toeplitz matrix by annihilating all entries below the k-th diagonal
function triu(A::Toeplitz, k = 0)
if k < 0
error("Second argument cannot be negative")
Expand All @@ -247,10 +257,10 @@ mutable struct SymmetricToeplitz{T<:BlasReal} <: AbstractToeplitz{T}
dft::Plan
end

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


Expand Down Expand Up @@ -285,6 +295,13 @@ ldiv!(A::SymmetricToeplitz, b::StridedVector) =
copyto!(b, IterativeLinearSolvers.cg(A, zeros(length(b)), b, strang(A), 1000, 100eps())[1])

# Circulant
"""
Circulant{T<:Number, S<:Number}
Subtype of `AbstractToeplitz{T}` for representing a circulant matrix
where the matrix elements have `eltype T`
and the underlying DFT has `eltype S`.
"""
mutable struct Circulant{T<:Number,S<:Number} <: AbstractToeplitz{T}
vc::Vector{T}
vcvr_dft::Vector{S}
Expand All @@ -298,7 +315,22 @@ function Circulant{T}(vc::Vector{T}) where T
end

Circulant{T}(vc::AbstractVector) where T = Circulant{T}(convert(Vector{T}, vc))
Circulant(vc::AbstractVector) = Circulant{promote_type(eltype(vc), Float32)}(vc)

"""
C = Circulant(vc::AbstractVector)
Create a circulant matrix `C` from its first column `vc`.
"""
function Circulant(vc::AbstractVector{T}) where T
V = promote_type(eltype(vc), Float32)
return Circulant{V}(convert(Vector{V}, vc))
end

"""
C = Circulant(A::AbstractMatrix)
Create a circulant matrix `C` from the first column of matrix `A`.
"""
Circulant{T}(A::AbstractMatrix) where T = Circulant{T}(A[:,1])
Circulant(A::AbstractMatrix) = Circulant(A[:,1])

Expand Down Expand Up @@ -686,7 +718,7 @@ getindex(A::Hankel, i::Integer, j::Integer) = A.T[i,end-j+1]
*(A::Hankel, B::AbstractMatrix) = A.T * flipdim(B, 1)
## BigFloat support

(*)(A::Toeplitz{T}, b::Vector) where {T<:BigFloat} = irfft(
(*)(A::Toeplitz{T}, b::AbstractVector) where {T<:BigFloat} = irfft(
rfft([
A.vc;
reverse(A.vr[2:end])]
Expand Down
73 changes: 45 additions & 28 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,33 +8,39 @@ nl = 2000

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

vc = LinRange(1,3,3) # for testing with AbstractVector
vv = Vector(vc)
vr = [1, 5.]

cases = [
(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"),
]

for (As, Al, st) in cases
@testset "Toeplitz: $st" begin
@test As * xs ≈ Matrix(As) * xs
@test As'* xs ≈ Matrix(As)' * xs
@test Al * xl ≈ Matrix(Al) * xl
Expand All @@ -46,7 +52,16 @@ xl = randn(nl, 5)
@test Matrix(As') == Matrix(As)'
@test Matrix(transpose(As)) == transpose(Matrix(As))
end
)
end

@testset "Mixed types" begin
@test eltype(Toeplitz([1, 2], [1, 2])) == Float32 # !!
@test Toeplitz([1, 2], [1, 2]) * ones(2) == fill(3, 2)
@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))
end

@testset "Real general rectangular" begin
Ar1 = Toeplitz(0.9.^(0:nl-1), 0.4.^(0:ns-1))
Expand Down Expand Up @@ -74,6 +89,7 @@ end
@test ldiv!(Al, copy(xl)) ≈ Matrix(Al) \ xl
@test StatsBase.levinson(As, xs) ≈ Matrix(As) \ xs
@test StatsBase.levinson(Ab, xs) ≈ Matrix(Ab) \ xs
@test Matrix(SymmetricToeplitz(vc)) == Matrix(SymmetricToeplitz(vv))
end

@testset "Hankel" begin
Expand Down Expand Up @@ -101,6 +117,7 @@ end
@test Hs * xs[:,1] ≈ Matrix(Hs) * xs[:,1]
@test Hs * xs ≈ Matrix(Hs) * xs
@test Hl * xl ≈ Matrix(Hl) * xl
@test Matrix(Hankel(reverse(vc),vr)) == Matrix(Hankel(reverse(vv),vr))
end

@testset "Complex square" begin
Expand Down

2 comments on commit 3c2ae6e

@dlfivefifty
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/19551

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.6.2 -m "<description of version>" 3c2ae6e71d042ebedb403d294c2b8028d23e9ff4
git push origin v0.6.2

Please sign in to comment.