Vico-Greengard-Ferrando quadratures in the tensor solver for integral equations

TL;DR

This paper introduces a tensor-based method using Vico-Greengard-Ferrando quadratures for efficient, spectrally accurate convolution computations with Green's functions, applicable to integral equations like Lippmann-Schwinger.

Contribution

It develops a low-rank tensor implementation using QTT decomposition for convolution with Green's functions, enhancing computational efficiency and accuracy.

Findings

Spectral accuracy achieved with FFT-based quadratures.

Efficient low-rank tensor representation developed.

Applicable to various integral equations.

Abstract

Convolution with Green's function of a differential operator appears in a lot of applications e.g. Lippmann-Schwinger integral equation. Algorithms for computing such are usually non-trivial and require non-uniform mesh. However, recently Vico, Greengard and Ferrando developed method for computing convolution with smooth functions with compact support with spectral accuracy, requiring nothing more than Fast Fourier Transform (FFT). Their approach is very suitable for the low-rank tensor implementation which we develop using Quantized Tensor Train (QTT) decomposition.

Tables

Table 1. Table 1 : Relative error (err) and effective rank (erank) for different grid sizes and rounding errors ϵ italic-ϵ \epsilon . Results are presented for two types of function q 𝑞 q : Gaussian ( 11 ) and smoothed cube ( 12 ).

	3D grid size	${32}^3$		${64}^3$		${128}^3$
		err	erank	err	erank	err	erank
Gaussian	$ϵ={10}^{-3}$	8e-1	14	8e-2	15	7e-2	14
	$ϵ={10}^{-5}$	8e-1	24	4e-4	25	4e-4	26
	$ϵ={10}^{-7}$	8e-1	30	5e-6	34	4e-6	38
Smoothed cube	$ϵ={10}^{-3}$	8e-1	22	8e-2	35	6e-2	34
	$ϵ={10}^{-5}$	8e-1	29	5e-4	57	5e-4	59
	$ϵ={10}^{-7}$	8e-1	29	2e-5	73	6e-6	80

Equations

$$h({\text{r}}^{\prime })={\int }_{{\mathbb{R}}^3}g(\text{r}-{\text{r}}^{\prime })\rho (\text{r})d\text{r},$$

$$h(\text{r})={\int }_{{\mathbb{R}}^3}g(\text{r}-{\text{r}}^{\prime })\rho ({\text{r}}^{\prime })d{\text{r}}^{\prime },$$

$${\int }_{{\mathbb{R}}^3}g(\text{r}-{\text{r}}^{\prime })\rho ({\text{r}}^{\prime })d{\text{r}}^{\prime }={\int }_{{\mathbb{R}}^3}{g}_L(\text{r}-{\text{r}}^{\prime })\rho ({\text{r}}^{\prime })d{\text{r}}^{\prime }.$$

$$g(\text{r})=\frac{e^{\mathrm{i}kr}}{4\pi r},$$

$$G_L(\text{s})={\int }_{{\mathbb{R}}^3}g(\text{r})\mathrm{rect}\left(\frac{r}{2L}\right)e^{\mathrm{i}\text{s}\cdot \text{r}}d\text{r}$$

$$=\frac{-1+e^{\mathrm{i}Lk}(\cos Ls-\mathrm{i}\frac{k}{s}\sin Ls)}{(k-s)(k+s)}.$$

$$h(\textbf{r})=\Big{(}\frac{1}{2\pi}\Big{)}^{3}\int_{\mathbb{R}^{3}}e^{\mathrm{i}\textbf{s}\cdot\textbf{r}}\widehat{\rho}(\textbf{s})G_{L}(\textbf{s})d\textbf{s}.$$

$${\rho }_{ijk}:=\rho (i{h}_1,j{h}_2,k{h}_3),h_i=\frac{L_i}{N_i}.$$

