A fast direct solver for boundary value problems on locally perturbed geometries

TL;DR

This paper introduces a fast direct solver for boundary value problems on geometries with local perturbations, leveraging low rank updates and precomputed solutions to significantly reduce computation time.

Contribution

It presents a novel method that efficiently updates solutions for perturbed geometries using low rank updates and the Sherman-Morrison formula, saving computational effort.

Findings

Solver is three times faster than building from scratch for fixed perturbations.

Method effectively handles local boundary modifications and refined discretizations.

Numerical results demonstrate significant computational savings.

Abstract

Many applications involve solving several boundary value problems on geometries that are local perturbations of an original geometry. The boundary integral equation for a problem on a locally perturbed geometry can be expressed as a low rank update to the original system. A fast direct solver for the new linear system is presented in this paper. The solution technique utilizes a precomputed fast direct solver for the original geometry to efficiently create the low rank factorization of the update matrix and to accelerate the application of the Sherman-Morrison formula. The method is ideally suited for problems where the local perturbation is the same but its placement on the boundary changes and problems where the local perturbation is a refined discretization on the same geometry. Numerical results illustrate that for fixed local perturbation the method is three times faster than building a new fast direct solver from scratch.

Tables

Table 1. Table 1. The rank numbers of 𝗔 k ​ c subscript 𝗔 𝑘 𝑐 \bm{\mathsf{A}}_{kc} of the test geometry shown in Figure 4 . k 0 superscript 𝑘 0 k^{0} denotes the length of the index vector J 𝐽 J produced by the factorization technique in Algorithm LABEL:alg:A_kc prior to the recompression step, and k 𝑘 k denotes the length of the index vector J 𝐽 J after the recompression step. k opt subscript 𝑘 opt k_{\,\mbox{\footnotesize opt}} is the number of singular values of 𝗔 k ​ c subscript 𝗔 𝑘 𝑐 \bm{\mathsf{A}}_{kc} that are greater than ϵ = 1 × 10 − 10 italic-ϵ 1 superscript 10 10 \epsilon=1\times 10^{-10} and is considered as the optimal rank numbers for compression.

$N_k$	$N_c$	$k^0$	$k$	$k_{\text{opt}}$
1200	80	245	17	15
4800	320	319	17	15
19200	1280	361	17	15

Table 2. Table 2. Timing results for the square with thinning nose geometry in Section 5.1 .

$N_o$	$T_{new,p}$	$T_{hbs,p}$	$r_p$	$T_{new,s}$	$T_{hbs,s}$	$r_s$
1168	1.68e-01	5.12e-01	3.28e-01	7.38e-03	7.54e-03	9.79e-01
2320	1.69e-01	6.13e-01	2.75e-01	1.13e-02	7.18e-03	1.57e+00
4624	2.36e-01	9.24e-01	2.56e-01	1.50e-02	1.11e-02	1.35e+00
9232	3.34e-01	1.37e+00	2.44e-01	2.06e-02	1.65e-02	1.25e+00
18448	5.47e-01	2.20e+00	2.49e-01	3.46e-02	2.82e-02	1.23e+00
36880	1.10e+00	3.76e+00	2.93e-01	6.18e-02	4.63e-02	1.34e+00
73744	1.98e+00	6.88e+00	2.87e-01	1.26e-01	8.96e-02	1.41e+00
147472	3.95e+00	1.32e+01	2.99e-01	2.37e-01	1.71e-01	1.39e+00

Table 3. Table 3. Timing results for the square with fixed nose geometry in Section 5.1 .

$N_o$	$T_{new,p}$	$T_{hbs,p}$	$r_p$	$T_{new,s}$	$T_{hbs,s}$	$r_s$
1168	1.55e-01	3.82e-01	4.05e-01	8.82e-03	7.36e-03	1.20e+00
2336	2.05e-01	5.27e-01	3.89e-01	1.34e-02	6.37e-03	2.10e+00
4672	3.09e-01	7.83e-01	3.95e-01	1.89e-02	1.00e-02	1.89e+00
9344	5.54e-01	1.15e+00	4.82e-01	2.40e-02	1.55e-02	1.55e+00
18688	1.00e+00	1.88e+00	5.34e-01	3.96e-02	2.54e-02	1.56e+00
37376	2.56e+00	3.53e+00	7.24e-01	7.46e-02	4.90e-02	1.52e+00
74752	7.51e+00	6.55e+00	1.15e+00	1.68e-01	9.20e-02	1.82e+00
149504	3.277e+01	1.290e+01	2.541e+00	7.038e-01	1.858e-01	3.787e+00

Table 4. Table 4. Timing results for the circle with shrinking bump geometry in Section 5.2 .

$N_o$	$T_{new,p}$	$T_{hbs,p}$	$r_p$	$T_{new,s}$	$T_{hbs,s}$	$r_s$
2000	8.61e-02	1.61e-01	5.34e-01	7.77e-03	8.26e-03	9.42e-01
4000	1.26e-01	2.50e-01	5.04e-01	1.25e-02	9.81e-03	1.27e+00
8000	2.19e-01	4.80e-01	4.55e-01	2.31e-02	1.93e-02	1.20e+00
16000	4.11e-01	9.41e-01	4.37e-01	4.38e-02	3.80e-02	1.15e+00
32000	8.31e-01	1.89e+00	4.39e-01	8.43e-02	7.71e-02	1.09e+00
64000	1.67e+00	3.78e+00	4.42e-01	1.71e-01	1.59e-01	1.08e+00
128000	3.43e+00	7.54e+00	4.55e-01	3.51e-01	3.11e-01	1.13e+00
256000	7.10e+00	1.51e+01	4.69e-01	6.80e-01	6.36e-01	1.07e+00

Table 5. Table 5. Timing results for the circle with fixed bump geometry in Section 5.2 .

$N_o$	$N_p$	$T_{new,p}$	$T_{hbs,p}$	$r_p$	$T_{new,s}$	$T_{hbs,s}$	$r_s$
1000	49	8.46e-02	1.48e-01	5.71e-01	3.96e-03	6.14e-03	6.46e-01
2000	99	5.74e-02	1.42e-01	4.04e-01	6.37e-03	4.83e-03	1.32e+00
4000	199	1.27e-01	2.60e-01	4.90e-01	1.16e-02	1.00e-02	1.16e+00
8000	399	3.40e-01	4.81e-01	7.08e-01	2.50e-02	1.90e-02	1.32e+00
16000	799	1.15e+00	9.05e-01	1.27e+00	6.04e-02	3.86e-02	1.57e+00
32000	1599	4.16e+00	1.79e+00	2.32e+00	1.61e-01	7.68e-02	2.10e+00
64000	3199	1.93e+01	3.64e+00	5.31e+00	6.38e-01	1.56e-01	4.10e+00
128000	6399	1.63e+02	7.43e+00	2.20e+01	1.69e+01	3.18e-01	5.31e+01

Table 6. Table 6. Timing results for the star with refined panels geometry in Section 5.3 .

