The Wiener–Hopf technique, its generalizations and applications: constructive and approximate methods

This paper reviews the modern state of the Wiener–Hopf factorization method and its generalizations. The main constructive results for matrix Wiener–Hopf problems are presented, approximate methods are outlined and the main areas of applications are mentioned. The aim of the paper is to offer an overview of the development of this method, and demonstrate the importance of bringing together pure and applied analysis to effectively employ the Wiener–Hopf technique.


Introduction
The Wiener-Hopf method has been motivated by interdisciplinary interests ever since its inception. It resulted from a collaboration between Norbert Wiener, who worked on stochastic processes, and Eberhard Hopf, who worked on partial differential equations (PDEs). The method was first described in their joint article [1] where they study a convolution-type integral equation for f (x) k(x − t)f (t)dt = g(x), x > 0, (1.1) with k(x) and g(x) given. The study of this equation was principally motivated by an interest of one of the authors, in a differential equation governing the radiation equilibrium of stars [2]. The Wiener-Hopf equation arises from extending the integral equation ( for some function h(x). Note that, although h(x) is unknown, it is not independent as it is uniquely defined by the left hand side of (1. 2) once f (x) is determined. Applying the Fourier transform to (1.2) results in an equation of type (2.1) (the interaction of a Fourier transform with its convolution has to be used, and analytic properties of half-range Fourier transforms employed 1 ). This laid the foundation for the study of scalar Wiener-Hopf equations. Independently, the theory of the more general Riemann-Hilbert boundary value problem has been developed. The Riemann-Hilbert boundary value problem for a particular case was first formulated by Riemann as a part of the problem of construction of complex differential equation of the Fuchs type with a given monodromy. In the latter form, the problem was presented by Hilbert as his 21st mathematical problem for the twentieth century [4]. Plemelj [5] proposed a method of solution to the Riemann-Hilbert boundary value problem based on the reduction to a case of the homogeneous Riemann boundary value problem proving its solvability by means of the Fredholm alternative. 2 An intensive study of the Riemann-Hilbert boundary value problem began in the 1930s. Properties of the factorization problem and its role in the study of the Wiener-Hopf integral equation were outlined in the seminal paper [6]. Further development of the factorization theory was described in the book by Litvinchuk & Spitkovskii [7]. Several known facts on factorization were presented in a recent survey [8] (see also [9]), where one can find, in particular, an extended list of references on the subject.
Simultaneously, the general theory of integral equations was developed intensively at the beginning of twentieth century. Of note are the contributions by Hilbert [10], Fredholm [11], Volterra [12] and Carleman [13]. The latter two authors also considered integral equations with kernels depending on the difference of arguments, but containing integrals with variable limits or along a finite interval, respectively. Nevertheless, the Wiener-Hopf integral equation (1.1) became an important independent subject of research due to its wide range of applicability and many connections. The Wiener-Hopf method has enabled researchers to solve analytically numerous previously intractable integral and partial differential equations, as well as many boundary value problems. To date, it still remains the standard method for solving a wide class of canonical physical problems. The method is an important cornerstone (among other methods) in areas of pure mathematics, for example in development of the abstract theory of singular linear operators, pseudo-differential operators, etc. Additionally, the exact solutions form the basis of approximate techniques applicable to more complicated problems. Moreover, even in case when the constructive (analytical) construction is available, it is not always possible to realize numerically with a confidence [14]. 1 In fact, in the original paper [1] the Laplace transform was applied and the corresponding functional equation is valid in a vertical strip. Interestingly the original paper contained a mistake; it claimed that γ n (z) = (z 2 − 1) n/2 → z n as z → ±∞ where n was the number of zeros of the kernel in the strip of analyticity (this is only true for n even). This was pointed out in ( [3], ch. 16), and resulted from an incorrect choice of branch cut. However, it did not influence the actual results they obtained in the second part of the paper, since they assumed that the kernel is even and hence n is even. For the sake of completeness, we are going to outline the usual procedure used to solve the Wiener-Hopf problem in the case of canonical scalar factorization. The key step in the Wiener-Hopf procedure is to find K + and K − (the Wiener-Hopf factors) such that K(z) = K + (z)K − (z), z ∈ H. (2.2) with K ± and their inverses analytic in H ± . This product is referred to as a factorization of the kernel and is the most important part of the procedure; it is straightforward to accomplish for scalar kernels.
Multiplying through by K −1 − we obtain Now performing an additive splitting we have (2.4) where C ± (z) are analytic in their indicated half-planes H ± , respectively. Finally, rearranging the equation yields The right-hand side of (2.5) is a function analytic in H + and the left-hand side of equation (2.5) is analytic in H − . Hence, together they offer analytic continuation from H into the whole of the z-plane, and so each side is equal to an entire function J(z), say. We also require at most algebraic growth of each side of equation (2.5) as |z| → ∞ in the respective half-planes of analyticity, i.e.
for some constant n. Applying the extended form of Liouville's theorem implies that J(z) is a polynomial of degree less than or equal to n. Hence Φ + (z) and Ψ − (z) are determined up to n + 1 unknown constants (typically n ≤ 0 or 1, and the constants can be fixed using extra information about the behaviour of the solution). Then the application of the inverse Fourier transform to Φ + (z) or Ψ − (z) would yield the solution to (1.1). If Φ + (z), Ψ − (z) and C(z) in (2.1) are vectors of functions with a matrix of functions K(z), then (2.1) is called matrix/vectorial Wiener-Hopf equations. All the above steps in the solution are analogous. The difficulty lies in performing the multiplicative factorization (2.2), which becomes the main challenge to overcome and is addressed in § §3 and 4.

(b) Riemann-Hilbert boundary value problem
The Wiener-Hopf equation is closely related to the Riemann-Hilbert boundary value problem 3 [18,19]. The latter problem is to determine two functions or vectors of functions Φ ± (z), which are analytic in the complementary domains D ± (i.e. D + ∪ L ∪ D − = C, with a given simple closed curve L = ∂D + = ∂D − ) and satisfy on L the following linear condition: with the function/matrix G(t) and the function/vector g(t) given on L.
In the scalar case, the complete solution of the Riemann problem (2.7) was given by Gakhov [18]. He has shown that solvability of the problem depends on the Cauchy index (winding number) . (2.8) In the case of Hölder-continuous functions G and g, a simple smooth contour L and in the case of non-negative index κ ≥ 0 the solution of (2.7) has the form 4 where P κ (z) is a polynomial of order κ, and In the case of negative index κ, the unique solution has the form Φ ± (z) = X ± (z)Ψ ± (z) provided that a finite number of solvability conditions are satisfied (otherwise there is no solution). The scalar and matrix function factorization, related to the Riemann-Hilbert boundary value problem, was introduced by Birkhoff [20] and also in a more straightforward form by Gakhov [21]. Factorization means obtaining the following representation of an n × n matrix function G, given on a bounded curve L (curve such that 0 ∈ D + , ∞ ∈ D − ), as where G ± (t) are boundary values of bounded analytic zero-free matrices G ± (z), z ∈ D ± , and where κ j ∈ Z, j = 1, 2, . . . , n, are integer numbers called partial indices. Representation (2.13) is called a left-sided factorization. Interchanging G + and G − we arrive at the right-sided factorization.
Respectively, the partial indices are called left (right) partial indices. In general, they are not the same; however, their sum is always equal to the index of det G (this can sometimes fail for functions with severe discontinuities worse than the piece-wise continuous case [9], theorem 3.21). In the case L = R the variable t in (2.14) is replaced by the ratio (t − i)/(t + i) (similar changes are need in (2.10) and (2.11)). In this case one needs to assume the continuity of G at infinity. Unfortunately, in the general case, there exists no method of determining the partial indices, while such information is important to reconstruct the factors, G ± , and is instrumental in proving the stability of any respective algorithm. Those issues are discussed further in §4. Factorization, i.e. determination of the partial indices κ j ∈ Z, j = 1, 2, . . . , n, and factors G + , G − is an independent problem playing a crucial role in many theoretical and applied problems.
We would like to make here a few comments on formal differences between the Wiener-Hopf and Riemann-Hilbert problem. Note that (2.13) is different to (2.2) since equation (2.2) does not contain Λ(t). It is because it is common in the Wiener-Hopf literature [17] to implicitly assume that κ = 0 in formulae (2.10) and (2.11), which is often the case in applications. On the other hand, (2.13) can be always presented in the form (2.2) by further factorizing the matrix Λ(t) = Λ + (t)Λ − (t), and then redistributing the respective diagonal terms. However, this destroys some of the assumptions on the behaviour of the factors G ± : the behaviour at infinity, satisfaction of the invertibility property (i.e. it introduces zeros into the half-planes of analyticity), etc. Thus, this form of the factorization is tightly linked to the conditions imposed on the factors (which space of analytic function is chosen) [7]. Note also that the factors K ± themselves, ( general have non-zero partial indices. Another important difference between the approaches is that (2.2) is valid on a strip (and thus, knowledge of the possible poles and zeros inside the strip provides useful information for solving the problem under consideration), while the factorization (2.13) is defined on a line (that can be analytically extended if possible). These differences are also discussed in [22]. Finally, we would like to stress that there are authors associating the Wiener-Hopf problem with the Riemann-Hilbert one and vice versa that only underlines the fact that the division is rather artificial and depends on the research communities. To avoid any unintentional misleading, we discuss particular factorization techniques below as they have appeared in the original presentations, with minor comments where necessary.

(c) Canonical matrix
The factorization problem (2.13) is almost equivalent to the homogeneous (i.e. with g(t) ≡ 0) Riemann-Hilbert boundary value problem (2.7): (2.15) In [23], the conditions for the existence of a solution to (2.15) are described. If the solution to (2.15) exists one can find 5 a so-called fundamental system of solutions Φ 1 (z), . . . , Φ n (z) and the corresponding fundamental matrix Φ(z) = (Φ 1 (z), . . . , Φ n (z)) with vectors Φ j (z) being their columns. Fundamentality of the system/matrix means that the determinant of the fundamental matrix cannot be identically zero. The next step (which is more constructive [23, p. 520]) is to transform the fundamental system of solutions to a normal form. By definition the normal system of solutions Ψ 1 (z), . . . , Ψ n (z) to (2.15) is the fundamental system such that the determinant of the fundamental matrix does not vanish anywhere on C (including the curve L). The matrix Ψ (z) = (Ψ 1 (z) . . . Ψ n (z)) of the normal system is called the normal matrix.
By elementary transformations, the normal system of solutions can be transformed to the so called canonical system. The canonical system is usually denoted X 1 (z), . . . , X n (z) and its matrix is X(z) = (X 1 (z) . . . X n (z)). The properties of the canonical system are (i) The canonical system is the normal system of solutions to (2.15), i.e. the determinant of the canonical matrix det X(z) is nowhere vanishing in C. (ii) The sum of the orders κ j at infinity of the columns of the matrix X(z) (i.e. solutions X j (z)) is equal to the Cauchy index of the determinant det X(z).
They are called partial indices of the Riemann-Hilbert boundary value problem (2.15). The integer numbers κ j are the same as partial indices of the factorization (2.13) of the matrix function G(t) (the coefficient in (2.15)).
Note that for a canonical matrix it is said that the matrix X − (z) has the normal form at infinity for a factorization (2.13). An example of the use of this method is presented in §3(c).

(d) Other convolution-type integral equations
Related to the integral equations of the convolution type (1.1) on the half-axis are paired equations and the so-called transpose to the dual equations: They can be both solved in a very similar manner to (1.1) by extending the range, taking the Fourier transform of each of the equations, and eliminating one unknown, to give rise to the usual Wiener-Hopf equation where K j (s) are Fourier transforms of the kernels k j (x), j = 1, 2, and F ± , T(s) are different for each equation (details can be found in [18,24]).

(e) Carleman-type boundary value problem
A closely related equation, but much less known is the so-called smooth-transition equation for −∞ ≤ x ≤ ∞. The interesting feature of this equation is that for x large and positive it is close to And for x large and negative it is close to which is similar to (2.17). Finally, for values of x = 0 close to zero, the equation takes the form: This analysis explains the name of the equation; it behaves as the dual equation when |x| 1 and as the transpose to the dual when |x| 1.
In order to analyse this equation, define the following [3] φ Then the original equation (2.20) can then be written as After application of Fourier transform to (2.24) and (2.25) (again, with a capital letter indicating the transform of a lower-case letter) we arrive at the following equations Eliminating F(z) we arrive at This is not a standard form for a Riemann-Hilbert boundary value problem but it can be transformed to one by defining a new function , (2.28) then the jump in the Riemann-Hilbert boundary value problem is across the cut arg ψ = 0. The smooth transition equation (2.20) is a special case of a more general Carleman-type boundary value problem: to find Φ(z) analytical in a domain D and satisfy on part of the boundary of D denoted S the following condition: (2.29) with the functions E(t), m(t) and g(t) given on S ([3], §15.1). Also it is required that m(t) is continuous and maps to the rest of the boundary of D. Additionally we require that t and m(t) transverse the boundary in different directions. Sometimes its solution can be reduced to a Riemann-Hilbert boundary value problem, but in general no closed-form solution is known [25]. Many applied problems are reduced to certain types of functional-difference equation (related to Wiener-Hopf and Riemann-Hilbert problems). We mention here the problems of applied mechanics [26][27][28][29][30] and diffraction theory [31][32][33].

(f) Discrete Wiener-Hopf equation
In this section, we will look at an infinite system of equations of Wiener-Hopf type, a discrete analogue of integral equation (1.1). For all infinite sequence a n considered here we will assume the following holds for some constant M We will also assume that the Fourier series A(θ ) = ∞ n=−∞ a n e inθ is Hölder continuous. The discrete Wiener-Hopf equations for x n is ∞ k=0 a n−k x k = c n , n = 0, 1, 2 . . . , (2.31) where the infinite sequences a n and c n are given. This is called an infinite system with a Toeplitz operator due to a presence of n − k in the index of a. Similar to the integral equation case, we first need to extend this equation for n negative ∞ k=0 a n−k x k = c n for n = 0, 1, 2, . . .
for some d n , which are unknown for n = −1, −2 . . . . Note that the d n are uniquely defined by the left-hand side of (2.32) once x k is determined. Next we apply the discrete analogue of Fourier transform, the Z-transform. The Z-transform of a sequence a n is defined as and as before is denoted by a capital letter. We multiply the nth equation in (2.32) by t n and sum over all equations: The main property is the interaction of the Z-transform with a Töplitz operator which holds on the unit circle, and the superscript +/− indicates the function is analytic inside/outside the unit circle. This Riemann-Hilbert boundary value equation can be solved as before to find X + (t) and then the inverse of the Z-transform applied to recover x n . A slightly more general infinite system which can be solved in an identical fashion is the discrete version of equation (2.17) ∞ k=−∞ a n−k x k = c n , for n = 0, 1, 2, . . .
where a k , b k , c k and d k are known. Another version of the discrete equation is the analogue to the transpose to the dual equation

Another interesting discrete system which has an analytic solution via a Wiener-Hopf technique is
where |r| = 1. With the help of a Z-transform, this can be reduced to a boundary value problem of Carleman type [3].
There are other methods for solving equation of type (2.31) without the use of Wiener-Hopf method. The discrete systems of equations with a Toeplitz operators have effective numerical methods of solutions, some which use Cauchy-like matrices, which have small numerical ranks [34]. Also there are some results about perturbation of Toeplitz operators for example solution of infinite systems like Au + Bu = f where A is Toeplitz operator and B is in some sense small [35].
Equations of type (2.31) naturally appear in applications when the boundary conditions are discrete and periodic in semi-infinite configurations. For example, this happens in the Sommerfeld half plane problem for discrete Helmholtz equation [36]. This is also the case when an semi-infinite discrete array of point scatterers is considered [37][38][39][40]. Additionally, this type of problems is common in crack propagation problems [41].  accounting for the index of the function K(z), where c is some non-zero complex constant. This can be seen by applying analytic continuation to K + /P + and P − /K − and then using the extended Liouville's theorem [42]. Here we deal with factorization of the matrix defined on a simple closed curve L in the complex plane. The corresponding result in the 2 × 2 matrix case is rather more complicated: ). Let the matrix G(t) admit the Wiener-Hopf factorization (2.13). Then: where is also a factorization of G(t) where H is a constant invertible matrix function if κ 1 = κ 2 and otherwise is of the form: This means that there is more freedom to choose the Wiener-Hopf factors when κ 1 = κ 2 . In the case of higher order matrix functions, similar results can be shown, while the matrix structure has a bit more complicated form [43].

(h) General class of PDEs that can be reduced to a Wiener-Hopf
As an illustration of the concept, we consider a class of partial differential equations (PDEs) that can be reduced to a Wiener-Hopf equations or a Riemann-Hilbert boundary value problem ([3], §20.3). This method is suitable for linear PDEs of the form where h pq and g are given and u is to be determined; note that h pq has no x dependence. In the simplest case, the boundary conditions are given on the lines y =const. or half lines y = const. for x > 0 or x < 0. Let these n lines be positioned at constants y 1 , . . . , y n . The boundary conditions have the form and There could also be some conditions for y → ∞ or y → −∞. In order to obtain a Riemann-Hilbert boundary value problem a list of systematic steps is performed in detail in ( [3], §20.3), and which is translated into English in [44]. We also refer an interested reader to the classic monographs, for example, [
The Wiener-Hopf method is extensively used for canonical scattering problems both in acoustics [17] and electromagnetism [49]. In both of these applications, the governing equation is Helmholtz's equation. There are fewer methods available to tackle Helmholtz's equation compared to Laplace's equation. For the Laplace equation, other techniques from complex analysis such as conformal mapping are applicable [70]. Although the Helmholtz equation is selfadjoint addition of boundary conditions, including the Sommerfeld radiation condition, means that the boundary value problem is no longer self-adjoint. This causes problems in demonstrating completeness of series expansions [71]. The Wiener-Hopf method on the other hand very naturally incorporates the radiation condition and complicated edge conditions. Simple models in aeroacoustics can also be reduced to acoustic problems via Lighthill's analogy and the reciprocal theorem [72].
Water waves is an area where the Wiener-Hopf techniques is accepted as one of key analytic tools [73]. The ice cover on water can be modelled as an elastic material for which flexural motions are able to be described by thin-plate or Timoshenko-Mindlin theory [74][75][76]. Eigenfunctionmatching method is also a commonly used method and in some cases this has been shown to be equivalent to the Wiener-Hopf method [77,78]. Recently, the problem of waves generated in a fluid and an ice sheet by a pressure region moving on the free surface of the fluid, along the edge of the semi-infinite ice sheet, has been solved using the Wiener-Hopf technique [79].
One of the new applications of the Wiener-Hopf technique has been nanophotonics. Graphene and other two-dimensional (2D) materials may sustain evanescent, fine-scale electromagnetic waves that are tightly confined to the boundary, called plasmonpolariton [66]. This is also closely related to edge magnetoplasmons, the theory of which was also derived using the Wiener-Hopf method [80] and the more realistic extension of which lead to more complicated Wiener-Hopf systems [81]. Recently, composites and metamaterials have also been investigated using the Wiener-Hopf method [37,67,82,83]. This is by no means an exhaustive list of applications; unfortunately space constrains us in this short survey. We now list a number of approaches to particular constructive exact or approximate factorization methods in the sections below.

A list of constructive procedures
By constructive factorization, we mean that there is an algorithmic way to obtain the factorization (2.13) or (2.2). Recall, for scalar functions, factorization is always possible via integrals of Cauchy type (2.11) and (2.12). Currently, there are no constructive factorization procedures for matrix functions in general, and only very specific classes can be exactly solved. What complicates the development of an algorithmic approach is that for a given matrix function it is difficult to determine if any known constructive procedure apply. Pre-and post-multiplication can sometimes transform the matrix to a known form, but this is undertaken mainly using ad hoc techniques. A systematic treatment of some transformations is given in [84].
The main classes of matrix kernels which have a constructive factorization are reviewed below. They are the rational matrix functions, commutative matrices, particular matrices with exponential factors, as well as several miscellaneous classes. We also refer the reader to several other important classes of piece-wise constant matrix functions not included in the review [43,85,86], the Meister and Speck matrix as discussed in [87] and for applications [88,89]. We also do not consider the case where there are singularities on the curve L [90]; for more information, the reader is refereed to [

(a) Rational matrix function
This is the simplest class of matrix functions for which an exact constructive factorization is known. This is because there are only a finite number of isolated singularities which need to be considered. A procedure to construct the Wiener-Hopf factorization is known in this case and is briefly described immediately below, followed by the standard 'pole removal' technique.

(i) Factorization of rational matrix functions
This subsection describes an algorithm to factorize a rational matrix function ( [9], §1.2). Given a matrix function M it can be expressed as P/q, where P is a polynomial matrix function and q is a scalar polynomial. It is enough to factor them separately. The scalar function q can be factorized easily, see §2(b). We will assume that det P has no roots on the real axis. What prevents a polynomial matrix P being a plus factor is the presence of zeros of det P in the upper halfplane. (Recall that, for factorization, a plus (minus) factor requires the matrix and its inverse to be analytic in the upper (lower) half plane.) A way of eliminating one zero at a time is outlined in the following lemma.
where z = z 0 is one of the aforementioned poles, while c k are computed constructively and uniquely.
We can repeat this procedure, eliminating all the roots of det P in the upper half-plane and hence construct the plus factor. Note that the matrix R is invertible. The minus factor can be found either by pre-multiplying P by the inverse of the plus factor, or by taking the inverse of the product of all of the R matrices. This completes the factorization.

(ii) 'Pole removal' for rational matrix functions
In applications the above rational factorization algorithm is rarely used; instead 'pole removal' or 'singularity matching' is employed [49,93] and is recapped below.
Suppose that A(z) is rational (every entry of the matrix function is a rational function), then we can solve (2.1) without multiplicative factorization. Perform an additive split A(z)Φ + (z) = (A(z)Φ + (z)) − + (A(z)Φ + (z)) + and C(z) = C − (z) + C + (z), where the first term is analytic in H − and second in H + . Since A(z) is rational we can write (A(z)Φ + (z)) − = n i=0 (A i /(z − z i )), where z i are the poles in the upper half-plane, and A i are constants (residues at z i ). Thus, separates terms on the left-and right-hand sides into functions analytic in the upper and lower half planes respectively. To find A i we need to solve a linear system of equations (by setting z = z i ). The pole removal method has an advantage over the rational factorization procedure. It enables one to remove undesired poles without changing the rest of the equations. This is one of the reason that it can not only be applied to rational Wiener-Hopf kernels but also to more general kernels which have one or more rational parts. Pole removal is closely linked to the generalized Liouville's theorem (which allows the unknown function to have poles). It can also be naturally extended to meromorphic functions by allowing an infinite sum, which would lead to an infinite linear system that would need to be truncated in order to be solved. The additional advantage of 'pole removal' is that the multiplicative factorization is reduced to an additive splitting. It is also worth noting that many alternative formulations, such as Fredholm factorization [49,94,95], approximate techniques (see §4) and numerical methods [96], are also based on this principle or replacing a multiplicative factorization by additive splittings.

Remark 3.2.
Note that the instability issue discussed in §4 still applies here. To recognize the challenge it is enough to refer to the fact that, generally speaking, all poles or zeros are computed approximately, while the system of linear equations mentioned above (which have to be solved to compute the constants) can be ill-conditioned. In the case of a stable set of the partial indices the system is not ill-conditioned, otherwise, at least for some of the poles, this would be expected [97].

(b) Analytic and meromorphic matrix functions
In [98] the factorization problem for the meromorphic matrix function was solved. The corresponding algorithm is based on the reduction of the problem to a finite number of systems of linear algebraic equations, whose coefficients' matrices are block Toeplitz matrices: with C j being the power moments of the inverse G −1 (t) to the given matrix with respect to the contour L The partial indices of the left and right factorization are represented in terms of ranks of the above Toeplitz matrices. An analogous method is applied in [99] for the case of analytic matrix functions which, being a special case of the latter, need less computation. A part of this algorithm related to factorization of the scalar polynomial is implemented in the Maple computer system as module PolynomialFactorization (see its description and illustrative examples in [100]).

(c) Triangular matrix functions
Upper or lower triangular matrix functions (with factorizable diagonal elements) can sometimes be reduced to a set of scalar equations. If the equations can be solved in turn and the previous solution used to make the next equation scalar, then the system decouples. However, even in the case of triangular matrix functions there may be complications which mean that the system does not decouple in this way. This happens when the partial indices are non-zero, for which case we describe a general procedure below. It relies on the idea of a canonical matrix, introduced in §2(c) and this was one of the first attempts to realize this construction explicitly (i.e. to determine canonical matrix factorization for (2.15)).
Denote κ j = ind Γ ζ j (t) and let x ± j (z) be canonical functions for the scalar homogeneous Riemann-Hilbert boundary value problems with coefficients ζ j (t), respectively [18], j = 1, 2. Then the piecewise analytic matrix , z ∈ D ± , satisfies the following boundary condition: X + (t) = G(t)X − (t). Let μ ≥ 1 be the order of the function φ − (t) at infinity. The orders of non-zero elements of the matrix X − (z) at infinity can be characterized by the following table: If κ 1 ≤ κ 2 + μ, then the matrix X − (z) has the normal form at infinity, i.e. it is the canonical matrix. Thus, partial indices of the matrix G(t) are equal (κ 1 , κ 2 ). In the case κ 1 > κ 2 + μ, Chebotarev proposed to use an expansion of 1/φ − (z) as a continued fraction where q γ i (z) are polynomials of orders γ i , respectively (γ 0 = μ). Using these polynomials in the scheme of the elementary transformations it was shown in [101] that in a finite number of steps the 'minus'-matrix can be rebuilt to the normal form at infinity and thus the partial indices will be found and factors in (2.13) will be determined. In [102], this approach was applied to obtain an explicit solution of the factorization problem for a case of the Khrapkov-Daniele matrix functions.
In [103], this approach was applied to solve R-linear conjugation problem. In [104] Chebotarev's approach was realized for triangular matrix function of arbitrary order with factorizable diagonal entries. The basis for the inductive steps is the following statement: let L be a simple smooth closed contour 0 ∈ L, and let B(t), for t ∈ L, be a non-singular Hölder continuous square matrix-function of order n having the following form:

Suppose that the non-singular square matrix-function A(t) of the order n − 1 admits the factorization
where Λ(t) = diag {t κ 1 , . . . , t κ n−1 }. Then the matrix-function B(t) possesses a factorization if the following matrix does: .
(3. Note that for a matrix with a specific structure, the general stability criterion for the partial indices can be extended (in other words, within a specific sub-algebra, there are cases when the perturbation preserving the structure remains stable, while there exists a general (small) perturbation of other structure that changes the set of the partial indices. For the 2 × 2 triangular matrices such analysis has been conducted in [97].

(d) Functionally commutative matrix functions
There has been significant progress in the case of the matrices possessing commutative factorizations [7]. One of the reasons is that one can apply some of the techniques that have been developed in the scalar setting. We present the corresponding result in the case of Hölder continuous coefficients G and g in (2.7).
In [105], the following question was considered: whether there exist a class of matrix coefficients G(t) such that formula (2.9) still gives the bounded solution of the vector-matrix problem (2.7) (i.e. for n > 1)? Such a question is related to the purely algebraic question of the fulfilment of the identity for certain matrices A(t) and B(t) connected to the coefficient G(t). It was shown in [105] that identity (3.7) holds (and thus Gakhov's formula (2.9) can be used to solve problem (2.7) in the matrix case), whenever the matrix function G(t) is functionally commutative. The latter condition means The question of solvability of the vector-matrix problem (2.7) in the case of a functionally commutative matrix G(t) is solved in [105]. In [106], the following necessary and sufficient condition for a matrix to be functionally commutative is reported: the matrix G(t) is functionally commutative if and only if it can be represented in the form where ϕ 1 (t), ϕ 2 (t), . . . , ϕ m (t) are linear scalar independent functions, and G 1 , G 2 , . . . , G m are constant linear independent mutually commutative matrices. The structure of lower-dimensional (up to n = 4) functionally commutative matrix functions was described in [106] too. In [107] a solution to the Riemann-Hilbert boundary value problem (2.7) and the corresponding factorization problem is illustrated by examples [108].

(e) Commutative factorization
Close to the above case is so called commutative factorization. Let us formulate it in the Wiener-Hopf setting. Let the matrix K(z) be defined on a strip H = {z ∈ C : a < Im z < b} around the real line. The commutative factorization of K(z) are factors K − (z) and K + (z) satisfying (2.2) and the factors K + , K − commute in H. The idea of commutative factorization stems from the work of Heins [109]. For criteria when commutative factorization is possible, the reader is referred to [110]. Several special cases of matrices admitting commutative factorizations are described below.

(i) Khrapkov-Daniele matrices
Matrices of this kind represent one of the simplest classes of matrices with commutative factorization. In 2 × 2 case, these matrices were discovered from certain diffraction problems [111][112][113]. They appear in the study of diffraction by geometries which had a right angle such as wedges or gratings (a grating can be thought of as a perpendicular strip in a waveguide) [114] as well as other configurations [ form K(α) = k 0 (α)I + k 1 (α)J(α), (3.10) where J(α) is a polynomial matrix-function with the property J 2 (α) = 2 (α)I, (3.11) with k 0 , k 1 being arbitrary functions analytic in H and having algebraic growth at infinity and 2 is a polynomial in α, Commutative factorization under the above conditions has the form [116] K ± (α) = r ± (α) exp{θ ± (α)J}, (3.12) where scalar functions r ± , θ ± satisfy the relations and (3.14) An interesting approximate technique, based on the Khrapkov-Daniele matrix form, which allows it to be extended to certain non-commutative matrices, has also been developed [119].

(ii) Jones class
Jones [120] established the commutative factorization of the following n × n class of matrices where a 1 (α), a 2 (α), . . . , a n (α) are functions analytic in the strip H, E(α) is an entire matrix-function, E n = q n (α)I. It is supposed also that det C(α) = 0 and tr E r = 0, r = 1, 2, . . . , n − 1. Under these conditions, it was shown that C(α) possesses the commutative factorization Here ω = e (2πi)/n and the functions δ + s , δ − s are analytic in the corresponding half-planes. In [116], the Wiener-Hopf factorization is constructed for matrices of the Jones class under assumption that E is an entire matrix-function with polynomial elements having distinct eigenvalues λ 1 , λ 2 , . . . , λ n such that the following condition holds where μ j are polynomials and n 1 + n 2 + . . . n p = n.

(iii) Moiseev class
Moiseev's method [121,122] is applied to class of matrices similar to that of Jones type which allow non-algebraic growth at infinity where a m (α) are scalar functions and E(α) is a polynomial matrix. In the simplest case of matrix E having distinct eigenvalues almost everywhere, E can be decomposed as where T is the matrix of the eigenvectors and Λ is a diagonal matrix composed of the eigenvalues. Both T and Λ are algebraic matrices (matrices with entries being algebraic functions), and in [121] the Riemann surface R was introduced on which T and Λ are single valued. Further, the matrix factorization problem reduced to a scalar Riemann-Hilbert boundary value problem on R. This problem can be solved in terms of Abelian integrals with the help of Jacobi's inversion problem as stated by Zverovich in [123]. So, the solution to the problem of factorization of (3.20) is known at least formally, and it possibly can be used for practical purposes. In some special cases, this scheme was realized in a series of articles by Antipov and Silvestrov, and the solution to the factorization problem (as well as to the corresponding problems arising in application) was constructed in an explicit form in terms of special functions (see [124,125] and references therein). The idea of reducing a matrix factorization to a scalar one on a Riemann surface is also developed in [126]. In [127], a simplified method for treating matrices of Moiseev's type is proposed. An additional restriction is that a m (α) are algebraic functions. The factorization problem is transformed by Hurd's procedure [49,128,129] into a Riemann-Hilbert boundary value problem on a set of cuts. The method relies on formulation of ordinary differential equation with an unknown coefficient. It is shown that the proposed procedure in the Khrapkov case is equivalent to the standard solution method.

(f) Factorization of matrices with exponential terms in the entries (i) Integral equation on a finite interval
Consider the system of the convolution equations on a finite interval (3.22) note the similarity and differences with the basic integral Wiener-Hopf equation (1.1). The application of the Fourier transform leads to the necessity to factorize the 2n × 2n matrix function of type where the n × n matrix K(z) is the Fourier transform of any extension k onto the whole real line R of the kernel k 0 . The factorization of matrix (3.23) is studied in [130] under condition that the operator (Bϕ)(t) := ϕ(t) + a 0 k 0 (t − s)ϕ(s)ds, (3.24) is invertible in L n×1 (0, a) and it is a priori known that partial indices of the matrix A(λ) are equal to zero. In this case, the corresponding factors are found explicitly.

(ii) A system of integral equations with an exponential kernel
In [131], the following system of two convolution-type equations on the real semi-axes is considered: (3.25) where the kernel k is the following matrix: Fourier transformation of this system leads to a functional equation with matrix kernel Factorization of matrices of such type, as well as some other matrices involving exponential elements, is given in [131].

(iii) AKM approach
The name of this approach is simply an abbreviation of the authors' names in [132] which considers factorization of 2 × 2 non-rational matrix functions of the form for an arbitrary real parameter k. Matrices of this type arise as (modified) scattering matrices for the one-dimensional Schrödinger equation and some related Schrödinger-type equations. For the above discussion, the so-called scattering matrix

S(t) = T(t) R(t) L(t) T(t)
, (3.29) plays an important role. Here T is called the transmission coefficient, and R, L are reflection coefficients from the right and from the left, respectively. It is supposed that the following conditions are satisfied  (t; k), as the function of t ∈ R, belongs to a suitable Banach algebra of 2 × 2 matrix functions within which factorization is possible. This may be the Wiener algebra W = W 2×2 or algebra of functions f (t) such that f * (ξ ) = f (i((1 + ξ )/(1 − ξ ))) ∈ H α (T), α ∈ (0, 1) (which we denote H 2×2 α ).
First the following statement was proved: let us consider where q(∞) = 0, and either q ∈ W or q ∈ H α (T), α ∈ (0, 1). Then W(t) has a unique (right) canonical factorization where either W ± ∈ W 2×2 or W ± ∈ H 2×2 α , and W ± (∞) = E 2 . It leads to the following result: let matrix (3.28) be a unitary matrix for all t ∈ R satisfying the following conditions: (i) T(t) = 0 for all t ∈ R, (ii) T(t) can be continued to a meromorphic function on Π + with continuous boundary values on R, and T(∞) = 1. (iii) T(t), R(t) and L(t) belong to either W or to H α , 0 < α < 1.
Then G(t; k) has a (right) factorization with equal partial indices κ 1 = κ 2 = κ, where κ is the number of zeros minus the number of poles of T(z) in Π + .
The factorization has the form where W − , W + are factors of the canonical decomposition of the matrix G(t; k)/T(t) and T − , T + are product factors of the scalar function T(t): Some approximate methods for matrices with exponential entries are considered in §4(c).

(g) Miscellaneous classes (i) EF algorithm
An algorithm based on the successive solution of the scalar boundary value problems was proposed by Feldman et al. [24]. They construct the solution of the factorization problem (2.13) when the curve L is either the unit circle T or the real line R.
In the case L = T the 2 × 2 matrix investigated is a non-singular matrix of the following type , t ∈ T, (3.31) satisfying the following conditions: (i) b, c, d belong to the Wiener algebra W(T) of absolutely converging Fourier series on T; (ii) b(t) = p(t)/q(t), where p ∈ W + (T) (i.e. is a Fourier series with vanishing coefficients of negative indices), and q is a polynomial q(t) = k j=1 (t − a j ) m j with distinct zeros a j lying in the unit disc, |a j | < 1.
In the case L = R, the Wiener algebra W(R) (and its subalgebras W + (R), W − (R)) of the functions φ(λ), −∞ < λ < +∞, are defined in the standard way where c is an arbitrary complex constant and ϕ ∈ L 1 on its domain of integration.
The matrix under consideration is the following: where p ∈ W + (R), and q is a rational function q(t) = k j=1 ((t − a j )/(t + iγ j )) m j with distinct zeros a j laying in the upper half-plane, i.e. Im a j > 0, and γ j > 0.
In [133], this algorithm is applied to the solution of two problems: (a) the factorization of matrix (3.31) with ratio f (t) = p(t)/q(t) being meromorphically extended either inside the unit disc or outside of it; (b) the factorization of the matrix (3.30) with q ∈ W(T) and q 2 being a rational function free of zeros and poles on T.

(ii) Linear-fractional problem approach
In [134], a constructive method of factorization was proposed, based on the relation of the homogeneous Riemann-Hilbert boundary value problem (2.15) and the so-called linear-fractional problem.
Consider a homogeneous Riemann-Hilbert boundary value problem with a non-singular Hölder continuous 2 × 2 matrix coefficient ) T the components of the solution to (3.34) and introduce the following functions: . (3.36) These functions are the components of a piece-wise meromorphic solution to a so-called linearfractional boundary value problem The approach in [134] (see also [135]) is based on different possibilities arising from the study of solvability and representation of the solution to (3.37). These situations are listed in [134,135]. The representation of the solution to the corresponding factorization problem is given for each special case. In [136], this approach is applied to the study of the Riemann-Hilbert boundary value problem with a 3 × 3 matrix coefficient G(t).

(h) Rawlins-Williams class of matrix Wiener-Hopf problems
In acoustics and electromagnetism, the matrix kernel A(α) sometimes is a function only of γ (α) = (k 2 − α 2 ) 1/2 . This motivates the class considered by Rawlins & Williams in [137] A , (3.38) where F, G and H are analytic functions (except possibly when γ = 0). The way that this class of matrices is solved is by first transforming them to a matrix Riemann-Hilbert problem on a half line (along a cut emanating from one of the branch-points of γ (α)). It then transpires that this system can be decoupled into two scalar Riemann-Hilbert problems by Hurd's method [128], which can be solved explicitly. Another class of matrices was considered by Jones in [138]: where f (s), g(s) are analytic in the strip and e 1 , e 2 and e 3 are analytic functions in the whole complex plane. A commutative factorization is constructed. There are also extensions to non-commutative factorization by pre-multiplying by analytic matrices. Although it is not obvious, the class of matrices considered by Rawlins and Williams can, in fact, be reduced to Jones' class [138]. For a review of other classes which originate from applications [49].

Approximate procedures
Owing to the lack of a general exact constructive factorization procedure, there has been considerable interest in approximate methods for Wiener-Hopf matrix factorization. Numerous approaches have been proposed in the literature [41,68,93,114,127,131,[139][140][141][142][143][144][145][146][147][148]. We describe some of the more widely applicable classes here. Most approximate constructive methods rely on approximating a Wiener-Hopf problem by a class where exact constructive methods exist. There are also methods which reduce the Wiener-Hopf problem to a different equation such as Fredholm integral equation [94,95,140] but this is outside the scope of this survey. Suppose there are two functions Question 4.1. Is it true that ifK(α) approximates K(α) thanK ± (α) approximates K ± (α)? As stated the question is too vague to have a definitive answer, but nevertheless this kind of question motivates the search for approximate solutions. For scalar functions and L p norm if |K(α) −K(α)| p ≤ p there are bounds for |K ± (α) −K ± (α)| p in terms of the p and computable quantities of K(α) [141]. The general answer for a matrix valued function K(α) is negative. This can be the case due to a choice of norm or it can be due to an instability.
The simplest example of instability is obtained by mapping an example [9] from the unit circle to the real line. Consider a diagonal matrix function with partial indices [1, Perturbing the matrix we have ⎛ This example demonstrates that a small perturbation can change the factors by an arbitrary amount and can also change the partial indices (from {1, −1} to {0, 0}). This is significant because the partial indices are uniquely defined. Note that the sum of the partial indices remains the same; this is true in general, by considering the determinant of both sides it is reduced to scalar factorization. For scalar factorization of function f the index is the winding number of the curve (Re f (t), Im f (t)) t ∈ R hence is stable under perturbations. The partial indices are linked to the growth at infinity of the Wiener-Hopf factorization [149].
The following surprising theorem provides the necessary and sufficient conditions for the partial indexes to remain the same under small enough perturbations. In fact, this condition is enough to ensure the stability of factors in the Wiener norm.

Theorem 4.2 (Shubin [7]).
Assume the matrix function G to have a Wiener-Hopf factorization and the tuple of its partial indices to be stable. Then, for every > 0 there exists a δ > 0 such that, for ||F − G|| < δ, the matrix function F admits a factorization in which ||F ± − G ± || < .
An obstacle in using this result in applications is that one cannot in general determine the partial indices without constructing the factorization. Thus, in general it is impossible to prove that some approximate factorization of a matrix function is a good approximation. There are some specific results for Khrapkov-Daniele matrix functions [148] and for small perturbations to the identity matrix [145]. This will be discussed in more detail in the sections below. Interestingly, all three methods below reduce the matrix factorization problem to a series of scalar additive Wiener-Hopf splittings. This, in some sense, avoids problems with instabilities of matrix factorization with unstable indices.

(a) Rational approximation
It is attractive to consider a rational matrixK(α) in (4.1) since its factorization can be computed exactly ( §3(a)), without the need for a Cauchy-type integral (2.12) representation. Approximate rational solutions were considered early on [150] but they were mainly constructed using ad hoc observations ([17], ch. 4.5). In 2000, a systematic way of approximating the Wiener-Hopf equations was developed using a two point Padé approximation (the two points being 0 and ∞ on the real axis) [151]. The power of the method came from the fact that not the whole matrix was approximated but only specific parts and then pole removal method was applied ( §3(a)). Since then it proved popular and found applications in different branches of mathematics [93,152,153] including finance [51].
From the last paragraph, it is clear that it is desirable to approximate a matrix function, the whole or various parts, by rational functionsK(α). A natural question to ask: which K(α) can be approximated by rational matrix functionsK(α)? For this question, we have to specify in what sense are K(α) andK(α) close. It is not an easy task to decide what the correct norm is to consider, but good candidates are L p norms or Sobolev spaces [154]. We will review some known facts about approximation by rational functions, mostly by Padé approximants. Functions with branch cuts (multi-valued) cannot be 'fully' approximated by rational functions which are single-valued. But for a function analytic at infinity the next best result is true, the maximal domain of K(α) where the function has a single-valued branch is the domain of convergence of the diagonal Padé approximants for K(α). In fact, most of the poles tend to the boundary of the domain of convergence and lie on the system of cuts that makes the function single-valued [155,156]. Numerical experimentation has confirmed this holds true for even small degree of rational approximation and, what is more, often the interlacing of poles and zeros occurs on the branch cut [141,157]. This simulates the behaviour of the branch cut; the values on either side of the cut made by poles and zeros has a jump that tends to the exact value as the rational approximation degree is increased [158].
The practical implementation of the rational approximations has now been well developed. In particular there are numerous algorithms on Chebfun (Matlab). One of the latest additions in Chebfun is the very fast AAA rational approximation [159]. It is also useful to construct a mapping of the real line to the unit interval since most existing algorithms are for bounded intervals [141]. The main difficulty with rational approximations arise when a branch cut goes though infinity as in γ (α) = α 2 − k 2 (which is frequently encountered in diffraction problems). The problem arises because the real line crosses the branch cut at infinity. Hence rational approximation can no longer be accurate on the whole contour (the real line). It can still be very accurate on a bounded domain and hence rational approximations of γ (α) are used in acoustics far-field calculations [144].

(b) Asymptotic Wiener-Hopf factorization
Let M be a matrix function analytic around the real axis R. Let I be the identity matrix and ε ∈ [0, 1]. We can write M(x) = I + εG(x), (4.5) and then the asymptotic algorithm produces a good approximate Wiener-Hopf factorization if εG(x) is small compared to I [145,147]. Asymptotic Wiener-Hopf factorization looks for an approximate factorization of the form We define the error j at a step j as the term neglected in the approximation: For the first step, we are looking for G − 0 and G + 0 , 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 G − 1 , G + 1 and thus 2 : In the same manner, arbitrary order approximations are constructed and, under some conditions, its convergence can be established [145] and used in [62] for the factorization on the unit circle. Furthermore, this technique can be extended for the case of a set of stable partial indices [41] or, in some cases, even for unstable sets [160].

(c) Neglecting coupling and iteration
The main difficulty in solving matrix Wiener-Hopf equation arises from the fact that there are coupled scalar equations. In some situations the coupling is small and a reasonable approximation can be obtained by neglecting it. In other words the original matrix is close to being diagonal or triangular. More often the coupling is too strong to give a desired approximation see Section 4.4 in [17] or [161], nevertheless the same ideas can still be used. One way is through the use of iterative methods, which in the context of Wiener-Hopf equations were proposed early on [162,163]. One of the reasons that they did not gain popularity in the past is that they are computationally intensive. But recent advances in computing Cauchy-type transforms, and significant computational resource, have now made it a practical approach [96,164,165]. In particular, it has been applied to a class of triangular matrix functions containing exponential factors [166]. The aim is to find functions Φ on the strip H. The functions A(α), B(α) and C(α) are known and L is a positive constant. The method has found applications in aeroacoustics [167], crack propagation problems [168]. It has been extended to n-dimensional matrix problems and applied to scattering from multiple colinear plates [169]. We would also like to note that there is a similarity of the above iteration methods with a concept in diffraction called the Schwarzschild series [170][171][172][173]. This illustrates a useful viewpoint: many problems which give rise to a matrix Wiener-Hopf equation can be thought of as an interaction of a number of distinct problems, each of which would lead to a scalar Wiener-Hopf equation [174].

Related methods and open problems
To solve Wiener-Hopf-type problems which arise in applications, the following are desirable For problems of scalar Wiener-Hopf-type, the above points can be considered relatively complete. In the case of a classical matrix Wiener-Hopf equations the first one is straightforward while the remaining points continue to offer future avenues of research. It is possible to develop alternative methods for solving equations of Wiener-Hopf-type, for example based on Fredholm integral equation theory [140,175,176], and also more theoretical operator based approaches to convolution integral equations [154,177].
There are equations of Wiener-Hopf type which have not been reviewed here. For example, PDEs on wedge domains are not immediately in a standard class, but can be stated as generalized Wiener-Hopf equations if appropriate transformations/mapping are considered [178,179]. The difficulty with wedge domains is that the Fourier transform is not the natural operator to use, hence a mapping is required. Alternatively, the Mellin transform [180] could be employed if the respective differential operators are of the same homogeneity, and a combination of the Mellin and Fourier transform could be used for multilayered-multiwedged domains [28,181,182]. There is also a link to the unified transform method which provides a systematic way of deriving a natural transform for some polygonal domains [71]. A natural but difficult generalization is to consider (1.1) in more variables (a double integral) [183,184] leading to interesting Wiener-Hopf related methods in multi-variate complex analysis [185,186]. A different direction is when the kernel (1.1) has support on the half-line but is not of convolution type (the equation is referred to as being of Hammerstein type) [187]. There is also a wealth of literature on spectral properties of Toeplitz matrices/operators [188][189][190][191] which is a related topic but outside the scope of this review.
Interest in Wiener-Hopf techniques has been steady ever since the 1970s, judging by the number of papers containing 'Wiener-Hopf' on mathscinet, with a total around 3000 to date (it is likely to be a significant underestimate since many of the references for this review either do not contain 'Wiener-Hopf' in the title, or are published in engineering and physics journals). The related Riemann-Hilbert boundary value method has also been popular with around 4,600 publications; with a small but steady increase in number of publication each year overall. The Wiener-Hopf method could be viewed through many lenses, has many related problems, arises in numerous applications and still has many interesting open problems. It can be expected that further work in coming years will resolve some of the long-standing questions related to, and limitations of, the Wiener-Hopf method, and hence extend its relevance and applicability still further.
Data accessibility. This article has no additional data. Authors' contributions. All authors contributed to the manuscript. Competing interests. We declare we have no competing interests.