Optimal design of optical analog solvers of linear systems

In this paper, given a linear system of equations A x = b, we are finding locations in the plane to place objects such that sending waves from the source points and gathering them at the receiving points solves that linear system of equations. The ultimate goal is to have a fast physical method for solving linear systems. The issue discussed in this paper is to apply a fast and accurate algorithm to find the optimal locations of the scattering objects. We tackle this issue by using asymptotic expansions for the solution of the underlyingpartial differential equation. This also yields a potentially faster algorithm than the classical BEM for finding solutions to the Helmholtz equation.


Introduction
In physical problems such as reflection of light in a three dimensional environment, aircraft simulations, or image recognition, we are searching for methods to numerically solve linear systems of equations of the form A x = b. Recently, optical analog computing has been introduced as an alternative paradigm to classical computational linear algebra in order to contribute to computing technology. In [9], the authors design a two dimensional structure with a physical property, which allows for solving a predetermined linear system of equations. In general, structures with such favourable physical properties are called meta-structures or meta-surfaces and are under very active research [3,7,8]. To be precise in their set-up, sending specific waves across the meta-structure modifies those waves, such that they represent the solution x to the problem.
One issue with this method is to quickly find the accurate structure for a given matrix A. In [9] the authors use physics software to gradually form such a structure. Here we demonstrate another method, which relies on using asymptotic expansions of solutions to partial differential equation. Such expansions have already been studied in different papers [2,1,4,5]. With that tool we can position and scale objects, on which the waves scatter, such that the resulting structure satisfies the desired requirements.
In the process of developing these asymptotic formulas, we have realized that we can compute the scattered wave using a method which is similar to the explicit Euler scheme. There we can numerically compute the solution to an ordinary differential equation by successively progressing in time with small time-steps until we reach the desired time. With our method we solve the partial differential equation around an obstacle by progressing the object-radius with small radius increments until we reach the full extent. We present the numerical application of that method on a circular domain.
This paper is organized as follows. In Section 2 we model the mathematical foundation for the underlying physical problem and define the fundamental partial differential equations for the asymptotic expansions. This leads us to the definition for the Neumann function on the outside. We then explain the connection of the Neumann function and the linear system of equations. In Section 3 we prove the asymptotic formulas concerning the Neumann function. There we discover special singular behaviours, which are essential to prove the asymptotic expansions. In Section 4, we first show the method to solve for the wave by increasing the radius by small steps and discuss the numerical error. Afterwards, we explain how we numerically build the meta-structure to solve the linear system of equations and discuss how well it operates. In Section 5, we conclude the paper with final considerations, open questions and possible future research directions. In the appendix we provide an interesting proof of a technical result and a modification of the trapezoidal rule, when we apply a logarithmic singularity.

Preliminaries
Let k ∈ (0, ∞) and let Ω be a finite union of disjoint, non-touching, simply connected, bounded and open C 2 -domains in R 2 be N ∈ N source points on the horizontal axis Λ = {x ∈ R 2 | x 1 ∈ (0, 1) and x 2 = 0}. We have N functions u j : R 2 → C, for j = 1, . . . , N, which solve the following partial differential equation: (2.1) Figure 1: This is the geometric set-up.
In black we have the domain Ω, in green the line segment Λ and in blue the source point z j . In violet, the absorbing layer is depicted.
where · | + denotes the limit from the outside of Ω to ∂Ω, and ∂ ν x denotes the outside normal derivative on ∂Ω. δ z j denotes the Dirac delta function at point z j and I j ∈ C denotes the intensity at the source z j . The first condition is the Helmholtz equation which represents the time-independent wave equation, and arises from the wave equation using the Fourier transform in time on the wave equation.
The second condition is known as the Neumann condition and models a material with a high electrical impedance. The third and fourth conditions model an absorbing layer on ∂R 2 + \ Λ and a Neumann condition on Λ. The fifth condition is known as the Sommerfeld radiation condition, and originates from a physical constraint for the behaviour of an outgoing wave.
We define Γ k to be the fundamental solution to the Helmholtz equation, that is Γ k solves PDE (2.1) without the Neumann boundary condition and a source at the origin. Furthermore we define Then we define the Neumann function N k Ω to be the solution to In contrast with u j , N k Ω is not only defined on the upper half. We recall that N k Ω (z, x) = N k Ω (x, z), which we can readily see using a Green's identity.
We can express N k Ω as a sum of Γ k and a smooth remainder, which satisfies PDE (2.2) with a vanishing right-hand side in the first equation. The same holds true for u j .
Using Green's identity on the convolution of u j with + k 2 N k Ω , we can infer for i = j that Using the trapezoidal rule we can approximate the integral in the last equation up to an error in O(N −1 ). Here we note that u j and N k Ω have a logarithmic singularity at z j , hence u j (z j ) is not well defined. Thus we use a slight modification in the trapezoidal rule, which is elaborated in Appendix A. After such modification, we define the complex column vector . Then we have that Our objective is to solve a linear system of equation A x = b using a physical procedure, in which an electrical signal is applied at z j , for every j = 1, . . . , N, and then is measured again at those points, for A ∈ C N×N and x, b ∈ C N , where A and b are given. The scattered wave, originated at z j , is in its Fourier space the function u j . Thus we are especially looking for a domain Ω, which yields and in searching so we keep track of the vector N j with the intention of rapidly determining the intensities I s such that ∑ N j=1 I j N j = b. In this paper we primarily consider rapidly finding the domain Ω such that Equation (2.3) holds.