$N_p$	$\frac{N_p}{N_o}$	$T_{new,p}$	$T_{hbs,p}$	$r_p$	$T_{new,s}$	$T_{hbs,s}$	$r_s$
96	0.03	1.34e-01	7.44e-01	1.80e-01	2.33e-02	2.49e-02	9.35e-01
192	0.06	1.35e-01	6.65e-01	2.02e-01	2.34e-02	2.45e-02	9.55e-01
384	0.12	1.67e-01	6.73e-01	2.48e-01	2.42e-02	2.48e-02	9.77e-01
768	0.24	2.34e-01	7.18e-01	3.25e-01	2.76e-02	2.55e-02	1.08e+00
1536	0.48	4.80e-01	8.33e-01	5.76e-01	5.87e-02	3.29e-02	1.78e+00
3072	0.96	1.29e+00	1.06e+00	1.22e+00	2.44e-01	3.94e-02	6.19e+00
6144	1.92	4.26e+00	1.63e+00	2.62e+00	1.70e+00	6.25e-02	2.72e+01
12288	3.84	2.79e+01	2.38e+00	1.17e+01	1.41e+01	9.81e-02	1.44e+02

Equations

$$\begin{array}{rl}-\Delta u(\mathsf{x})=& 0,\mathsf{x}\in \Omega ,\\ u(\mathsf{x})=& g(\mathsf{x}),\mathsf{x}\in \Gamma .\end{array}$$

$$u(\mathsf{x})={\int }_{\Gamma }D(\mathsf{x},\mathsf{y})\sigma (\mathsf{y})ds(\mathsf{y}),\mathsf{x}\in \Omega ,$$

$$-\frac{1}{2}\sigma (\mathsf{x})+{\int }_{\Gamma }D(\mathsf{x},\mathsf{y})\sigma (\mathsf{y})ds(\mathsf{y})=g(\mathsf{x}).$$

$$\mathsf{A}\sigma =\mathsf{g}$$

$$\left(\underbrace{\left[\begin{array}[]{cc}\bm{\mathsf{A}}_{oo}&\bm{\mathsf{0}}\\ \bm{\mathsf{0}}&\bm{\mathsf{A}}_{pp}\end{array}\right]}_{\tilde{\bm{\mathsf{A}}}}+\underbrace{\begin{bmatrix}0&\begin{pmatrix}-\bm{\mathsf{A}}_{kc}\\ -\bm{\mathsf{B}}_{cc}\end{pmatrix}&\bm{\mathsf{A}}_{op}\\ \bm{\mathsf{A}}_{pk}&0&0\end{bmatrix}}_{\bm{\mathsf{Q}}}\right)\begin{pmatrix}\bm{\mathsf{\sigma}}_{k}\\ \bm{\mathsf{\sigma}}_{c}\\ \bm{\mathsf{\sigma}}_{p}\end{pmatrix}=\underbrace{\begin{pmatrix}\bm{\mathsf{f}}_{k}\\ \bm{\mathsf{0}}\\ \bm{\mathsf{f}}_{p}\end{pmatrix}}_{\bm{\mathsf{f}}_{\rm ext}}$$

$$\bm{\mathsf{M}}\sim\left[\begin{array}[]{ccccc}\bm{\mathsf{D}}_{1}&\bm{\mathsf{M}}_{1,2}&\bm{\mathsf{M}}_{1,3}&\cdots&\bm{\mathsf{M}}_{1,p}\\ \bm{\mathsf{M}}_{2,1}&\bm{\mathsf{D}}_{2}&\bm{\mathsf{M}}_{2,3}&\cdots&\bm{\mathsf{M}}_{2,p}\\ \vdots&\vdots&\vdots&&\vdots\\ \bm{\mathsf{M}}_{p,1}&\bm{\mathsf{M}}_{p,2}&\bm{\mathsf{M}}_{p,3}&\cdots&\bm{\mathsf{D}}_{p}\end{array}\right]$$

$$\Vert {\mathsf{U}}_{\sigma }{\tilde{\mathsf{M}}}_{\sigma ,\tau }{\mathsf{V}}_{\tau }^{\ast }-{\mathsf{M}}_{\sigma ,\tau }\Vert \le ϵ\sigma ,\tau \in \{1,2,\dots ,p\},\sigma \ne \tau .$$

$$\mathsf{M}\sim \mathsf{U}\tilde{\mathsf{M}}{\mathsf{V}}^{\ast }+\mathsf{D},$$

$$\mathsf{U}=\text{mbox}diag({\mathsf{U}}_1,{\mathsf{U}}_2,\dots ,{\mathsf{U}}_p),\mathsf{V}=\text{mbox}diag({\mathsf{V}}_1,{\mathsf{V}}_2,\dots ,{\mathsf{V}}_p),\mathsf{D}=\text{mbox}diag({\mathsf{D}}_1,{\mathsf{D}}_2,\dots ,{\mathsf{D}}_p),$$

$$\tilde{\bm{\mathsf{M}}}=\left[\begin{array}[]{cccc}0&\tilde{\bm{\mathsf{M}}}_{12}&\tilde{\bm{\mathsf{M}}}_{13}&\cdots\\ \tilde{\bm{\mathsf{M}}}_{21}&0&\tilde{\bm{\mathsf{M}}}_{23}&\cdots\\ \tilde{\bm{\mathsf{M}}}_{31}&\tilde{\bm{\mathsf{M}}}_{32}&0&\cdots\\ \vdots&\vdots&\vdots\end{array}\right].$$

$$\bm{\mathsf{M}}\sim\bm{\mathsf{U}}^{(3)}\bigl{(}\bm{\mathsf{U}}^{(2)}\bigl{(}\bm{\mathsf{U}}^{(1)}\,\bm{\mathsf{B}}^{(0)}\,(\bm{\mathsf{V}}^{(1)})^{*}+\bm{\mathsf{B}}^{(1)}\bigr{)}(\bm{\mathsf{V}}^{(2)})^{*}+\bm{\mathsf{B}}^{(2)}\bigr{)}(\bm{\mathsf{V}}^{(3)})^{*}+\bm{\mathsf{D}}^{(3)},$$

$$\begin{array}[]{cccccccccccccccccc}\bm{\mathsf{U}}^{(3)}&\bm{\mathsf{U}}^{(2)}&\bm{\mathsf{U}}^{(1)}&\bm{\mathsf{B}}^{(0)}&(\bm{\mathsf{V}}^{(1)})^{*}&\bm{\mathsf{B}}^{(1)}&(\bm{\mathsf{V}}^{(2)})^{*}&\bm{\mathsf{B}}^{(2)}&(\bm{\mathsf{V}}^{(3)})^{*}&\bm{\mathsf{D}}^{(3)}.\\ \includegraphics[scale={0.3}]{ffig_U3.eps}&\includegraphics[scale={0.3}]{ffig_U2.eps}&\includegraphics[scale={0.3}]{ffig_U1.eps}&\includegraphics[scale={0.3}]{ffig_B0.eps}&\includegraphics[scale={0.3}]{ffig_V1.eps}&\includegraphics[scale={0.3}]{ffig_B1.eps}&\includegraphics[scale={0.3}]{ffig_V2.eps}&\includegraphics[scale={0.3}]{ffig_B2.eps}&\includegraphics[scale={0.3}]{ffig_V3.eps}&\includegraphics[scale={0.3}]{ffig_D3.eps}&\end{array}$$

$$[\mathsf{P},J]=\mathrm{id}(\hat{\mathsf{A}})$$

$$\mathsf{W}=\mathsf{P}\mathsf{W}(J(1:l),:)$$

