Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Pull request - extended circulant functionality + DFT representation. #13

Closed
wants to merge 10 commits into from
Closed

Pull request - extended circulant functionality + DFT representation. #13

wants to merge 10 commits into from

Conversation

willtebbutt
Copy link

Hi.

I have extended the package in a couple of ways. Firstly, the Discrete Fourier Transform (DFT) matrix has been added in src/DFT.jl. This is interesting semantically as it allows for (in certain circumstances) a more natural way of expressing the Fourier Transform of vectors and matrices (where a matrix is just viewed as a collection of vectors to be transformed). Note that there are primarily two implementations of this matrix, Unitary DFT (UDFT) uses the matrix U defined such that ctranspose(U) = inv(U), whereas this does not hold for DFT.

This is particularly useful for cleaning up the implementation of the Circulant matrix, as we can now describe all of the operations one might expect to be able to perform with Circulant matrices in terms of multiplication by the DFT matrix or it's conjugate transpose (or some re-scaling thereof). Furthermore, I have added a significant amount of extra functionality involving Circulant matrices, the vast majority of which is tested in test/circulant.jl. (Note that I have separated out the functionality involving Toeplitz and Circulant matrices). Furthermore, I have changed the behaviour of the multiplication of real and integer Circulant matrices. Multiplication now always returns a complex-valued float, but helper functions have been added (see the end of src/circulant.jl) to extend the casting functionality from Base necessary to enable the user to explicitly round and cast. It is my feeling that this is the correct thing to do as it is unwise to mask type conversions that may result in some loss of accuracy from the user. Commenting in src/circulant.jl is fairly sparse as most code should be self-documenting. Let me know if more comments are wanted.

I appreciate that this is a fairly substantial modification / addition, so please review and let me know what you think!

Will Tebbutt added 5 commits January 12, 2017 16:47
Basic DFT functionality added. Coverage of the expected operations is…
9406a8f 
…, however, somewhat patchy.
Re-implmented existing Circulant multiplication functionality to corr…
ff5320a 
…ectly use dispatch to determine what to cast back to. Previous implementation used dynamic type checking and missed some edge cases. Not re-implemented the in-place multiplication functionality in this commit, nor have the tests been updated. Also added typealias for Circulant to reduce verbosity.
Completed + tested basic functionality for circulant matrices, includ…
20ed4b4 
…ing multiplication in terms of DFT matrices, (conjugate) transposition, addition, subtraction, eigenvalues + eigenvectors. dft.jl has been further extended to cover more functionality and a few bugs have been fixed.
Commit to bring remote up to speed.
f449a93
Circulant functionality sufficiently complete for a first vesion.
da86f31
@ivanslapnicar
Copy link
Contributor

ivanslapnicar commented Jan 20, 2017 via email

@willtebbutt
Stupid mistake in circulant
93b7673 
Removed two lines of redundant test code from circulant.jl.
@willtebbutt
Copy link
Author

Hi Ivan,

Thanks for the response. With regards to your concerns about the underlying Circulant implementation, I have been careful to ensure that (U)DFT matrices are never explicitly evaluated - as you point out doing so would require O(n^2) operations + memory which would indeed defeat the point of working with Circulant matrices to get the complexity down to O(n logn).

As evidence for this, I would first highlight lines 10-14 of src/dft.jl. The UDFT (Unitary DFT, eigenvectors of the Circulant matrix of a given size) stores only a single Int and Float64 (the latter being for convenience). Furthermore, there is no indexing or 'full' functionality defined for the (U)DFT matrix, so it is not possible for the JIT to unexpectedly compile code which explicitly evaluates the full (U)DFT matrix and default to a naive approach to performing multiplication operations. Furthermore, if you look at, for example, line 45 of src/dft.jl you will see that the left-multiplication of a StridedVecOrMat (SVM) by the DFT matrix is implemented in terms of the fft, as one would hope. The rest of the file is more of the same - it is all convenience functionality to mean that you can write things like

X * U

