Pull request - extended circulant functionality + DFT representation.

src/SpecialMatrices.jl

@@ -1,16 +1,26 @@
 module SpecialMatrices
 
-import Base: A_mul_B!, full, getindex, inv, isassigned, size, *
+import Base: A_mul_B!, full, getindex, inv, isassigned, size, *, length, +, -,
+    Ac_mul_B, A_mul_Bc, At_mul_B, A_mul_Bt, eig, eigfact, Eigen, \, /,
+    transpose, ctranspose, copy, conj, conj!, inv, det, logdet, real, convert,
+    round
 
-include("cauchy.jl") #Cauchy matrix
-include("companion.jl") #Companion matrix
-include("frobenius.jl") #Frobenius matrix
-include("hankel.jl") #Hankel matrix
-include("hilbert.jl") #Hilbert matrix
-include("kahan.jl") #Kahan matrix
-include("riemann.jl") #Riemann matrix
-include("strang.jl") #Strang matrix
-include("toeplitz.jl") #Toeplitz matrix, Circulant matrix
-include("vandermonde.jl") #Vandermonde matrix
+typealias SV StridedVector
+typealias SM StridedMatrix
+typealias SVM StridedVecOrMat
+
+include("dft.jl") # Discrete Fourier Transform matrices.
+
+include("cauchy.jl") # Cauchy matrix
+include("circulant.jl") # Circulant matrix.
+include("companion.jl") # Companion matrix
+include("frobenius.jl") # Frobenius matrix
+include("hankel.jl") # Hankel matrix
+include("hilbert.jl") # Hilbert matrix
+include("kahan.jl") # Kahan matrix
+include("riemann.jl") # Riemann matrix
+include("strang.jl") # Strang matrix
+include("toeplitz.jl") # Toeplitz matrix
+include("vandermonde.jl") # Vandermonde matrix
 
 end # module

src/circulant.jl

