Van-Phu Ha
email: van-phu.ha@inria.fr

Olivier Sentieys
email: olivier.sentieys@inria.fr

Leveraging Bayesian Optimization to Speed Up Automatic Precision Tuning

Using just the right amount of numerical precision is an important aspect for guaranteeing performance and energy efficiency requirements. Word-Length Optimization (WLO) is the automatic process for tuning the precision, i.e., bit-width, of variables and operations represented using fixed-point arithmetic. However, state-of-the-art precision tuning approaches do not scale well in large applications where many variables are involved. In this paper, we propose a hybrid algorithm combining Bayesian optimization (BO) and a fast local search to speed up the WLO procedure. Through experiments, we first show some evidence on how this combination can improve exploration time. Then, we propose an algorithm to automatically determine a reasonable transition point between the two algorithms. By statistically analyzing the convergence of the probabilistic models constructed during BO, we derive a stopping condition that determines when to switch to the local search phase. Experimental results indicate that our algorithm can reduce exploration time by up to 50%-80% for large benchmarks.

I. INTRODUCTION

The rapid development in scientific and technological innovations during the last decade opens a new era for intelligent and sophisticated systems. The demand for integrating many large applications in a limited silicon area poses new challenges for energy-efficient computing. Recently, Approximate Computing (AC) is considered as a good solution to address energy efficiency issues. The primary objective of approximation techniques is to trade quality of service for cost saving. One of the popular AC techniques is to use Fixed-Point arithmetic for lowprecision computation in Digital Signal Processing or Machine Learning systems. This technique always requires a floatingpoint to fixed-point conversion that optimizes the fixed-point word-lengths for a good compromise between cost and quality requirement. This procedure, called Word-Length Optimization (WLO), accounts for 25-50% of design time [START_REF] Clark | Accelerating fixed-point design for mb-ofdm uwb systems[END_REF] and is still considered as a problem of interest to reduce time-to-market.

Approaches to address WLO can be classified in two groups: analytical and simulation-based approaches. Analytical approaches relax the WLO problem for convexity and then apply some convex optimization algorithms to directly obtain the optimal solution [START_REF] Chan | Wordlength optimization of linear timeinvariant systems with multiple outputs using geometric programming[END_REF], [START_REF] Fiore | Efficient approximate wordlength optimization[END_REF]. Despite handling WLO quickly, these approaches require the accuracy to be modeled as convex functions, which cannot be analytically constructed in general. Simulation-based approaches solve WLO by iterative search using simulations [START_REF] Cantin | An automatic word length determination method[END_REF]- [START_REF] Nguyen | Novel algorithms for word-length optimization[END_REF]. They are thus generic with all systems and quality metrics. However, these approaches do not scale well in large applications where many variables are involved since the number of simulations increases dramatically with the number of variables. In addition, most of these iterative searches are based on local search that moves with a short distance in the discrete domain. Hence, convergence speed is slow if the initial point is far from a local minima.

Bayesian Optimization (BO) is a popular approach for tuning hyper-parameters in machine learning algorithms [START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF]. This approach aims to optimize problems where the mathematical expression of the target function is unknown and without derivatives. BO constructs a probabilistic model based on past samples to suggest new points [START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF]. Thanks to this model, at a certain state, BO can ignore neighboring points if they produce a low probability of being good solutions to search more quickly. This feature might overcome weakness of the simulation-based approaches. However, there is no study yet indicating the performance of Bayesian optimization for WLO.

In this paper, we propose a Hybrid algorithm combining Bayesian optimization and a local search to improve the scalability of simulation-based approaches. We first show how this combination can lead to large improvements in exploration time. Then, we design a reasonable transition point between the two approaches to leverage the efficiency of the combination, which is also the core of our approach. Using design points sampled by BO, we derive a statistical metric to evaluate the convergence of the models during BO. Experimental results show that our Hybrid algorithm outperforms latest simulationbased approaches, reducing exploration time by up to 50%-80%, while leading to similar cost solutions.

The rest of the paper is organized as follows. We introduce necessary background and discuss related work in Section II. A motivation for our work is presented in Section III. Then, we describe our proposed Hybrid approach in Section IV. We show the performance of our approach in a comparison with different benchmarks and latest approaches in Section V followed by a conclusion in Section VI.

II. BACKGROUND AND RELATED WORK

In this section, we introduce the WLO problem and discuss earlier work before presenting Bayesian optimization.

A. Word-Length Optimization and Classical Approaches

A number represented in fixed-point arithmetic contains integer and fractional word-lengths (WLs), represented on I and F bits, respectively. The integer WL covers the dynamic range whereas the fractional WL controls the precision. In this work, we focus on the WLO for fractional WL, which is the time consuming part of the exploration. Let the vector W = [W 0 , W 1 , . . . , W N -1 ] denote a word-length configuration with N effective variables to be explored for fixed-point conversion.

The main objective of WLO is to determine a good-enough word-length configuration that minimizes a cost function under a quality constraint:

min C(W ) Subject to λ(W ) ≥ λ obj (1)
where C and λ are functions that express cost and quality, respectively, and the quality target is given as λ obj . How the functions C and λ are realized varies across work, ranging from analytical models to simulation-based approaches. Some approaches [START_REF] Chan | Wordlength optimization of linear timeinvariant systems with multiple outputs using geometric programming[END_REF], [START_REF] Fiore | Efficient approximate wordlength optimization[END_REF], [START_REF] Parashar | A polynomial time algorithm for solving the word-length optimization problem[END_REF] construct analytical models for noise power through mathematical expression to avoid costly simulations during the exploration process. These analytical approaches take advantage of a property of errors, under the hypothesis of linear and time-invariant (LTI) systems, that their propagation do not interfere with each other. Hence, the error introduced at a noise source may be propagated through the system independently and aggregated afterwards. Thus, these approaches cannot be directly extended to handle general non-LTI programs. Moreover, complex quality metrics, such as Structural Similarity (SSIM), which are not directly related to noise power, are hard to model analytically.

Many simulation-based approaches [START_REF] Cantin | An automatic word length determination method[END_REF]- [START_REF] Nguyen | Novel algorithms for word-length optimization[END_REF] were proposed based on iterative search using heuristics to address WLO. These approaches use variants of gradient descent algorithms that evaluate neighboring solutions which differ by one or few bits from the current solution at each iteration. Since the number of neighboring solutions increases with the number of variables, the number of solutions that must be evaluated at each iteration quickly increases with the complexity of the application. Additionally, due to short movements at each iteration (one or few bits), these approaches require many iterations to converge if the starting point is far from a local minima. These are the main drawbacks causing a huge number of simulations especially in large applications and hence leading to an extremely long exploration time. Min+1, Min+b [START_REF] Cantin | An automatic word length determination method[END_REF], [START_REF] Han | Numerical word-length optimization for cdma demodulator[END_REF] and Max-1 [START_REF] Cantin | An automatic word length determination method[END_REF] are typical approaches of this group. The Min+1 (resp. +b) algorithm begins with a design where each variable is assigned by its minimum WL (MWL), i.e., the WL satisfying the quality target when other variables are set to the highest precision. At each iteration, the algorithm moves towards neighboring candidates by increasing by 1 (resp. b) bit(s) the best variable among the candidates until obtaining a solution satisfying the quality constraint. By contrast, Max-1 initializes variables with the highest possible WL and decreases variables by 1 bit until to reach the final solution. Afterwards, many variants of these two procedures were proposed. Heuristic approaches [START_REF] Sung | Simulation-based word-length optimization method for fixed-point digital signal processing systems[END_REF] or Hybrid combinations of Min+b and Max-1 outperform the original algorithms [START_REF] Cantin | A comparison of automatic word length optimization procedures[END_REF]. Likewise, the Greedy Randomised Adaptive Search Procedure (GRASP) [START_REF] Nguyen | Novel algorithms for word-length optimization[END_REF] combines local search and stochastic optimisation. GRASP is an iterative two-phase procedure. In the construction phase, the search algorithm (similar to Min+1) randomly selects one of the best neighboring candidates during gradient descent. Then, a Tabu search is applied to refine the solution found by the first phase. Tabu allows movements in both directions (+1/-1) and uses a Tabu list to skip some explored variables. These two phases are iterated and the randomization of the construction phase avoids to stay in local minima. Some noise budgeting techniques [START_REF] Ha | Towards generic and scalable wordlength optimization[END_REF], [START_REF] Novo | Cracking the complexity of fixed-point refinement in complex wireless systems[END_REF] decompose large applications into smaller kernels to break down the exponential complexity of WLO. However, these methods still face the limitation of the classical above-mentioned approaches, which are used for each kernel to solve local WLO problems.

B. Bayesian Optimization

Bayesian Optimization (BO) is a machine-learning-based optimization method [START_REF] Shahriari | Taking the human out of the loop: A review of bayesian optimization[END_REF] aiming to optimize functions which usually have no mathematical expression and/or derivatives. Despite being widely applied in many real world problems, there is no study of BO for WLO. Generally, BO has two key elements: i) a probabilistic surrogate model for modeling the unknown objective function based on already observed samples and ii) an acquisition function that optimizes over the surrogate model to suggest next samples. One main difference in BO methods is the process of selecting the surrogate models. While Gaussian Processes (GP) are often used for continuous-domain moderate-size problems, tree-based models like Random Forest and Tree-structured Parzen Estimator (TPE) facilitate discretedomain large-size problems [START_REF] Eggensperger | Towards an empirical foundation for assessing bayesian optimization of hyperparameters[END_REF]. TPE works by identifying points that could have been drawn, and that appear promising on the basis of the evaluation of a loss function at other points. This loss function is very important to guide the optimizer, as detailed in Section IV-A.

Unlike gradient-based algorithms, BO inspects past iterations to construct the interplay among variables to evaluate a loss function via a probabilistic model. Based on the knowledge from this model, BO is able to evaluate and ignore neighboring designs which have low probability of being good solutions to search more globally. As a result, it can speedup the search and better avoid local minima if current state is on a plateau or a bad local minima. Thus, BO is a promising candidate to overcome limitation of classical simulation-based WLO approaches. However, the computational complexity of BO is quadratic O(i 2 ), whereas most heuristic approaches have a linear complexity O(i), i being the number of iterations. Thus, BO still faces scalability issues, especially for large applications in which BO requires more evaluations to obtain a good solution. Based on their features and properties, we propose a new method that combines Bayesian optimization and local search with a reasonable transition condition to improve the exploration efficiency. TPE is selected for BO because it is suited to discrete problems as WLO, and Tabu, a state-of-theart WLO method, is used for local search.

III. MOTIVATIONS

In this section, we compare the performance of TPE with Tabu combined with Min+1 gradient descent for WLO. Then, we indicate benefits of combining TPE and Tabu via empirical evidence, which motivates the proposed Hybrid approach presented in Section IV. We used an Infinite Impulse Response (IIR) filter with Peak Signal to Noise Ratio (PSNR) equal to 40 dB as a quality constraint. Note that the trends in this section are consistent with the other benchmarks used in Section V-B.

A. Performance Analysis: TPE vs. Tabu

Figure 1 shows the search process of Tabu and TPE during WLO of the IIR. Tabu searches in a narrow range constituted by near neighboring solutions which are some bits different to each other. Thus, it moves slowly towards optimal region. At some points in the initialization procedure, the search takes a significant period to surpass the plateau region which has no quality improvement. As a result, local search must wait for a long time and is only started when a feasible solution is found in the initialization procedure. Meanwhile, TPE first searches randomly in the solution space (very first solutions scattered in around 35-55 dB). Afterwards, it explores solutions using probabilistic models in a wide range and quickly finds feasible solutions around the quality target. Note that by using models, TPE can stick around the quality target to enhance the search efficiency instead of whole solution space as random search.