$$h_{ijk}=\frac{1}{4N_1}\frac{1}{4N_2}\frac{1}{4N_3}\sum _{s_1=-2N_1}^{2N_1-1}\sum _{s_2=-2N_2}^{2N_2-1}\sum _{s_3=-2N_3}^{2N_3-1}{G}_L(\text{s})\hat{\rho }(\text{s})e^{2\pi \mathrm{i}\frac{s_{1}i}{4N_1}}{e}^{2\pi \mathrm{i}\frac{s_{2}j}{4N_2}}{e}^{2\pi \mathrm{i}\frac{s_{3}k}{4N_3}},$$

$$\hat{\rho }(\text{s})=\sum _{i^{\prime }=0}^{4N_1-1}\sum _{j^{\prime }=0}^{4N_2-1}\sum _{k^{\prime }=0}^{4N_3-1}{\rho }_{i^{\prime }{j}^{\prime }{k}^{\prime }}{e}^{-2\pi \mathrm{i}\frac{s_1{i}^{\prime }}{4N_1}}{e}^{-2\pi \mathrm{i}\frac{s_2{j}^{\prime }}{4N_2}}{e}^{-2\pi \mathrm{i}\frac{s_3{k}^{\prime }}{4N_3}}.$$

$$h_{ijk}=\sum _{i^{\prime }=0}^{4N_1-1}\sum _{j^{\prime }=0}^{4N_2-1}\sum _{k^{\prime }=0}^{4N_3-1}{G}_{i-i^{\prime },j-j^{\prime },k-k^{\prime }}^M{\rho }_{i^{\prime }{j}^{\prime }{k}^{\prime }},$$

$$G_{i-i^{\prime },j-j^{\prime },k-k^{\prime }}^M=\frac{1}{4N_1}\frac{1}{4N_2}\frac{1}{4N_3}\sum _{s_1=-2N_1}^{2N_1-1}\sum _{s_2=-2N_2}^{2N_2-1}\sum _{s_3=-2N_3}^{2N_3-1}{G}_L(\text{s})e^{2\pi \mathrm{i}\frac{s_1(i-i^{\prime })}{4N_1}}{e}^{2\pi \mathrm{i}\frac{s_2(j-j^{\prime })}{4N_2}}{e}^{2\pi \mathrm{i}\frac{s_3(k-k^{\prime })}{4N_3}}.$$

$$h_{ijk}=\sum _{i^{\prime }=0}^{N_1-1}\sum _{j^{\prime }=0}^{N_2-1}\sum _{k^{\prime }=0}^{N_3-1}{G}_{i-i^{\prime },j-j^{\prime },k-k^{\prime }}^M{\rho }_{i^{\prime }{j}^{\prime }{k}^{\prime }}.$$

$${\mathcal{X}}_{i_1{i}_2\dots i_d}=X^{(1)}(i_1)X^{(2)}(i_2)\dots X^{(d)}(i_d),$$

$$i=i_1+2^1{i}_2+\cdots +2^{l_i-1}{i}_d,i_m=0,1,m=1,\dots ,d$$

$${\rho }_{ijk}\equiv {\tilde{\rho }}_{i_1\dots i_{l_1}{j}_1\dots j_{l_2}{k}_1\dots k_{l_3}}.$$

$${\mathcal{Y}}_{i_1\dots i_d}=\sum _{j_1\dots j_d}{\mathcal{A}}_{i_1\dots i_d{j}_1\dots j_d}{\mathcal{X}}_{j_1\dots j_d}.$$

$${\mathcal{A}}_{i_1\dots i_d{j}_1\dots j_d}=A^{(1)}(i_1,j_1)A^{(2)}(i_2,j_2)\dots A^{(d)}(i_d,j_d),$$

