Hello everyone,
I want to calculate the explicit DFT matrix (and its inverse). I have been looking at related Julia packages to see if this functionality is somewhere available. While I could calculate it myself, I am always wary of implementing this kind of numerical tools myself (sensitive numerics, performance), I was hoping to find a “standard” implementation somewhere. Is somebody perhaps aware of a package that returns the DFT matrix? thanks.
You could calculate it by taking the FFT of unit vectors, using a trustworthy algorithm like FFTW.
2 Likes
For the benefit of the reader, DFT probably stands for “Discrete Fourier Transform” (not e.g. “Density Functional Theory” 
)
3 Likes
I’ve updated the title to make this clear.
3 Likes
Something easy which probably does the job:
julia> function DFT_matrix(N, norm=1 / √(N))
ω = exp(-2π * 1im / N)
W = [norm * ω^((i-1) * (j-1)) for i ∈ 1:N, j ∈ 1:N]
end
DFT_matrix (generic function with 2 methods)
julia> DFT_matrix(4, 1)
4×4 Matrix{ComplexF64}:
1.0+0.0im 1.0+0.0im 1.0+0.0im 1.0+0.0im
1.0+0.0im 6.12323e-17-1.0im -1.0-1.22465e-16im -1.83697e-16+1.0im
1.0+0.0im -1.0-1.22465e-16im 1.0+2.44929e-16im -1.0-3.67394e-16im
1.0+0.0im -1.83697e-16+1.0im -1.0-3.67394e-16im 5.51091e-16-1.0im
julia> return DFT_matrix(2, 1)
2×2 Matrix{ComplexF64}:
1.0+0.0im 1.0+0.0im
1.0+0.0im -1.0-1.22465e-16im
1 Like
That could even be done slightly easier as W = [norm * ω^(i+j) for i ∈ 0:(N-1), j ∈ 0:(N-1)] (typo: I of course meant i*j).
(or stick to the original range and write . (typo, see below)(i+j-2))
Are you sure?
It’s i * j for i, j = 0, 1, ..., N-1.
Also of course (i + j - 2) can’t be derived.
Oh – yeah that’s a typo, sorry its Friday after and exhausting week. Of course I meant *.
2 Likes
Surprisingly (?), the FFT approach suggested by @John_Gibson is faster than building the matrix in a comprehension for larger N:
using LinearAlgebra
using FFTW
DFT_mat_fft(N) = stack(fft.(eachcol(I(N))))
julia> DFT_matrix(10, 1) ≈ DFT_mat_fft(10)
true
julia> DFT_matrix(2_500, 1) ≈ DFT_mat_fft(2_500)
true
julia> @btime DFT_matrix(2_500, 1);
785.036 ms (2 allocations: 95.37 MiB)
julia> @btime DFT_mat_fft(2_500);
160.114 ms (20005 allocations: 287.08 MiB)
2 Likes
Even more simply, fft(I(n), dims=1)fft(I(n), 1)
What’s slowing down the comprehension approach is recomputing the exponential for every element. There are various ways to speed it up, but I suspect that brevity is the priority here.
Of course, computing the matrix explicitly, as opposed to using an FFT as a linear operator, is usually unnecessary. I’m curious to know what the OP’s application is.
6 Likes
The dims keyword is not used for AbstractFFTs.jl
julia> fft(I(10), dims=1)
ERROR: MethodError: no method matching fft(::Diagonal{Bool, Vector{Bool}}; dims::Int64)
Closest candidates are:
fft(::AbstractArray{<:Real}, ::Any) got unsupported keyword argument "dims"
@ AbstractFFTs ~/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:214
fft(::AbstractArray, ::Any) got unsupported keyword argument "dims"
@ AbstractFFTs ~/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:67
fft(::AbstractArray) got unsupported keyword argument "dims"
@ AbstractFFTs ~/.julia/packages/AbstractFFTs/4iQz5/src/definitions.jl:66
...
Stacktrace:
[1] top-level scope
@ REPL[20]:1
julia> fft(I(10), (1,))
10×10 Matrix{ComplexF64}:
1.0+0.0im 1.0+0.0im 1.0+0.0im … 1.0+0.0im
1.0+0.0im 0.809017-0.587785im 0.309017-0.951057im 0.809017+0.587785im
1.0+0.0im 0.309017-0.951057im -0.809017-0.587785im 0.309017+0.951057im
1.0+0.0im -0.309017-0.951057im -0.809017+0.587785im -0.309017+0.951057im
1.0+0.0im -0.809017-0.587785im 0.309017+0.951057im -0.809017+0.587785im
1.0+0.0im -1.0+0.0im 1.0+0.0im … -1.0+0.0im
1.0+0.0im -0.809017+0.587785im 0.309017-0.951057im -0.809017-0.587785im
1.0+0.0im -0.309017+0.951057im -0.809017-0.587785im -0.309017-0.951057im
1.0+0.0im 0.309017+0.951057im -0.809017+0.587785im 0.309017-0.951057im
1.0+0.0im 0.809017+0.587785im 0.309017+0.951057im 0.809017-0.587785im
One simple way to get the matrix representation of an FFT or in fact any function that represents a linear operator is to get some AD package to compute its Jacobian for you:
julia> using Diffractor: DiffractorForwardBackend
julia> using AbstractDifferentiation: jacobian
julia> using FFTW
julia> only(jacobian(DiffractorForwardBackend(), fft, rand(5)))
5×5 Matrix{ComplexF64}:
1.0+0.0im 1.0+0.0im 1.0+0.0im 1.0+0.0im 1.0+0.0im
1.0+0.0im 0.309017-0.951057im -0.809017-0.587785im -0.809017+0.587785im 0.309017+0.951057im
1.0+0.0im -0.809017-0.587785im 0.309017+0.951057im 0.309017-0.951057im -0.809017+0.587785im
1.0+0.0im -0.809017+0.587785im 0.309017-0.951057im 0.309017+0.951057im -0.809017-0.587785im
1.0+0.0im 0.309017+0.951057im -0.809017+0.587785im -0.809017-0.587785im 0.309017-0.951057im
# some other function
julia> only(jacobian(DiffractorForwardBackend(), x -> 3circshift(x,2), rand(5)))
5×5 Matrix{Float64}:
0.0 0.0 0.0 3.0 0.0
0.0 0.0 0.0 0.0 3.0
3.0 0.0 0.0 0.0 0.0
0.0 3.0 0.0 0.0 0.0
0.0 0.0 3.0 0.0 0.0
(will depend on the AD library being able to handle your function, here eg Diffractor works but ForwardDiff didn’t)
2 Likes
Sorry, it’s just fft(I(n), 1) — it’s positional, not a keyword. (This API pre-dates the existence of keyword arguments.)
1 Like
Just applying it to the columns of I is a lot clearer to me (and is a lot more likely to work … lots of functions aren’t AD-able).
I will be needing to calculate derivatives with respect to the vector x in DFT{x}. I thought that if had the DFT matrix F such that DFT{x} = F*x, then I could calculate all the derivatives I might need without relying on AD packages being able to differentiate through the DFT routines.
But this doesn’t really answer the question — why do you need the Jacobian matrix explicitly (i.e. the explicit DFT matrix F)?
In most applications, you should only need a linear operator that acts the Jacobian (or its inverse or transpose) on a vector. In that case you are much better off using an FFT (or an inverse FFT).
(That being said, there are ChainRules defined for AbstractFFTs.jl, so many AD packages should be able to handle them for you. But it’s still good to know how to do it manually in an efficient way.)
1 Like
I have a long computation (compositions of many functions) that depends on some parameters. One step in these computations, involves the DFT of the result so far. I will be needing to calculate the gradient of this long computation with respect to the parameters. I still haven’t figured out how I will do this. What’s for sure is that, it personally helps me to write out the computations on paper using the DFT matrix F so that I can reason about applying the chain rule. I can then test my ideas numerically, by replicating them in code using again the DFT matrix F.
Also: typically I use ForwardDiff.jl as it appears to be more stable than other options, even though it can be slower. ForwardDiff.jl doesn’t seem to work with AbstractFFTs.jl, but I might be wrong…
In most applications, you should only need a linear operator that acts the Jacobian (or its inverse or transpose) on a vector. In that case you are much better off using an FFT (or an inverse FFT).
I think your advice should hold in my case too. Given my little experience in using the DFT, having the explicit matrix will help me verify my calculations and strengthen my understanding of it.
If you are computing a gradient, then definitely you will only need a vector-Jacobian product. This is the heart of reverse-mode differentiation, also called backpropagation or adjoint methods, and it can be done manually as well as by AD.
In the case of a DFT, the transpose or adjoint are both DFTs as well that can be computed with FFTs. So you will not need an explicit matrix in the end if you do things properly.
4 Likes
Thank you for your input! Thanks everyone for their comments. My indirect question lead to an answer to the direct problem I am facing…