# FFT and IntervalArithmetic

**URL:** <https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756>\
**Category:** General Usage\
**Created:** [November 19, 2018, 10:49pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756 "2018-11-19T22:49:18Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [November 19, 2018, 10:49pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/1 "2018-11-19T22:49:18Z")

</div>

Good evening,

I keep getting the following error message

`ERROR: LoadError: MethodError: no method matching plan_bfft(::Array{Complex{IntervalArithmetic.Interval{Float64}},2}, ::UnitRange{Int64})`

Does anyone know if `IntervalArithmetic` is compatible with `fft` ? I could not find any documentation about it.  
Or perhaps this is due to a different bug in my code.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 19, 2018, 11:24pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/2 "2018-11-19T23:24:06Z")

</div>

I’ve spent some time thinking about this myself, mostly in the context of convolution. My tentative conclusion is that it’s not a straightforward operation, and at least for the purposes of convolution it plays quite badly with respect to interval arithmetic’s tendency to widen the interval for each arithmetic operation. If your goal is to compute convolutions with small kernels, I’ve had more luck with real-space approaches.

---

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [November 19, 2018, 11:35pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/3 "2018-11-19T23:35:54Z")

</div>

Unfortunately, in general, I need to do convolutions on arrays that are not small. Using loops, computing these convolutions (a priori any number of them) becomes a nightmare.

Do you think this would eventually be implemented in the package ?

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 20, 2018, 12:24am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/4 "2018-11-20T00:24:08Z")

</div>

As a C-library wrapper, FFTW probably can’t support this properly. But a pure-Julia implementation of the fast DFT should do it straightforwardly. I seem to remember the existence of such a package but my search skills were not up to finding it quickly. Maybe someone else knows?

Whether it’s _useful_ is another question, I’d be interested in your experience.

---

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [November 20, 2018, 12:42am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/5 "2018-11-20T00:42:25Z")

</div>

About the use:  
I work on computer-assisted proofs in differential equations. Along the steps, a key part of the work is to write the nonlinear vector field as a polynomial.  
Using the regularity/property of the solutions of the differential operator, we pick a basis and wish to compute the coefficients (e.g Taylor/Fourier series).  
Since the vector field is a polynomial, discrete convolutions occurs and this is where `FFT` comes in.

So, computing convolutions efficiently with interval arithmetics is a crucial part of the work. Until now, the work was done using Matlab and IntLab. Hopefully there is a solution, my experience with Julia and `IntervalArithmetic` has been great 😁

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 20, 2018, 1:56am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/6 "2018-11-20T01:56:26Z")

</div>

Do you actually need the FFT? Isn’t the DFT enough? I.e. just write out the sum that occurs in the definition. In this case, it should “just work” with IntervalArithmetic.jl, except possibly that complex number support is not so great still.

I am not an expert, but it seems like something like the following should work:

```julia

julia> using IntervalArithmetic

julia> mypi = pi_interval(Float64)
[3.14159, 3.1416]

julia> ii = 0 + im
0 + 1im

julia> exp(-2*mypi*ii / 10)
[0.809016, 0.809017] - [0.587785, 0.587786]*im

julia> xx = [i..i for i in 1:1024];

julia> k = 3
3

julia> N = length(xx)
1024
julia> Xk = sum(xx[n] * exp(-2*mypi*ii / 1024 * k * n) for n in 1:N)
[511.999, 512.001] + [55627.1, 55627.2]*im

```

---

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [November 20, 2018, 2:16am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/7 "2018-11-20T02:16:20Z")

</div>

Isn’t FFT just a fast implementation of the DFT ?  
Indeed, as for now the DFT will have to do.  
I am curious, could you elaborate on “complex number support is not so great still” ? Or perhaps provide a link ?

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 20, 2018, 2:18am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/8 "2018-11-20T02:18:13Z")

</div>

Yes, exactly, FFT is just a fast implementation of DFT. But it’s quite possible / probable that you don’t really need it to be “fast”?

I was wondering whether the complex intervals in that equation would work, but apparently they do. At some point we are expecting to no longer have an `Interval` be a subtype of `Real`, in which case we will have to write our own complex routines (some of which already exist).

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [November 20, 2018, 2:19am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/9 "2018-11-20T02:19:06Z")

</div>

If you want to do anything rigorous with this stuff, then normally you will need to estimate the error when using a finite number of points in your DFT. But I’m guessing that you already have this under control.

---

<div class="post-metadata">