$$\begin{array}{rl}{G}_{i-i^{\prime },j-j^{\prime },k-k^{\prime }}^M=& G^{(1)}(i_1,i_1^{\prime })\dots G^{(l_1+2)}(i_{l_1+2},i_{l_1+2}^{\prime })\\ & G^{(l_1+3)}(j_1,j_1^{\prime })\dots G^{(l_1+l_2+4)}(j_{l_2+2},j_{l_2+2}^{\prime })\\ & G^{(l_1+l_2+5)}(k_1,k_1^{\prime })\dots G^{(l_1+l_2+l_3+6)}(k_{l_3+2},i_{l_3+2}^{\prime }).\end{array}$$

$$\mathrm{DOF}(\mathcal{X})=2(r_1+r_d)+\sum _{i=2}^{d}2r_{i-1}{r}_i.$$

$$\sigma (\text{r})+k^{2}q(\text{r}){\int }_{{\mathbb{R}}^3}\frac{e^{\mathrm{i}k|\text{r}-{\text{r}}^{\prime }|}}{|\text{r}-{\mathbf{r}}^{\prime }|}\sigma ({\text{r}}^{\prime })d{\text{r}}^{\prime }=-k^{2}q(\text{r}){\varphi }^{\mathrm{inc}},$$

$${\varphi }^{\mathrm{scat}}={\int }_{{\mathbb{R}}^3}\frac{e^{ik|\text{r}-{\text{r}}^{\prime }|}}{|\text{r}-{\text{r}}^{\prime }|}\sigma ({\text{r}}^{\prime })d{\text{r}}^{\prime }.$$

$$q(\text{r})=e^{\frac{-|\text{r}-{\text{r}}^{\prime }{|}^2}{2a^2}},$$

$${\varphi }^{\mathrm{inc}}(x,y,z)=e^{\mathrm{i}x}.$$

$$q(\text{r})=e^{-0.5{\left(\frac{|\text{r}-{\text{r}}^{\prime }|}{a}\right)}^8},$$

$$q(\text{r})=\sum _{i=1}^{19}\sum _{j=0}^3\sum _{k=0}^1{q}_{(0.25+0.5j,0.25+0.5k,0.5i),0.1}(\text{r})$$

$${\varphi }^{\mathrm{inc}}(x,y,z)=e^{4\pi \mathrm{i}z}.$$

Full text

Vico-Greengard-Ferrando quadratures in the tensor solver for integral equations

V. Khrulkov1, M. Rakhuba1, **I. Oseledets1,2

1**Skolkovo Institute of Science and Technology, Russia

2Institute of Numerical Mathematics, Russia

Abstract— Convolution with Green’s function of a differential operator appears in a lot of applications e.g. Lippmann-Schwinger integral equation. Algorithms for computing such are usually non-trivial and require non-uniform mesh. However, recently Vico, Greengard and Ferrando developed method for computing convolution with smooth functions with compact support with spectral accuracy, requiring nothing more than Fast Fourier Transform (FFT). Their approach is very suitable for the low-rank tensor implementation which we develop using Quantized Tensor Train (QTT) decomposition.

1 Introduction

In this paper we propose an algorithm for computing an approximation to a convolution

[TABLE]

where $G$ is a continuous Green’s function of some PDE. It has been shown recently in [1] that discrete approximation to $h$ yieding spectral accuracy for compactly supported smooth functions $\rho$ can be constructed using relatively simple idea. Idea is based on the fact that continuous Green’s function can be replaced by the truncated one, Fourier transform of which is in $\mathcal{C}^{\infty}(\mathbb{R}^{3})$. Applying convolution theorem one gets answer as a Fourier type integral of a rapidly decreasing function from $\mathcal{C}^{\infty}(\mathbb{R}^{3})$ which is then discretized on finer grid and evaluated using Discrete Fourier Transform (DFT). This scheme achieves spectral order of accuracy due to the superalgebraic approximation of the continuous Fourier transform by DFT and has been shown to also give second order with small constant for continuos functions. We implemented this scheme using data compression via Quantized Tensor Train (QTT). Approach is straightforward – we compute discrete truncated Green’s function using formulas from [1] and convert it to QTT-format. Then we use standard algorithms in QTT which allow for logarithmic storage with respect to the grid size and accuracy of approximation. Moreover, discrete convolution and solution of integral equations are then performed also with logarithmic complexity by the algorithms described in [2],[3]. Paper is organized as follows:

