# Discrete Fourier Transform (DFT) matrix and inverse

**URL:** <https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952>\
**Category:** Numerics\
**Tags:** fourier\
**Created:** [September 27, 2024, 11:33am UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952 "2024-09-27T11:33:57Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Nikos\_Gianniotis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikos_gianniotis/32/11487_2.png) [@Nikos\_Gianniotis](https://discourse.julialang.org/u/Nikos_Gianniotis)\
**Post date:** [September 27, 2024, 11:33am UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/1 "2024-09-27T11:33:57Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![John\_Gibson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/john_gibson/32/5321_2.png) [@John\_Gibson](https://discourse.julialang.org/u/John_Gibson)\
**Post date:** [September 27, 2024, 11:53am UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/2 "2024-09-27T11:53:38Z")

</div>

You could calculate it by taking the FFT of unit vectors, using a trustworthy algorithm like FFTW.

---

<div class="post-metadata">

**Author:** ![pablosanjose](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pablosanjose/32/7006_2.png) [@pablosanjose](https://discourse.julialang.org/u/pablosanjose)\
**Post date:** [September 27, 2024, 12:22pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/3 "2024-09-27T12:22:59Z")

</div>

For the benefit of the reader, DFT probably stands for “Discrete Fourier Transform” (not e.g. “Density Functional Theory” 😆)

---

<div class="post-metadata">

**Author:** ![carstenbauer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/carstenbauer/32/4981_2.png) [@carstenbauer](https://discourse.julialang.org/u/carstenbauer)\
**Post date:** [September 27, 2024, 1:03pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/4 "2024-09-27T13:03:23Z")

</div>

I’ve updated the title to make this clear.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [September 27, 2024, 1:18pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/5 "2024-09-27T13:18:00Z")

</div>

Something easy which probably does the job:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [September 27, 2024, 1:34pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/6 "2024-09-27T13:34:42Z")

</div>

> [@roflmaostc](#):
>
> `W = [norm * ω^((i-1) * (j-1)) for i ∈ 1:N, j ∈ 1:N]`

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 `(i+j-2)`)~~. (typo, see below)

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [September 27, 2024, 1:37pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/7 "2024-09-27T13:37:13Z")

</div>

> [@kellertuer](#):
>
> (i+j)

Are you sure?  
It’s `i * j` for `i, j = 0, 1, ..., N-1`.

Also of course `(i + j - 2)` can’t be derived.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [September 27, 2024, 1:42pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/8 "2024-09-27T13:42:46Z")

</div>

Oh – yeah that’s a typo, sorry its Friday after and exhausting week. Of course I meant `*`.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [September 27, 2024, 4:15pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/9 "2024-09-27T16:15:07Z")

</div>

Surprisingly (?), the FFT approach suggested by @John_Gibson is faster than building the matrix in a comprehension for larger N:

```julia
using LinearAlgebra
using FFTW

DFT_mat_fft(N) = stack(fft.(eachcol(I(N))))

```

```julia-repl
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)

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 27, 2024, 4:31pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/10 "2024-09-27T16:31:23Z")

</div>

> [@bertschi](#):
>
> `DFT_mat_fft(N) = stack(fft.(eachcol(I(N))))`

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.

---

<div class="post-metadata">

**Author:** ![roflmaostc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/roflmaostc/32/30123_2.png) [@roflmaostc](https://discourse.julialang.org/u/roflmaostc)\
**Post date:** [September 27, 2024, 4:48pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/11 "2024-09-27T16:48:19Z")

</div>

> [@stevengj](#):
>
> `fft(I(n), dims=1)`

The `dims` keyword is not used for AbstractFFTs.jl

```julia
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

```

---

<div class="post-metadata">

**Author:** ![marius311](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marius311/32/3953_2.png) [@marius311](https://discourse.julialang.org/u/marius311)\
**Post date:** [September 27, 2024, 4:51pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/12 "2024-09-27T16:51:41Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 27, 2024, 5:11pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/13 "2024-09-27T17:11:14Z")

</div>

> [@roflmaostc](#):
>
> > [@stevengj](#):
> >
> > `fft(I(n), dims=1)`
> 
> The `dims` keyword is not used for [AbstractFFTs.jl](https://juliahub.com/ui/Packages/General/AbstractFFTs)

Sorry, it’s just `fft(I(n), 1)` — it’s positional, not a keyword. (This API pre-dates the existence of keyword arguments.)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 27, 2024, 5:13pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/14 "2024-09-27T17:13:49Z")

</div>

> [@marius311](#):
>
> 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:

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

---

<div class="post-metadata">

**Author:** ![Nikos\_Gianniotis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikos_gianniotis/32/11487_2.png) [@Nikos\_Gianniotis](https://discourse.julialang.org/u/Nikos_Gianniotis)\
**Post date:** [September 27, 2024, 8:52pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/15 "2024-09-27T20:52:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 27, 2024, 9:04pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/16 "2024-09-27T21:04:07Z")

</div>

> [@Nikos\_Gianniotis](#):
>
> 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.)

---

<div class="post-metadata">

**Author:** ![Nikos\_Gianniotis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikos_gianniotis/32/11487_2.png) [@Nikos\_Gianniotis](https://discourse.julialang.org/u/Nikos_Gianniotis)\
**Post date:** [September 27, 2024, 9:13pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/17 "2024-09-27T21:13:26Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [September 27, 2024, 11:58pm UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/18 "2024-09-27T23:58:42Z")

</div>

> [@Nikos\_Gianniotis](#):
>
> 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.

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.

---

<div class="post-metadata">

**Author:** ![Nikos\_Gianniotis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nikos_gianniotis/32/11487_2.png) [@Nikos\_Gianniotis](https://discourse.julialang.org/u/Nikos_Gianniotis)\
**Post date:** [September 28, 2024, 7:02am UTC](https://discourse.julialang.org/t/discrete-fourier-transform-dft-matrix-and-inverse/119952/19 "2024-09-28T07:02:40Z")

</div>

Thank you for your input! Thanks everyone for their comments. My indirect question lead to an answer to the direct problem I am facing…