@@ -0,0 +1,122 @@
+export Circulant, CircEig, tocirc
+
+immutable Circulant{T} <: AbstractMatrix{T}
+    c::Vector{T}
+end
+typealias Circ Circulant
+
+function full{T}(C::Circ{T})
+    n = size(C, 1)
+    M = Array(T, n, n) 
+    for i=1:n
+        M[i:n,i] = C.c[1:n-i+1]
+        M[1:i-1,i] = C.c[n-i+2:n]
+    end
+    M
+end
+
+# getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
+isassigned(C::Circ, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
+size(C::Circ, r::Int) = (r==1 || r==2) ? length(C.c) :
+    throw(ArgumentError("Invalid dimension $r"))
+size(C::Circ) = size(C,1), size(C,2)
+
+copy(C::Circ) = Circ(C.c)
+conj(C::Circ) = Circ(conj(C.c))
+conj!(C::Circ) = (conj!(C.c); C)
+
+transpose(C::Circ) = Circ(circshift(reverse(C.c), 1))
+ctranspose(C::Circ) = conj!(transpose(C))
+
++(C::Circ, a::Number) = Circ(C.c + a)
++(a::Number, C::Circ) = Circ(a + C.c)
++(C::Circ, D::Circ) = Circ(C.c + D.c)
+-(C::Circ, a::Number) = Circ(C.c - a)
+-(a::Number, C::Circ) = Circ(a - C.c)
+-(C::Circ, D::Circ) = Circ(C.c - D.c)
+-(C::Circ) = Circ(-C.c)
+
+# Return the eigen-factorisation of C. Constituents are efficient to compute
+# and the eigenvectors are represented efficiently using a UnitaryDFT object.
+eig{T<:Real}(C::Circ{T}) = (N = length(C.c); (DFT{Complex{T}}(N) * C.c, UDFT{Complex{T}}(N)))
+eig{T<:Complex}(C::Circ{T}) = (N = length(C.c); (DFT{T}(N) * C.c, UDFT{T}(N)))
+eigfact(C::Circ) = Eigen(eig(C)...)
+typealias CircEig{T<:Complex} Eigen{T, T, UDFT{T}, Vector{T}}
+
+# Utility function to convert back to usual circulant formulation.
+tocirc{T}(C::CircEig{T}) = (N = length(C.values); Circ(DFT{T}(N)'C.values / N))
+
+# Check for conformal arguments.
+function conf{T}(C::CircEig{T}, X::SVM)
+    if size(X, 1) != size(C.vectors, 2)
+        throw(ArgumentError("C and X are not conformal."))
+    end
+end
+function conf{T}(X::SVM, C::CircEig{T})
+    if size(X, 2) != size(C.vectors, 1)
+        throw(ArgumentError("X and C are not conformal."))
+    end
+end
+
+# Always compute the eigen-factorisation to compute the multiplication. If
+# this operation is to be performed multiple times, then we should just cache
+# the factorisation to save re-computing it for each multiplication operation.
+function *{T}(C::CircEig{T}, X::SVM)
+    conf(C, X)
+    Γ, U = Diagonal(C.values), C.vectors
+    Y = U * X
+    Ac_mul_B!(U, Γ * Y)
+end
+function *{T}(X::SVM, C::CircEig{T})
+    conf(X, C)
+    Γ, U = Diagonal(C.values), C.vectors
+    Y = A_mul_Bc(X, U)
+    A_mul_B!(A_mul_B!(Y, Y, Γ), U)
+end
+*(C::Circ, X::SVM) = eigfact(C) * X
+*(X::SVM, C::Circ) = X * eigfact(C)
+
+# "In-place" multiplication.
+function A_mul_B!{T<:Number, V<:Complex}(Y::SVM{V}, C::CircEig{T}, X::SVM{V})
+    conf(C, X) 
+    Γ, U = Diagonal(C.values), C.vectors
+    Ac_mul_B!(U, A_mul_B!(Y, Γ, A_mul_B!(U, X)))
+end
+function A_mul_B!{T<:Number, V<:Complex}(Y::SVM{V}, X::SVM{V}, C::CircEig{T})
+    conf(X, C)
+    Γ, U = Diagonal(C.values), C.vectors
+    A_mul_B!(A_mul_B!(Y, A_mul_Bc!(X, U), Γ), U)
+end
+A_mul_B!(Y::SVM, C::Circ, X::SVM) = A_mul_B!(Y, eigfact(C), X)
+A_mul_B!(Y::SVM, X::SVM, C::Circ)  = A_mul_B!(Y, X, eigfact(C))
+
+inv{T}(C::CircEig{T}) = Eigen(1 ./ C.values, UDFT{T}(length(C.values)))
+inv(C::Circ) = tocirc(inv(eigfact(C)))
+
+det{T}(C::CircEig{T}) = prod(C.values)
+det(C::Circ) = det(eigfact(C))
+
+logdet{T}(C::CircEig{T}) = sum(log(C.values))
+logdet(C::Circ) = logdet(eigfact(C))
+
+# Use the eigen-factorisation to compute the inv(C) * X. The inverse
+# factorisation should be cached if the operation is required multiple times.
+function \{T}(C::CircEig{T}, X::SVM)
+    conf(C, X)
+    Γ, U = Diagonal(C.values), C.vectors
+    Ac_mul_B(U, (Γ \ (U * X)))
+end
+function /{T}(X::SVM, C::CircEig{T})
+    conf(X, C)
+    Γ, U = Diagonal(C.values), C.vectors
+    (A_mul_Bc(X, U) / Γ) * U
+end
+\(C::Circ, X::SVM) = eigfact(C) \ X
+/(X::SVM, C::Circ) = X / eigfact(C)
+
+# Helper functionality for casting.
+real(C::Circ) = Circ(real(C.c))
+convert{T}(::Type{Circ{T}}, C::Circ) = Circ(convert(Vector{T}, C.c))
+
+round(C::Circ) = Circ(round(C.c))
+round{T}(::Type{Circ{T}}, C::Circ) = Circ(round(Vector{T}, C.c))

src/dft.jl

@@ -0,0 +1,86 @@
+export DFT, UnitaryDFT, UDFT
+
+# Non-unitary DFT matrix.
+abstract AbstractDFT{T} <: AbstractMatrix{T}
+immutable DFT{T<:Complex} <: AbstractDFT{T}
+    N::Int
+end
+
+# Unitary DFT Matrix.
+immutable UnitaryDFT{T<:Complex} <: AbstractDFT{T}
+    N::Int
+    sqrtN::Float64
+    UnitaryDFT(N) = new(N, sqrt(N))
+end
+typealias UDFT UnitaryDFT
+
+Base.show(io::IO, U::DFT) = print(io, "DFT matrix of size $U.N.")
+Base.show(io::IO, U::UDFT) = print(io, "Unitary DFT matrix of size $U.N.")
+
+size(M::AbstractDFT) = M.N, M.N
+length(M::AbstractDFT) = M.N^2
+
+# Don't bother to implement an in-place transpose as U is sufficiently
+# light-weight not to care in the slightest.
+transpose(U::AbstractDFT) = copy(U)
+
+# Indexing is explicitely disabled as accidentally using the naive version of
+# the DFT is a really bad idea. Basically any functionality that would allow
+# you to do so has been disabled.
+# getindex(M::DFT, m::Int, n::Int) =
+#     exp(-2π * im * (m-1) * (n-1) / M.N)
+# getindex(M::UnitaryDFT, m::Int, n::Int) =1
+#     exp(-2π * im * (m-1) * (n-1) / M.N) / M.sqrtN
+# getindex(M::AbstractDFT, m::Int) = getindex(M, rem(m, M.N), div(m, M.N))
+
+# isassigned(M::AbstractDFT, i::Int) = 1 <= i <= M.N ? true : false
+
+# Helper functions to check whether stuff is conformal or not.
+conf(U::AbstractDFT, x::SV) = assert(U.N == length(x))
+conf(x::SV, U::AbstractDFT) = assert(false)
+conf(U::AbstractDFT, X::SM) = assert(U.N == size(X, 1))
+conf(X::SM, U::AbstractDFT) = assert(size(X, 2) == U.N)
+
+# Multiplication from the left and right by the DFT matrices.
+*(U::DFT, X::SVM) = (conf(U, X); fft(X, 1))
+*(X::SVM, U::DFT) = (conf(X, U); fft(X, 2))
+*(U::UDFT, X::SVM) = (conf(U, X); fft(X, 1) / U.sqrtN)
+*(X::SVM, U::UDFT) = (conf(X, U); fft(X, 2) / U.sqrtN)
+
+# In-place functionality matrix multiplication.
+A_mul_B!(Y::SV, U::DFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y))
+A_mul_B!(Y::SV, X::SV, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y))
+A_mul_B!(Y::SV, U::UDFT, X::SV) = (conf(U, X); copy!(Y, X); fft!(Y); Y ./= U.sqrtN)
+A_mul_B!(Y::SV, X::SV, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y); Y ./= U.sqrtN)
+A_mul_B!(Y::SM, U::DFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1))
+A_mul_B!(Y::SM, X::SM, U::DFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2))
+A_mul_B!(Y::SM, U::UDFT, X::SM) = (conf(U, X); copy!(Y, X); fft!(Y, 1); Y ./= U.sqrtN)
+A_mul_B!(Y::SM, X::SM, U::UDFT) = (conf(X, U); copy!(Y, X); fft!(Y, 2); Y ./= U.sqrtN)
+A_mul_B!(U::DFT, X::SVM) = (conf(U, X); fft!(X, 1))
+A_mul_B!(X::SVM, U::DFT) = (conf(X, U); fft!(X, 2))
+A_mul_B!(U::UDFT, X::SVM) = (conf(U, X); fft!(X, 1); X ./= U.sqrtN)
+A_mul_B!(X::SVM, U::UDFT) = (conf(X, U); fft!(X, 2); X ./= U.sqrtN)
+
+# Conjugate transpose U'X and X*U'.
+Ac_mul_B(U::DFT, X::SVM) = (conf(U, X); bfft(X, 1))
+A_mul_Bc(X::SVM, U::DFT) = (conf(X, U); bfft(X, 2))
+Ac_mul_B(U::UDFT, X::SVM) = (conf(U, X); bfft(X, 1) / U.sqrtN)
+A_mul_Bc(X::SVM, U::UDFT) = (conf(X, U); bfft(X, 2) / U.sqrtN)
+
+# Conjugate transpose operations in-place.
+Ac_mul_B!{T<:SVM}(Y::T, U::DFT, X::T) = (conf(U, X); copy!(Y, X); bfft!(Y, 1))
+A_mul_Bc!{T<:SVM}(Y::T, X::T, U::DFT) = (conf(X, U); copy!(Y, X); bfft!(Y, 2))
+Ac_mul_B!{T<:SVM}(Y::T, U::UDFT, X::T) = (conf(U, X); copy!(Y, X); bfft!(Y, 1); Y ./= U.sqrtN)
+A_mul_Bc!{T<:SVM}(Y::T, X::T, U::UDFT) = (conf(X, U); copy!(Y, X); bfft!(Y, 2); Y ./= U.sqrtN)
+Ac_mul_B!(U::DFT, X::SVM) = (conf(U, X); bfft!(X, 1))
+A_mul_Bc!(X::SVM, U::DFT) = (conf(X, U); bfft!(X, 2))
+Ac_mul_B!(U::UDFT, X::SVM) = (conf(U, X); bfft!(X, 1); X ./= U.sqrtN)
+A_mul_Bc!(X::SVM, U::UDFT) = (conf(X, U); bfft!(X, 2); X ./= U.sqrtN)
+
+# Explicit implementation of X'U and U*X' to avoid fallback in Base being called.
+Ac_mul_B(X::SVM, U::AbstractDFT) = (Xc = ctranspose(X); conf(Xc, U); A_mul_B!(Xc, U))
+A_mul_Bc(U::AbstractDFT, X::SVM) = (Xc = ctranspose(X); conf(U, Xc); A_mul_B!(U, Xc))
+
+# All DFT matrices are symmetric.
+At_mul_B(U::AbstractDFT, X::SVM) = U * X
+A_mul_Bt(X::SVM, U::AbstractDFT) = X * U