$${\left(\mathsf{A}+\mathsf{L}\mathsf{R}\right)}^{-1}={\mathsf{A}}^{-1}+{\mathsf{A}}^{-1}\mathsf{L}{\left(\mathsf{I}+\mathsf{R}{\mathsf{A}}^{-1}\mathsf{L}\right)}^{-1}\mathsf{R}{\mathsf{A}}^{-1}$$

$$\begin{array}[]{cccccccccc}{\mathbf{A}}_{kc}&\approx&{\mathbf{L}}_{kc}&{\mathbf{R}}_{kc},&&&{\mathbf{A}}_{op}&\approx&{\mathbf{L}}_{op}&{\mathbf{R}}_{op},\mbox{ and }\\ N_{k}\times N_{c}&&N_{k}\times k_{kc}&k_{kc}\times N_{c}&&&N_{o}\times N_{p}&&N_{o}\times k_{op}&k_{op}\times N_{p}\\ &&&{\mathbf{A}}_{pk}&\approx&{\mathbf{L}}_{pk}&{\mathbf{R}}_{pk}&&&\\ &&&N_{p}\times N_{k}&&N_{p}\times k_{pk}&k_{pk}\times N_{k}&&&\\ \end{array}$$

$$\begin{array}[]{cccc}\mathbf{Q}&=&\mathbf{L}&\mathbf{R}\\ N_{\rm ext}\times N_{\rm ext}&&N_{\rm ext}\times k&k\times N_{\rm ext}\end{array}$$

$$\mathbf{L}=\left[\begin{array}{cc}\left(\begin{array}{cc}-{\mathbf{L}}_{kc}& 0\\ 0& -{\mathbf{B}}_{cc}\end{array}\right)& {\mathbf{L}}_{op}\\ {\mathbf{L}}_{pk}& 0\end{array}\right]\text{mbox}and\mathbf{R}=\left[\begin{array}{ccc}{\mathbf{R}}_{pk}& 0& 0\\ 0& \left(\begin{array}{c}{\mathbf{R}}_{kc}\\ {\mathbf{I}}_{cc}\end{array}\right)& 0\\ 0& 0& {\mathbf{R}}_{op}\end{array}\right].$$

$$\bm{\mathsf{U}}_{\tau}=\left[\begin{array}[]{cc}\bm{\mathsf{U}}_{\sigma_{1}}&\bm{\mathsf{0}}\\ \bm{\mathsf{0}}&\bm{\mathsf{U}}_{\sigma_{2}}\end{array}\right]\bm{\mathsf{P}}_{\tau}.$$

$${\mathsf{A}}^{-1}\mathsf{L}$$

$$=\left[\begin{array}{cc}-{\mathbf{A}}_{oo}^{-1}\left(\begin{array}{cc}{\mathbf{L}}_{kc}& 0\\ 0& {\mathbf{B}}_{cc}\end{array}\right)& {\mathbf{A}}_{oo}^{-1}{\mathbf{L}}_{op}\\ {\mathbf{A}}_{pp}^{-1}{\mathbf{L}}_{pk}& 0\end{array}\right].$$

$$E=\frac{\Vert {\mathbf{u}}_{exact}-{\mathbf{u}}_{new}{\Vert }_2}{\Vert {\mathbf{u}}_{exact}{\Vert }_2},$$

Full text

A fast direct solver for boundary value problems on locally perturbed geometries

*Y. Zhang, and A. Gillman

Department of Computational and Applied Mathematics, Rice University

Yabin.Zhang@rice.edu, adrianna.gillman@rice.edu*

Abstract: Many applications involve solving several boundary value problems on geometries that are local perturbations of an original geometry. The boundary integral equation for a problem on a locally perturbed geometry can be expressed as a low rank update to the original system. A fast direct solver for the new linear system is presented in this paper. The solution technique utilizes a precomputed fast direct solver for the original geometry to efficiently create the low rank factorization of the update matrix and to accelerate the application of the Sherman-Morrison formula. The method is ideally suited for problems where the local perturbation is the same but its placement on the boundary changes and problems where the local perturbation is a refined discretization on the same geometry. Numerical results illustrate that for fixed local perturbation the method is three times faster than building a new fast direct solver from scratch.

1. Introduction

This manuscript presents a fast direct solver for boundary integral equations where the geometry for each problem corresponds to a local perturbation of the original geometry. In particular, we are interested in problems where the local perturbation to the geometry is much smaller than the original geometry. Since a direct solver is constructed, the technique is ideal for problems with many right hand sides and/or suffer from ill-conditioning due to geometric complexity. Boundary value problems involving locally perturbed geometries arise in a variety of applications such as optimal design [19], and adaptive discretization techniques [5]. For example, finding the optimal placement of an attachment to a large geometry which minimizes the radar cross section involves solving many problems where the local perturbation is the same but the placement on the boundary changes.

For many boundary integral equations, the linear system that results from the discretization of an integral equation is amenable to fast direct solvers such as those built from hierarchically semiseparable (HSS) [20, 4], $\mathcal{H}$-matrix [12], hierarchically block separable (HBS) [7], hierarchical interpolatory factorization (HIF) [14] and hierarchical off-diagonal low rank (HODLR) [1] representations. These direct solvers utilize the fact that the off-diagonal blocks of the dense matrix are low rank. The different variants correspond to different ways of exploiting this property. Let $\bm{\mathsf{A}}$ denote the matrix resulting from discretization of the boundary integral equation. The factored approximation of the matrix $\bm{\mathsf{A}}$, denoted by $\tilde{\bm{\mathsf{A}}}$, is constructed so that $\|\bm{\mathsf{A}}-\tilde{\bm{\mathsf{A}}}\|\leq\epsilon$ for a user defined tolerance $\epsilon$. $\tilde{\bm{\mathsf{A}}}$ is called the compressed representation of $\bm{\mathsf{A}}$. The inverse of the compressed matrix is then constructed via a variant of a Woodbury formula or by expanding the matrix out to a larger sparse system and using a sparse direct solver. We refer the reader to the references for further details.

Building from the approach in [9], the solution technique presented in this paper casts the linear system for problem on the perturbed geometry as an extended linear system which consists of a two-by-two block diagonal matrix plus a low rank update. The block diagonal matrix has a block equal to the matrix for the original geometry. By using the Sherman-Morrison formula, the approximate inverse of the original system can be exploited and the approximate inverse of the extended system can be applied rapidly. The compressed representation matrix in the original system is utilized to reduce the cost of computing the low rank factorization of the update matrix. The techniques used in this work draw from earlier work in [16, 2, 17, 21].

The method presented in this paper is ideally suited for problems where the local perturbation is the same over many placements on the boundary or the number of removed points on the boundary is not large. Since the solution technique does not modify the original compressed representation, it can be combined with any fast direct solver.

1.1. Related work

The paper [18] presents a technique for updating the Hierarchical interpolative factorization (HIF) of the matrix $\bm{\mathsf{A}}$. This task involves locating and updating the relevant low rank factors and (potentially) modifying the underlying hierarchical tree structure. An approximate inverse is then constructed for the updated compressed representation. The inversion step is one of the most expensive steps in the precomputation of a fast direct solver.

1.2. Outline of paper

