Combining Rational Approximation with Asymptotic Wiener–Hopf Factorization Algorithm for Matrix Functions: Implementation and Testing

This paper discusses how an asymptotic Wiener–Hopf factorization can be implemented for a wide class of functions. Asymptotic Wiener–Hopf factorization was discussed [19] and the convergence for matrices sufficiently close to the identity matrix is shown. We demonstrate how the algorithm can be successfully implemented with the help of the rational approximations. The idea is to simplify the matrix first by rationally approximating and then perform the approximate factorisation. There is no compromise in accuracy since the factorisation is approximate anyway and rational approximations are very precise usually. There is also a mapping of real line discussed and implemented to make the rational approximation more optimal. The code is tested against some easy examples which are calculated by hand. The use of this code is illustrated with some more complicated matrix functions motivated by applications. The method has been implemented for 2 × 2 and 4 × 4 matrices but can easy be adapted for any size matrix. The code will be made available with the publication of this paper. We note that to date there are very few implemented Wiener–Hopf factorisation available due to instabilities, so this paper will make an important contribution to this area. Citation: Rougerie T, Kisil A (2017) Combining Rational Approximation with Asymptotic Wiener–Hopf Factorization Algorithm for Matrix Functions: Implementation and Testing. J Appl Computat Math 6: 374. doi: 10.4172/2168-9679.1000374

The most challenging step of the method is Wiener-Hopf factorizations. There are fundamentally two different cases to consider: the scalar and the matrix Wiener-Hopf factorizations. For the scalar case the solution of Wiener-Hopf factorizations can be expressed in terms of a Cauchy type integral. Matrix Wiener-Hopf factorizations is a natural extension of the traditional scalar Wiener-Hopf factorizations. It allows to solve more advanced and realistic problems and is the field of current intensive research. Approximating matrix Wiener-Hopf equations is complex due to instabilities. This paper concerns reliable approximate matrix Wiener-Hopf factorizations and its implementation. The existence of matrix factorizations under general assumptions has been known for a long time [12,22]. Nevertheless, to date constructive matrix factorizations remains a formidable challenge. By constructive factorizations it is meant that there is an algorithmic way to obtain the factorizations. Due to its complexity, different classes of matrix functions have to be treated separately. Below we outline some of the few classes of matrix functions with known constructive factorizations and give references that an interested reader can follow. Constructive factorizations algorithms exist for classes of Veitch and Peake [25].
• Rational matrix functions. An algorithm for factorization has been known to exist for a long time. Nevertheless it has not been successfully implemented (as far as the authors know) due to instabilities.
• Upper or lower triangular matrix functions (with factorable diagonal elements). They can be reduced to a scalar equation.
• Functionally commuting matrix functions [18]. One of the reasons is that one could apply some of the techniques that have been developed in the scalar setting. For some conditions under which the commutative factorisation is possible [4].
What complicates the application of algorithms further is that, given a matrix, it is difficult to determine which class it belongs to, if any. Pre-and post-multiplication can transform the matrix to a right form but this is done mainly using ad hoc techniques. A systematic treatment of some transformations is given [8].
The next best option is approximate constructive factorisation which is what is mainly used in applications [6,17]. There are a number of approaches but they are usually tailored to a particular application and hard to generalize and make algorithmic. In other cases such as rational matrix functions there does exist an algorithm for Wiener-Hopf factorizations but it is not implemented due to instabilities. The authors are not aware of any implemented, freely available algorithms for Wiener-Hopf factorisation of matrix functions. This paper will be concerned with implementation of one type of approximate constructive factorisation called asymptotic factorisation [19], which is outlined below.
The structure of the paper is as follows. We outline the mathematical ideas in the proposed procedure. The next section comments on the implementation of the algorithm in MATLAB. The use of the code to compute examples starting with easy cases which can be verified by hand then by examples which arise in acoustics. We end with a few conclusive remarks.

Procedure for Approximate Wiener-Hopf Factorizations
There are a few steps in the implemented approximate Wiener-Hopf factorizations, which are described below.