• •

We briefly describe the algorithm developed in Section 2

• •

We present necessary definitions and algorithms from TT theory in Section 3.1

• •

We describe TT implementation of the algorithm in Section 3.2

• •

We present numerical results in Section 4

The main reference for Vico-Greengard-Ferrando quadrature is the original paper [1], in the presentation of the algorithm we follow [4].

2 Algorithm

We wil present the scheme described in [1] in a way which is most suitable for transitioning to TT-format. Let $\rho(\textbf{r})$ be a smooth function such that $\operatorname{supp}\rho\subset D=[0,L_{1}]\times[0,L_{2}]\times[0,L_{3}]$ and we are interested in computing

[TABLE]

where $g(\textbf{r})$ is a Green’s function of some differential operator. Let us assume that $g(\textbf{r})$ depends only on $r=|\textbf{r}|$. If we seek restriction of the solution $h(\textbf{r})$ to $D$ then $|r-r^{\prime}|$ in $\eqref{eq1}$ doesn’t exceed $L=\sqrt{L_{1}^{2}+L_{2}^{2}+L_{3}^{2}}$. Thus if we replace $g(\textbf{r})$ by $g_{L}(\textbf{r})=g(\textbf{r})\operatorname{rect}{\frac{r}{2L}}$,

[TABLE]

Advantage of using $g_{L}(\textbf{r})$ is that since it has compact support it’s Fourier transform $G_{L}(\textbf{s})$ is in $C^{\infty}(\mathbb{R}^{3})$, and is straightforward to compute for many differential operators. We mainly focus on Helmholtz differential operator $\nabla^{2}+k^{2}$ for which the following formula holds:

[TABLE]

It is easy to check that $G_{L}(\textbf{s})$ is indeed smooth and nonsingular. If we denote Fourier transform of $\rho(\textbf{r})$ by $\widehat{\rho}(\textbf{s})$ we obtain the final formula:

[TABLE]

This integral is then discretized using trapezoidal rule and computed using DFT. However to cancel slightly oscillatory behavior of $\widehat{G_{L}}(s)$ zero padding by a factor of at least $3$ is required (we will use factor $4$ to keep grid size being a power of $2$ - for the analysis see [4]). Suppose that domain $D$ is discretized using uniform grid with $N_{i}$ nodes in corresponding dimensions and function $\rho$ is sampled on this grid yielding an array $\rho_{ijk}$ i.e.

[TABLE]

Algorithm then is summarized as following

2.1 Convolution form of the algorithm

Let us write Algorithm 1 more explicitly. All steps together can be represented as follows

[TABLE]

where

[TABLE]

By plugging (5) into (4) and by changing the order of summation it is easy to see that

[TABLE]

where

[TABLE]

Moreover, since $\rho_{i^{\prime}j^{\prime}k^{\prime}}$ is [math] for $i^{\prime}\geq N_{1},j^{\prime}\geq N_{2},k^{\prime}\geq N_{3}$ (see step $1$ of Algorithm 1) and we truncate the result, formula (6) simplifies and finally:

[TABLE]

We see that (8) takes the form of a discrete aperiodic convolution with discrete Green’s function $G^{M}$ (which we will call mollified Green’s function). One can notice that to fully determine $G^{M}$ it is sufficient to run Algorithm 1 once for a special right hand side $\rho_{ijk}=\delta_{i0}\delta_{j0}\delta_{k0}$. Formula (8) plays essential role in the further analysis. Multiplication by multilevel Toeplitz matrix generated by $G^{M}$ can be performed with logarithmic complexity in QTT format as described in [2], and we discuss neccessary definitions and algorithms in the next section.