This manuscript begins by presenting the boundary integral formulation, the discretized linear system and the extended linear system for a model problem in section 2. The fast direct solver presented in this paper for a locally perturbed geometry utilizes the factors computed in the fast direct solver for the original geometry. While the method can be utilized in conjunction with any fast direct solver, we chose to review the HBS method and its physical interpretation in section 3 for simplicity of presentation. Next, the construction of the new fast direct solver is presented in section 4. Finally numerical experiments report on the performance of the solver in section 5. Section 6 reviews the methods and highlights the potential of the solution technique.

2. Model problem

This section begins by reviewing the boundary integral approach for solving a Laplace boundary value problem. Then the technique for writing the linear system corresponding to the discretized boundary value problem on a geometry that is a local perturbation of the original is presented in section 2.2.

2.1. Boundary integral equations

For simplicity of presentation, we consider the Laplace boundary value problem

[TABLE]

Figure 1(a) illustrates a sample geometry. The vector $\bm{\mathsf{\nu}}_{\bm{\mathsf{y}}}$ denotes the outward normal vector at the point $\bm{\mathsf{y}}\in\Gamma$. For $\bm{\mathsf{x}}\in\Omega$, we represent the solution to (1) as a double layer potential

[TABLE]

where $D(\bm{\mathsf{x}},\bm{\mathsf{y}})=\partial_{\bm{\mathsf{\nu}}_{\bm{\mathsf{y}}}}G(\bm{\mathsf{x}},\bm{\mathsf{y}})$ is double layer kernel, $G(\bm{\mathsf{x}},\bm{\mathsf{y}})=-\frac{1}{2\pi}\log|\bm{\mathsf{x}}-\bm{\mathsf{y}}|$ is the fundamental solution and $\sigma(\bm{\mathsf{x}})$ is the unknown boundary charge distribution. By taking the limit of $u(\bm{\mathsf{x}})$ as $\bm{\mathsf{x}}$ goes to the boundary and setting it equal to $g(\bm{\mathsf{x}})$, we find the boundary charge distribution $\sigma(\bm{\mathsf{x}})$ satisfies the following boundary integral equation

[TABLE]

Discretization of the boundary integral equation (3) with either a Nyström or a boundary element method results in a linear system of the form

[TABLE]

where the solution $\bm{\mathsf{\sigma}}$ is the approximation of $\sigma(\bm{\mathsf{x}})$ at the discretization points on $\Gamma$.

2.2. An extended linear system for local perturbations

In this section, we consider the Laplace boundary value problem (1) on the geometry illustrated in Figure 1(b) where $\Gamma=\Gamma_{k}\cup\Gamma_{p}$, $\Gamma_{k}$ is the remaining portion of the original boundary, $\Gamma_{c}$ is the removed portion of the original boundary and $\Gamma_{p}$ is the newly added portion of the boundary. We define the boundary of the original geometry by $\Gamma_{o}=\Gamma_{k}\cup\Gamma_{c}$.

The discretized problem on $\Gamma$ can be expressed as an extended linear system [9] by

[TABLE]

where $\bm{\mathsf{A}}_{kc}$ denotes the submatrix of $\bm{\mathsf{A}}_{oo}$ corresponding to the interaction between $\Gamma_{k}$ and $\Gamma_{c}$, $\bm{\mathsf{A}}_{op}$ denotes the discretization of the double layer integral operator on $\Gamma_{p}$ evaluated on $\Gamma_{o}$, $\bm{\mathsf{A}}_{pk}$ denotes the discretization of the double layer integral operator on $\Gamma_{k}$ evaluated on $\Gamma_{p}$, and $\bm{\mathsf{B}}_{cc}$ denotes the sub-matrix of $\bm{\mathsf{A}}_{oo}$ corresponding to the interaction of $\Gamma_{c}$ with itself but the diagonal entries are set to zero. The matrix $\bm{\mathsf{Q}}$ is called the update matrix.

For many choices of $\Gamma_{p}$ and $\Gamma_{c}$, the update matrix $\bm{\mathsf{Q}}$ is low rank. These are the cases where section 4 presents a technique for efficiently creating the low rank factorization of $\bm{\mathsf{Q}}$ and how to exploit the factorization to make a fast direct solver.

3. A fast direct solver for boundary integral equations

In order to construct the low rank factorization of $\bm{\mathsf{Q}}$ as efficiently as possible, we reuse factors in the compressed representation of $\bm{\mathsf{A}}_{oo}$. To give the new work context and introduce vocabulary, this section presents a brief review of the construction of HBS representation of the matrix $\bm{\mathsf{A}}$ in equation (4). More details and the inversion technique are presented in [7] and [3]. Other direct solution techniques, such as $\mathcal{H}$-matrix, HSS, etc., use similar factorization techniques and can be coupled to the new solver.

Fast direct solvers for the linear system in equation (4) begin by creating a data-sparse representation which approximates $\bm{\mathsf{A}}$. Roughly speaking, a data-sparse representation of a matrix is a factorization which requires $O(N)$ memory to store where $N$ is the number of discretization points. For many fast direct solvers including the HBS method, the reduction in memory is achieved by exploiting low rank approximations of off-diagonal blocks. These data-sparse representations also yield fast matrix vector multiplication and fast inversion schemes.

This section begins by briefly reviewing the Hierarchically block separable (HBS) representation of a dense matrix in section 3.1. Section 3.2 reviews the physical interpretation of the hierarchical method. Then, section 3.3 presents a fast technique for creating the low rank approximations of off-diagonal blocks.

3.1. HBS representation

This section reviews the HBS representation of a matrix $\bm{\mathsf{M}}$

Consider the $np\times np$ block partition of a matrix $\bm{\mathsf{M}}$ into $p\times p$ blocks each of size $n\times n$:

[TABLE]

Given a desired tolerance $\epsilon$, for each $\tau=1,\,2,\,\dots,\,p$, there exists a constant $k_{\tau}$ and $n\times k_{\tau}$ matrices $\bm{\mathsf{U}}_{\tau}$ and $\bm{\mathsf{V}}_{\tau}$ such that each off-diagonal block $\bm{\mathsf{M}}_{\sigma,\tau}$ of $\bm{\mathsf{M}}$ admits an approximate factorization

[TABLE]

The columns of $\bm{\mathsf{U}}_{\sigma}$ form a column basis for the columns of all off-diagonal blocks in row $\sigma$. Likewise, the columns of $\bm{\mathsf{V}}_{\tau}$ form a row basis for all the rows of all off-diagonal blocks in columns $\tau$.

This factorization allows $\bm{\mathsf{M}}$ to be approximated in the following factored form

[TABLE]

where

[TABLE]

and

[TABLE]

This is a one level factorization of $\bm{\mathsf{M}}$ which can be inverted via a variation of the Sherman-Morrison-Woodbury formula (see Lemma 3.1 in [7] or [3]).

When the matrix $\tilde{\bm{\mathsf{M}}}$ can be factored in the same manner, the matrix is called Hierarchically block separable (HBS). A three level factorization of this kind is expressed as

[TABLE]

where the block structure of the factors is

[TABLE]

3.2. Discretized boundary integral equation to HBS form

The HBS representation is based on a binary tree partitioning of an index vector $I=[1,\ldots,N]$. For simplicity, we present the technique with a uniform binary tree. The root of the tree is $I_{1}=I$. In the next level of the tree, the index vector $I$ is split into two equilength index vectors $I_{2}$ and $I_{3}$. This process is repeated until each index vector has less than some preset number $n$ of entries. A leaf node in the tree is an index vector that is not split. A non-leaf node $\tau$ has children $\sigma_{1}$ and $\sigma_{2}$. The nodes $\tau$ is the parent of $\sigma_{1}$ and $\sigma_{2}$ if $I_{\tau}=I_{\sigma_{1}}\cup I_{\sigma_{2}}$. Figure 2 illustrates a three level binary tree where $N=400$ and $n=50$.