Asymptotic Wiener-Hopf factorization
Let M be a matrix of the following form: where the m ij (x) are functions analytic around the real axis R. Let I be the identity matrix and ε ∈ [0,1]. We can write the asymptotic algorithm produces a good approximate Wiener-Hopf factorization if εG(x) is small compared to I. This will be discussed in numerical examples.
We are going to find step by step approximate Wiener-Hopf multiplicative factorizations (we will only consider the case of zero partial indices), we aim to find M − (x) M + (x) where: on the real line and in a superscript+(or −) indicates that the function is analytic in the upper (or lower) half-plane. We will look for an approximate factorization of the form [19]: We define the error ∆ j at a step j as the term neglected in the approximation: We will describe here the first four steps that will be used in the code. It would be trivial to include more steps in the code if there was a need for it.
For the first step, we are looking for 0 0 G and G − + respectively analytic on the lower half-plane and the upper half-plane, so that Hence, at the first order of ε, we have the Wiener-Hopf additive splitting: The error is For the second order of ε: where the term in ε 2 should be equal to zero.
Therefore, we can deduce 1 1 G G − + and thus ∆ 2 : ( ) Recursively, for the third and fourth orders of ε, we have There are two issues which need to be addressed for the approximate procedure to be successful. The first is the question of convergence i.e., if the error decreases at every step. This has been discussed [19] and the general criterion formulated, but this criteria is not easy to verify. In practice, when using this procedure this is easily verified by displaying the ∆ i . Secondly, the above procedure is designed for matrices which have zero partial indices. It has been modified for partial indices differing by at most one [20]. So in order to apply it we first want to determine the partial indices. This is in general a difficult issue. Fortunately, (I+εG(α)) and εG(α) is small enough the partial indices are zero.

Rational matrix functions
One could find 0 0 G and G − + using theorem B, but the repeated evaluation of Cauchy type integral would cause computational problems. Instead we first implement the rational approximation and then split the matrices at each step. The splitting/factorization of the rational matrix function is then straightforward as is described below. When G is in rational form, we split the poles of the lower and the upper half-planes of each rational function in order to find 0 0 G and G − + . This is implemented in split add.m.
Let M be a rational matrix function, represented by 11 where p ij and q ij are polynomial functions, q ij not identically zero on the domain, for any i, j ∈ {1,2}.
Using a partial fraction decomposition, we can write (in the case of no repeated roots, otherwise can be easily adapted), for any i,j ∈ {1,2}: where m ij d are the m poles of this fraction, ( ) m ij r x the corresponding residues, and k(x) is the direct term.
We assume that none of the poles are on the real axis. If at least one of them is, we can rotate the x-axis using x↦exp(iαx), with α ∈ (0,π/4] and find the factorisation with respect to the new x-axis. We can split the sum between the poles of the lower and the upper half-planes: where and We finally have All the G i , i>0 will be rational matrix functions, as sums and products of rational matrix functions G j , j<i. Thus the additive splitting of G i will be done as explained above.

Mapping of the real axis for even functions
The rational approximation of even function over R can be improved by employing transformations. We can map the domain by a succession of two conformal mapping in order to obtain better rational approximations. These choices of maps are based on [16]. This is implemented in map.m which is used in rat approx. m. If a non-even function was to be studied in our algorithm, the rational approximation would be performed without the mapping in rat approx. noteven.m.
In this section, for any even function G of the variable y, analytic on ℝ, we will write G(z) and G(x) respectively for its mapped version on the unit circle and on the interval [−1,1].
The first step is the mapping of the real line into the unit circle, using the conformal M obvious map: ( ) 1 1 1 : , This map does the transformation 1↦ i, −1 ↦ −i, ∞↦ 1.
We then apply a Joukowski map to project the unit circle to the interval [−1,1]: Here, only half of the circle is projected on the interval [−1,1]. The function being even, the other half is the same.
The two maps are composed together: , Hence, G(y) can be mapped into G(x) with and unmapped using G(y)=G(JM(y)).
This transformation from G to G by unmapping [−1,1] to the realaxis keeps the rational property of the function (performed in map roots).

MATLAB Implementation of the Algorithm File main.m.
We implemented an algorithm in MATLAB in order to compute the different steps of the approximation. The file main.m has six phases: • Importing data X.m, and computing a rational approximation if the data is composed of non-rational functions, Step 1: Computing Step 2: Computing ( ) Step 3: Computing Step 4: Computing ( ) A schematics of the code is given in Figure 1. For more details see read me.txt file in the code.