src/hankel.jl

@@ -3,7 +3,8 @@ export Hankel
 immutable Hankel{T} <: AbstractArray{T, 2}
     c :: Vector{T}
 end
-Hankel{T}(c::Vector{T}) = length(c) % 2 == 1 ? Hankel{T}(c) : throw(ArgumentError(""))
+(::Hankel{T}){T}(c::Vector{T}) = length(c) % 2 == 1 ? Hankel{T}(c) : throw(ArgumentError(""))
+
 
 getindex(H::Hankel, i::Int, j::Int) = H.c[i+j-1]
 isassigned(H::Hankel, i::Int, j::Int) = isassigned(H.c, i+j-1)

src/toeplitz.jl

@@ -1,4 +1,4 @@
-export Circulant, Toeplitz
+export Toeplitz
 
 immutable Toeplitz{T} <: AbstractArray{T, 2}
 	c :: Vector{T}
@@ -38,48 +38,3 @@ function full{T}(To::Toeplitz{T})
 	end
 	M
 end
-
-
-
-immutable Circulant{T} <: AbstractArray{T, 2}
-	c :: Vector{T}
-end
-
-getindex(C::Circulant, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
-isassigned(C::Circulant, i::Int, j::Int) = isassigned(C.c,mod(i-j,length(C.c))+1)
-size(C::Circulant, r::Int) = (r==1 || r==2) ? length(C.c) :
-    throw(ArgumentError("Invalid dimension $r"))
-size(C::Circulant) = size(C,1), size(C,2)
-
-# Fast matrix x vector via fft()
-# see Golub, van Loan, Matrix Computations, John Hopkins, Baltimore, 1996, p. 202 
-function *{T}(C::Circulant{T},x::Vector{T})
-    xt=fft(x)
-    vt=fft(C.c)
-    yt=vt.*xt
-    typeof(x[1])==Int ? map(Int,round(real(ifft(yt)))): ( (T <: Real) ? map(T,real(ifft(yt))) : ifft(yt))
-end
-
-function A_mul_B!{T}(y::StridedVector{T},C::Circulant{T},x::StridedVector{T})
-    xt=fft(x)
-    vt=fft(C.c)
-    yt=ifft(vt.*xt)
-    if T<: Int
-        map!(round,y,yt) 
-    elseif T<: Real
-        map!(real,y,yt)
-    else
-        copy!(y,yt)
-    end
-    return y
-end
-
-function full{T}(C::Circulant{T})
-	n=size(C, 1)
-	M=Array(T, n, n)
-	for i=1:n
-		M[i:n,i] = C.c[1:n-i+1]
-		M[1:i-1,i] = C.c[n-i+2:n]
-	end
-	M
-end