3.3. Efficient construction of low rank factorizations of off-diagonal blocks

Constructing the HBS representation of $\bm{\mathsf{M}}$ via pure linear algebraic techniques would result in an $O(N^{2})$ compression scheme. When the matrix $\bm{\mathsf{M}}$ results from the discretization of a boundary integral equation, physics can be exploited to reduce the computational cost of compression to linear.

First, we note that the index vectors in the binary tree structure have a physical interpretation. Since each index corresponds to a point on $\Gamma$, an index vector corresponds to a collection of discretization points on $\Gamma$ which can be thought of as segment(s) on $\Gamma$.

For presentation purposes, consider the task of creating the low rank factorization of the sub-matrix of the discretized linear system (4) corresponding to the interaction of $\Gamma_{\tau}$ and $\Gamma_{\tau}^{c}=\Gamma/\Gamma_{\tau}$ (see Figure 3(a)). Let $\bm{\mathsf{A}}_{\tau,c}$ denote this matrix.

Instead of factoring $\bm{\mathsf{A}}_{\tau,c}$, we partition $\Gamma_{\tau}^{c}$ into the “near” and “far” portions. The near portion of $\Gamma_{\tau}^{c}$, denoted by $\Gamma_{\tau}^{\rm near}$, lies inside a proxy surface denoted by $\Gamma^{\rm proxy}_{\tau}$. In this work, we take $\Gamma^{\rm proxy}_{\tau}$ to be a circle with radius $r\sim 1.5r_{\tau}$, where $r_{\tau}$ is the radius of $\Gamma_{\tau}$, concentric with $\Gamma_{\tau}$. Figure 3(b) illustrates the proxy surface $\Gamma_{\tau}^{\rm proxy}$, proxy points, and near points $\Gamma_{\tau}^{\rm near}$ for $\Gamma_{\tau}$. The portion of $\Gamma_{\tau}^{c}$ outside of $\Gamma_{\tau}^{\rm proxy}$ is the far portion of the boundary, denoted by $\Gamma_{\tau}^{\rm far}$.

From potential theory ideas similar to those employed in the fast multipole method [10], it is known that the evaluation of a smooth kernel such as (2) for points that are far from each other can be expressed with a small number of basis functions to arbitrary accuracy. The basis functions we use are a collection of single poles lying on $\Gamma^{\rm proxy}_{\tau}$. This collection of points placed on the proxy surface are called proxy points. For the examples considered in this paper, we found that it is enough to have $75$ proxy points. Let $\bm{\mathsf{A}}_{\rm proxy}$ denote a matrix characterizing the interaction between the discretization points on $\Gamma_{\tau}$ and the proxy points. Let $\bm{\mathsf{A}}_{\tau,\rm near}$ denote the sub-matrix of $\bm{\mathsf{A}}_{\tau,c}$ corresponding to the interaction between $\Gamma_{\tau}$ and the near points. Then, we compute the low rank factorization of $\hat{\bm{\mathsf{A}}}=\left[\bm{\mathsf{A}}_{\tau,\rm near}|\bm{\mathsf{A}}_{\rm proxy}\right]$ which has dimension $N_{\tau}\times(N_{\rm near}+N_{\rm proxy})$ where $N_{\rm near}$ denotes the number of near points, and $N_{\rm proxy}$ denotes the number of proxy points.

Instead of computing a QR or SVD, we use an interpolatory decomposition defined in definition 3.1. The underlying algorithms for computing such a factorization efficiently include rank-revealing QR [11] and randomized sampling [15, 13]. Since one of the factors is a sub-matrix of the matrix being factored, the potential theory associated with the boundary value problem extends to the factor. Thus the proxy surface technique can be applied recursively to create factorizations corresponding to unions of intervals. Let

[TABLE]

denote the process of computing the interpolatory decomposition where $\bm{\mathsf{P}}$ is the interpolation matrix, and $J$ is the corresponding index vector. Algorithm LABEL:alg:compress provides a pseudocode for computing the low rank factorization of $\bm{\mathsf{A}}_{\tau,c}$ described in this section.

Definition 3.1**.**

The interpolatory decomposition of a $m\times n$ matrix $\bm{\mathsf{W}}$ that has rank $l$ is the factorization

[TABLE]

where $J$ is a vector of integers $j_{i}$ such $1\leq j_{i}\leq m$, and $\bm{\mathsf{P}}$ is a $m\times l$ matrix that contains a $l\times l$ identity matrix. Namely, $\bm{\mathsf{P}}(J(1:l),:)=\bm{\mathsf{I}}_{l}$.

4. The fast direct solver for the extended linear system

This section presents the fast direct solver for extended linear system. Recall the matrix $\bm{\mathsf{Q}}$ has subblocks $\bm{\mathsf{A}}_{kc}$, $\bm{\mathsf{A}}_{op}$ and $\bm{\mathsf{A}}_{pk}$ which correspond to interactions between small portions of the boundary and the remainder or portions of the remainder of the boundary. From section 3.3, we know these matrices are low rank. Let $\bm{\mathsf{Q}}\sim\bm{\mathsf{L}}\bm{\mathsf{R}}$, where $\bm{\mathsf{L}}\in\mathbb{R}^{n\times k}$ and $\bm{\mathsf{R}}\in\mathbb{R}^{k\times n}$, be the rank $k$ approximation of $\bm{\mathsf{Q}}$. Then the inverse of $\bm{\mathsf{A}}+\bm{\mathsf{Q}}$ can be approximated by the following Sherman-Morrison formula [8]

[TABLE]

Recall $\bm{\mathsf{A}}$ is a block diagonal matrix with block $\bm{\mathsf{A}}_{oo}$ and $\bm{\mathsf{A}}_{pp}$. Since $\bm{\mathsf{A}}_{oo}^{-1}$ is already approximated by a fast direct solver and the size of $\bm{\mathsf{A}}_{pp}$ is small for applications of interest, the application of the approximate inverse of $\bm{\mathsf{A}}$ is fast.

Remark 1**.**

If $\Gamma_{p}$ remains the same for multiple perturbed geometries, the cost of constructing $\bm{\mathsf{A}}^{-1}_{pp}$ is not included in the precomputation. In applications where $N_{p}$ is large, an approximation of $\bm{\mathsf{A}}^{-1}_{pp}$ can be constructed via a fast direct solver.

For any local perturbation, the matrices $\bm{\mathsf{L}}$, $\bm{\mathsf{R}}$, $\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}$ and $\left(\bm{\mathsf{I}}+\bm{\mathsf{R}}\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}\right)^{-1}$ need only be computed once. The construction of the factorizations and $\left(\bm{\mathsf{I}}+\bm{\mathsf{R}}\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}\right)^{-1}$ comprise the precomputation of the new solver. Once constructed, the application of the Woodbury formula (10) can be evaluated for linear cost with small constant.