Ω
Let Ω be a finite union of disjoint, non-touching, simply connected, bounded and open C 2 -domains in R 2 . Let k ∈ (0, ∞), then for z ∈ Ω, where Ω is the topological closure of the open set Ω, we define the outside Neumann function N k Ω (z, x), for x ∈ Ω, as the solution to the partial differential equation (2.2). Let B r (ζ) be a ball centred at ζ ∈ R 2 with radius r > 0.
We define Ω r as the union of some set Ω as defined above and of the ball B r (ζ), where B r (ζ) does not intersect Ω. Let N k Ω r and N k Ω be the outside Neumann function to Ω r and Ω, respectively. Theorem 3.1 For kr small enough and for all z, x ∈ Ω r , x = z, we have that where ∇ denotes the gradient to the second input in N k Ω , '·' denotes the dotproduct and O(·) denotes the limiting behaviour for r → 0. Additionally, for y ∈ ∂(Ω r \ Ω), we have that where (∇∇ T ) denotes the Hessian matrix, γ ≈ 0.57721 denotes the Euler-Mascheroni constant, and R k Ω (z, w) := N k Ω (z, w) − Γ k (z, w) has a removable singularity at w = z. For k · r > 0 small enough and for all y, w ∈ ∂Ω r , y = w, we have that Proof 1 Using Green's identity and the PDE (2.2) we readily see that

