ToeplitzMatrices.jl

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

Multiplication of large matrices and sqrt, inv, LinearAlgebra.eigvals, LinearAlgebra.ldiv!, and LinearAlgebra.pinv for circulant matrices are computed with FFTs. To be able to use these methods, you have to install and load a package that implements the AbstractFFTs.jl interface such as FFTW.jl:

using FFTW

If you perform multiple calculations with FFTs, it can be more efficient to 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.

A Toeplitz matrix has constant diagonals. It can be constructed using

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.

Toeplitz(1:3, [1.,4.,5.])

is a sparse representation of the matrix

[ 1.0  4.0  5.0
  2.0  1.0  4.0
  3.0  2.0  1.0 ]

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}.

A Hankel matrix has constant anti-diagonals, for example,

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

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

The reverse can transform between Hankel and Toeplitz. It is used to achieve fast Hankel algorithms.

Note that scalar multiplication, conj, + and - could be removed once broadcast is implemented.

reverse(Hankel) returns a Toeplitz, while reverse(AbstractToeplitz) returns a Hankel.

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,

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.

Symmetric, UpperTriangular and LowerTriangular from LinearAlgebra are also overloaded for convenience.

Symmetric(T::Toeplitz) = SymmetricToeplitz(T)
UpperTriangular(T::Toeplitz) = UpperTriangularToeplitz(T)
LowerTriangular(T::Toeplitz) = LowerTriangularToeplitz(T)

TriangularToeplitz is reserved for backward compatibility.

TriangularToeplitz = Union{UpperTriangularToeplitz,LowerTriangularToeplitz}

The old interface is implemented by

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.

Methods in this section are not exported.

_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.

_circulate(v::AbstractVector) converts between the vr and vc of a Circulant.

isconcrete(A::Union{AbstractToeplitz,Hankel}) decides whether the stored vector(s) are concrete. It calls Base.isconcretetype.