In order for this to be a fast direct solver, the low rank factorization of $\bm{\mathsf{Q}}$ must scale linearly with the number of discretization points on the original geometry. Before detailing how to efficiently factorize $\bm{\mathsf{Q}}$, we introduce some notation. Let $N_{k}$ denote the number of discretization points on $\Gamma_{k}$, $N_{c}$ denote the number of discretization points on $\Gamma_{c}$, and $N_{p}$ denote the number of discretization points on $\Gamma_{p}$. Then the number of discretization points on the original boundary $\Gamma_{o}$ is $N_{o}=N_{c}+N_{k}$, the number of discretization points on the new geometry is $N_{n}=N_{k}+N_{p}$, and the dimension of $\bm{\mathsf{Q}}$ is $N_{\rm ext}\times N_{\rm ext}$ where $N_{\rm ext}=N_{o}+N_{p}$.

The remainder of this section describes how to efficiently construct the low rank factorization of $\bm{\mathsf{Q}}$. The factorization is achieved by constructing low rank factorizations of the submatrices. Let

[TABLE]

denote the low rank factorizations of the submatrices. Then

[TABLE]

where

[TABLE]

Section 4.1 and section 4.2 present a linearly scaling technique for factorizing $\bm{\mathsf{A}}_{kc}$ and $\bm{\mathsf{A}}_{op}$ respectively. The low rank factorization of $\bm{\mathsf{A}}_{pk}$ is achieved via a technique similar to that as factoring $\bm{\mathsf{A}}_{op}$ and thus is not presented.

4.1. The efficient factorization of $\bm{\mathsf{A}}_{kc}$

The matrix $\bm{\mathsf{A}}_{kc}$ is a sub-matrix of $\bm{\mathsf{A}}_{oo}$ and thus much of the low rank factorization can be extracted from the compressed representation of $\bm{\mathsf{A}}_{oo}$. This section presents a technique for efficiently creating a low rank factorization of $\bm{\mathsf{A}}_{kc}$ by reusing as much information as possible from the compressed representation of $\bm{\mathsf{A}}_{oo}$. Roughly speaking, many of the factors can be found by collecting the $\bm{\mathsf{U}}$ or $\bm{\mathsf{V}}$ factors in the HBS factorization (such as (8)). For simplicity of presentation, this section is limited to collecting the $\bm{\mathsf{U}}$ factors.

The algorithm begins by first constructing the low rank factorizations that cannot be extracted from the HBS representation. This consists of leaf boxes $\tau$ that have points in $\Gamma_{k}$ and $\Gamma_{c}$. Let $W$ denote the set of all such leaf boxes. For each $\tau\in W$, the points on $\tau\in\Gamma_{k}$ are identified and labeled $J_{\tau,k}$ and the matrix $\bm{\mathsf{A}}_{kc}(J_{\tau,k},:)$ is compressed via the method presented in section 3.3. The result is the interpolation matrix $\bm{\mathsf{P}}_{\tau}$ and index vector $J_{\tau}$.

Next we extract as much of the information from the HBS representation of the matrix $\bm{\mathsf{A}}_{oo}$ as possible. Let $T$ denote the set of boxes $\tau\subset\Gamma_{k}$. Information is extracted by going through the boxes in $T$ from smallest to the largest. For a leaf box $\tau\in T$, the interpolation matrix $\bm{\mathsf{U}}_{\tau}$ is extracted from the factorization. For a non-leaf box $\tau\in T$, let $\sigma_{1}$ and $\sigma_{2}$ denote the children boxes. Then let $\bm{\mathsf{P}}_{\tau}$ denote the interpolation matrix extracted from the HBS representation. Then

[TABLE]

Let $V=W\cup T=\{\tau_{1},\ldots,\tau_{m}\}$. A low rank factorization can result from letting $\bm{\mathsf{L}}_{kc}$ denote the block diagonal matrix where the subblocks are the matrices $\bm{\mathsf{U}}_{\tau_{j}}$ for $j=1,\ldots,m$ and $\bm{\mathsf{R}}_{kc}=\bm{\mathsf{A}}_{kc}(J,:)$ where ${J}=[{J}_{\tau_{1}},\ldots,{J}_{\tau_{m}}]$ denotes the corresponding list of indices. Unfortunately, the size of these factors is significantly larger than optimal. To prevent this from hampering the performance of the solver, we must do an additional compression step. We call this extra step recompression. Table 1 illustrates the approximate ranks computed via the different compression schemes for the geometry illustrated in Figure 4, the number of points $N_{k}$ on $\Gamma_{k}$, the number of points $N_{c}$ on $\Gamma_{c}$, the length $k^{0}$ of the index vector ${J}$, the size $k$ of the second dimension of $\bm{\mathsf{L}}_{kc}$ after recompression, and the optimal rank $k^{\rm opt}$ of $\bm{\mathsf{A}}_{kc}$. The size of the factors is close to optimal after the recompression step.

To recompress the factorization, we focus our attention on the submatrices of $\bm{\mathsf{A}}_{kc}$ which contain all the rank information. To construct the low rank factorization of $\bm{\mathsf{A}}_{kc}(J,:)$, we start by constructing the low rank factorization of $\bm{\mathsf{A}}_{kc}([{J}_{\tau_{1}},{J}_{\tau_{2}}],:)$ via the method in section 3.3. The result is an interpolation matrix $\bm{\mathsf{P}}_{12}$ and index vector ${J}_{12}$. We proceed by constructing the low rank factorization of $\bm{\mathsf{A}}_{kc}([{J}_{12},{J}_{\tau_{3}}],:)$ via the method in section 3.3. This continues until all $m$ boxes have been processed. The new matrix $\bm{\mathsf{L}}_{kc}$ is the old matrix multiplying by the interpolation matrices and the final index vector defines $\bm{\mathsf{R}}_{kc}$ as $\bm{\mathsf{R}}_{kc}=\bm{\mathsf{A}}_{kc}(J,:)$. Let $k_{kc}$ denote the length of the index vector $J$. Table 1 illustrates the near optimal rank $k_{kc}$ resulting from this procedure. Algorithm LABEL:alg:recompress provides a pseudocode for the recompression scheme.

Algorithm LABEL:alg:A_kc provides a pseudocode for the algorithm presented in this section.

4.2. The efficient compression of $\bm{\mathsf{A}}_{op}$

This section presents an efficient technique for constructing the low rank factorization of $\bm{\mathsf{A}}_{op}$. While this matrix is not a sub-matrix of $\bm{\mathsf{A}}_{oo}$, the far field information from the HBS representation of $\bm{\mathsf{A}}_{oo}$ can be reused in constructing the low rank factorization of $\bm{\mathsf{A}}_{op}$. Recall from section 3.3 that all far field interactions can be captured via the interaction with a proxy surface. This means we can reuse the interpolation matrices for all boxes that are far from $\Gamma_{p}$.

The technique for creating this factorization is based on partitioning $\Gamma_{o}$ into the portions near $\Gamma^{\rm near}$ and far $\Gamma^{\rm far}$ from $\Gamma_{p}$. For all boxes $\tau$ contained in $\Gamma^{\rm far}$, the technique for extracting the factors from the HBS representation of $\bm{\mathsf{A}}_{oo}$ is the same as in section 4.1. To compress near field interactions, the technique is similar to creating an HBS factorization from scratch.