3 Low-rank tensor approach

3.1 TT and QTT formats

To understand the QTT format let us start with describing the TT-format, which is a nonlinear low-paramentric representation of multidimensional arrays, called tensors. Tensor $\mathcal{X}\in\mathbb{C}^{n_{1}\times n_{2}\times\dots\times n_{d}}$ is said to be in the TT-format if it represents as

[TABLE]

where $X^{(k)}(i_{k})\in\mathbb{C}^{r_{k-1}\times r_{k}}$, $r_{0}=r_{d}=1$, $i_{k}=1,\dots,n_{k}$. Matrices $X^{(k)}$ are called TT-cores and $r_{k}$ are called TT-ranks. Notice that if $r=\max_{k}r_{k}$ is small, then there is a significant compression to store $\mathcal{X}$. Indeed, initial tensor requires storing $n^{d}$ parameters, while to store its TT-representation only $\mathcal{O}(dnr^{2})$ parameters are needed.

In fact, one could use TT representation to store and to work with arising in Algorithm 1 3-dimensional arrays. However we will use a more sophisticated approach called QTT format, which allows for additional storage reduction compared to TT. QTT format is the following modification of the TT format. First we assume that $d=3$, $n_{i}=2^{l_{i}}$, $i=1,2,3$. Then each “physical” index $i$, $j$, $k$ is represented in the binary format, i.e.

[TABLE]

and we have initial tensor $\rho_{ijk}$ encoded as a $(l_{1}+l_{2}+l_{3})$-dimensional array $\tilde{\rho}$:

[TABLE]

TT decomposition of $\tilde{\rho}$ is called the QTT decomposition. The storage of the QTT decomposition is $\mathcal{O}(r^{2}(l_{1}+l_{2}+l_{3}))=\mathcal{O}(r^{2}\log n)$, so if $r=\max_{i}r_{i}$ is bounded, the total storage scales logarithmically. In practice tensors of exact low-rank rarely occur. Typically one fixes accuracy $\epsilon$ and tries to find best approximation with this accuracy. It has been shown that in some applications ranks grow as $r=\mathcal{O}(\log^{\alpha}\epsilon^{-1})$, $\alpha>1$ [5, 6].

3.2 Translation of the algorithm to the QTT format

Computation of kernel $G^{M}$.

To use (8) we first need to find $G^{M}$ (7). For this purpose we run Algorithm 1 for $\rho_{ijk}=\delta_{i0}\delta_{j0}\delta_{k0}$. Precomputations are done in the full format, in other words we form the whole dense tensor $G^{M}$ and utilize TT-SVD algorithm [7] to find its QTT representation. TT-SVD algorithm is based on the computation of SVD decompositions of tensors reshaped into full 2D matrices and therefore is quite expensive. In principle one could use DFT in the QTT format [8] to avoid forming full tensors. Unfortunately, we found that intermediate tensors arising in Algorithm 1 are of large rank. We will address this problem in our future work.

Computation of $\rho$.

Tensor $\rho$ can be already given in the QTT representation. This can happen, e.g. if we are running a certain iterative process involving computation of convolution (1) and all operations in this process are done within the QTT format. Otherwise, $\rho_{ijk}$ can be approximated with logarithmic complexity by using the cross approximation method [9], which adaptively samples elements of a tensor. In this case we just need $\rho_{ijk}$ be given as a function which returns value by given 3 indices $i,j,k$.

Computation of convolution $G^{M}*\rho$.

Next goal is to find convolution of tensors $G^{M}$ and $\rho$ (8). The convolution can be considered as a multiplication of multilevel Toeplitz matrix generated by $G^{M}$ and vector $\rho$. Matrices can also be represented in the TT and by analogy in the QTT format. The definition is similar to that of TT-tensor: given matrix (operator) $\mathcal{A}_{i_{1}\dots i_{d}j_{1}\dots j_{d}}$, which acts on vector $\mathcal{X}_{j_{1}\dots j_{d}}$ such that