**Author:** ![OlivierHnt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivierhnt/32/6227_2.png) [@OlivierHnt](https://discourse.julialang.org/u/OlivierHnt)\
**Post date:** [November 20, 2018, 2:34am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/10 "2018-11-20T02:34:21Z")

</div>

Well I guess you are right, some problems do not require a “fast” DFT.  
Nonetheless, the faster the better 🙂. Especially if we have convolutions on a 2d / 3d solution. To give some idea about the order, I am dealing with an array containing about 1000 elements for which I do a cubic convolution, so I would have liked to get the runtime O(n\*log(n)) given by fft.

Yes indeed, I have to be more careful with the error if I go with the DFT.  
Thanks for your help!

---

<div class="post-metadata">

**Author:** ![baggepinnen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/baggepinnen/32/693_2.png) [@baggepinnen](https://discourse.julialang.org/u/baggepinnen)\
**Post date:** [November 20, 2018, 5:00am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/11 "2018-11-20T05:00:43Z")

</div>

I have successfully used the implementation here  
[https://github.com/MikaelSlevinsky/FastTransforms.jl/blob/master/src/fftBigFloat.jl#L16](https://github.com/MikaelSlevinsky/FastTransforms.jl/blob/master/src/fftBigFloat.jl#L16)  
The functions there are type restricted to BigFloat etc., but these restrictions can be removed and they work equally well on DualNumbers etc.

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [November 20, 2018, 7:57am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/12 "2018-11-20T07:57:06Z")

</div>

> [@tim.holy](#):
>
> at least for the purposes of convolution it plays quite badly with respect to interval arithmetic’s tendency to widen the interval for each arithmetic operation

This is surprising since the FFT is both an orthogonal operation and stable. It seems like one would be able to find sharp L2 bounds on the error in the FFT in the literature, and then one could embed a vector of intervals inside a ball.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 20, 2018, 8:55am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/13 "2018-11-20T08:55:36Z")

</div>

It’s not necessarily “error in the FFT” that you’re trying to estimate here. Suppose you want to take the DFT of `fill(a, n)` where `a = -1..1`: among the concrete instantiations of this signal are `fill(0, n)` (which has a DFT of all-zeros) and `sin.(α*x)` (which has all power concentrated at single frequencies). As a consequence, to me it seems the DFT will have to span something like `fill((-n..n) + (-n..n)*im, n)`. Then when you take the IFFT of this you get something considerably wider than the original signal.

What you’re running up against here is conceptually equivalent to the fact that while `1 - 1 == 0`, `(0..1) - (0..1) == (-1..1)`. You might be able to do better using affine arithmetic (because it keeps track of where each interval comes from and thus can permit more cancellation) but it’s not obvious to me that this can be performed efficiently.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 20, 2018, 9:17am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/14 "2018-11-20T09:17:54Z")

</div>

> To give some idea about the order, I am dealing with an array containing about 1000 elements for which I do a cubic convolution, so I would have liked to get the runtime O(n\*log(n)) given by fft.

If you’re doing naive FFTs, the `O(n^2)` nature of the problem is the same as if you did it by direct (real-space) convolution. Moreover, with interval arithmetic the two are _not_ guaranteed to come to the same answer. My intuition is that you will likely get a much better answer (i.e., with tighter bounds) doing it in real space, since the FFT “mixes” the widths of all the intervals whereas the real-space version may (for kernels with many zeros) do less mixing.

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [November 20, 2018, 11:09am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/15 "2018-11-20T11:09:15Z")

</div>

@dlfivefifty well there is a “problem” (we consider this a beautiful feature) with complex multiplication that it rotates. So a simple multiplication by `e^(iπ/2)` of a square box will ask You to widen the box by `sqrt(2)`…

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [November 20, 2018, 1:41pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/16 "2018-11-20T13:41:40Z")

</div>

Yes if you rotate 2 entries at a time it widens a lot. But if you have a bound on the L^2 norm you should be able to bound the total error ( that is, embed a product of intervals in an n-ball) without bounding each operation at a time.

---

<div class="post-metadata">

**Author:** ![abulak](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abulak/32/28314_2.png) [@abulak](https://discourse.julialang.org/u/abulak)\
**Post date:** [November 21, 2018, 4:50pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/17 "2018-11-21T16:50:05Z")

</div>

That’s true, but with IntervalArithmetic You only get L-∞ norm. The one You’re probably thinking about is ball arithmetic?where You store best approximation of `mid` and minimal `radius` of ball in L^p

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [November 21, 2018, 5:10pm UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/18 "2018-11-21T17:10:26Z")

</div>

I guess what I’m thinking is if you know the worst case L^infty norm you don’t need to do each step of the algorithm with interval arithmetic: do it with floating point arithmetic using the midpoints and then decide the output intervals using the norm. Using L2 norms would just be an intermediary step, as it’s an orthogonal operation in L2.

Unless I’m mistaken the intervals should only grow like O(n), which should be pretty sharp.

---

<div class="post-metadata">

**Author:** ![tim.holy](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tim.holy/32/52_2.png) [@tim.holy](https://discourse.julialang.org/u/tim.holy)\
**Post date:** [November 22, 2018, 11:42am UTC](https://discourse.julialang.org/t/fft-and-intervalarithmetic/17756/19 "2018-11-22T11:42:22Z")

</div>

If you come up with something concrete I’d certainly be interested in hearing about it.