ToeplitzMatrices.jl

Fast matrix multiplication and division for Toeplitz matrices in Julia
Popularity
59 Stars
Updated Last
12 Months Ago
Started In
September 2013

ToeplitzMatrices.jl

Build Status Coverage Coverage

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

Note

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.

Introduction

Toeplitz

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 ]

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,

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

Implemented interface

Operations

  • ✓ implemented
  • ✗ error
  • _ fall back to Matrix
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!

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

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

LinearAlgebra

Constructors and conversions

Toeplitz Symmetric~ Circulant UpperTriangular~ LowerTriangular~ Hankel
from AbstractVector
from AbstractMatrix
from AbstractToeplitz
to supertype
to Toeplitz -
to another eltype

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

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.

Unexported interface

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.