4)
where the normal vector still points outwards. Using an analogous argument by integrating N k Ω with itself, we see that N k We let x go to ∂B r (ζ) and apply the gradient on both sides and apply the normal at y. Then we obtain the equation where y as well as w are elements of ∂B r (ζ). We remark here that we cannot pull the normal derivative inside the integral. Let us consider next the decomposition is the fundamental solution to the Helmholtz equation, that means that ( w + k 2 )Γ k (y, w) = δ y (w), and R k Ω (y, w) is the remaining part of the PDE (2.2). Γ k can be expressed through 0 is the Hankel function of first kind of order zero and R k Ω is smooth [6]. From the decomposition in Equation (3.5) to arrive at by using the fact that the integral over ∂B r (ζ) decays linearly for r → 0. Transforming the normal derivative in the integral using polar coordinates, where we use that we have an integral over the boundary of a circle, we can infer that With the Taylor series for the Hankel function H 1 , for k r small enough, and considering the asymptotic behaviour of the forthcoming terms and applying some trigonometric identities we readily see that Using integration by parts, where we consider that N k Ω r (z, w(·)) is a periodic function, we obtain that Before we can proceed, we have to study a linear operator we callH, which takes a 2π periodic C 2 function ϕ and maps it to We can readily show that where 'p.v.' stands for the 'principle value'. This equation follows by integration by parts on both sides of the equation, and by using the integrability of the logarithm function. Now, we can state thatH is an invertible operator up to a constraint, according to [10, § 28], and that the solution toH[ϕ] = ψ is given , where the constraint is that 2π 0 ϕ = 0. Then we can infer thatH where we used thatH[O(r log(r))] = O(r log(r)). Thus it follows that for a constant function C in t. Next, we approximate the known function N k Using thatH[sin(·)] = cos(·) , we see that Analogously to Equation (3.4), we can formulate the statement that where z ∈ Ω and y ∈ ∂B r (ζ), and where we use that the Dirac measure located at y, which is at the boundary of the integration domain, which is a C 2 boundary, yields only half of the evaluation of the integrand at y. We then apply Equation (3.6) to the last equation and see that C = N k Ω (z, ζ). Applying it again for the second order term, while using Taylor expansions and comparing coefficients of the same order in r, we readily obtain the second equation in Theorem 3.1. For the first equation in Theorem 3.1, we apply the formula for N k Ω r (z, y(t)), and the Taylor expansion up to second order for ∂ ν y N k Ω (z, y(t)) to Equation (3.4) to obtain where (∇∇ T ) denotes the Hessian matrix which emerges from the Taylor expansion. We evaluate the two integrals explicitly, use that N k Ω = −k 2 N k Ω and obtain Equation (3.1). For Equation (3.3) we use Green's identity and obtain Solving for 1 2π π −π N k Ω r (y, u(t))dt and substituting we obtain Equation (3.3).
Let R > r > 0 and let Ω R be defined in the way that Ω r was introduced, that is Ω R is a ball of radius R at ζ ∈ R 2 adjoined to the domain Ω, hence Ω r Ω R . Then for any z r ∈ ∂B r (ζ), we define z R ∈ ∂B R (ζ) to be the projection of z r along the normal vector to B R (ζ). Thus we have that Proof 2 Equations (3.8) and (3.9) follow by readily using Green's identity on the convolution of N k Ω r with Γ k , and PDE (2.2), where we have to consider that an integral whose integration-boundary is over the singularity of the Dirac measure leads to half of the evaluation of the integrand.
For Equation (3.10), its proof is a simplification of the derivation of Equation (3.11). For Equation (3.11) we have with Green's identity that Splitting N k Ω in its singular part Γ k and its smooth remainder and subsequently extracting the singularity in Γ k , and doing so for N k Ω r as well, where we use Equation (3.9), we obtain that Using the technical derivation shown in Appendix A we prove Equation (3.11).
We decompose N k Ω r (z r (t z ), x r (t x )), for z r , x r ∈ ∂B r (ζ), into its singular part and a smooth enough part, that is, and furthermore we express N k Ω r through a Fourier series as Theorem 3.3 For kr > 0 small enough and for all z, x ∈ Ω r , x = z, we have that where the O L 2 term is a function with a L 2 (∂Ω r ) norm, which is in O (R−r) 2 r 2 , in the x r variable. Moreover, , The idea of proving this theorem is to extract the singularities developed in Lemma 3.2 in the integral expression for N k Ω R . Then any explicitly appearing integrals are solved in a similar way as described in Appendix A by using Fourier series.
Proof 3 Assuming z r = x r , we can use Taylor's theorem to obtain that for some w R,r ∈ R 2 between z R and z r . We note that ∂ ν zr N k Ω r (z r , x r ) = 0. We need the term 1 2 (R − r) 2 ∂ 2 ν zr N k Ω r (w R,r , x r ) to be in O L 2 , but that is not the case due to the singular term in N k Ω r . Hence we extract the singular term from N k Ω r (z R , x r ) and then we can infer that Extracting the term log( R r ) from the logarithm term and using the Taylor approximation for R → r, on that extraction, we then obtain Equation (3.13). Considering Green's identity we have that Next we apply ∂ x R on both sides and then interchange the integral and ∂ x R . This leads to the term ∂ ν x R ∂ ν yr N k Ω (x R , y r ), whose singular part we extract from ∂ ν x R ∂ ν yr N k Ω (x R , y r ). The equation then reads for all n ∈ N 0 , which we can readily show from the 2π-periodicity by using trigonometric formulas and applying an induction on n ≥ 1. Furthermore, we have that