[TABLE]

its TT-decomposition is defined as

[TABLE]

where $A^{(k)}(i_{k},j_{k})\in\mathbb{C}^{R_{k-1}\times R_{k}}$, $R_{0}=R_{d}=1$, $i_{k}=1,\dots,n_{k}$. QTT decomposition of 3D operator $G^{M}_{i-i^{\prime},j-j^{\prime},k-k^{\prime}}$ is defined by analogy with the QTT decomposition of a tensor – we quantize indices $i,j,k$ and $i^{\prime},j^{\prime},k^{\prime}$, group them pairwise and then compute TT decomposition:

[TABLE]

We use approach from [2] and analytically construct QTT representation of the induced multilevel Toeplitz matrix $G^{M}_{i-i^{\prime},j-j^{\prime},k-k^{\prime}}$ given QTT representation of $G^{M}_{i,j,k}$. Then matrix-vector product (8) can be done in different ways. We used optimization procedure AMEn (alternating minimal energy method) [10, 11] which allows for rank adaptation compared to standard ALS (alternating least squares) [12] optimization which works with the representation of a given size.

4 Numerical experiments

Approximating $G^{M}$ using QTT.

Firstly we show that using QTT representation greatly reduces number of degrees of freedom (DOF) of $G^{M}$. Suppose that tensor $\mathcal{X}$ is given in the QTT format with ranks $\{r_{1},\dots,r_{d}\}$. Then it is easy to count total number of DOF of $\mathcal{X}$:

[TABLE]

By applying this formula to $G^{M}$ computed as described in 3.2 for various values of $k$ (while keeping $L=1$) we obtained the following results (see Figure 1). This shows advantages of using QTT.

Solving scattering problems.

To further test our approach we solve the Lippmann-Schwinger equation which is used for solving scattering problems:

[TABLE]

and then we find

[TABLE]

We used rounding by $\epsilon=10^{-7}$ in our computations. Firstly we fixed $k=1$ and $L=32\pi$ and took $q(\textbf{r})$ to be a $3D$ gaussian:

[TABLE]

with $\textbf{r}^{\prime}=\left(\frac{L}{2},\frac{L}{2},\frac{L}{2}\right),$ and $a=\frac{L}{10},$ and

[TABLE]

To solve arising systems in the QTT format we used AMEn [10, 11] , which allows for rank adaptation.

Taking the solution computed on a grid with size $256^{3}$ as a reference we measured relative error of the solutions computed on smaller grids for various rounding parameters $\epsilon$. Results are given in the Table 1. We see that error is roughly equal to the $\epsilon$ even for modest number of grid nodes per wavelength.

As a next experiment we performed the same computations for $q(r)$ representing smoothed cube:

[TABLE]

for $\textbf{r}^{\prime}=\left(\frac{L}{2},\frac{L}{2},\frac{L}{2}\right)$ and $a=\frac{L}{4}$. For the results see Figure 3.

Scattering problems on quasiperiodic structures.

As a final test we took $q(\textbf{r})$ to be a periodic grid of smoothed cubes (structures like this are extremely suitable for QTT computations). Namely, denoting $q(\textbf{r})$ defined by formula (12) by $q_{\textbf{r}^{\prime},a}(\textbf{r})$ we solve the equation (10) in the domain $[0,2]\times[0,1]\times[0,10]$ on the grid $64\times 64\times 1024$ for $q(\textbf{r})$ defined as

[TABLE]

and for $k=4\pi$. $\phi^{\mathrm{inc}}$ in this case is a plane wave propagating in $z$-direction

[TABLE]

Figure 4 demonstrates our results.

Acknowledgements

This study was supported by the Ministry of Education and Science of the Russian Federation (grant 14.756.31.0001), by RFBR grants 16-31-60095-mol-a-dk, 16-31-00372-mol-a and by Skoltech NGP program.