For cost comparison (see Section V-A for more details on cost), Figure 2 shows the best solution obtained so far by TPE as a function of exploration time and the final solution obtained by Tabu. TPE converges quickly in the beginning. However, the cost is improving slowly for latest iterations, taking up a large portion of the total execution time. With this behavior, we emphasize the interest of stopping TPE in pruning phase to switch to a more efficient local search, such as Tabu.

B. Initial Combinations of TPE and Tabu

From the performance analysis of TPE and Tabu, we empirically perform initial combinations of the two algorithms. TPE is stopped at a certain moment returning a temporary configuration, defined as a transition point, which is then served as starting point for Tabu. Among the designs explored by TPE, we use only the best solution obtained so far. These transition points (orange crosses in Figure 2) are selected if they have a significant improvement in terms of cost compared to the previous adjacent selected one. This ensures a high chance of having different WL configurations.

In Figure 3, we compare the performance of Tabu with TPE+Tabu for different transition points. The cost of the final solutions depends on the moment when TPE is stopped. Stopping TPE too early (H-8 to H-283) does not lead to good solutions compared to Tabu alone, even though their

TPE Tabu

Fig. 3: A comparison in energy cost (top), total exploration time (middle) and separated exploration time (bottom) of Tabu and different combinations between TPE and Tabu given the selected transition points. The x-axis separates the data for different cases; H-X represents for the different combinations where X is the index of the stopping point of TPE. convergence speed is much faster. Stopping TPE after a certain amount of time (H-578, H-1368) can lead to an equivalent or even better solution cost, still in a shorter time. The benefit of the superior convergence speed is gradually decreased if TPE is stopped too late (H-1629 to H-2538). It is important to notice that, except for the early ones, the choice of the exact transition point does not impact much the cost of the final solution.

From this analysis, we highlight that TPE can quickly prune the design space and find a solution from which Tabu will be able to fine-tune the cost, resulting in a significant reduction in exploration time. This reduction depends on choosing the right transition point, which is not an obvious task especially in an automated way, and is the aim of the next section.

IV. PROPOSED HYBRID APPROACH

In this section, we construct the loss function used in our Bayesian optimization and define a method to automatically find the transition point between TPE and Tabu.

A. Loss Function

Derived from the original WLO problem of Eq. 1, we construct the loss function using the Lagrangian form as

f (W ) = C(W ) -α(λ(W ) -λ obj ), (2)
with a positive and big enough α. TPE will then follow f (W ) to find solutions that minimize the original problem.

Based on many WLO experiments, we experienced that the cost and quality functions tend to be proportional to W . A high quality solution mostly corresponds to a high cost. Thus, the best solutions to satisfy the constraint with a small cost are likely to be around the quality target λ obj . Therefore, for faster convergence, we force the loss function to cover only a narrow range [q l , q h ] around λ obj . The loss function is then defined as

f (W ) = C(W ) -α(λ(W ) -λ obj ) if λ(W ) ∈ [q l , q h ] +∞ otherwise (3)
Moreover, to penalize solutions below λ obj with regards to solutions above the target, we choose α = 0.5 for solutions below λ obj , and α = 0 for other solutions in the range [q l , q h ].

B. Transition Point for Hybrid Approach

In Bayesian optimization relying on TPE [START_REF] Bergstra | Algorithms for hyper-parameter optimization[END_REF], hyperparameters -WL configurations W i in our case-are chosen uniformly over the search space and evaluated by the loss function f (W ). Then, obtained samples {W i , f (W i )} stored in an observation history H are divided into two groups. The first group contains good samples where the loss f i is less than a threshold γ * , whereas the second group consists of the remaining, considered as bad samples. TPE uses these two sample groups to model two likelihood probability density functions l(W ) = p(W |f (W ) < γ * ) and g(W ) = p(W |f (W ) ≥ γ * ), respectively. Then, it decides which hyper-parameter to try in the next iteration by maximizing the ratio

l(W ) g(W ) = p(W |f (W ) < γ * ) p(W |f (W ) ≥ γ * ) .
Algorithm 1 describes TPE algorithm with a stopping condition Algorithm 1 TPE with stopping condition

W * = argmax l(W ) g(W )
4:

Evaluate f (W * )

5:

H ← H ∪ (W * , f (W * )) 6:
if stopping condition is satisfied then Update l(W ) and g(W ) given H 10: end for checked at each iteration. The algorithm stops if the stopping condition is satisfied and then returns a WL configuration W stop which serves as the starting point of Tabu. How the condition is constructed is described as follows. At t th iteration, we consider the vector S t L comprising the L top-performing solutions from the observation history H. Our main objective is to evaluate if the next iteration of TPE can lead to a new top-performing solution, that is significantly different in terms of WL configuration compared with the solutions in S t L . Then, we decide to stop TPE if the solutions in S t L share similar WL configurations. Otherwise, we still wait for further exploration iterations. The similarity of WL configurations in S t L reflects that TPE is likely to converge to a local region formed by a number of similar WL solutions, which only differ by few bits. Thus, continuing to spend time with TPE is likely to be inefficient. Instead, a local search like Tabu is more reasonable to converge quickly in that local region. This is the key idea in our method for the combination. We use TPE for narrowing the solution space down to a good region, from which Tabu will continue the search for fine-tuning to provide the final solution following its greedy gradient.

Each solution in S t L is represented by a WL configuration with N effective variables. Let the L × N matrix W t contain WL configurations of the L top performing solutions at iteration t. Elements w t i,j of W t represent the number of bits of variable j ∈ [1, N ] from solution i ∈ [1, L]. We use the standard deviation to evaluate the distribution of WL values for each variable (i.e., each column of W t ) as

Σ t = [σ t 0 , σ t 1 , . . . , σ t N -1 ], (4)
where σ t j is the individual standard deviation of WL values of the j th variable and is calculated as

σ t j = 1 L -1 L-1 i=0 (w t i,j -W t j ), W t j = 1 L L-1 i=0 w t i,j . (5)
The smaller the standard deviation, the more similar the solutions. To evaluate the similarity on all WL configurations at a certain iteration, we aggregate standard deviations by the average of all σ t j with j ∈ [0, N -1] as

Σ t = 1 N N -1 i=0 σ t i . (6)
Σ t represents the average difference on the number of bits for the WL configurations of the L top-performing solutions. By evaluating the convergence of TPE through Σ t at each iteration, we obtain the same behavior on all our benchmarks. In Figure 4, we show an illustration of this behavior for the IIR example.

In the beginning, the curve of Σ t varies dynamically with a high variance, around 4.5 bits. It constantly finds new good solutions which are relatively different to each other, resulting in dissimilar WL configurations in S t L and then producing a high value of Σ t . Afterwards, it enters a more stable area with a value of Σ t less than 1 bit. This indicates that TPE is likely to converge to a local region. So, we should stop TPE and switch to a greedy algorithm like Tabu to search locally around this region. Based on that behavior, we derive two indicators to detect the stable region to decide to stop TPE. At the t th iteration, we evaluate the stability of Σ t in the M latest samples via mean and standard deviation as

µ(Σ, M ) = 1 M t i=t-M -1 Σ i (7) σ(Σ, M ) = 1 M -1 t i=t-M -1 (Σ i -µ(Σ, M )) (8)
µ(Σ, M ) indicates the average magnitude of Σ. σ(Σ, M ) is chosen with a small value to ensure the stability of Σ. Then, two corresponding thresholds are chosen to decide the adequate stopping moment for TPE. In our evaluation, we monitor L = 25 top-performing solutions at each iteration and choose two conditions, µ(Σ, M ) ≤ 1 and σ(Σ, M ) ≤ 0.1, with M = 100 to stop TPE. Figure 4 indicates with a dashed line the iteration where TPE is stopped.

V. EVALUATION A. Experiment setup

Experiments are performed on an Intel Xeon E5640 2.67GHz with 4GB memory running Linux. We used Adaptive TPE (ATPE), an extension on top of TPE implemented on Hyperopt [START_REF] Bergstra | Hyperopt: a python library for model selection and hyperparameter optimization[END_REF], which provides a machine learning model to automatically tune the hyper-parameters, e.g., γ * , of TPE. Tabu search is described in [START_REF] Nguyen | Novel algorithms for word-length optimization[END_REF]. We use energy as our cost model. The energy model counts the number of operations performed by each operator, and calculates the total cost based on the energy consumption of an operation empirically gathered from several ASIC synthesis/simulation for different WLs. An operator is characterized by the WLs of the operands, the WL of the result, and the arithmetic operation performed. Characterization is performed using Synopsys Design Compiler and Prime Time using a 28nm FDSOI technology.

We compare our Hybrid approach with reference algorithms (Max-1, Min+1, Tabu search and GRASP, as in Section II-A) in terms of solution cost and exploration time. GRASP is configured with T RCL = 3 running in 10 iterations as in [START_REF] Nguyen | Novel algorithms for word-length optimization[END_REF]. During the execution of Min+1 and Tabu, we experienced that the search process usually gets stuck in situations that keep increasing the number of bits and thus cost, without a significant quality improvement. Actually, in such situations there are more distant neighbors that gain a significant amount of quality. As a result, more quality can be achieved without increasing the number of bits too much, which reduces redundant cost. The bias versions of Min+1 and Tabu search add a bias to the selection criterion to reduce the priority of the candidates for which a significant number of bits has already been added, which leads to optimized versions of the original algorithms.

We evaluate cost and exploration time on five applications. FIR and IIR filters are implemented with 5-stage cascaded structure of 33 th and 2 nd order filters, respectively. The number of effective variables to optimize is 17 for FIR and 33 for IIR. Block-Matching and 3D Filtering (BM3D) [START_REF] Lebrun | An analysis and implementation of the bm3d image denoising method[END_REF] is a stateof-the-art image processing noise reduction algorithm with 45 effective variables. ISP is an Image Signal Processor chain which applies a sequence of processing kernels on raw data from camera sensors to produce an enhanced color image. ISP comprises four kernels (Non-Local Means (NLM) denoising filter [START_REF] Buades | A non-local algorithm for image denoising[END_REF], Demosaicing, Gamma Correction, and Unsharp filtering) and has 74 effective variables. Stereo Matching (SM) is widely used in computer vision to extract a depth information of a scene from two images taken from two cameras. SM includes 85 effective variables to optimize.

We use PSNR as the quality metric for filters and Structural Similarity (SSIM) for image processing applications. Two quality constraints λ obj are evaluated for PSNR (40 dB and 50 dB) and SSIM (0.9 and 0.99). We use the narrow range [q l , q h ] as ±2 dB around λ obj for PSNR, ±0.09 for SSIM=0.9, and ±0.009 for SSIM=0.99.

B. Performance Evaluation

Figure 5 presents the comparison between our Hybrid approach and the references in terms of solution cost (top) and exploration time (bottom) for our five benchmarks and two accuracy constraints. Max-1 converges fast but always leads to the worst solution cost. Conversely, GRASP obtains relatively competitive solutions but with a long exploration time. The runtime of GRASP is proportional to the number of iterations used to obtained different randomized starting points to avoid local minima. Most solutions found by GRASP are slightly better than those obtained by Min+1, Min+1[bias] and Tabu. However, GRASP does not outperform Tabu[bias] and is therefore not that efficient in spending longer time through randomization (see e.g., ISP-0.9 example where GRASP still gets stuck in bad local minima).