Then we use Equation (3.13) and this leads us to the equation
With that identity we can determine all integrals in Equation (3.16) and show Equation (3.14). For an elaborated calculation of the third integral, see Appendix A.
Using Green's identity on N k Ω R (z R , x R ), we can infer that Similar to the derivation of Equation (3.14), we can compute that Using that π −π 1 2π we readily see that where the logarithm term is derived similarly as in Appendix A. We consider the integral term in Equation (3.18). To this end we will apply Equation (3.14) and consider the singular parts of ∂ ν x R N k Ω r (z R , x R ). Hence we define , because using Taylor series we have that Then, applying the singular decomposition to the integral in Equation (3.18), and using the same techniques as are those used in Appendix A, we have that

Then we can apply Equation (3.20) and obtain
We can further simplify this approximation by using Green's identity on N k Ω r (z r , x r ), with N k Ω (z r , x r ), and using Taylor series on N k Ω (z R , x R ) and on N k Ω . This leads us to the equation Using that ∂ ν y R (t) N k Ω r (x R , y R (t)) = O(R − r), 1 − ( r R ) n = O(R − r) and that the logarithm in the second integral is in O L 2 (R − r), we can infer that and thus make further simplifications, which lead to Equation (3.15) and finishes the proof. Remark 3. 4 We note here, that Equation (3.22) is numerically more stable than Equation (3.15) in Theorem 3.3. We expect the reason to be that the constant error in the first step is lowered by the factor (R − r) and the factor 1/2 in every subsequent step.

Applying Theorem 3.3 -Gradually Increasing the Radius
With Theorem 3.3 we are able to evaluate the Neumann function N k Ω r (z r , x r ) while we increase the radius of the circular sub-domain B r (ζ) in Ω r by ∆r, where z r , x r ∈ ∂B r (ζ), with an error in O((R − r) 2 ). Similar to how we numerically evaluate the solution to an ordinary differential equation y(t) = f (t, y(t)), y(t 0 ) = y 0 , using the explicit Euler scheme, where we start at t 0 and then evaluate the function y at t 0 + ∆t with an error in O((∆t) 2 ), we can now evaluate the function N k Ω r (z r , x r ) at radius R = r + ∆r. For the Euler scheme, we can show using Grönwall's inequality that the global error is O(∆t). Thus we expect the global error of N k Ω r (z r , x r ) to be O(∆r). The domain Ω r for the numerical evaluation is set to be Ω r = B r ([0, 0] T ) ∪ B 1 ([1, 2.5] T ). We increase the radius of B r in Ω r by ∆r successively until the radius reaches 1. In every step we compute the first N f Fourier coefficients of the smooth part of N k Ω r (z r , · ), see (3.12), using Theorem 3.3 with Remark 3.4, where we also have to discretize z r in such a way that we have N f equidistant points on ∂B r (0), where one point is set at [−r, 0] T . For the first step we use Equation (3.3) in Theorem 3.1.
In Figure 2 we have depicted the error, which we calculated using MATLAB, between the actual Neumann function and the approximation given through the algorithm corresponding to ∆r. To be more precise, we computed all possible N 2 f discretized values of the smooth enough part of N k Ω 1 (z 1 , x 1 ) for z 1 , x 1 ∈ ∂B 1 (ζ) and averaged them in the numerical approximation. The actual Neumann function was numerically computed using the BEM with a very large number of boundary points. The Figure shows that we indeed achieve an error in O((∆r) 1 ). It seems that we even achieve a higher order than only a linear one, but this is not further investigated here.
Comparing this numerical approximation with the BEM, we see that for this approximation we have have a runtime complexity of O(N 2 f ) × Figure 2: We have depicted the average error between the smooth enough Neumann function N k Ω (z 1 , x 1 ) and the numerical approximation with respect to the radius increase ∆r, at the N f ∈ N points z 1 , x 1 = 1 · [cos(t n ), sin(t n )] T , t n = −π, ..., π. We see that the error is at least linear in ∆r. We set here k = 1, N f = 2 8 . O((∆r) −1 ) and an error in O(∆r) multiplied to an error with respect to N f , which in the above numerical experiments had no influence. For the BEM we have to invert a N c × N c matrix, where N c is the amount of discretisation points used on the boundary, which yields an error in O(N −1 c ) and has a complexity runtime of O(N 3 c ) in simple algorithms.