The Chebfun library
We used the open-source Chebfun MATLAB library to compute a rational approximation of non-rational functions. The approximation of the g ij (x), i,j ∈ {1,2} as rational functions is done by the chebpade function of Chebfun which consists in a Chebyshev-Pade approximation: chebpade (g ij (x),d 1 ,d 2 ) returns two polynomials p ij (x) and q ij (x) of degrees respectively d 1 and d 2 such that . Choosing d 1= d 2 ≤ 8 usually gives some good results. When unmapped, the half of the even function becomes entire so that the resulting polynomial functions will have a degree of ≤ 16.
In practice, we first create a chebfun function on a [−1,1] interval for each G ij (x), i,j ∈ {1,2} where G is the mapped version of G. In our algorithm, we chose a [−1,0.9] interval rather than [−1,1], in order to map the function without including all the values of the function to +∞.

Rational matrix function with complex poles
The first example (data A.m in the code) is a rational function given by ( ) We derive by hand the approximate factorisation where ε=1 in order to compare it to the code output. The partial fraction decomposition gives: At this point, the part neglected equals ( ) ( ) Then, at the second step, Another partial fraction decomposition on G 1 leads to: According to eqn. (7), the error is the third order of ε gives: and ( ) ( ) Here, the term neglected was, as in (8) For the fourth step: In the end, the error equals: The results given by the algorithm code match these derivations with an error of 10 −16 due to MATLAB approximations.
We expect the error ε↦∆ i (ε) for a given step i to be decreasing when εA is small compared to I ([?art1] Mishuris & Rogosin 2014). We check this property for ε=1, ε=0.5, ε=0.05 and ε=0.005. We find that as long as εA is less than 0.8−0.9 in all elements, convergence occurs. Otherwise ∆ i (ε) does not decrease and can even increase. Figure 2 shows that the method of our algorithm works much better for εA small compared to I: the error is decreasing much faster, hence less steps are required to obtain a good approximation.

Rational matrix functions witTh poles on tThe real axis
The second example (data B.m in the code) is a rational function defined by: The poles are y=2 and y=−2 so the x-axis can be rotated, using y → e ik y (rot.m), where k ∈ [0,2π] is a parameter of the algorithm. This has the effect of making sure that y=2 is included in the poles in the lower half plane and y=−2 with the poles in the upper half plane. Hence and y=2e −ik and y=−2e −ik are the new poles. The complex plane can now be split in two halves, depending of the sign of the imaginary part of each poles. Again we calculate the decomposition by hand and compare them to the code output. The agreement is again of order of magnitude 10 -16 .
The study of I+εB, ε=1 gives some poor results because B is larger in magnitude than I. Even at the fourth step, |∆ 4 | has a magnitude of more than 0.7. This huge error gives a bad approximation of I+εB (Figure 3).
For ε=0.005, the algorithm gives a very good approximation because εB is now small compared to I. At the fourth step, |∆ 4 | has a very small magnitude of 10 −7 . The approximation of I+εB is now very close to the actual I+εB (Figure 4).

A more general rational matrix functions
We tested the algorithm on a third example (data C.m in the code) without any symmetry: In this case again, the error i ↦ ∆ i is not always decreasing ( Figure  5), for example for ε=1, but once the perturbation matrix is small enough the expected decrease is achieved.

Non rational matrix functions
This example (data G.m in the code) is motivated by acoustics, the Wiener-Hopf problem arising when the plane wave is incident at an angle π/2 onto a periodic grating of length 2a which is spaced 2b apart [9]. The matrix to factorize is ( ) . We will find the rational approximation ( )  Considering the rational approximation directly gives poor results. Instead, we will calculate a rational approximation after a conformal mapping.
For parameters a=2.5, b=1.5 and rotate the line by exp(−0.16i), the functions can be approximated and the error is less than 10 −6 . The poles lie where it is desired on the rotated imaginary axis and rotated real axis, Figure 6. The approximations are shown for the first step in Figure  7 and the error is large in the non-diagonal elements. The forth step Figure 8, shows that the error is very small for all entries of the matrix.

Conclusions
This paper presents an implementation of asymptotic factorization for matrix functions of the form I+G. The implementation relies on rational approximation.
It is shown that this procedure converges for a number of examples as long as the elements of G are smaller than 0.8−1. Typically values around 0.7 or smaller give good results in all cases. In order to make this procedure applicable to a wider class it is possible to consider suitable pre-and post-multiplications [8]. For example, S can be far from I but once we transform F − SK + , it can be close to I matrix. The code is currently written for 2 × 2 and 4 × 4 matrices but can easily be extended to deal with other dimensions.

Authors Contributions
This paper resulted from a summer project of Tiphaine Rougerie under the supervision of Anastasia Kisil. The code was mostly written by Tiphaine based on the earlier code of Anastasia.

Acknowledgements
Anastasia Kisil benefited from useful discussions with Professors Abrahams,