Thanks to the local bi-directional search, Tabu always find better solutions compared to Min+1, and sometimes the bias version is better, sometimes not. However, for a fair comparison, we always use the best reference as a competitor to our Hybrid approach. We also provide the solutions obtained by TPE only iterating for the same time as Tabu[bias]. Solutions obtained by TPE are relatively competitive to Tabu[bias]. For BM3D-0.99, SM-0.9, and SM-0.99, TPE found better solutions than Tabu, which indicates that TPE is able to avoid bad local minima better than Tabu.

The results show that our Hybrid approach always find competitive solutions in shorter time than other methods. For small problems solved in minutes such as FIR-40 and FIR-50, it is sufficient to use classical approaches because they are able to converge in short time. For larger problems solved in a few hours to a few days, our approach can reduce by 50% to 80% (66% on average) exploration time compared to the best reference with similar cost. Our approach even defeats TPE in most benchmarks with slightly better solutions found in a much shorter time (2x-5x faster than TPE, 3.7x on average). Apart from IIR-50 and SM-0.99 where solutions obtained by our approach is 0.4% and 14.5% worse than TPE, respectively, the other cases record a 0.4%-6.3% solution improvement (3.3% on average) compared to TPE. It is clear that the speedup of our approach comes from the benefits of i) TPE to speedup the first phase and ii) our transition condition to stop TPE early. The solution cost improvement takes advantage of the good convergence of Tabu in local regions. In absolute execution time, our approach reduces exploration from 7.2 days to 1.5 days compared to TPE for SM-0.9, from 4 days to 18.8 hours compared to Tabu[bias] for BM3D-0.9, and from 2.5 days to 13.3 hours compared to Tabu[bias] for ISP-0.99.

VI. CONCLUSION

In this paper, we propose a Hybrid method leveraging Bayesian Optimization and local search to improve the scalability issue that WLO faces for large applications. We first show how TPE is efficient at narrowing down the search on good regions thanks to its randomness and Tabu for fine-tuning because of its bi-directional local search. Then, the core of our approach is to derive a statistical metric to evaluate model convergence of TPE, to automatically find the right transition moment between TPE and Tabu, and to maintain the search efficiency. For large applications, our approach significantly outperforms state-ofthe-art approaches, with a 50%-80% reduction in exploration time, while still providing similar or better cost solutions. The exploration of the parameters used in our methods is important but left as future work.

Fig. 1 :

 1 Fig. 1: Search process of Tabu and TPE. Random search is used as a reference. Points correspond to quality of different designs evaluated in the search process (for Tabu, all neighboring solutions in an iteration are plotted at a corresponding time). Tabu has two procedures: i) Initialization uses Min+1 to search in infeasible region until a solution satisfying the quality target is found and ii) Local search, a combination of Min+1 and Max-1, optimizes locally around the quality target.

Fig. 4 :

 4 Fig. 4: Stopping point for IIR, quality target PSNR = 40dB.

Fig. 5 :

 5 Fig. 5: Energy cost of solutions (top) and exploration time (bottom) normalized to Tabu[bias]. Numbers in parentheses are the execution time in minutes of Tabu[bias] used for the normalization. Results for Max-1 are removed for better readability when leading to bad solutions.

  Manually selected stopping points on the best solution obtained so far by TPE. The curve was constructed by best designs, satisfying the quality target with minimum cost, at a certain time. The numbers indicate index of the selected points.

	Energy Cost (nJ)	0.0009 0.0010 0.0011 0.0012 0.0013 0.0014	Minutes 0 25 50 75 100 125 150 175 200 575 1368 1629 1980 2341 2538 8 283 113 92 36 TPE (the best solution found so far) Tabu solution Stopping point
	Fig. 2:

ACKNOWLEDGMENT

This work was supported by the French National Research Agency under the ARTEFaCT (ANR-RF-2015-01) and Ade-quateDL (ANR-18-CE23-0012-01) projects, and by Inria. The authors also want to thank Tomofumi Yuki for his support and valuable feedback in applying TPE for WLO and designing experiments to construct the Hybrid approach.


Martin J Gander

Lukáš Jakabčin

Michal Outrata

The Schur complement on a bounded domain is a spectral Padé approximation about infinity of the Schur complement on the unbounded domain

We show for a specific model problem that the truncation of an unbounded domain by an artificial Dirichlet boundary condition placed far away from the domain of interest is equivalent to a specific absorbing boundary condition at the boundary of the domain of interest. In particular, using Schur complement techniques, we prove that the absorbing boundary condition obtained is a spectral Padé approximation about infinity of the transparent boundary condition. We also study numerically two improvements for this boundary condition -the truncation with an artificial Robin boundary condition placed far away from the domain of interest, and a Padé approximation about different point than infinity. Both of these give new and substantially better results compared to the artificial Dirichlet boundary condition.

Seen through the optic of linear algebra, we show that the Schur complement of our model problem written with respect to the eigenbasis can be identified with a truncation of a certain continued fraction. We use the theory of continued fractions to establish an approximation result of this truncation and hence interpreting the Schur complement as the Padé approximation of the optimal boundary operator in the eigenbasis. We then look to further improve the approximation qualities by changing some of the structure of the continued fraction so that the approximation is more accurate around a point of our choice and propose two different ways of achieving this result.

Introduction