The compression of the near field begins by creating a new binary tree (similar to the one in Figure 2) for the points $I_{\rm near}$ on $\Gamma^{\rm near}$. Let $\mathcal{T}_{near}$ denote the binary tree, then the factorization $\bm{\mathsf{A}}(I_{near},I_{p})$ can be constructed via a nested factorization. If a box $\tau$ is far from $\Gamma_{p}$, a proxy surface is used to create the low rank factorization of $\bm{\mathsf{A}}_{\tau,p}$. If a box $\tau$ is near $\Gamma_{p}$, a proxy surface $\Gamma^{\rm proxy}$ is placed around $\tau$ and any points in $\Gamma_{p}$ that are inside $\Gamma^{\rm proxy}$ are near $\tau$. Let $I_{p}$ denote the indices of those points. The interpolatory decomposition is then applied to $[\bm{\mathsf{A}}_{op}(\tau,I_{p}),\bm{\mathsf{A}}^{\rm proxy}]$ where $\bm{\mathsf{A}}^{\rm proxy}$ denotes the interaction between $\tau$ and the proxy surface. Let $\bm{\mathsf{P}}$ and $J$ denote the interpolation matrix and index vector respectively resulting from applying the interpolatory decomposition to $[\bm{\mathsf{A}}_{op}(\tau,I_{p}),\bm{\mathsf{A}}^{\rm proxy}]$. Algorithm LABEL:alg:geo provides a pseudocode for factorizing a near field interaction.

Let $\bm{\mathsf{L}}_{\rm far}$ and $\bm{\mathsf{R}}_{\rm far}=\bm{\mathsf{A}}_{op}(J_{\rm far},:)$ denote the low rank factorization of the far field where $J_{\rm far}$ denotes the index vector resulting from the far field compression. Likewise let $\bm{\mathsf{L}}_{\rm near}$ and $\bm{\mathsf{R}}_{\rm near}=\bm{\mathsf{A}}_{op}(J_{\rm near},:)$ denote the low rank factorization of the near field. A low rank factorization of $\bm{\mathsf{A}}_{op}$ can be $\bm{\mathsf{L}}_{op}\bm{\mathsf{R}}_{op}$ where \bm{\mathsf{L}}_{op}=\left[\begin{array}[]{cc}\bm{\mathsf{L}}_{\rm far}&0\\ 0&\bm{\mathsf{L}}_{\rm near}\end{array}\right] and \bm{\mathsf{R}}_{op}=\left[\begin{array}[]{cc}\bm{\mathsf{R}}_{\rm far}\\ \bm{\mathsf{R}}_{\rm near}\end{array}\right]. As seen in creating the low rank factorization of $\bm{\mathsf{A}}_{kc}$, the approximate rank given by this factorization is likely far from optimal. Thus another compression is necessary. Let $J_{\rm tot}=[J_{\rm far},J_{\rm near}]$. Applying the interpolatory decomposition to $\bm{\mathsf{A}}_{op}(J_{\rm tot},:)$ results in another interpolation matrix $\bm{\mathsf{P}}$ and index vector $J$. Thus the final factors are $\bm{\mathsf{L}}_{op}=\bm{\mathsf{L}}_{op}\bm{\mathsf{P}}$ and $\bm{\mathsf{R}}_{op}=\bm{\mathsf{A}}_{op}(J,:)$ where the approximate rank $k_{op}$ is the length of the index vector $J$.

Algorithm LABEL:alg:A_op gives a pseudocode for the technique presented in this section.

In the case where $N_{c}<<N_{k}$, the total compression cost of $\bm{\mathsf{A}}_{op}$ and $\bm{\mathsf{A}}_{kc}$ can be further reduced by combining the two and only performing one transversal of the binary tree on $\Gamma_{o}$.

4.3. Computational cost for the precomputation

Recall that the proposed direct solver is comprised of two steps: precomputation and solve. The precomputation step is more expensive than the solve step but it needs to only be computed once. As stated in the beginning of this section, the precomputation consists of computing the low rank factorization of the update matrix $\bm{\mathsf{Q}}$, $\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}$ and inverting $\left(\bm{\mathsf{I}}+\bm{\mathsf{R}}\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}\right)$. This section details the computational cost of the precomputation step.

The computational cost of constructing the low rank factorization of $\bm{\mathsf{Q}}$ is $O((N_{o}+N_{p})k_{Q})$, where $k_{Q}=\max\{k_{\rm HBS},k_{p}\}$, $k_{\rm HBS}$ denotes the rank of the HBS factors and $k_{p}$ denotes the rank of the $\bm{\mathsf{A}}_{op}$ factorization. The value of $k_{p}$ depends on the relationship of $\Gamma_{p}$ and $\Gamma_{o}$.

Let the low rank factors $\bm{\mathsf{L}}$ and $\bm{\mathsf{R}}^{T}$ of $\bm{\mathsf{Q}}$ have size $N\times k$ where $N=N_{o}+N_{p}$ and $k=k_{kc}+N_{c}+k_{pk}+k_{op}$. Constructing $\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}$ can be done in a block fashion by

[TABLE]

The computational cost of computing the upper left and right blocks is $O\left(N_{o}(k_{kc}+N_{c})\right)$ and $O(N_{o}k_{op})$, respectively. If $N_{p}$ is small, it is efficient to compute $\bm{\mathsf{A}}_{pp}^{-1}$ via dense linear algebra for $O(N_{p}^{3})$ cost. For $N_{p}$ large, an approximate inverse of $\bm{\mathsf{A}}_{pp}$ can be constructed via fast direct solver for $O(N_{p})$ computational cost. The computational cost of evaluating $\bm{\mathsf{A}}_{pp}^{-1}\bm{\mathsf{L}}_{pk}$ is $O(N_{p}^{2}k_{pk})$ via dense linear algebra and $O(N_{p}k_{pk})$ via fast linear algebra. Thus for $k_{op}$ small and $N_{c}$ constant, the cost of constructing $\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}$ is linear with respect to $N_{o}$.

For problems of interest, the matrix $\left(\bm{\mathsf{I}}+\bm{\mathsf{R}}\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}\right)$ is small enough to be inverted rapidly via dense linear algebra for $O(k^{3})$ computational cost.

4.4. Computational cost of the solve step

The solve step consists of applying the approximate inverse of $\bm{\mathsf{A}}$ to the the vector $\bm{\mathsf{f}}_{\rm ext}$ in equation (5), matrix vector multiplications and vector addition. As discussed in section 4.3, $\bm{\mathsf{A}}^{-1}$ is applied by blocks for a cost $O(N_{o}+N_{p}^{2})$. Applying $\bm{\mathsf{R}}$, $\left(\bm{\mathsf{I}}+\bm{\mathsf{R}}\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}\right)^{-1}$, and $\bm{\mathsf{A}}^{-1}\bm{\mathsf{L}}$ have computational cost $O((N_{o}+N_{p})k)$, $O(k^{2})$, and $O((N_{o}+N_{p})k)$. Thus the total computation cost of the solve step is $O((N_{o}+N_{p})k+k^{2}+N_{o}+N_{p}^{2})$. Thus the solve step is linear with respect to $N_{o}$.

5. Numerical experiments

This section illustrates the performance of the new fast direct solver for three types of locally-perturbed geometries. The geometries under consideration are

Square with a nose: The original boundary $\Gamma_{o}$ is a square with corners rounded via the method in [6]. The local perturbation $\Gamma_{p}$ is a rectangle with rounded corners that is attached to $\Gamma_{o}$ again using the method in [6] to smooth. Figure 6 illustrates the geometry. The length of the local perturbation rectangle is fixed but the height $d$ can vary depending on experiment. Composite Gaussian quadrature is used to discretize the integral equation.

