Lecture 15. Low-rank methods and oscillations

Previous lecture

• Boundary integral equations for Helmholtz equation
• A word about multiscale
• Fast direct solvers for sparse matrices

Todays lecture

• Low-rank strikes back: butterfly & directional low-rank property
• A few words about multiscale
• High-dimensional problems and tensors(?)

Fourier transform

Fast Fourier Transform relies on the uniform grid!

$$g(x) = \int_a^b e^{i x y} f(y) dy $$

Then, the FFT is just

$$g_k = \sum_l w^{kl} f_l, \quad w^n = 1.$$

In fact, it should be FFFT (Fast Finite Fourier Transform), but FFT is better.

Quiz: Can you find low-rank submatrices in the Fourier matrix?

Reading

http://statweb.stanford.edu/~candes/papers/ButterflyFIO.pdf

Non-uniform case

Consider the evaluation of Fourier integral operators

$u(x) = \sum_{k \in \Omega} e^{2 \pi i \Phi(x,k)} f(k), \quad x \in X,$

where $\Omega = \{ -\frac{N}{2} \leq k_1, k_2 \leq \frac{N}{2} \}$.

Both $x$ and $k$ can be two (three) dimensional.

In two dimensions with $n \times n$ grid evaluation will take $\mathcal{O}(n^4)$

Assumption: $\Phi(x, k)$ is smooth and homogenious, $\Phi(x, \lambda k) = \lambda \Phi(x, k)$

For simplicity, consider 2D case.

Step 1

First step: represent frequency variables in polar coordinates

$k = (k_1, k_2) = \frac{1}{\sqrt{2}} Nnp_1 e^{2 \pi i p_2}, \quad e^{2 \pi i p_2} = (\cos 2 \pi i p_2, \sin 2 \pi i p_2)$

This guarantees that each point $(p_1, p_2)$ belongs to the unit square $[0, 1]^2$.

Due to homogenuity, $\Phi(x, k) = N \Psi(x, p)$

Step 2

Consider two square boxes $A$ and $B$ in $[0, 1]^2$, centered at $x_0$ and $p_0$, such that their side-lengths satisfy

$w(A) w(B) \leq \frac{1}{n}$

Introduce

$R^{AB}(x, p) = \Psi(x, p) - \Psi(x_0(A), p) - \Psi(x, p_0(B)) + \Psi(x_0(A), p_0(B))$

Step 3

Then, $ e^{2 \pi i n \Psi(x, p)} = e^{2 \pi i n \Psi(x_0(A), p)} e^{2 \pi i n \Psi(x, p_0(B))} \cdot e^{-2 \pi i n \Psi(x_0(A), p_0(B))} e^{2 \pi i n R^{AB}} $

Applying mean-value theorem to $R^{AB}$ we get $R^{AB} = \mathcal{O}(\frac{1}{n})$

It means, that the phase is bounded, thus can be approximated by low-rank

Butterfly algorithm

We need blocks (boxes $A$ and $B$) such that $w(A) w(B) = \frac{1}{N}$

How can we do that?

Compute interactions of large boxes with small boxes

Algorithm

1. Compute two dyadic trees $T_X$, and $T_P$ (as usual)
2. Initialization: Construct expansion coeficients $\delta^{AB}_t$ for root-leaf interaction
3. Recursion: Visit level $l$ in $T_X$ and $L - l$ in $T_P$ We can recompute from parent of $A$ and all children of $B$

Directional low-rank property

Finally, what to do for the Helmholtz kernel? Let $G = \exp(i k r)/r$.

What to do?

Use directional low-rank property

Directional low-rank property

Let us scale the problem such that $k = 1$ Important!!!

The algorithm

The final algorithm is equivalent to the adapted separation criteria.

Summary

Thank you.

Questions & App Period!