Reconstructing a Matrix
In this section we use the approximation shown in the last section to determine a specific scattering matrix S, as it is elaborated in Section 2. Different than in Equation (2.3) we search here for a matrix S, which is as close as possible in average value to all entries to a predetermined Matrix, which we call in this subsection matrix A ∈ C N×N . Thus we try to minimise the value e(S) = 1 N 2 ∑ i,j |A i,j − S i,j |. The procedure to form such a matrix S is as follows. We have N source points (z i ) N i=1 equidistantly distributed in (0, 1) × {0}. When there are no scattering objects placed in R 2 + , then the Neumann function N k Ω is simply the Γ k function, and hence S i,j = ∂ (z j ) 2 N k Ω (z i , z j ) = 0. Next we place a small ball within R 2 + , where we place the center so that the error e(S) is minimised, which we in turn calculate using Theorem 3.1. We did this minimization classically using a grid of points, but can in general be realized with more sophisticated methods as for example with the gradient descend method. Given the initial ball, we increase its radius using Theorem 3.3 as it is shown in the previous section. After every increase we compute the Neumann function at the source points using the associated integral formulation, that is, Thus we can compute e(S), with the objective to see whether the error decreases or increases and whether we should increase the radius further or not. As soon as an increase in the radius does not yield a lower error, we search for a place to add another small ball. We again use Theorem 3.1 to determine the next best place to center the ball. In addition, we need to calculate N k in order to apply Theorem 3.3, where ζ i , ζ j are values in R 2 + \ Ω r and where (∇ ζ i ∇ T ζ j ) denotes a 2 × 2 matrix in which the entries are the respective coordinate differentiations. To this end, we use the integral formulation above, in which we can interchange integration and differentiation. In practice, we used a linear interpolation to speed up the calculation. After we established a new place for the small ball, we can also increase it until the error e(S) does not decrease any further. And then we search for a place for a third ball, and then a fourth and so forth until we cannot decrease e(S) any further. This algorithm is explicit and does not use the inversion of any matrix as it is commonly done using a BEM.
For our first numerical experiment, we set our predetermined matrix A to be the scattering matrix of a predetermined domain, which is given on the left-hand side in Figure 3. Using the algorithm described above, we obtain the domain on the right-hand side. On the left-hand side in Figure  4 we see a heat-map of the real part of the matrix A and on the right-hand side we see a heat-map of the real part of the scattering matrix S.
For more general matrices A we need more sources than given by the size of A. To have such a more extended matrix we have to cast A to a integral of the form Λ K(z i , y)u(y)dσ y , and finally discretize that integral, and then apply the algorithm to the discretization.

Concluding Remarks
We considered the physical experiment presented in [9], in which scattering objects were placed in front of signal sources. Those sources send waves which reflect at the object and then receiving points collect the wave intensity. The registered intensity is the solution to a predetermined linear system of equations. Hence, instead of solving the linear system with mathematical means, we can solve it using a physical set-up, which is substantially faster. The complication arises in finding the exact configuration of the scattering objects. Using a mathematical model for the underlying physical problem we were able to describe the PDE using the Neumann function. Studying its asymptotic behaviour when we place tiny scattering objects and when we increase the extent of those objects successively, we were able to develop an explicit algorithm to place and enlarge objects such that the scattering matrix approaches the predetermined matrix, which is needed to solve the linear system of equation. In Section 4 we showed that the numerical implementation for calculating the Neumann function when we enlarge an object works better then intended, in regard of the explicit Euler scheme. With such an algorithm we have a new and faster numerical method to calculate the Neumann function than using the ordinary BEM. We then applied that process to approach a desired matrix.
In this paper we considered circular scattering objects. It would be interesting to have more complicated domains such as ellipses, which would allow for one more easily accessible degree of freedom to control the waves. We think that the mathematical proofs in Section 3 can be readily extended to more complicated C 2 -boundaries. To this end, we need to consider a function ϕ : (−π, π) → R 2 , which described the boundary, and consider it in the integration formulae.
In the last section we mentioned that reconstructing a more general matrix A in a linear system of equation does not work well. We need more options in our algorithm, or a bigger matrix, which has similar properties to A, and additionally can be described as a kernel of an integration operator. In [9], the authors set the matrix to be the lower left quadrant of their scattering matrix.
We are looking forward to see these asymptotic formulae being used in other physical problems concerning scattering problems. We are also very curious to see improvements in the object reconstruction of general linear systems and hope that our research will lead to an improvement of mathematical and technological tools for numerical computing.