Circle with a bump: The original geometry is a circle. The perturbed geometry replaces an arc of the circle by a smooth bump with central angle $\theta$. Figure 12 illustrates the geometry. This geometry is considered in [18]. Trapezoidal rule is used to discretize the integral equation.

Star with refined panels: The original geometry is a star geometry discretized with composite Gaussian quadrature (see Figure 9(a)). To create the perturbed boundary, three Gaussian panels illustrated in Figure 9(b) are replaced with more panels. Figure 9(b) illustrates a local perturbation $\Gamma_{p}$ consisting of six panels.

To ensure the solution technique preserves accuracy, we test it on problems with a known exact solution given by $u_{exact}(\vec{x})=\sum_{j=1}^{10}q_{j}G(\vec{x},\vec{s}_{j})$, where $\{\vec{s}_{j}\}_{j=1}^{10}$ are point charges placed outside of $\Omega$ and $\{q_{j}\}_{j=1}^{10}$ are the charge values. We define the the relative error to be

[TABLE]

where the vector $\bm{\mathsf{u}}_{new}$ and $\bm{\mathsf{u}}_{exact}$ contain the approximate and exact solution, respectively, at ten points $\{\vec{t}\}_{j=1}^{10}\in\Omega$. For all problems under consideration, the geometries are fully resolved. Thus, with the tolerance of the compression schemes set to $\epsilon=10^{-10}$, the relative error $E$ is approximately $10^{-9}$ for all choices of $N_{o}$ and $N_{p}$.

For each geometry, we report the following:

$N_{o}$ the number of discretization points on the original geometry.

$N_{p}$ the number of discretization points on the added geometry.

$T_{new,\,p}$ the time in seconds for the precomputation step for the new solver.

$T_{hbs,\,p}$ the time in seconds to construct a new HBS solver.

$r_{p}:=\frac{T_{new,\,p}}{T_{hbs,\,p}}$.

$T_{new,\,s}$ the time in seconds to apply the new solver.

$T_{hbs,\,s}$ the time in seconds to apply the HBS approximation of the inverse.

$r_{s}:=\frac{T_{new,\,s}}{T_{hbs,\,s}}$.

The ratios $r_{p}$ and $r_{s}$ illustrate the performance of the new solver relative to building a new HBS solver from scratch for the boundary value problem on the perturbed geometry.

5.1. Square with a nose

This section reports the performance of the new solver on the square with nose geometry illustrated in Figure 6. Two choices of “nose” height $d$ are considered: thinning and fixed. For the thinning nose geometry, the height $d$ is decreased as $N_{o}$ increases allowing $N_{c}$ to remain constant. For the fixed nose geometry, $N_{c}$ is thus increasing at the same rate as $N_{o}$. For both geometries, $N_{p}$ varies between 700 and 900 due to the corner-smoothing procedure.

Figure 7 presents log-log plots of the time in seconds versus $N_{o}$ for the (a) precomputation and (b) solve steps for the new solver and HBS solver for the thinning nose geometry. This figure and the timings reported in Table 2 illustrate that the new solver does scale linearly with respect to $N_{o}$ for this problem. The entry $r_{p}$ in Table 2 reports that the precomputation step of the new solver is approximately three times faster than the precomputation of the HBS solver. The solve step of the new solver is slower than the solve step of the HBS solver. However given the much larger constant associated with the precomputation step, it would take 100 to 260 solves to make the new solver slower than building a new HBS solver from scratch.

Figure 8 presents log-log plots of the time in seconds versus $N_{o}$ for the (a) precomputation and (b) solve steps for the new solver and HBS solver for the fixed nose geometry. As predicted in section 4.3, the computational cost for both steps does not scale linearly with respect to $N_{o}$. Table 3 reinforces this statement. Specifically, notice that the time to build the HBS solver ($T_{hbs,p}$) scales linearly with $N_{o}$ while the time for the new solver does not, and thus $r_{p}$ does not approach a constant. The lack of linear scaling is a direct result of the fact that the rank of the update matrix ${\mathbf{Q}}$ is dependent on $N_{o}$.

5.2. Circle with a bump

This section reports the performance of the new solver on the circle with bump geometry illustrated in Figure 12. Two choices of “bump” are considered: shrinking and fixed. For the shrinking bump geometry, the angle $\theta$ decreases as $N_{o}$ increases so that $N_{p}=N_{c}=199$ independent of $N_{o}$. For the fixed bump geometry, the angle $\theta$ remains fixed as $N_{o}$ increases. Thus both $N_{c}$ and $N_{p}$ increase at the same rate as $N_{o}$. These are the same geometries considered in [18].

Figure 10 presents log-log plots of the time in seconds versus $N_{o}$ for the (a) precomputation and (b) solve steps for the new solver and HBS solver for the shrinking bump geometry. As reported in Table 4, the precomputation of the new solver is linear and two times faster than the precomputation of the HBS solver. Since $r_{s}\sim 1$, the cost of the solve step is nearly the same. Thus, the new solver is the more efficient choice for this geometry.

Figure 11 presents log-log plots of the time in seconds versus $N_{o}$ for the (a) precomputation and (b) solve steps for the new solver and HBS solver for the fixed bump geometry. For this geometry, linear scaling of the new solver is not expected since the size of the update matrix ${\mathbf{Q}}$ grows with $N_{o}$ and $N_{p}$. The timings reported in Table 5 and Figure 10 support this statement.

5.3. Star with refined panels

This section considers the star with refined panels geometry. Such a geometry occurs in the construction of an adaptive boundary integral equation discretization technique. For this experiment $N_{o}$ is fixed, $N_{o}=3200$, while $N_{p}$ increases. Figure 13 illustrates that even though the direct solver does not scale linearly with respect to $N_{p}$, there is a range of $N_{p}$ such that the new solver is faster than constructing a fast direct solver from scratch. In addition to the information listed in the beginning of this section, Table 6 also reports the ratio of $N_{p}$ to $N_{o}$. This ratio says that for $N_{p}$ less than half $N_{o}$ the new solver is at worst two times faster than building a new direct solver from scratch. The addition of a fast direct solver applied $\bm{\mathsf{A}}_{pp}$ will keep the speed up factor large for larger $N_{p}$.

6. Summary

This manuscript presented a fast direct solver for boundary value problems on locally perturbed geometries. The solution technique is ideal for problems where the local perturbation involves removing a small number of points. Thus making the solver useful for optimal design problems where the perturbed geometry is placed in different portions of the original geometry and for improving the efficiency of adaptive refinement strategies. For problems where the number of cut points $N_{c}$ is constant (corresponding to the optimal type problems) the method is three times faster than building a new direct solver from scratch. For the adaptive refinement approach, using the new solver is faster when the number of new points is less than fifty percent of the number of the original points on the geometry.

Future work will include the non-trivial extension of the solver presented in this paper to three dimensional boundary value problems. This will involve careful management and processing of the geometry to make the best use of memory and limit communication. The integration of the new solver to an adaptive integral equation discretization technique is also a future project.

7. Acknowledgements

This research is supported by the Alfred P. Sloan foundation, the NFS (DMS-1522631) and the Ken Kennedy Institute for Information Technology.