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
.
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]
, whileSymmetricToeplitz(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
.