for some StridedVecOrMat X and (U)DFT U, and the correct (i)fft operation will be used to compute the requested operations efficiently under the hood.

This propagates through to the Circulant matrices in src/circulant.jl. First observe line 41; this line 'computes' the eigen-factorisation of a circulant matrix without actually computing the eigenvectors explicitly, by using the previously described UDFT matrix. This means that multiplication routine on line 64 works efficiently. I'll explain what is going on there: conf(U,X) checks that U and X are conformal. Line 66 simply gets pointers to the eigenvalues and stores them in a Diagonal matrix (which is implemented efficiently under the hood to do the things that you would expect), and get a pointer to the UDFT matrix which represents the eigenvectors. Consequently, line 67 is computed efficiently using the multiplication function on line 47 of src/dft.jl. Γ * Y on line 68 is computed efficiently by exploiting the efficient representation of Γ and the Ac_mul_B! is computed efficiently using code from line 77 of src/dft.jl.

Note that the discussed multiplication function accepts the eigen-factorisation of some Circulant matrix. This is done so that if one is, for example, left-multiplying multiple vectors/matrices by a Circulant matrix, the initial fft required to obtain the eigenvalues can be done once and cached (this is relevant to some of the work that I am doing). However, line 76 of src/circulant.jl provides a wrapper which computes the eigen-factorisation and passes it on such that a 'raw' Circulant matrix can be provided if the described caching behaviour is not required.

The rest of the functionality in src/circulant.jl is more of the same - care has been taken to ensure that all operations have at most complexity O(n logn).

The only notable changes to your API are that the indexing on line 18 of src/circulant.jl has been commented out - depending upon your use case it might make sense for this not to be the case.

Best,
Will

iamed2
iamed2 suggested changes Jan 27, 2017
View reviewed changes
src/circulant.jl

function full{T}(C::Circ{T})
n = size(C, 1)
M = Array(T, n, n)
Copy link

Choose a reason for hiding this comment

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

This syntax is deprecated and Array{T}(n, n) is now preferred

Copy link

Choose a reason for hiding this comment

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

In general when creating an array with a known number of dimensions you should call a specific constructor (in this case Array{T, 2}(n, n) or Matrix{T}(n, n).

test/runtests.jl
include("toeplitz.jl") #Toeplitz matrix
include("vandermonde.jl") #Vandermonde matrix
include("vandermonde.jl") #Vandermonde matrix
Copy link

Choose a reason for hiding this comment

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

Add back newlines at the end of your files

src/circulant.jl
M
end

# getindex(C::Circ, i::Int, j::Int) = C.c[mod(i-j,length(C.c))+1]
Copy link

Choose a reason for hiding this comment

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

Since show relies on getindex, you'll have to override it for Circulant so it doesn't error when trying to display Circulants

src/circulant.jl
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))
Copy link

Choose a reason for hiding this comment

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

With the two round methods you should pass-through trailing args... to allow the other options to round to be used.

src/dft.jl
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
Copy link

Choose a reason for hiding this comment

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

Since the DFT types are immutable and the contents are immutable, you can just return the item itself, e.g. transpose(U::AbstractDFT) = U and transpose!(U::AbstractDFT) = U

src/circulant.jl
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}}
Copy link

Choose a reason for hiding this comment

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

Calling show on an instance of this type causes a StackOverflowError.

show(STDOUT, eigfact(Circulant(collect(1.0:10.0)))); to replicate

src/circulant.jl
# 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)...)
Copy link

Choose a reason for hiding this comment

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

This function should convert/promote to the right types such that the answer is a CircEig. Currently that is not the case for Circ{T} where T <: Integer. That should either be fixed or disallowed.

@willtebbutt willtebbutt mentioned this pull request Aug 24, 2017
Closed
@willtebbutt willtebbutt closed this Jun 6, 2019
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers

@iamed2 iamed2 iamed2 requested changes

Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.
3 participants
@willtebbutt @ivanslapnicar @iamed2