In order to numerically solve a problem on an unbounded domain, we need to truncate the domain to a finite size to perform computations. This domain truncation problem was first studied in [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF], where the authors introduced the so called absorbing boundary conditions (ABC) for wave propagation phenomena, see also [START_REF] Bayliss | Radiation boundary conditions for wave-like equations[END_REF]. A second major technique for the truncation of infinite domains are the perfectly matched layers (PML), presented in the seminal paper [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. The pole condition, introduced in [START_REF] Schmidt | Pole condition: A numerical method for Helmholtz-type scattering problems with inhomogeneous exterior domain[END_REF], is yet another new way to construct and study domain truncations and in [START_REF] Gander | The Pole condition: A Padé approximation of the Dirichlet to Neumann operator[END_REF][START_REF] Gander | On the relationship between the pole condition, absorbing boundary conditions and perfectly matched layers[END_REF] the authors relate the pole condition technique, ABC and PML, obtaining an important understanding of all of these techniques.

At the discrete level, the ABC and PML techniques can be identified with techniques approximating the Schur complement in some sense. A number of iterative solvers has been derived based on this connection, see, e.g., [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods[END_REF] and the references therein. Our approach builds upon the eigendecomposition of the Schur complement, which for our model problem is very closely linked with the Fourier analysis of the Schur complement or, equivalently, the frequency domain analysis. Notably, the question of the optimal PML for problems with finite difference grids has been discussed in [START_REF] Ingerman | Optimal finite difference grids and rational approximations of the square root : I. Elliptic problems[END_REF][START_REF] Asvadurov | On optimal finitedifference approximation of PML[END_REF] for the Laplace equation and then also extended to the Helmholtz equation in [START_REF] Druskin | Near-optimal perfectly matched layers for indefinite Helmholtz problems[END_REF]. Our results go in a similar direction but are qualitatively different.

Domain truncation is also important in domain decomposition where a given computational domain is decomposed into many smaller subdomains, and then subdomain solutions are computed independently in parallel. An iteration process is used to obtain better and better approximations of the true solution on the entire domain. The solutions on the smaller subdomains can naturally be interpreted as solutions on truncated domains, and thus it is of interest to use ABC or PML techniques at the interfaces between the subdomains to enhance the convergence. Based on the insight from [START_REF] Nataf | Optimal interface conditions for domain decomposition methods[END_REF] this led to a new class of optimal and optimized Schwarz methods, see [START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | Schwarz methods over the course of time[END_REF] and references therein, and also the review on Schwarz methods [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods[END_REF].

The classical Schwarz method [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF] uses Dirichlet transmission conditions between subdomains and an overlap to achieve convergence [START_REF] Toselli | Domain Decomposition Methods -Algorithms and Theory[END_REF]. In what follows the goal is to interpret the overlap as a specific ABC once the unknowns of the overlap are folded onto the interface, similarly to the patch method in [START_REF] Magoulès | Algebraic approximation of Dirichlet-to-Neumann maps for the equations of linear elasticity[END_REF][START_REF] Gander | Analysis of patch substructuring methods[END_REF]. In this sense, the overlapping domain decomposition method is shown to be equivalent to a non-overlapping domain decomposition method with a particular ABC transmission condition. Although the Schwarz method is not explicitly mentioned in what follows, it is one of the main applications for our results.

We start in Section 2 by defining the model problem and its discrete counterpart. We recall the notion of the Schur complement for this problem in Section 3. Section 4 contains the main theoretical results of this paper, obtained by spectral analysis: we show that there exists a limit of the Schur complement as the width of the overlap goes to infinity, and that the Schur complement of a finite width truncation with a Dirichlet condition is a spectral Padé approximation around infinity of the unbounded width limit. We next explore numerically how the spectral approximation changes when the Dirichlet condition is replaced by a Robin condition in Section 5, and present an optimized choice for the Robin parameter, and also a generalization in volume. We give concluding remarks and discuss possible extensions in Section 7.

Model Problem

We use as our model problem the partial differential equation (PDE)

(η -∆)u = f in Ω, η > 0, u = 0 on ∂Ω, (1)
where the domain Ω := (0, +∞) × (0, 1), see Figure 1. We assume that the support of the right-hand side function f is localized, namely suppf ⊂ Ω a := (0, a) × (0, 1) ⊂ Ω.

Let b ≥ a and set Ω b := (0, b) × (0, 1) ⊂ Ω as a larger region, containing Ω a . Solving problem [START_REF] Asvadurov | On optimal finitedifference approximation of PML[END_REF] on Ω b and with homogeneous Dirichlet boundary condition at x = b (i.e., u(b, y) = 0) gives a discrete approximation u b of the true solution u on the unbounded domain Ω. Increasing b, the error of u b as an approximation of u is decreasing. Using a finite difference discretization of [START_REF] Asvadurov | On optimal finitedifference approximation of PML[END_REF] on Ω b with an equidistant mesh with mesh size h := 1 N +1 , the unknowns approximate the function values at the points of the grid. We can simplify the notation by gathering the unknowns with identical x coordinate into a vector of unknowns of length N . Assuming that the boundaries coincide with the mesh points, i.e., a = (N a + 1)h and b = (N b + 1)h,

we obtain the system of linear equations

A b u b = f b , (3)
where

f b = [f T 1 , • • • , f T N a , 0 T , • • • , 0 T ] T , with f i ∈ R N for i = 1, 2, .
. . , N a and 0 ∈ R N . The system matrix A b can be written in the classical block tridiagonal form

A b = 1 h 2             D 1 -I N -I N . . . . . . . . . D N a -I N -I N D N a +1 . . . . . . . . . -I N -I N D N b             ∈ R (N •N b )×(N •N b ) , (4)
with I N the identity matrix of size N × N and

D i = D =    ηh 2 + 4 -1 -1 . . . -1 -1 ηh 2 + 4    ∈ R N ×N , i = 1, 2, . . . , N b . (5)
In this notation, we can also identify the problem only on Ω a , i.e.,

A a u a = f a , (6)
with

f a = [f T 1 , . . . , f T N a ] T and A a = 1 h 2       D 1 -I N -I N . . . . . . . . . . . . -I N -I N D N a       . (7)
Finally, we introduce also the discretization of the unbounded problem without domain truncation, i.e.,

Au = f , (8)
where the solution as well as the right-hand side are vectors of infinite size, i.e.,

f = [f T 1 , • • • , f T N a , 0 T , • • • , 0 T , . . . ] T , u = [u T 1 , • • • , u T N a , u T N a +1 , • • • , u T N b , .
. . ] T , and the system matrix A is also of infinite dimension, i.e.,

A = 1 h 2                 D 1 -I N -I N . . . . . . . . . D N a -I N -I N D N a +1 . . . . . . . . . -I N -I N D N b . . . . . . . . .                 .
In what follows, it is enough to understand this infinite system simply as the limit of (4) as b → +∞; for more details on infinite matrices, see, e.g., the historical overview [START_REF] Bernkopf | A history of infinite matrices[END_REF].

Because the data (the function f ) vanishes outside of Ω a we can formulate two new problems on Ω a so that their solutions coincide with the solution of the problem (3) on Ω b and with the solution of the problem [START_REF] Gander | Optimized Schwarz methods[END_REF] on Ω. This can be done both at the continuous level, using the Dirichlet-to-Neumann map (see, e.g., [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods[END_REF]Section 5.2]), and at the discrete level using the Schur complement. Here we choose the latter and introduce the Schur complement in the following section in more detail.

The Schur complement operator

To reduce the system (3) to a smaller one of the dimension of ( 6), we eliminate the variables (u b N a +1 , . . . , u b N b ). Since those satisfy the equations in (3), i.e.,

- u b N b -1 h 2 + D N b u b N b h 2 = 0, - u b i-1 h 2 + D i u b i h 2 - u b i+1 h 2 = 0, for i = N b -1, . . . , N a + 1, (9)
we can recursively eliminate them. This is the well-known process of block Gaussian elimination and in the matrix formulation is known as the Schur complement approach.

Definition 3.1 (Schur complement) The Schur complement T b N a is defined recursively by

T b N b = D N b h 2 = D h 2 , T b i = D i h 2 - (T b i+1 ) -1 h 4 = D h 2 - (T b i+1 ) -1 h 4 , for i = N b -1, . . . , N a . ( 10
)
Using the Schur complement T b N a and the particular zero structure of the right-hand side in (3), we can reduce the problem (3) to a problem on Ω a , namely

Ãa ũa = f a , (11)
where

Ãa =       D 1 h 2 -I N h 2 -I N h 2 . . . . . . . . . D N a -1 h 2 -I N h 2 -I N h 2 T b N a       . ( 12
)
The only difference between ( 6) and ( 11) is the last block in [START_REF] Gander | On the relationship between the pole condition, absorbing boundary conditions and perfectly matched layers[END_REF] where the Dirichlet boundary condition block has been replaced by the Schur complement T b N a , representing the "far-field" domain (or the overlap) unknowns in Ω b \Ω a that have been folded in. Since the solution ũa of the modified system [START_REF] Gander | The Pole condition: A Padé approximation of the Dirichlet to Neumann operator[END_REF] coincides with the solution ũb of (3) restricted to the points in Ω a , the smaller system gives a better approximation of the underlying solution u on the infinite domain compared to the solution of [START_REF] Druskin | Near-optimal perfectly matched layers for indefinite Helmholtz problems[END_REF].

To further improve the solution, we can move the boundary point b to the right to some b > b (and thus having an N b > N b ) and hence, the Schur complement matrix T b N a will be defined by a longer recurrence (see Definition 3.1). If b goes to infinity, the corresponding Schur complement matrix T ∞ N a will be governed by the limit of ( 10), namely

T ∞ N a = D h 2 - (T ∞ N a ) -1 h 4 , (13)
which does not depend on N a . Thus, T ∞ N a ≡ T ∞ for any N a , where T ∞ satisfies the same equation ( 13), i.e., one has

T ∞ = D h 2 - (T ∞ ) -1 h 4 . ( 14
)
Using T ∞ instead of T b N a in [START_REF] Gander | On the relationship between the pole condition, absorbing boundary conditions and perfectly matched layers[END_REF], the corresponding solution ũa will coincide with u on the points in Ω a . Equation ( 14) can be reformulated into a quadratic matrix equation for T ∞ given by

(T ∞ ) 2 - D h 2 T ∞ + I h 4 = 0. ( 15
)
This equation has two solutions, each of them corresponding to a different underlying solution u. But only one of these corresponds to a bounded solution u, which is of interest. In order to solve [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF], it is convenient to change the basis we work in to the eigenbasis of D, which, in this case, corresponds to the discrete sine basis1 . We present the change of the basis and its implications for the computation in the following section.

Spectral analysis

The eigenpairs of the matrix D in (5) can be evaluated by writing D = D yy + 2I, where D yy is the 3-point finite difference stencil discretization of η -∂ yy multiplied by h 2 . The eigenpairs of D yy are known in closed form, i.e., D yy = Q T ZQ with

Z := diag(z 1 , . . . , z N ), z k := ηh 2 + 4 sin 2 kπ 2(N + 1) , (16)
and Q unitary and symmetric, with the eigenvectors q k in its columns,

q k := 2 N sin kπ N +1 j N j=1 = 2 N + 1    sin kπ N +1 . . . sin kπ N +1 N    ∈ R N . (17)
Thus, we can diagonalize D with the same basis, i.e.,

D = Q T ΛQ
with eigenvalues given by Λ = diag(λ 1 , . . . , λ N ) and

λ k = 2 + z k . (18)
For the rest of the text we fix the index notation so that the index k = 1, . . . , N corresponds to a particular eigenmode (or frequency or Fourier mode in different terminologies).

Remark 1 Changing the basis so that the diagonal blocks of the matrices A, A b , A a and also Ãa diagonalize will on the one hand make the Schur complement analysis significantly easier -we will be able to treat each eigenmode separately -but on the other hand more complicated -treating each eigenmode separately adds yet another index to the already loaded notation.

Throughout the text we keep the index k reserved for the eigenmode notation.

Although it is in our eyes not possible to simplify the notation while avoiding confusion, we recognize the difficulties associated with the notation of the following subsection.

Diagonalization of the Schur Complement

Changing the basis for the Schur complement definition in [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods[END_REF] gives

T b N b = QT b N b Q T = Q D h 2 Q T = Λ h 2 , T b i = QT b i Q T = Q D h 2 Q T -Q (T b i+1 ) -1 h 4 Q T = Λ h 2 - ( T b i+1 ) -1 h 4 , (19)
where i = N b + 1, . . . , N a and all of the matrices T b i are diagonal. Working with the diagonal entries only, each mode (frequency) also follows the analogous recurrence, namely

tb N b ,k = λ k h 2 , tb i,k = λ k h 2 - 1 h 4 tb i+1,k , for i = N b + 1, . . . , N a . (20)
Moreover, we can also write an analogue of the above recurrence for the recurrence we establish for the solution u b in [START_REF] Gander | Schwarz methods over the course of time[END_REF]. In particular, having

- u b N b -1 h 2 + Q T ΛQu b N b h 2 = - u b N b -1 h 2 + Q T T b N b Qu b N b = 0, - u b i-1 h 2 + Q T ΛQu b i h 2 - u b i+1 h 2 = - u b i-1 h 2 + Q T T b i Qu b i h 2 = 0, (21)
with i = N b -1, . . . , N a , setting ûb i := Qu b i and multiplying [START_REF] Toselli | Domain Decomposition Methods -Algorithms and Theory[END_REF] on the left by Q, we obtain

- ûb N b -1 h 2 + Λû b N b h 2 = - ûb N b -1 h 2 + T b N b ûb N b = 0, - ûb i-1 h 2 + Λû b i h 2 - ûb i+1 h 2 = - ûb i-1 h 2 + T b i ûb i h 2 = 0, (22)
where i = N b + 1, . . . , N a .

In the following section we use the above recurrences to analyze the the Schur complement T b N a and its limit as b → +∞. Notice that in the context of Section 1, obtaining the limit lim b→+∞ T b N a allows us to construct ũa such that it coincides with u (restricted to Ω a unknowns) and thus to obtain the best possible approximation in this framework.

Convergence of the Schur Complement

Focusing on the limit case lim b→+∞ T b N a , (19 -20) implies that we can compute the limit for each eigenmode independently, meaning we can focus on N scalar limits

lim b→+∞ tb N a ,k =: t∞ N a ,k , (23)
for k = 1, . . . , N . Since [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF] holds for any b fixed, the limit t∞ N a ,k satisfies the same equation, meaning we have

( t∞ N a ,k ) 2 - λ k h 2 t∞ N a ,k + 1 h 4 = 0. ( 24
)
Since the equation does not depend on N a , neither will the solutions, i.e., the limits, and thus we can omit the N a subscript in (23), obtaining t∞ N a ,k ≡ t∞ k . Returning to equation (24), there are two solutions,

τ ∞,1 k = λ k + λ 2 k -4 2h 2 and τ ∞,2 k = λ k -λ 2 k -4 2h 2 , (25)
and using the definition of λ k in [START_REF] Nataf | Optimal interface conditions for domain decomposition methods[END_REF], both solutions are real and positive for all k = 1, . . . , N .

Taking any k and using the Vieta formulas for (24), we have

(h 2 τ ∞,1 k )(h 2 τ ∞,2 k ) = 1, (26)
and, moreover,

0 < h 2 τ ∞,2 k < 1 < h 2 τ ∞,1 k . ( 27
)
The next step is to show that one of the solutions τ ∞,1 k , τ ∞,2 k indeed acts as the limit Schur complement for our solution vector ũa .

The key observation is that the characteristic polynomial of the recurrence relation in ( 22) is preserved through the limit process and thus the solutions τ ∞,1 k , τ ∞,2 k of the limit equation (24) coincide with the roots of the characteristic polynomial

p k (r) = -r 2 + λ k r -1.
This together with the explicit formula for the solution of the recurrence relation ( 22) is enough to establish the following result. N a defined in (10) converges to T ∞,1 , i.e., the solution of the formal limit equation [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF], as b → +∞. To be more specific, the eigenvectors of those matrices are equal and the eigenvalues of T b N a converge to the ones of T ∞,1 , i.e.,

t∞ k ≡ lim N b →∞ tb N a ,k = τ ∞,1 k = λ k + λ 2 k -4 2h 2 , ( 28
)
where

λ k = ηh 2 + 2 + 4 sin 2 hkπ 2 , see (18)
.

Proof Fixing b, we take the solution subvector ûb i of length N for any i = N a , . . . , N b and denote its scalar entries by ûb i,k for k = 1, . . . , N . These entries follow the system of difference equations in ( 22) and as such they can be written as a linear combination of powers of the roots of the characteristic polynomial of the difference equation. In other words, there exist two constants ν b k and µ b k such that the solution ûb i,k of ( 22) is given by

ûb i,k = µ b k h 2 τ ∞,1 k i-N a + ν b k h 2 τ ∞,2 k i-N a , ( 29
)
where the constants

µ b k , ν b k depend on N b and N a . Recalling the inequality of the roots τ ∞,1 k , τ ∞,2 k in (27) it follows that h 2 τ ∞,1 k N b -N a → +∞.
Since we look for the limit of solutions satisfying a homogeneous Dirichlet boundary condition at x = b, the limit has to decay at +∞. Thus the bounded solution of the form (29) has to satisfy µ b k → 0 as b → +∞. Hence, setting ν ∞ k as the limit of ν b k as b → +∞ we obtain

û∞ i,k = ν ∞ k h 2 τ ∞,2 k i-N a . ( 30
)
Returning to the system of difference equations [START_REF] Wall | Analytic Theory of Continued Fractions[END_REF], the scalar equations become

-û b N a ,k + λ k ûb N a +1,k -ûb N a +2,k = -û b N a ,k + h 2 tb N a +1,k ûb N a +1,k = 0,
and solving for h 2 tb N a +1,k , we get

h 2 tb N a +1,k = ûb N a ,k ûb N a +1,k
.

Inserting the values of ûb i,k from (30), we find lim

N b →∞ h 2 tb N a +1,k = lim N b →∞ ûb N a ,k lim N b →∞ ûb N a +1,k = û∞ N a ,k û∞ N a +1,k = ν ∞ k h 2 τ ∞,2 k i-N a ν ∞ k h 2 τ ∞,2 k i-N a +1 = 1 h 2 τ ∞,2 k = h 2 τ ∞,1 k ,
where the last step follows from (26), finishing the proof.

Remark 2 All of the computations starting at the beginning of this subsection with (23) and finishing by Theorem 4.1 were focusing on computations for the full Schur complement, i.e., the full recurrence in [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF]. However, starting in (23) with N a + 1, N a + 2 or any fixed N a +n, all of the computations remain the same. This confirms that indeed the limit result is independent of N a as the same result holds also if we stop the recurrence elimination process one, two or any other fixed number of steps earlier (that is before reaching the domain of interest Ω a ).

The limit t∞

k can be written in a more convenient form, using the notation λ k = z k + 2 from [START_REF] Nataf | Optimal interface conditions for domain decomposition methods[END_REF]. A direct computation gives

t∞ k = λ k + λ 2 k -4 2h 2 = 2 + z k + (2 + z k ) 2 -4 2h 2 = 1 h 2 + z k 2h 2 + z 2 k + 4z k 2h 2 = 1 h 2 + z k 2h 2 + z 2 k (1 + 4 zk ) 2h 2 = 1 h 2 + z k 2h 2 + z k 2h 2 1 + 4 z k , (31)
and thus, for any k = 1, .., N , the limit t∞ k can be written in the form

t∞ k = 1 h 2 1 + z k 2 1 + 1 + 4 z k , (32)
with z k := ηh 2 + 4 sin 2 (hkπ/2). Each limit mode t∞ k can thus be seen as a function of an argument z,

t∞ (z) = 1 h 2 1 + z 2 1 + 1 + 4 z , (33)
evaluated at z = z k ∈ [ηh 2 + 4 sin 2 (hπ/2), ηh 2 + 4 sin 2 (hN π/2)].
In the same way, the recurrence relation [START_REF] Schmidt | Pole condition: A numerical method for Helmholtz-type scattering problems with inhomogeneous exterior domain[END_REF] for tb i,k can be interpreted as a function evaluation at the points

z = z k , tb i,k = tb i (z k ) for i = N a , . . . , N b . The first three functions are tb N b (z) = 2 + z h 2 , tb N b -1 (z) = 2 + z h 2 - 1 h 4 2+z h 2 = 1 h 2 2 + z - 1 2 + z , tb N b -2 (z) = 2 + z h 2 - 1 h 4 tb N b -1 (z) = 2 + z h 2 - 1 h 4 1 h 2 2 + z -1 2+z = 1 h 2 2 + z - 1 2 + z -1 2+z ,
and by the recursive definition in [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF], the general term is

tb i (z) = 2 + z h 2 - 1 h 2 2 + z - 1 2+z- . . . 2+z-1 2+z . ( 34
)
Objects of this form are called continued fractions. Their theory links various areas of mathematics, e.g., Padé approximations, orthogonal polynomials, Vorobyev's moment matching problem, Gauss quadrature and the method of conjugate gradients (see [START_REF] Lorentzen | Continued Fractions with Applications[END_REF] and also [START_REF] Liesen | Krylov Subspace Methods: Principles and Analysis[END_REF] for further references) and we will use some of these links to establish our main result in Section 4.4 later on. In the light of this observation we adjust our notation, clarifying everything in Remark 3 below.

Remark 3 Notice that in the continued fraction representation of tb i (z) in (34), the continued fraction has exactly N b -i levels. In order to simplify the notation, we will from now on change the subscript i to correspond to the "number of levels" or "depth" of the continued fraction. Hence, for the rest of the text we will write

tb 0 (z) = 2 + z h 2 , tb 1 (z) = 2 + z h 2 - 1 h 4 2+z h 2 = 1 h 2 2 + z - 1 2 + z , tb 2 (z) = 2 + z h 2 - 1 h 4 tb N b -1 (z) = 2 + z h 2 - 1 h 4 1 h 2 2 + z -1 2+z = 1 h 2 2 + z - 1 2 + z -1 2+z
, and so on. This means that the index i changes the meaning from the number of grid columns in the domain Ω b to the number of grid columns in the domain Ω b \ Ω a . Whenever i increases, we understand it as extension of the domain Ω b , rather than the mesh size changing.

We continue by a simple observation regarding the functions t∞ and tb i .

Remark 4 By a direct computation we obtain

t∞ (z) = 2 + z - 1 t∞ (z) ,
and by re-insertion we also get

t∞ (z) = 2 + z - 1 2 + z -1 t∞ (z)
, and so on. This suggests that the function t∞ (z) is equal to the infinite continued fraction

t∞ (z) = 2 + z - 1 2 + z -1 2+z-...

and the functions tb

i (z) in (34) are approximations in the sense of a truncation of this infinite continued fraction after i levels.

Since the continued fractions are not the primary focus, which is to say that we will use the continued fractions only as a tool to arrive at our main result, we choose [START_REF] Baker | Padé Approximants Part I: Basic theory[END_REF] as the main reference, which is a book written as an overview for Padé approximations and the continued fractions are approached mainly from that perspective. We refer the interested readers to [START_REF] Lorentzen | Continued Fractions with Applications[END_REF] and [START_REF] Wall | Analytic Theory of Continued Fractions[END_REF] for more detailed expositions of the theory of continued fractions. We continue with Section 4.3 where we recall some standard terminology and results of the field and give some auxiliary lemmas in the following subsection.

Padé Approximation and Continued Fractions

We start by recalling the Padé approximation theory. We follow the notation from [START_REF] Baker | Padé Approximants Part I: Basic theory[END_REF], i.e., the

[M/L]-Padé approximant of f (z) is denoted by [M/L] f ≡ [M/L] f (z).
1. Let α, β ∈ R. Then α + β[M/L] f = [M/L] α+βf .
2. Let m ≥ 1 be fixed and assume f (z) = +∞ j=0 c j z j to be a formal power series (we do not consider the convergence question here). Setting

g(z) = 1 z m f (z) - m-1 j=0 c j z j and assuming M -m ≥ L -1 we have [M -m/L] g (z) = 1 z m [M/L] f (z) - m-1 j=0 c j z j . 3. Assume f (0) = 0 and set g(z) = 1/f (z). Then [M/L] g (z) = 1/[L/M ] f (z).
We continue by introducing the basic terminology of continued fractions following [2, Chapter 4]. Definition 4.3 A continued fraction is given by sequences of real numbers {a j } j , {b j } jthe numerator and the denominator sequence of the continued fraction -and has the general form

b 0 + a 1 b 1 + a 2 b 2 + a 3 . . . =: b 0 + +∞ j=1 a j b j + ≡ b 0 + a 1 b 1 + a 2 b 2 + . . . ,
where the sum is to be understood only formally. The continued fraction is called infinite as long as a j , b j = 0 for all j.

The n-th truncation (or convergent) of a continued fraction is given by

A n B n = b 0 + n j=1 a j b j + = b 0 + a 1 b 1 + a 2 b n-2 + . . . b n-1 + an bn
, where A n and B n are the n-th truncation (or convergent) numerator and denominator.

Replacing the scalars a j and/or b j by linear (or affine) functions of a real variable z, A n and B n become polynomials in z and the n-th truncation of the continued fraction becomes a rational function in z. Different settings of this framework lead to different types of continued fractions. Most notably, a continued fraction is called regular C-fraction (short for regular classical continued fraction), provided it has the form

b 0 + a 1 z 1 + a 2 z 1+ a 3 z . . . ≡ b 0 + a 1 z 1 + a 2 z 1 + . . . ,
with a j = 0 for all j. If, moreover, a j > 0 for all j, then it is called S-fraction (short for Stieltjes continued fraction). If the continued fraction takes the form

b 0 + r 1 z + s 1 - r 2 z+s 2 - r 3 . . . ≡ b 0 + r 1 z + s 1 - r 2 z + s 2 - . . . ,
with k j = 0 for all j then it is called J-fraction (short for Jacobi continued fraction).

Next, we give some remarks on the definition above. First, let us emphasize that we have ignored the questions of convergence of infinite continued fractions and we refer the reader to [START_REF] Lorentzen | Continued Fractions with Applications[END_REF] and [START_REF] Wall | Analytic Theory of Continued Fractions[END_REF]. Also, notice that one function can be represented by two seemingly different continued fractions (different possibly in type and/or in the coefficient values) and one way to recognize the equality of two continued fractions is via the three-term recurrence relation of the numerators and denominators of the continued fraction truncations (convergents), see [2, Theorem 4.1.1, pp.106]. We have that

A -1 = 1, A 0 = b 0 , A n = b n A n-1 + a n A n-2 , (35) B -1 = 0, B 0 = 1, B n = b n B n-1 + a n B n-2 .
and assuming the n-th truncation (convergent) of two continued fractions are equal for any n, the infinite continued fractions are equal as well. For more details on the introduced types of continued fractions as well as other types of continued fractions (e.g., non-regular C-fraction, T-fraction, P-fraction,. . . ) we refer to [START_REF] Lorentzen | Continued Fractions with Applications[END_REF] and [START_REF] Wall | Analytic Theory of Continued Fractions[END_REF] and references therein. Last but not least, we want to note that some authors will call a continued fraction an S-fraction even though the fraction itself does not meet the definition above but can be transformed into a continued fraction that does. We next recall a basic transformation rule of continued fractions. For purposes of this text we present immediately the continued fraction result for the square root function, which is of interest to us 2 . We state the result in Theorem 4.5, referencing to the book of Baker but the original result is due to Gauss, who showed a much more general result for the hypergeometric function 2 F 1 ; for more details derivation we refer the reader to [22, Chapter XVIII] or [START_REF] Lorentzen | Continued Fractions with Applications[END_REF]Chapter VI].

Theorem 4.5 ([2, Section 4.6, Theorem 4.4.3 and formula (6.4) on pp. 139]) For any α ∈ (-1, +∞) 3 we have

√ 1 + α = 1 + α 2 1 + α 2 2+ α 2 1+ α 2 b n-2 + . . . b n-1 + an bn+ a n+1 ... = 1 + α 2 1 + α 2 2 + α 2 1 + . . . + a n b n + . . . , (36)
where the denominator sequence is given by b

0 = 1, b j = 3+(-1) j 2
and the numerator sequence is given by a j = α 2 , j ≥ 1. Moreover, for any n the [n, n]-Padé approximation of √ 1 + α expanded about α = 0 is given by the (2n)-th truncation of the continued fraction in (36) and the [n + 1, n]-Padé approximation of √ 1 + α expanded about α = 0 is given by the (2n + 1)-st truncation of the continued fraction in (36).

Remark 5 By a direct computation we see that

√ 1 + α = 1 + α 2 + α 2 + α 2 + . . . ,
and thus the above representation in (36) can be equivalently written as a cyclic S-fraction 4 with a j = 1/2 for all j.

We finish this subsection by proving some preparatory results, the first of which will be useful in linking a truncation of the S-fraction introduced in Theorem 4.5 and a truncation of the J-fraction from Remark 4. However, notice that the continued fraction considered there is not identical with the one in (36) -they obey the same recurrence but they differ at the beginning. Lemma 4.6 Let α be real and consider the continued fraction 2 Although we do not use them, there are many ways to create a continued fraction representation of a function based on its formal power series. We do not consider these here in more detail but rather refer the interested reader to [2, Section 4. 3 There is a misprint in [2, equation (6.4), page 139]. The authors state the convergence "for all z except -∞ < z ≤ 1)" but the the result also holds for z ∈ (-1, 1]. 4 Infinite continued fractions such that the sequences {a j }, {b j } are periodic are called cyclic continued fractions.

τ (α) := α 2 2 + α 2 1+ α 2 2+ α 2 1+...
and denote its n-th truncation by A n (α)/B n (α). Considering the J-fraction

σ(α) := 1 1 + 4 α - 1 2+ 4 α - 1 2+ 4 α - 1 2+ 4 α -...
with n-th truncation C n (α)/D n (α) we have

A 2n (α)/B 2n (α) = C n (α)/D n (α)
for any n = 1, 2, . . . .

Proof

We start by transforming the continued fraction τ by the rules of Lemma 4.4 and without further relabeling we obtain

τ (α) := 1 4 α + 1 1+ 1 4 α + 1 1+ 1 4 α + 1 1+... . (37)
First, by a direct computation, we confirm the equality for n = 1, obtaining

1 4 α + 1 1 = 1 4 α + 1
, and next, we notice that the continued fraction (37) can be written in a cyclic form with the core R given by

R = 4 α + 1 1 + 1 R . ( 38
)
That is, the continued fraction can be obtained by a successive re-insertion of the core equality (38) into itself, e.g.,

1 4 α + 1 1 = A 2 (α) B 2 (α) , 1 4 α + 1 1+ 1 4 α + 1 1 = A 4 (α) B 4 (α)
,

1 4 α + 1 1+ 1 4 α + 1 1+ 1 4 α + 1 1 = A 6 (α) B 6 (α)
, . . . .

Notice that in this way every re-insertion adds two elements of the numerator and denominator sequences as highlighted by the indices of the truncations. Using the algebraic identity

1 1 + 1 R = 1 - 1 1 + R ,
we reformulate the core equality (38) to obtain

R = 4 α + 1 - 1 1 + R ,
or, more conveniently,

1 + R = 2 + 4 α - 1 1 + R . ( 39
)
Notice that using the algebraic identity above we contracted the two-level core that is to be re-inserted to only a one-level core, while not changing the resulting value of the fraction. Moreover, the core equality in (39) is the one that generates the J-fraction σ(α).

Hence we have shown that for any n ≥ 2 the 2n re-insertions of the core R in the core equality (38) is equal to only n re-insertions of the core 1 + R in the equality (39). It follows that the 2n-th convergent A 2n (z)/B 2n (z) is equal to the n-th convergent C n (z)/D n (z), yielding the result.

We build upon Lemma 4.6 with the following result that contracts the S-fraction in (36) into a J-fraction. Proposition 4.7 Let α be real and set the continued fractions τ (α) and σ(α) as in Lemma 4.6. Moreover, we define the continued fractions

τ (α) := 1 1 + τ (α)
and φ(α

) := 1 - 1 2 + 4 α - 1 2+ 4 α - 1 2+ 4 α -...
with n-th truncations Ãn (α)/ Bn (α) and E n (α)/F n (α) respectively. Then

A 2n+1 (α)/B 2n+1 (α) = E n (α)/F n (α)
for any n = 0, 1, 2, . . . with E 0 = F 0 = 1.

Proof The equality for n = 0 holds by inspection. Taking n ≥ 1, we use Lemma 4.6 for the continued fraction τ (α) and we obtain

Ã2n+1 (α)/ B2n+1 (α) = 1 1 + A 2n (α)/B 2n (α) = 1 1 + C n (α)/D n (α)
and it remains to show that

1 1 + C n (α)/D n (α) = 1 -E n (α)/F n (α), (40)
where C n (α), D n (α) are the truncations of the J-fraction σ(α) from Lemma 4.6. We first notice that the cyclic parts of both J-fractions σ(α) and φ(α) coincide and we denote them by σ(α),

σ(α) := 1 2 + 4 α -1 2+ 4 α -... (41)
We then notice that

σ(α) = 1 1 + 1 1+ 4 α -σ(α)
, and

φ(α) = 1 - 1 2 + 4 α -σ(α)
.

Therefore, in order to show (40) it is enough to show that

1 1 + 1 1+ 4 α -σ = 1 - 1 2 + 4 α - σ ,
as σ contains the common part. By a direct computation we obtain

1 1 + 1 1+ 4 α -σ = 1 + 4 α -σ(α) 2 + 4 α -σ(α) , and
1 - 1 2 + 4 α - σ = 1 + 4 α -σ(α) 2 + 4 α -σ(α)
, finishing the proof.

Approximation Properties of the Schur Complement

We now show that the function tb i (z) representing the Schur complements T b i is a Padé approximation about z = +∞ of the function t∞ (z) representing the infinite Schur complement T ∞ . To obtain this result in Theorem 4.8 we use a similar technique as in [START_REF] Gander | The Pole condition: A Padé approximation of the Dirichlet to Neumann operator[END_REF] where the authors compute a Padé approximation of the Dirichlet to Neumann operator. This similarity is not a coincidence: the Schur complement and the Dirichlet-to-Neumann map have a very deep connection, see, e.g., [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized schwarz methods[END_REF]Section 5.2].

Also, in the proof and in the rest of the document we change the variables back to the ones from our application, i.e.,

α = 4 z . (42)
Notice that the term change of variables here is somewhat misleading. Usually, change of variables in approximation theory then requires a re-computation of the approximationbecause of the way derivation of a composite function works. However, our result considers expansion about +∞, which is defined by considering the expansion of the same function about zero but of the reciprocal argument, e.g., not of z but of 1/z. Recalling this convention, our change of variables in fact consists only of multiplying by 4. Also, notice that such change of variables does not require a re-computation of the derivatives.

Theorem 4.8 The function tb i (z) defined in (34) as

tb i (z) = 1 h 2      2 + z - 1 2 + z - 1 2+z- . . . 2+z-1 2+z     
is the [i, i]-Padé approximation about the expansion point z = +∞ of the Schur complement function t∞ (z) defined in (33) as

t∞ (z) = 1 h 2 1 + z 2 + z 2 1 + 4 z .
Proof First, we drop the 1/h 2 factor for both of the functions and transpose the expansion point z = +∞ to α = 0 as in (42) and without a further relabeling of the functions we

obtain tb i (α) = 2 + 4 α - 1 2 + 4 α - 1 2+ 4 α - . . . 2+ 4 α -1 2+ 4 α , t∞ (α) = 1 + 2 α + 2 α √ 1 + α.
Recalling (26), we have

{ t∞ } -1 (α) := 1 t∞ = 1 + 2 α - 2 α √ 1 + α (43)
and thus by Theorem 4.2 part 3., we have

[i/i] t∞ (α) = 1 [i/i] { t∞ } -1 (α)
for any i ≥ 1. By a direct computation we obtain

{ t∞ } -1 (α) = 1 + 2 α - 2 α √ 1 + α = 1 -2 1 α √ 1 + α -1 ,
and by application of Theorem 4.2 parts 1. and 2. we obtain

[i/i] { t∞ } -1 (α) = 1 -2 1 α [i + 1/i] √ 1+α (α) -1 .
Using the continued fraction representation of the Padé approximant from Theorem 4.5, we obtain

[i/i] { t∞ } -1 (α) = 1 - 2 α 1 + 1 + A 2i+1 (α) B 2i+1 (α) -1 = 1 - 2 α            α 2 1 + α 2 2+ α 2 1+ α 2 b 2i-1 + . . . b 2i + a 2i+1 b 2i+1           
, where A 2i+1 (α), B 2i+1 (α) are the (2i + 1)-st truncation numerator and denominator of the continued fraction τ (α) in (36) and the sequences {a j } j , {b j } j are taken as in (36). Hence we have

[i/i] { t∞ } -1 (α) = 1 - 1 1 + α 2 2+ α 2 1+ α 2 b 2i-1 + . . . b 2i + a 2i+1 b 2i+1
, and recalling Proposition 4.7 we that the continued fraction is the (2i+1)-st truncation of the continued fraction τ (α) defined there. Therefore, Proposition 4.7 gives that

[i/i] { t∞ } -1 (α) = 1 -1 - Cn (α) Dn (α) = 1 2 + 4 α - 1 2+ 4 α - 1 2+ 4 α - . . . 2+ 4 α -1 2+ 4 α i-1 "levels"
, where Ci (α), Di (α) are the (2i+1)-st truncation numerator and denominator of the J-fraction σ(α) in (41) in Proposition 4.7. As a result we have that

[i/i] t∞ (α) = 1 1 2+ 4 α - 1 2+ 4 α - 1 2+ 4 α - . . . 2+ 4 α -1 2+ 4 α = 2 + 4 α - 1 2 + 4 α - 1 2+ 4 α - 1 2+ 4 α - . . . 2+ 4 α -1 2+ 4 α i-1 "levels"
for any i ≥ 1, finishing the proof. where i denotes the number of grid columns that were folded into the Schur complement, i.e., the number of layers between a and b, and thus also the degree of the Padé approximation. We show plots of the function err D (z, i) for small i in Figure 2.

As expected, we see that the error err D (z, i) decreases when more grid columns are added, i.e., when i and thus also b increases. Also, for any fixed i, the error is decreasing when z is increasing. This behavior was established for z → +∞ by Theorem 4.8.

On the other hand, the error is still quite large for small z. Recalling the role of z, this can be identified with poor quality of the approximation for low frequency modes following from our spectral analysis. Reducing the maximum of the error would yield a better approximation of the true solution u by the approximation ũa , see [START_REF] Gander | The Pole condition: A Padé approximation of the Dirichlet to Neumann operator[END_REF]. We explore this direction in the next section.

ηh 2 + 4 sin 2 π 2 1 N + 1 , ηh 2 + 4 sin 2 π 2 N N + 1 ≈ ηh 2 , ηh 2 + 4 , (44)
the expansion point z = +∞ of the Padé approximation is far away from the interval where we want to approximate t∞ , i.e., this expansion point is not well-chosen and far from optimal. An idea is to replace the Dirichlet boundary condition used to truncate the recurrence relation by a homogeneous Robin boundary condition6 with Robin parameter p ≥ 0 at b, i.e.,

∂u ∂n + pu = 0 at x = b.
The Dirichlet condition is the limit of the Robin condition as p → +∞. Taking any finite positive p will change the approximation error function in Figure 2 and can give a better approximation of t∞ (z) in the interval (44), as we show next by introducing a discretization. Using a centered finite difference approximation as before, the Robin condition can be discretized with the so-called ghost point trick: first, we discretize the Robin condition,

u N b +1 -u N b -1 2h + pu N b = 0 =⇒ u N b +1 = u N b -1 -2phu N b .
Then in the discretized equation at b, -u

N b +1 + D N b u N b -u N b -1 = 0,
Ãb = 1 h 2             D 1 -I N -I N . . . . . . . . . D N a -I N -I N D N a +1 . . . . . . . . . -I N -I N DN b             (45) with DN b = 1 2 (D N b + (2ph)I N )
. This also modifies the Schur complement, yielding

T b N b = DN b 2h 2 , with the diagonal entries tb N b ,k = λ k 2h 2 + p h , Ti = D i h 2 - T -1 i+1 h 4 = D h 2 - T -1 i+1 h 4 , for i = N b -1, . . . , N a . ( 46
)
The first three functions representing the diagonal entries are

tb 0 (z) = 1 + ph + z 2 h 2 , tb 1 (z) = 2 + z h 2 - 1 h 4 1+ph+ z 2 h 2 , = 1 h 2 2 + z - 1 1 + ph + z 2 , tb 2 (z) = 2 + z h 2 - 1 h 4 tb N b -1 (z) = 2 + z h 2 - 1 h 4 1 h 2 2 + z - 1 1+ph+ z 2 = 1 h 2 2 + z - 1 2 + z - 1 1+ph+ z 2 ,
where we keep the subscript notation (see Remark 3) and by the recursive definition in [START_REF] Schwarz | Über einen Grenzübergang durch alternierendes Verfahren[END_REF], the general form is

tb i (z) = 2 + z h 2 - 1 h 2 2 + z - 1 2+z- . . . 2+z- 1 1+ph+ z 2 . ( 47
)
Notice that by letting p → +∞ we can recover the original Dirichlet boundary condition with one less level of the continued fraction. In other words, having h fixed and the Robin parameter p large enough, we can interpret this as a Dirichlet boundary condition in a geometry where Ω b was extended by one additional grid column after b, i.e., by the strip (b, b + h) × (0, 1).

We see that the Robin condition just changes the last denominator in the continued fraction representation. With (47), we can numerically explore the effect of the Robin err_R(z_k,i) parameter p on the behavior of the error 7 err R (z, i) := t∞ (z) -tb i (z), where i denotes again the number of grid columns that were folded into the Schur complement, i.e., the number of layers added after the point a. Similarly to Figure 2 we show the evolution of the error function err R (z, i) in Figure 3.

i = 1 i = 2 i = 3 i = 4 i = 5 i = 6 i = 7
We see that the behavior around the right endpoint of the interval is still present but the Robin condition introduced another point z p around which the approximation is accurate. In Figure 3 we see that for z p ≈ 0.75. To confirm the hypothesis, we use a finer mesh in z in Figure 4. This confirms that the Robin parameter p affects the approximation error by minimizing it around a particular point z p . Shifting z p towards 0, i.e., towards the area where both 7 The subscript R stands for the "Robin" boundary condition at the end point x = b. for z with a given p. Assuming that err R (z, i) is smooth except at a finite number of points, which is the case based on the numerical experiments above, equation (48) defines z p as an implicit function of p and also the other parameters of the problem. For example, for i = 1 we get

err R (z, i) = 1 h 2 1 + z 2 1 + 1 + 4 z - 1 + ph + z 2 h 2 = 1 h 2 z 2 1 + 4 z -ph ,
which gives z p as the positive root of the quadratic equation

z 2 p + 4z p -4p 2 h 2 = 0 =⇒ z p = -2 + 2 1 + p 2 h 2 . ( 49
)
Based on the observation in Figure 4, varying the number of layers i does not substantially change the value of z p . In order to optimize the parameter p in the Robin condition, one can minimize the infinity norm of the error in Figures 3 and/or 4. Using an optimization routine to minimize the infinity norm of the error over p, the local maxima equioscillate, see Figure 5.

Here we did not use the natural values of z k but rather chose to span the entire interval with logarithmically equidistant points. We used N = 200, i = 5, η = 2, and the optimal i p * (i) value of p was found to be p = 7.167. The Dirichlet error in the infinity norm is equal to 0.1371 while the Robin error in the infinity norm is equal to 0.0245, an almost 5.6 fold improvement. Running the optimization while varying i, i.e., the number of grid columns from a to b we obtain Table 1, again for N = 200 and η = 2.

err D ∞ err R ∞ 1 27
This shows that the improvement over the Dirichlet truncation increases with increasing number of layers.

The corresponding results over a larger range of i are shown graphically in Figure 6. Again we minimized the maximum norm err R (z, i) ∞ over the interval in (44). We used again N = 200, η = 2 and varied i as powers of 2 from 2 1 = 2 to 28 = 256 on the left and then up to 2 15 on the right. One can see a linear dependence in the log-log scale on the left, i.e., for values i ≤ 256 and fitting the line gives the law p * (i) ∼ C • i q (50) with C ≈ 11 and q ≈ -1

for i ≤ 256. For practical purposes having i > 256 is not desirable and thus the above law basically covers a sufficient domain of i. However, Figure 6 shows that in general p * (i) does not follow the proposed relation (50).

Although we have achieved considerable improvement over the Dirichlet boundary condition, the Robin boundary condition still could be improved. Observing Figure 4 the Robin boundary condition error function still decreases for large z and it seems that it inherited most of the approximation qualities around z = +∞. This is, however, not useful for our application as we simply want to minimize the error only over the interval spanned by the eigenvalues, i.e., over [ηh 2 , ηh 2 + 4]. A possible way to do this is to shift the expansion point of the Padé approximation -an approach discussed in the following section.

Shifting the Padé expansion point

We start with shifting the expansion point of √ 1 + α. Taking some α 0 > 0, a direct computation gives

√ 1 + α = √ 1 + α 0 1 + α -α 0 1 + α 0 ,
and denoting α := α -α 0 1 + α 0 , we can write

√ 1 + α = √ 1 + α 0 √ 1 + α,
where α is small for values of α around 8 α 0 . Thus, using Theorem 4.5 we obtain the continued fraction representation

√ 1 + α = √ 1 + α 0        1 + α 2 1 + α 2 2+ α 2 1+ α 2 2+...       
, and analogously to the proof of Theorem 4.8, we combine Lemma 4.6 and Proposition 4.7 to obtain the J-fraction representation

√ 1 + α = √ 1 + α 0 1 + α 2 1 - 1 2 + 4 α -1 2+ 4 α -... .
Let us emphasize here that the equality is valid only for the infinite continued fraction and once we truncate, the correspondence is again as established in Proposition 4.7.

Similarly, we have that

α = α • (1 + α 0 ) + α 0 , (51)
so that we can rewrite the function t∞ (z) as

t∞ (α) = 1 + 2 α + 2 α √ 1 + α = 1 + 2 α(1 + α 0 ) + α 0 + 2 α(1 + α 0 ) + α 0 √ 1 + α 0 • √ 1 + α = 1 + 2 α(1 + α 0 ) + α 0 1 + 1 + α 2 √ 1 + α 0 - 2 α(1 + α 0 ) + α 0 • α 2 • 1 2 + 4 α -1 2+ 4 α -... and set t∞ α 0 (α) := 1+ 2 α(1 + α 0 ) + α 0 1 + 1 + α 2 √ 1 + α 0 - 2 α(1 + α 0 ) + α 0 • α 2 • 1 2 + 4 α -1 2+ 4 α -...
.

If we want to approximate t∞ around α 0 it is intuitive to take some truncation of the continued fraction in t∞ α 0 and based on Theorem 4.8 the expectation is that taking i levels will result in an [i + 1, i + 1]-Padé approximant of t∞ around α 0 .

Following the development in Section 5, we define the approximation Without relabeling the function 9 we can write

ti α 0 (α) := 1 + 2 α(1 + α 0 ) + α 0 1 + 1 + α 2 √ 1 + α 0 - 2 α(1 + α 0 ) + α 0 • α 2 • 1 2 + 4 α - . . .
ti z 0 (z) = 1 + z 2 1 + 1 + 2 z 0 z -1 4 + z 0 1 + 4 z 0 - 1 2 + 4+z 0 z 0 z -1 - 1 . . . 2 + 4+z 0 z 0 z -1 - 1 2+ 4+z 0 z 0 z -1 i "levels" . ( 52
)
Hence we can write the error function err P (z, i) (P for Padé) as

err P (z, i) := | t∞ (z) -ti z 0 (z)|.
The expectation is that the error function err P (z, i) should have one root at z 0 = 4/α 0 , which should get more pronounced as i increases. Indeed, the numerical results shown in Figure 7 support this fully. However, as Remark below emphasizes, this doesn't constitute the result yet.

Remark 6

The results above do not prove that the function ti z 0 is a Padé approximant of either t∞ z 0 or t∞ about the point z 0 (and then analogously in the α variable) because the construction took into account only the approximant of √ 1 + α in the α domain but neglected the rest of the function t∞ α 0 . However, the numerical results suggest that in spite of that the function ti z 0 has good approximation qualities. We end this section by comparing the error function err P (z, i) to the other ones, i.e., to err D (z, i) and err R (z, i). First, we present the equi-oscillation result in Figure 8, analogously to Figure 5.

We again did not use the natural values of z k but rather chose to span the entire interval with logarithmically equidistant points on the left and equidistant points on the right. We used N = 200, i = 5, η = 2, obtaining the optimal value of z 0 to be z 0 ≈ 0.12 (the optimal value of α 0 amounts to α 0 ≈ 32.57). The Dirichlet error in the infinity norm is equal to 0.1371 the Robin error in the infinity norm is equal to 0.0245 while the error of the "Shifted Padé" (see Remark 6 is equal to 0.0052, an almost 4.7 fold improvement over the optimal Robin parameter p * approximation and an overall 26.2 fold improvement over the original Dirichlet approximation.

Running the optimization while varying i, i.e., the number of grid columns from after a we obtain Table 2 (where we keep N = 200 and η = 2).

This shows that the improvement over both of the previous options increases with increasing number of layers and does so quite rapidly.

The corresponding results over a larger range of i are shown graphically in Figure 9.

Again we minimized the maximum norm err P (z, i) ∞ over the interval in (44). We used again N = 200, η = 2 and varied i as powers of 2 from 2 0 = 1 to 2 8 = 256 on the left and then up to 2 15 on the right. We can see that for i ≤ 64 there seems to be a trajectory for the optimal choice of z 0 . The evolution afterwards is caused by the finite precision. For i ≥ 80 the optimization routine (we used the Nelder-Mead algorithm) reaches mach on the entire interval [ηh 2 , 4 + ηh 2 ] seemingly for any z 0 beneath a certain threshold. As we increase the value of i this threshold increases as well up to a point where the initial guess z 0 = 1 already for the so-called absolute term in (52). Already here we can notice that the structure of the absolute term seems quite similar to the exact formula in (33), especially compared with the formulas (34) and (47) corresponding to the approximations given by the Dirichlet and Robin boundary conditions at b. Recalling the elimination process ( 19) -( 22), we recover the last N b -N a -1 block rows in A b (i.e., the block rows governing the unknowns in the interior of Ω b ) from (53) and the (N a -1)-st block row (i.e., the block row governing the interface of Ω a and Ω b ) from (54). In particular, denoting the modified matrix by Āb we have

i optimal z 0 err D ∞ err P ∞ err R ∞ err P ∞ 1
Āb = 1 h 2             D 1 -I N -I N . . . . . . . . . DN a -I N -I N DN a +1 . . . . . . . . . -I N -I N DN b             (55)
where DN b = DN b -1 = . . . = DN a +1 = DN a . Reversing the process of diagonalization of the Schur complement from Section 4.1 and realizing that µ 0 , η and h are only real constants, we get for any i

= N b -1, . . . , N a + 1 Di = Q T    2 + z 1 4+z 0 µ 0 -µ 1 . . . 2 + z N 4+z 0 µ 0 -µ N    Q = Q T    2 + (4 + z 0 )Z    µ 0 -µ 1 . . . µ 0 -µ N    -1    Q = 2I + (4 + ηh 2 + µ 0 )(D -2I)(µ 0 I -D yy ) -1 , (56)
with Z defined as the diagonal matrix Z = diag(z 1 , . . . , z N ) (see ( 16)), Q being the discrete Fourier sine basis (see [START_REF] Magoulès | Algebraic approximation of Dirichlet-to-Neumann maps for the equations of linear elasticity[END_REF]), D being the diagonal block of the original problem (see ( 5)) and D yy being the three-point finite difference stencil discretization of the second derivative in y (see [START_REF] Lorentzen | Continued Fractions with Applications[END_REF] and above).

Focusing on the block DN a we obtain

DN a = Q T      1 + z 1 2 + z 1 2 1 + 4 z 0 + µ 0 -µ 1 z 0 +4 1 + 4 z 0 . . . 1 + z N 2 + z N 2 1 + 4 z 0 + µ 0 -µ N z 0 +4 1 + 4 z 0      Q = Q T    1 2 (2I + Z) + 1 + 4 z 0 2 Z + 1 + 4 z 0 z 0 + 4    µ 0 -µ 1 . . . µ 0 -µ N       Q = D + 1 + 4 z 0 2 Z + 1 + 4 z 0 z 0 + 4 (µ 0 I -D yy ),
with the notation as in (56). We finish this section with the following remark.

Remark 7 Notice that the formula (56) contains an explicit inverse. This raises the question whether it is more reasonable to perform the inverse operation (or, in practice, the solve operation) in the original Schur complement (and thereby obtain the exact result) rather than on this approximation. This and other practical challenges as well as an overall deeper understanding of Āb and its continuous counterpart are clearly of interest and will be discussed in future work.

Conclusion and future work

We proved for a model problem that truncation of the unbounded computational domain by a Dirichlet boundary conditions at a certain distance away from the domain of interest is a spectral Padé approximation about infinity of the transparent boundary condition at the boundary of the domain of interest, and that the degree of the Padé approximation increases with the distance. We then replaced the Dirichlet truncation condition by a Robin truncation condition and showed that this greatly improves the behavior around a different point in the spectrum. We showed how to optimize the Robin parameter leading to an equioscillation property, but this is not a Padé approximation of the transparent boundary condition any more.

Aiming to obtain the Padé approximation about a different point we have proposed a different approximant in the eigenspace (leading to a new PML method for this problem), which poses a significant improvement over the Robin truncation. However, the theoretical proof of the approximation property is an open problem, which needs to be addressed properly on its own. We showed numerical results on the optimal choice of the parameter z 0 , i.e., the shifted expansion point.

There are many further roads of exploration opened up by our approach: first, one could try to obtain an asymptotic formula for p * (i) as h → 0, which would require to obtain the first expansion terms in our closed form formula for the error function err R (z, i, p) for the Robin condition. Similarly one could try to obtain an asymptotic formula for the best parameter in volume. Both would need expansions of the finite continued fraction (47), which could be quite technical. Analogously, the same direction is worth exploring also for the approximant with shifted expansion point z 0 . Moreover, the question about the nature of the approximation produced by the Robin truncation is still open as well as the question whether the second method matches the Padé approximation about the point z 0 . Last but not least, putting the above in the context of the work on the Zolotarev approximation in [START_REF] Ingerman | Optimal finite difference grids and rational approximations of the square root : I. Elliptic problems[END_REF][START_REF] Asvadurov | On optimal finitedifference approximation of PML[END_REF] seems also beneficial. We intend to address these in a future work.

Recognizing that our results were developed for a very particular problem, namely for the η -∆ equation on an unbounded strip in R 2 , there are some straightforward generalizations as none of our computations required the particular choice of D in [START_REF] Bernkopf | A history of infinite matrices[END_REF]. As long as D is symmetric and positive-definite, all of the computations still work and the only change is in the interval of interest for the minimization of the Robin parameter p and the shifted expansion point z 0 in Section 5 and Section 6.

This even holds if D is only symmetric, non-singular and with eigenvalues outside the interval (-∞, -1]. If the spectrum intersected the interval (-∞, -1], the square root becomes a complex number and we would have to move to the complex domain with the continued fractions. The same holds in fact for any diagonalizable non-singular normal matrix D. If D is not normal, then the eigenvectors cannot be chosen to form an orthonormal basis of R N (or C N ). In that case, the formulas would follow (based on the spectrum) one of the above mentioned cases in the same way, but one could not use the results directly. For example, the improvement factor would not be of immediate interest as the condition number of the eigenbasis would play an important role in computing the optimal Robin parameter p. If the matrix is diagonalizable and singular, then the modes corresponding to the zero eigenvalues do not admit the formulation of the function ti k (z) as in (31) but the analysis would work for the rest of the modes, based on the normality and spectrum of the matrix. In the case that the matrix is not diagonalizable, it is not immediately clear how to generalize any of the results based on the available Jordan form.

Finally, as we mentioned in Section 4, the three term recurrence (and thus the continued fraction formulation) has a deep, non-trivial connection with many other areas of mathematics, such as orthogonal polynomials, Gauss quadrature and the conjugate gradient method. Investigating this further would certainly be a worthwhile effort.

Figure 1 :

 1 Figure 1: The unbounded strip domain in R 2 with Ω = (0, +∞) × (0, 1).

Theorem 4 . 1

 41 The Schur complement T b

  Theorem 4.2 ([2, Theorem 1.5.3, 1.5.4, 1.5.1]) Let f (z) be a real function of a real variable. Then the following holds provided the Padé approximants exist :

Lemma 4 . 4 ([ 2 ,+ a 1 b 1 + a 2 b 2 + a 3 b 3 +

 4423 Section 4.1, pp. 105-106]) Let {a j } j , {b j } j be two real sequences of the numerators and denominators of a continued fraction as in Definition 4.3. Let {e j } j be a sequence of real numbers different from zero. Then we have b 0 • • • = b 0 + e 1 a 1 e 1 b 1 + e 1 e 2 a 2 e 2 b 2 + e 2 e 3 a 3 e 3 b 3 + . . . ,

  2, 4.4, 4.5], [16, Chapter V] or [22, Part II].

Theorem 4 .

 4 8 allows us to quantify the error 5 err D (z, i) := t∞ (z) -tb N a +i (z),

Figure 2 :

 2 Figure 2: Plots of the function err D (z, i) at the points z k with the parameters set to N = 20 and η = 2. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

  we can eliminate the unknowns u N b +1 -the ghost points -to get (D N b + 2hpI)u N b -2u N b -1 = 0, and thus the modified matrix A b becomes

Figure 3 :Figure 4 :

 34 Figure 3: Plots of the function err R (z, i) at the points z k evaluated for different values of i, for p = 20, N = 20 and η = 2. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

Figure 5 :

 5 Figure 5: Left: minimization over p of the infinity norm of the Robin condition error, clearly showing the equioscillation. Right: optimized error compared with the corresponding Dirichlet condition error. We set N = 200, i = 5 and η = 2. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

Figure 6 :

 6 Figure 6: Dependence of the optimized Robin parameter p * (i) on the number of layers i added after a compared with the predicted behavior. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

,

  and continue by focusing on the formulation of ti α 0 as a function of z rather than α. Recalling the definition of α in (51) we have

Figure 7 :

 7 Figure 7: Plots of the function err P (z, i) at points equally spaced in the interval [0, 4] evaluated for different values of i, for α 0 = 4 (and thus z 0 = 1), N = 200 and η = 2. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

Figure 8 :

 8 Figure 8: Left: minimization over z 0 (or equivalently over α 0 ) of the infinity norm of the error function err P , clearly showing the equi-oscillation. Right: optimized error compared with the error functions err R for the optimal p * and the error function err D . We set N = 200, i = 5 and η = 2. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

Figure 9 :

 9 Figure 9: Dependence of the optimal choice of z 0 (and consequently α 0 = 4/z 0 ) on the number of layers i added after a.

Table 1 :

 1 Evolution of the optimized Robin parameter p * (i) depending on the number of layers i and the improvement ratio from the Dirichlet condition error to the Robin condition error in the infinity norm. The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

		.4013	2.569
	2	13.7783	3.924
	4	8.2295	5.167
	8	5.6016	6.598
	16 4.3271	8.940

Table 2 :

 2 Evolution of the optimized expansion point z 0 depending on the number of layers i and the improvement ratio from the Dirichlet and Robin boundary condition error to the error of the approximation ti z 0 . The value i corresponds to the number of grid columns in Ω b \Ω a , see Remark 3.

		0.4356	3.691	1.441
	2	0.2101	10.091	2.572
	4	0.1409	18.446	3.569
	8	0.0932	86.163	13.058
	16	0.0680	3595.822 402.186

Section de Mathématiques, Université de Genève,

2-4 rue du Lièvre, CP 64, CH-1211 Genève. Email: martin.gander@unige.ch, michal.outrata@unige.ch 2 Laboratoire de Mécanique et d'Acoustique, Université d'Aix-Marseille, Technopôle Château-Gombert, 4, impasse Nicola Tesla, CS 40006, 13453 MARSEILLE Cedex 13, Email: jakabcin@lma.cnrs-mrs.fr

It would be possible to do the analysis that follows for a much more general PDE, but for simplicity and clarity, we use our model problem (1) throughout; more comments can be found in the concluding Section 7.

The subscript D stands for the "Dirichlet" boundary condition at the end point x = b.

A Robin boundary condition is a simple approximation to the transparent boundary condition and works in general substantially better than a Dirichlet condition; for subdomain truncation in domain decomposition see[START_REF] Gander | Optimized Schwarz methods[END_REF].

Meaning we do not consider α, α 0 such that |α -α 0 | > 1 in order to ensure the argument of the square root being positive.

However we do change the expansion point α 0 to z 0 in the subscript.

Acknowledgement

We would like to thank prof. Zdeněk Strakoš for his very useful comments and references to the continued fraction literature.

suffices.

If the approximant ti z 0 would be indeed the [i + 1, i + 1]-Padé approximant around z 0 = 4/α 0 (see Remark 6), then this behavior was to be expected from some value of i onward. The fact that the optimal "expansion point" z 0 tends towards the left endpoint was also expected as we have seen that the error is largest around that point.

We conclude this section by linking the above proposed approximation back to the physical problem and its solution methods by introducing a new PML technique that stems from the above approximation in the following section.

A new PML technique

Although the function ti z 0 (z) in (52) may not be the Shifted Padé approximant the numerical results suggest that it still leads to a very good approximation and by the construction we will be able to propose good and very accurate PML compared to the previously proposed ones in this paper.

Recalling the definition of A b in (4), the definition of the Schur complement in Definition 3.1 and the process of its transformation into a continued fraction, we need to reverse this process (starting with the eigenvalue reccurences and working up to the block matrix A b ). This requires returning to z k from the artificial z and linking the formulas containing the variables z k to the blocks in A b and thus obtain the PML method.

First, recalling

, as in [START_REF] Lorentzen | Continued Fractions with Applications[END_REF] we denote now the eigenvalues of the matrix D yy by µ k , obtaining

having an artificial z 0 as we have above can now be translated as having an artificial µ 0 and the denominator of the cyclic part of the continued fraction in ti z 0 (z) in (52) can be now rewritten back in terms of the variables η, h, µ k that are of interest. Returning to z k instead of z, a direct computation gives

for the reoccurring term of the continued fraction in (52) and