G. Germano, C. E. Phelan, D. Marazzina, G. Fusai

# Solution of Wiener-Hopf and Fredholm integral equations by fast Hilbert and
Fourier transforms

Guido Germano Carolyn E. Phelan Daniele Marazzina Gianluca Fusai  Department
of Computer Science, University College London
Systemic Risk Centre, London School of Economics and Political Science
<EMAIL_ADDRESS>Department of Computer Science, University College London
<EMAIL_ADDRESS>(corresponding author) Dipartimento di Matematica,
Politecnico di Milano<EMAIL_ADDRESS>Business School (formerly
Cass), City, University of London
Dipartimento di Studi per l’Economia e l’Impresa, Università del Piemonte
Orientale Amedeo Avogadro, Novara<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present numerical methods based on the fast Fourier transform (FFT) to
solve convolution integral equations on a semi-infinite interval (Wiener-Hopf
equation) or on a finite interval (Fredholm equation). We extend and improve a
FFT-based method for the Wiener-Hopf equation due to Henery, expressing it in
terms of the Hilbert transform, and computing the latter in a more
sophisticated way with sinc functions. We then generalise our method to the
Fredholm equation reformulating it as two coupled Wiener-Hopf equations and
solving them iteratively. We provide numerical tests and open-source code.
Wiener-Hopf, Fredholm, integral equation, fast Fourier transform, fast Hilbert
transform.

##  1 Introduction

We consider the linear integral equation of convolution type with constant
limits of integration

$\lambda f(x)-\int_{a}^{b}k(x-x^{\prime})f(x^{\prime})dx^{\prime}=g(x),\quad
x\in(a,b),$ (1)

where $f(x)$ is the unknown function, $k(x)$ is a given kernel, and $g(x)$ is
a given so-called forcing function. The domain of $f(x)$ and $g(x)$ is
$(a,b)$, the domain of $k(x)$ is $(a-b,b-a)$; the endpoints can be included if
they are finite. If $a=-\infty$ or $b=+\infty$ Eq. (1) is called a Wiener-Hopf
equation (Wiener & Hopf,, 1931; Noble,, 1958; Krein,, 1962; Polyanin &
Manzhirov,, 1998; Lawrie & Abrahams,, 2007); if both integration limits are
finite, it is called a Fredholm equation (Fredholm,, 1903; Whittaker &
Watson,, 1927; Polyanin & Manzhirov,, 1998). The latter case is also called a
Wiener-Hopf equation on a finite interval (Voronin,, 2004) or, because of an
application in electrotechnics, a longitudinally modified Wiener-Hopf equation
(LMWHE), while the former case is also called a classical Wiener-Hopf equation
(CWHE) (Daniele & Zich,, 2014). If $\lambda=0$ it is an equation of the first
kind; if $\lambda\neq 0$ it is an equation of the second kind. In the latter
case it can be assumed that $\lambda=1$, dividing the kernel and the forcing
function by values of this parameter different from 1. Historically these
equations arose in physics, e.g. to describe diffraction in the presence of an
impenetrable wedge or of planar waveguides (Daniele & Lombardi,, 2007), but
also for problems in crystal growth, fracture mechanics, flow mechanics (Choi
et al.,, 2005), geophysics, and diffusion (Lawrie & Abrahams,, 2007). The
connection of the Wiener-Hopf equation with probabilistic problems was noticed
by Spitzer, (1957) and is discussed by Feller, (1971), together with the
application of Fourier transform methods to stochastic processes. More
recently these equations have become of interest in finance to price
discretely monitored path-dependent options like barrier, first-touch,
lookback (or hindsight), quantile and Bermudan options (Fusai et al.,, 2006;
Green et al.,, 2010; Fusai et al.,, 2012; Marazzina et al.,, 2012; Fusai et
al.,, 2016; Phelan et al.,, 2018, 2019, 2020). The Wiener-Hopf method is
employed also to solve a large collection of mixed boundary value problems
(Duffy,, 2008).

##  2 Mathematical tools

###  2.1 Fourier transform and projection operators

We define the Fourier transform of a function $f(x)$ as

$\widehat{f}(\xi)=\mathcal{F}_{x\to\xi}[f(x)](\xi)=\int_{-\infty}^{+\infty}f(x)e^{i\xi
x}dx$ (2)

and its inverse as

$f(x)=\mathcal{F}_{\xi\to
x}^{-1}[\widehat{f}(\xi)](x)=\frac{1}{2\pi}\int_{-\infty}^{+\infty}\widehat{f}(\xi)e^{-ix\xi}a\xi$
(3)

where $i$ is the imaginary unit. We are aware that it is advised to define
Fourier space in terms of frequency $\nu$ (see e.g Press et al.,, 2007) rather
than angular frequency or pulsation $\xi=2\pi\nu$ (this is physics terminology
if $x$ is interpreted as time). With $\nu$ the inverse transform lacks the
factor $1/(2\pi)$, making it more symmetric with respect to the forward
transform, and the Nyquist relation between grids in the normal and Fourier
spaces simplifies to $\Delta x\Delta\nu=1/N$, where $N$ is the number of grid
points, without a factor $2\pi$ on the right-hand side. We also know that
putting the minus sign in the exponent of the forward Fourier transform
stresses the relationship with the Laplace transform and is the more common
choice in fast Fourier transform (FFT) libraries, including the FFTW (Frigo &
Johnson,, 2005) used in Matlab. However, we chose the definition of the
Fourier transform used normally in the two main application fields of Eq. (1),
physics and finance. An inconvenience is that the Fourier transform
$\mathcal{F}$ of our equations translates into ifft in the Matlab code
examples given in the supplementary material, and the inverse transform
$\mathcal{F}^{-1}$ into fft.

We define the projection of a function $f(x)$ on the positive or on the
negative half-axis through the multiplication with the indicator function of
that set,

$\displaystyle f_{+}(x)$ $\displaystyle=$
$\displaystyle\mathcal{P}_{+,x}f(x)=1_{\mathbb{R}_{+}}(x)f(x)$ (4)
$\displaystyle f_{-}(x)$ $\displaystyle=$
$\displaystyle\mathcal{P}_{-,x}f(x)=1_{\mathbb{R}_{-}}(x)f(x).$ (5)

A function that, like $f_{+}(x)$, is 0 for $x<0$ and nonzero for $x>0$ is
called “causal”, because it can be used to describe the effect of something
that happens at $x=0$ and causes the function to become nonzero. The two half-
range Fourier transforms are

$\displaystyle\widehat{f_{+}}(\xi)$ $\displaystyle=$
$\displaystyle\mathcal{F}_{x\to\xi}[f_{+}(x)](\xi)=\int_{0}^{+\infty}f(x)e^{i\xi
x}dx$ (6) $\displaystyle\widehat{f_{-}}(\xi)$ $\displaystyle=$
$\displaystyle\mathcal{F}_{x\to\xi}[f_{-}(x)](\xi)=\int_{-\infty}^{0}\,f(x)e^{i\xi
x}dx.$ (7)

The order in which the operators are applied matters:

$\widehat{f_{+}}(\xi)=\mathcal{F}_{x\to\xi}[\mathcal{P}_{+,x}f(x)](\xi)\neq\widehat{f}_{+}(\xi)=\mathcal{P}_{+,\xi}\mathcal{F}_{x\to\xi}[f(x)](\xi).$
(8)

In other words, $\widehat{f_{+}}(\xi)$ is the Fourier transform of a function
$f(x)$ that vanishes for negative arguments $x$, but $\widehat{f_{+}}(\xi)$
does not vanish itself for negative arguments $\xi$, which instead happens
with $\widehat{f}_{+}(\xi)$; similarly for the $-$ case. The function
$\widehat{f_{+}}(\xi)$ is analytic (or holomorphic), i.e. locally given by a
convergent power series, in an upper complex half-plane that includes the real
line; the function $\widehat{f_{-}}(\xi)$ is analytic in a lower complex half-
plane that includes the real line. The half-range Fourier transforms can be
considered special cases of the Laplace transform,

$\widetilde{f}(s)=\mathcal{L}_{x\to
s}[f(x)](s)=\int_{0}^{+\infty}f(x)e^{-sx}dx,\quad s\in\mathbb{C},$ (9)

where $s=\pm i\xi$, while the Fourier transform can be considered a special
case of the bilateral or two-sided Laplace transform. Except possibly for
$x=0$, the indicator function $1_{\mathbb{R}_{+}}(x)$ coincides with the
Heaviside step function $\Theta(x)$, and $1_{\mathbb{R}_{-}}(x)$ with
$1-\Theta(x)$; $\Theta(x)$ is 1 if $x>0$ and 0 if $x<1$, while for $x=0$ it
can be assigned the value 0 (left-continuous choice), 1 (right-continuous
choice), or 1/2 (symmetric choice). When integrating as in Eqs. (6) and (7),
the value for $x=0$ matters only numerically and only if $x=0$ is a grid
point, as analytically the measure of a point is zero. Clearly the sum of the
two projections, Eqs. (4) and (5), is the full function,

$f_{+}(x)+f_{-}(x)=f(x),$ (10)

and the sum of the two half-range Fourier transforms, Eqs. (6) and (7), is the
full Fourier transform,

$\widehat{f_{+}}(\xi)+\widehat{f_{-}}(\xi)=\widehat{f}(\xi).$ (11)

###  2.2 Gibbs phenomenon

As explained in Section 2.1, we numerically implement the forward and inverse
Fourier transform using the FFTW library in Matlab. The ranges of $x$ and
$\xi$ cease to be infinite and continuous, and are approximated with grids of
size $N$. The other parameter which defines both grids, which we centre around
zero, is the truncation in the $x$ domain $x_{\max}$. The step is $\Delta
x=2x_{\max}/N$ and the $x$ grid is

$x_{n}=n\Delta x,\quad n=-\frac{N}{2},-\frac{N}{2}+1,\dots,\frac{N}{2}-1.$
(12)

The step of the $\xi$ grid is given by the Nyquist relation,
$\Delta\xi=2\pi/(N\Delta x)=\pi/x_{\max}$; the truncation in the $\xi$ domain
is $\xi_{\max}=\pi/\Delta x$ and the $\xi$ grid is

$\xi_{m}=m\Delta\xi,\quad m=-\frac{N}{2},-\frac{N}{2}+1,\dots,\frac{N}{2}-1.$
(13)

The discrete forward and inverse Fourier transforms are

$\displaystyle\widehat{f}(\xi_{m},\Delta x,N)$ $\displaystyle=\Delta
x\sum^{N/2-1}_{n=-N/2}f\left(x_{n}\right)e^{i\xi_{m}x_{n}}$ (14)
$\displaystyle f(x_{n},\Delta\xi,N)$
$\displaystyle=\frac{\Delta\xi}{2\pi}\sum^{N/2-1}_{m=-N/2}\widehat{f}\left(\xi_{m}\right)e^{-ix_{n}\xi_{m}}.$
(15)

The truncation of the sums in Eqs. (14) and (15) causes the Gibbs phenomenon.
For a detailed explanation of its effect on the solution to Wiener-Hopf type
equations see Phelan et al., (2019). In this case we must consider two main
issues: firstly, if the function $f(x)$ has a discontinuity, the truncation of
$\widehat{f}(\xi_{m},\Delta x,N)$ causes oscillations in
$f(x_{n},\Delta\xi,N)$ close to the discontinuity; secondly, the error away
from that discontinuity will decay with the grid size $N$ as
$|f(x_{n})-f(x_{n},\Delta\xi,N)|=O(1/N)$.

There have been many different approaches to solve or mitigate the Gibbs
phenomenon; see e.g. Vandeven, (1991); Gottlieb & Shu, (1997); Tadmor &
Tanner, (2005); Tadmor, (2007); Ruijter et al., (2015). As in Phelan et al.,
(2019), we apply a spectral filter in the Fourier domain, specifically the
exponential filter of Gottlieb & Shu, (1997)

$\displaystyle\sigma(\eta)=e^{-\vartheta\eta^{p}},$ (16)

where $p\in\mathbb{N}$ is even and $\eta=\xi/\xi_{\max}$. This function does
not strictly meet the usual filter requirements described e.g. by Vandeven,
(1991) as it does not go exactly to zero when $|\eta|=1$, nor do so its
derivatives. However, if we select $\vartheta>-\log\varepsilon$, where
$\varepsilon$ is machine precision, then the filter coefficients are within
computational accuracy of the requirements. Advantages of the exponential
filter are its good performance, its simple form, and the order of the filter
being equal to the parameter $p$ which is directly input into the filter
equation.

We also investigated the use of the Planck taper (McKechan et al.,, 2010),
which is defined piecewise as

$\displaystyle\sigma(\eta)=\begin{cases}0,&\eta\leq\eta_{1},\ \
\qquad\eta_{1}=-1\\\ \frac{1}{e^{\;z(\eta)}+1},\
z(\eta)=\frac{\eta_{2}-\eta_{1}}{\eta-\eta_{1}}+\frac{\eta_{2}-\eta_{1}}{\eta-\eta_{2}},&\eta_{1}<\eta<\eta_{2},\
\,\eta_{2}=\epsilon-1\\\ 1,&\eta_{2}\leq\eta\leq\eta_{3},\
\,\eta_{3}=1-\epsilon\\\ \frac{1}{e^{\;z(\eta)}+1},\
z(\eta)=\frac{\eta_{3}-\eta_{4}}{\eta-\eta_{3}}+\frac{\eta_{3}-\eta_{4}}{\eta-\eta_{4}},&\eta_{3}<\eta<\eta_{4},\
\,\eta_{4}=1\\\ 0,&\eta\geq\eta_{4}.\end{cases}$ (17)

Here, the value of $\epsilon$ gives the proportion of the range of $\eta$
which is used for the slope regions. In common with the findings by Phelan et
al., (2019), the Planck taper, whilst having some interesting characteristics
such as a flat central section and a filter order of $\infty$, when tested did
not offer any advantage over the exponential filter, so we did not pursue its
use any further.

###  2.3 Hilbert transform and Wiener-Hopf factorisation

The Hilbert transform (see e.g Pandey,, 1996; Vergara Caffarelli & Loreti,,
1999; King,, 2009) of $f(x)$ is defined as the Cauchy principal value of the
convolution with $1/(\pi x)$,

$\displaystyle\mathcal{H}_{x}f(x)$ $\displaystyle=$
$\displaystyle\mathrm{p.v.}\,\frac{1}{\pi
x}*f(x)=\mathrm{p.v.}\,\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{f(x^{\prime})}{x-x^{\prime}}dx^{\prime}$
(18) $\displaystyle=$ $\displaystyle\lim_{\epsilon\to
0^{+}}\frac{1}{\pi}\left(\int_{-1/\epsilon}^{-\epsilon}\frac{f(x^{\prime})}{x-x^{\prime}}dx^{\prime}+\int_{\epsilon}^{1/\epsilon}\frac{f(x^{\prime})}{x-x^{\prime}}dx^{\prime}\right).$
(19)

The principal value avoids that the improper integral evaluates to the
indefinite form $+\infty-\infty$. Since with this definition the Hilbert
transform often appears multiplied by the imaginary unit (see the following
equations), some authors such as Stenger, (1973) define the Hilbert transform
as the principal value of the convolution with $i/(\pi x)$. To stress that the
Hilbert transform is actually a functional like the Fourier and Laplace
transforms, one could write $\mathcal{H}_{x^{\prime}\to x}[f(x^{\prime})](x)$
instead than $\mathcal{H}_{x}f(x)$, but this is too cumbersome and also less
useful because, unlike $x$ and $\xi$ or $x$ and $s$, $x^{\prime}$ and $x$
belong to the same space. A subscript like $x\to\xi$ or $x$ can be omitted
when there is no misunderstanding about which variable the operators
$\mathcal{F,F}^{-1},\mathcal{P_{+},P_{-}}$ and $\mathcal{H}$ act on, notably
when the argument function depends on a single variable; this is mostly the
case here. The operator $i\mathcal{H}$ is its own inverse,

$(i\mathcal{H})^{2}f(x)=f(x);$ (20)

equivalently, $\mathcal{H}^{-1}=-\mathcal{H}$. The convolution theorem

$\mathcal{F}[f(x)*g(x)]=\widehat{f}(\xi)\widehat{g}(\xi),$ (21)

which maps the convolution product in real space to a simple product in
Fourier space, together with the Fourier transform (Weisstein,, 2021)

$\mathrm{p.v.}\,\mathcal{F}\frac{1}{\pi x}=i\,\mathrm{sgn}\,\xi$ (22)

enables to express the Hilbert transform through a forward and an inverse
Fourier transform,

$\mathcal{H}f(x)=\mathcal{F}^{-1}\big{[}i\,\mathrm{sgn}\,\xi\,\widehat{f}(\xi)\big{]}.$
(23)

Thus a fast method to compute the Hilbert transform numerically consists
simply in evaluating Eq. (23) through a forward and an inverse FFT. In the
next subsection, 2.4, we shall see more sophisticated numerical methods.
Swapping the forward Fourier transform with its inverse, Eq. (22) becomes

$\mathrm{p.v.}\,\mathcal{F}^{-1}\frac{1}{\pi\xi}=-i\,\mathrm{sgn}\,x,$ (24)

and Eq. (23) becomes

$i\mathcal{H}\widehat{f}(\xi)=\mathcal{F}\big{[}\mathrm{sgn}\,x\,f(x)\big{]};$
(25)

substituting $\mathrm{sgn}\,x=1_{\mathbb{R}_{+}}(x)-1_{\mathbb{R}_{-}}(x)$
(this is true also for $x=0$, while $\mathrm{sgn}\,x=2\Theta(x)-1$ is
fulfilled for $x=0$ only if $\Theta(0)=1/2$), and applying the definitions of
the half-range Fourier transforms, Eqs. (6) and (7), yields

$\widehat{f_{+}}(\xi)-\widehat{f_{-}}(\xi)=i\mathcal{H}\widehat{f}(\xi).$ (26)

This can be shown also evaluating the integral in Eq. (18) with contour
integration methods in the complex plane. Together, Eqs. (11) and (26) are
known as Plemelj-Sokhotsky relations (Pandey,, 1996; Vergara Caffarelli &
Loreti,, 1999; King,, 2009). They can be rearranged as

$\displaystyle\widehat{f_{+}}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}\widehat{f}(\xi)+i\mathcal{H}\widehat{f}(\xi)\big{]}$
(27) $\displaystyle\widehat{f_{-}}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}\widehat{f}(\xi)-i\mathcal{H}\widehat{f}(\xi)\big{]}$
(28)

or, with a different notation involving the Fourier-transform and projection
operators, as

$\displaystyle\mathcal{F}\mathcal{P_{+}}f(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}\mathcal{F}f(x)+i\mathcal{H}\mathcal{F}f(x)\big{]}$
(29) $\displaystyle\mathcal{F}\mathcal{P_{-}}f(x)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{[}\mathcal{F}f(x)-i\mathcal{H}\mathcal{F}f(x)\big{]}.$
(30)

Substituting $f(x)$ with $\mathcal{P}_{+}f(x)$ in Eq. (29) and $f(x)$ with
$\mathcal{P}_{-}f(x)$ in Eq. (30), and taking into account that projection
operators are idempotent, i.e., $\mathcal{PP}f(x)=\mathcal{P}f(x)$, shows that
the half-range Fourier transforms are eigenfunctions of the Hilbert transform
operator,

$\displaystyle i\mathcal{H}\widehat{f_{+}}(\xi)$ $\displaystyle=$
$\displaystyle\phantom{-}\widehat{f_{+}}(\xi)$ (31) $\displaystyle
i\mathcal{H}\widehat{f_{-}}(\xi)$ $\displaystyle=$
$\displaystyle-\widehat{f_{-}}(\xi).$ (32)

This is evident also substituting $f(x)$ with $f_{+}(x)$ or with $f_{-}(x)$ in
Eq. (25), or applying the operator $i\mathcal{H}$ to both sides of Eqs.
(27)–(28) and simplifying with Eq. (20). Eqs. (31) and (32) allow to obtain
Eq. (26) applying the operator $i\mathcal{H}$ to both sides of Eq. (11);
conversely, Eq. (11) can be reobtained applying $i\mathcal{H}$ to both sides
of Eq. (26). Eqs. (27) and (28) are invariant with respect to an application
of $i\mathcal{H}$ to both sides.

The key step in the Wiener-Hopf solution of Eq. (1) described in the following
section is the decomposition of a function $\widehat{f}$, i.e., the reverse of
Eq. (11),

$\widehat{f}(\xi)=\widehat{f_{+}}(\xi)+\widehat{f_{-}}(\xi).$ (33)

The factorisation of a function $\widehat{f}(\xi)$

$\widehat{f}(\xi)=\widehat{f_{+}}(\xi)\widehat{f_{-}}(\xi),$ (34)

which is required too, can be reduced to a decomposition by taking logarithms,

$\log\widehat{f}(\xi)=\log\widehat{f_{+}}(\xi)+\log\widehat{f_{-}}(\xi);$ (35)

this procedure is called logarithmic decomposition. As described by Rino,
(1970), Henery, (1974) and Bart et al., (2004), the decomposition can be
achieved by

$\displaystyle\widehat{f_{+}}(\xi)$ $\displaystyle=$
$\displaystyle\mathcal{FP_{+}F}^{-1}\widehat{f}(\xi)$ (36)
$\displaystyle\widehat{f_{-}}(\xi)$ $\displaystyle=$
$\displaystyle\mathcal{FP_{-}F}^{-1}\widehat{f}(\xi),$ (37)

as can also be seen from the definitions of the half-range Fourier transforms,
Eqs. (6) and (7). More generally by the Plemelj-Sokhotsky relations, Eqs. (27)
and (28) can be used (Stenger,, 1973). Eqs. (36) and (37) are a special case
of the latter if the Hilbert transform is computed through Eq. (25).

The definition of the two half-range Fourier transforms, Eqs. (6) and (7), can
be generalized splitting the $x$ axis around a constant $a\neq 0$. Feng &
Linetsky, (2008) showed how the shift theorem,

$\mathcal{F}_{x\to\xi}[f(x+a)](\xi)=\widehat{f}(\xi)e^{-ia\xi},$ (38)

can be exploited to generalise the Plemelj-Sokhotsky relations to

$\displaystyle\widehat{f_{+}}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{\\{}\widehat{f}(\xi)+e^{ia\xi}i\mathcal{H}_{\xi}\big{[}\widehat{f}(\xi)e^{-ia\xi}\big{]}\big{\\}}$
(39) $\displaystyle\widehat{f_{-}}(\xi)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\big{\\{}\widehat{f}(\xi)-e^{ia\xi}i\mathcal{H}_{\xi}\big{[}\widehat{f}(\xi)e^{-ia\xi}\big{]}\big{\\}}.$
(40)

It might be a good idea to write $\widehat{f_{+,a}}$ and $\widehat{f_{-,a}}$
on the left-hand side, but we will not do it to not overburden the notation,
since it will be clear from the context with respect to which parameter a
function is decomposed. In the above formulas Eq. (25) generalizes to

$\displaystyle
e^{ia\xi}i\mathcal{H}_{\xi}\big{[}\widehat{f}(\xi)e^{-ia\xi}\big{]}$
$\displaystyle=$ $\displaystyle
e^{ia\xi}\mathcal{F}_{x\to\xi}\big{[}\mathrm{sgn}(x)f(x+a)\big{]}$ (41)
$\displaystyle=$
$\displaystyle\mathcal{F}_{x\to\xi}\big{[}\mathrm{sgn}(x-a)f(x)\big{]}$ (42)
$\displaystyle=$
$\displaystyle\mathcal{F}_{x\to\xi}\big{[}(1_{(a,+\infty)}(x)-1_{(-\infty,a)}(x))f(x)\big{]},$
(43)

and thus it is easy to show that

$\displaystyle\lim_{a\to-\infty}e^{ia\xi}i\mathcal{H}_{\xi}\big{[}\widehat{f}(\xi)e^{-ia\xi}\big{]}=\widehat{f}(\xi)$
(44)
$\displaystyle\lim_{a\to+\infty}e^{ia\xi}i\mathcal{H}_{\xi}\big{[}\widehat{f}(\xi)e^{-ia\xi}\big{]}=-\widehat{f}(\xi)$
(45)

and that $\lim_{a\to-\infty}\widehat{f_{+}}(\xi)=\widehat{f}(\xi),\
\lim_{a\to+\infty}\widehat{f_{+}}(\xi)=0,\
\lim_{a\to-\infty}\widehat{f_{-}}(\xi)=0,\
\lim_{a\to+\infty}\widehat{f_{-}}(\xi)=\widehat{f}(\xi)$. These limits are
useful to retrieve the results for the classical Wiener-Hopf equation from
those for the Fredholm equation.

###  2.4 Fast Hilbert transform with sinc functions

Eq. (23) provides a straightforward method to evaluate numerically the Hilbert
transform. Since it is based on two FFTs, its computational cost is $O(N\log
N)$, where $N$ is the number of grid points, and thus is called fast. We
compared this method with the quadrature method described in King, (2009, Eqs.
(4.19)–(4.20)), where the summation is taken over every second point in order
to avoid the singularity which results when $x_{i}-x_{j}=0$. We tried various
quadrature weights, including Simpson’s rule and 3rd and 4th order quadrature
(Press et al.,, 2007). For our implementation, see the Matlab functions htq.m
and weights.m in the supplementary material. All weights give the same result
and have polynomial convergence with $N$. Therefore, as with quadrature the
computation speed is $O(N^{2})$, the FFT-based method is preferable because of
its higher speed.

An alternative, but equally fast $O(N\log N)$ approach to compute numerically
the Hilbert transform is based on the sinc expansion approximation of
analytical functions. Sinc functions are deeply studied in the books by
Stenger, (1993, 2011), where the author shows that a function $f$ analytical
in the whole complex plane and of exponential type with parameter $\pi/h$,
i.e.,

$|f(z)|\leq Ce^{\pi|z|/h},$ (46)

can be reconstructed exactly from the knowledge of its values on an equispaced
grid of step $h$. Defining the sinc functions

$S_{n}(z,h)=\frac{\sin(\pi(z-nh)/h)}{\pi(z-nh)/h},\quad n\in\mathbb{Z},$ (47)

the function $f$ admits the sinc expansion (Stenger,, 1993, Theorem 1.10.1)

$f(z)=\sum_{n=-\infty}^{+\infty}f(nh)S_{n}(z,h).$ (48)

Moreover, also its Fourier transform admits the sinc expansion

$\widehat{f}(\xi)=h\sum_{n=-\infty}^{+\infty}f(nh)e^{i\xi nh}\quad\textrm{if}\
|\xi|<\pi/h,$ (49)

while it is zero if $|\xi|\geq\pi/h$. Finally, also the integrals of $f$ and
$|f|^{2}$ can be written as sinc expansions,

$\int_{-\infty}^{+\infty}f(x)dx=h\sum_{n=-\infty}^{+\infty}f(nh),\quad\int_{-\infty}^{+\infty}|f(x)|^{2}dx=h\sum_{n=-\infty}^{+\infty}|f(nh)|^{2}.$
(50)

The above results show in particular that the trapezoidal quadrature rule with
step size $h$ is exact. Using the following result on the Hilbert transform of
the sinc functions (Feng & Linetsky,, 2008, Corollary 6.1)

$\mathcal{H}_{\xi}S_{n}(\xi,h)=\frac{1-\cos(\pi(\xi-nh)/h)}{\pi(\xi-nh)/h},$
(51)

also the Hilbert transform can be evaluated exactly as

$\mathcal{H}_{\xi}f(\xi)=\sum_{n=-\infty}^{+\infty}f(nh)\frac{1-\cos(\pi(\xi-
nh)/h)}{\pi(\xi-nh)/h}.$ (52)

This holds for a function $f$ that is analytic in the whole complex plane.
However, this can be used also to approximate a function that is analytic only
in a strip including the real axis, which is the case considered in this
article. As shown in Stenger, (1993, Theorems 3.1.3, 3.2.1 and 3.1.4), in this
case the trapezoidal approximation has a discretisation error that decays
exponentially with respect to $h$. Feng & Linetsky, (2008, Section 6.5) also
shows that the computation of the Hilbert transform via a sinc expansion can
be performed using the FFT. The idea is that to compute a discrete Hilbert
transform it is necessary to do matrix-vector multiplications involving
Toeplitz matrices. As is well known, this kind of multiplications can be
performed exploiting the FFT, once those matrices are embedded in a circulant
matrix as in Feng & Linetsky, (2008, Appendix B) and Fusai et al., (2012). In
particular, Feng & Linetsky, (2008, 2009) prove the following convergence
result: if a function is analytic in a suitable strip around the real axis,
then the discretisation error of its numerical factorisation or decomposition
decays exponentially with respect to the step size between discretisation
points $h$. Feng & Linetsky, (2009, Theorem 3.3) considers the computation of
the Hilbert transform and Feng & Linetsky, (2008, Theorems 6.5) and Feng &
Linetsky, (2009, Theorems 3.4) are concerned with the calculation of the
Plemelj-Sokhotsky relations from Eqs. (27) and (28) in particular. An
implementation of both the $O(N\log N)$ methods presented above can be found
in the fuction ifht.m in the supplementary material.

However, the error due to the infinite sum in Eq. (52) being truncated to the
number of FFT grid points depends on the shape of the function under
transform. This has been explored further in the error bounds developped in
Feng & Linetsky, (2008, Section 6.4.2) and Phelan et al., (2019, Section 3).
In the case of a function that decays exponentially as
$|\xi|\rightarrow\infty$, exponential convergence is obtained. However, if a
function decays polynomially then the truncation error convergence is only
polynomially decreasing. In their paper on lookback options (Feng & Linetsky,,
2009, Theorem 3.3) report a result by Stenger that proves the exponential
convergence of the discrete sinc-based fast Hilbert transform to the
continuous Hilbert transform. They also examine the truncation error,
specifically observing in a footnote that this converges exponentially only
under certain conditions, notably $f(x)\leq\kappa\exp(c|x|^{\nu})$ for some
$\kappa,c,\nu>0$. This algorithm can be obtained with an eigenfunction
expansion of $\mathcal{H}$ and is identical to the Kress and Martensen method,
which was introduced with a proof that its error converges exponentially
(Kress & Martensen,, 1970; Weideman,, 1995).

##  3 The classical Wiener-Hopf method for a convolution equation on a semi-
infinite interval

Consider Eq. (1) with $b=+\infty$. The lower integration limit $a$ can be set
to 0 shifting the $x$ scale horizontally by the constant $a$ to a new scale
$x^{\prime}=x-a$; the prime is dropped hereafter. The functions $f(x)$ and
$g(x)$, whose domain is $[0,+\infty)$, are extended to the whole real axis
defining $f_{0}(x)=0$ for $x<0$, $f_{0}(x)=f(x)$ for $\geq 0$ and $g_{0}(x)=0$
for $x<0$, $g_{0}(x)=g(x)$ for $x\geq 0$. Define moreover the auxiliary
function

$f_{1}(x)=\int_{0}^{\infty}k(x-x^{\prime})f(x^{\prime})dx^{\prime}=\int_{-\infty}^{+\infty}k(x-x^{\prime})f_{0}(x^{\prime})dx^{\prime},\quad
x<0,$ (53)

and $f_{1}(x)=0$ for $x\geq 0$, i.e.,
$\widehat{f_{1}}=\big{(}\widehat{k}\widehat{f_{0}}\big{)}_{-}$. Since $f_{0}$
and $g_{0}$ are $+$ functions and $f_{1}$ is a $-$ function, it is customary
to denote these functions $f_{+}$, $g_{+}$, and $f_{-}$ respectively. With
them Eq. (1) is extended to

$\lambda
f_{+}(x)-\int_{-\infty}^{+\infty}k(x-x^{\prime})f_{+}(x^{\prime})dx^{\prime}+f_{-}(x)=g_{+}(x),\quad
x\in\mathbb{R},$ (54)

or, with a more compact notation for the convolution,

$\lambda f_{+}(x)-k(x)*f_{+}(x)+f_{-}(x)=g_{+}(x),\quad x\in\mathbb{R}.$ (55)

The extension of the integration domain to the whole real axis does not affect
the equation and its solution on the positive half-axis, but allows to apply
the convolution theorem, Eq. (21), to obtain the equation in Fourier space,

$\widehat{l}(\xi)\widehat{f_{+}}(\xi)+\widehat{f_{-}}(\xi)=\widehat{g_{+}}(\xi),$
(56)

where $\widehat{l}(\xi)=\lambda-\widehat{k}(\xi)$, ($l(x)$ and
$\widehat{l}(\xi)$ are the functional derivatives of the equation with respect
to the solution in normal and Fourier space respectively). Dropping the
argument $\xi$ for brevity, factorising
$\widehat{l}=\widehat{l_{-}}\widehat{l_{+}}$ and dividing the equation by
$\widehat{l_{-}}$ gives

$\widehat{l_{+}}\widehat{f_{+}}+\widehat{l_{-}}^{-1}\widehat{f_{-}}=\widehat{l_{-}}^{-1}\widehat{g_{+}}.$
(57)

Defining

$\widehat{c}=\widehat{l_{-}}^{-1}\widehat{g_{+}}$ (58)

and decomposing it as $\widehat{c}=\widehat{c_{+}}+\widehat{c_{-}}$ yields
finally

$\widehat{l_{+}}\widehat{f_{+}}+\widehat{l_{-}}^{-1}\widehat{f_{-}}=\widehat{c_{+}}+\widehat{c_{-}}.$
(59)

The $+$ and $-$ components can be separated into

$\displaystyle\widehat{f_{+}}=\widehat{l_{+}}^{-1}\widehat{c_{+}}$ (60)
$\displaystyle\widehat{f_{-}}=\widehat{l_{-}}\widehat{c_{-}},$ (61)

which allows to obtain the sought solution from

$f_{+}=\mathcal{F}^{-1}\big{(}\widehat{l_{+}}^{-1}\widehat{c_{+}}\big{)},$
(62)

while $f_{-}$ was introduced as an auxiliary function and is not of further
interest.

The case with $a=-\infty$ is treated in a similar fashion. The upper
integration limit $b$ can be set to 0 shifting the $x$ scale horizontally by
the constant $b$ to a new scale $x^{\prime}=x-b$; the prime is dropped
hereafter. The functions $f(x)$ and $g(x)$, whose domain is $(-\infty,0]$, are
extended to the whole real axis defining $f_{0}(x)=f(x)$ for $x\leq 0$,
$f_{0}(x)=0$ for $x>0$ and $g_{0}(x)=g(x)$ for $x\leq 0$, $g_{0}(x)=0$ for
$x>0$. Define moreover the auxiliary function

$f_{2}(x)=\int_{-\infty}^{0}k(x-x^{\prime})f(x^{\prime})dx^{\prime}=\int_{-\infty}^{+\infty}k(x-x^{\prime})f_{0}(x^{\prime})dx^{\prime},\quad
x>0,$ (63)

and $f_{2}(x)=0$ for $x\leq 0$, i.e.,
$\widehat{f_{2}}=\big{(}\widehat{k}\widehat{f_{0}}\big{)}_{+}$. Now $f_{0}$
and $g_{0}$ are $-$ functions and $f_{2}$ is a $+$ function, so it is
customary to denote these functions $f_{-}$, $g_{-}$, and $f_{+}$
respectively. With them Eq. (1) is extended to equations identical to Eqs.
(54)–(62), except that the $+$ and $-$ indices are swapped. In particular, the
sought solution is obtained from

$\displaystyle\widehat{c}=\widehat{l_{+}}^{-1}\widehat{g_{-}}$ (64)
$\displaystyle
f_{-}=\mathcal{F}^{-1}\big{(}\widehat{l_{-}}^{-1}\widehat{c_{-}}\big{)}.$ (65)

A more elegant alternative to shifting the $x$ scale forth and back by the
constant $a$ or $b$ is to modulate the functions in Fourier space decomposing
$\widehat{c}$ with respect to this constant by the generalized Plemelj-
Sokhotsky relations, Eqs. (39) and (40). The function $\widehat{l}$ is always
factorized with respect to 0, while $\widehat{c}$ is decomposed with respect
to $a$ when $b=+\infty$ or to $b$ when $a=-\infty$. For details, see the
function whf_gmf_filt4.m in the supplementary material.

##  4 Generalisation of the Wiener-Hopf method to a convolution equation on a
finite interval: the Fredholm equation

###  4.1 Theory

In the Fredholm equation both integration limits $a$ and $b$ are finite;
either $a$ or, less usually, $b$ can be set to 0 shifting the $x$ scale, but,
unlike with the classical Wiener-Hopf equation described in the previous
section, we prefer not to modify any of the two integration limits; instead,
we will use the generalized Plemelj-Sokhotsky relations. The functions $f(x)$
and $g(x)$, whose domain is $[a,b]$, are extended to the whole real axis
defining $f_{0}(x)=f(x)$ for $x\in[a,b]$, $f_{0}(x)=0$ for $x\notin[a,b]$ and
$g_{0}(x)=g(x)$ for $x\in[a,b]$, $g_{0}(x)=0$ for $x\notin[a,b]$. The kernel
$k(x)$, whose domain is $[a-b,b-a]$, is extended to the whole real axis
defining $k_{0}(x)=k(x)$ for $x\in[a-b,b-a]$ and $k_{0}(x)=0$ for
$x\notin[a-b,b-a]$. Define moreover the two auxiliary functions

$f_{1}(x)=\int_{a}^{b}k(x-x^{\prime})f(x^{\prime})dx^{\prime}=\int_{-\infty}^{+\infty}k_{0}(x-x^{\prime})f_{0}(x^{\prime})dx^{\prime},\quad
x<a,$ (66)

$f_{1}(x)=0$ for $x\geq a$, i.e.,
$\widehat{f_{1}}=e^{ia\xi}\big{(}\widehat{k_{0}}\widehat{f_{0}}e^{-ia\xi}\big{)}_{-}=\widehat{f_{-}}$,

$f_{2}(x)=\int_{a}^{b}k(x-x^{\prime})f(x^{\prime})dx^{\prime}=\int_{-\infty}^{+\infty}k_{0}(x-x^{\prime})f_{0}(x^{\prime})dx^{\prime},\quad
x>b,$ (67)

and $f_{2}(x)=0$ for $x\leq b$, i.e.,
$\widehat{f_{2}}=e^{ib\xi}\big{(}\widehat{k}\widehat{f_{0}}e^{-ib\xi}\big{)}_{+}=\widehat{f_{+}}$.
Because $k_{0}(x)=0$ for $x\notin[a-b,b-a]$, $f_{-}(x)=0$ also for
$x<a-(b-a)=2a-b$ and $f_{+}(x)=0$ also for $x>b-(a-b)=2b-a$. Thus Eq. (1)
extends to

$\lambda
f_{0}(x)-\int_{-\infty}^{+\infty}k_{0}(x-x^{\prime})f_{0}(x^{\prime})dx^{\prime}+f_{-}(x)+f_{+}(x)=g_{0}(x)$
(68)

or, with a more compact notation for the convolution,

$\lambda f_{0}(x)-k_{0}(x)*f_{0}(x)+f_{-}(x)+f_{+}(x)=g_{0}(x),$ (69)

and upon Fourier transformation, setting
$\widehat{l}(\xi)=\lambda-\widehat{k_{0}}(\xi)$,

$\widehat{l}(\xi)\widehat{f_{0}}(\xi)+\widehat{f_{-}}(\xi)+\widehat{f_{+}}(\xi)=\widehat{g_{0}}(\xi).$
(70)

Eqs. (68)–(70) look similar to Eqs. (54)–(56), but now we have two auxiliary
functions, $\widehat{f_{-}}(\xi)$, which is $-$ with respect to any $c\geq a$,
and $\widehat{f_{+}}(\xi)$, which is $+$ with respect to any $d\leq b$, while
both the unknown function $\widehat{f_{0}}(\xi)$ and the forcing function
$\widehat{g_{0}}(\xi)$ are $+$ with respect to $a$ (or any number $\leq a$)
and $-$ with respect to $b$ (or any number $\geq b$). Therefore the usual
approach is to split Eq. (70) into two coupled Wiener-Hopf equations, one with
the origin shifted to $a$, the other with the origin shifted to $b$ (Green et
al.,, 2010). These functions typically involve the four redundant unknowns
$\widehat{f_{0}}(\xi)e^{-ia\xi}$, $\widehat{f_{+}}(\xi)e^{-ib\xi}$ (which are
$+$ functions), $\widehat{f_{0}}(\xi)e^{-ib\xi}$,
$\widehat{f_{-}}(\xi)e^{-ia\xi}$ (which are $-$ functions). In Section 4.2 the
functions $\widehat{f_{-}}(\xi)$ and $\widehat{f_{+}}(\xi)$ correspond to
$J_{-}$ and $J_{+}$ from Green et al., (2010), Eq. (2.51), while
$\widehat{c}_{1}$ and $\widehat{c}_{2}$ correspond to $P$ and $Q$ from Eqs.
(2.12) and (2.24) in that paper.

###  4.2 Iterative solution

We solved the system of integral equations described in Eqs. (66)–(68)
iteratively observing that, if we know $\widehat{f_{+}}(\xi)$ and subtract it
from both sides of Eq. (70) with the origin of the $x$ axis shifted to $a$,
the result looks like Eq. (56), so that we can use the method described in
Sec. 3 to obtain $\widehat{f_{-}}(\xi)$; similarly, if we know
$\widehat{f_{-}}(\xi)$ and subtract it from both sides of Eq. (70) with the
origin of the $x$ axis shifted to $b$, we can use the method described in Sec.
3 to obtain $\widehat{f_{+}}(\xi)$. Thus, once again dropping the argument
$\xi$ for brevity of notation, our procedure is to write Eq. (70) divided once
by $\widehat{l_{-}}$, as in Eq. (57), and once by $\widehat{l_{+}}$,

$\displaystyle\widehat{l_{+}}\widehat{f_{0}}+\widehat{l_{-}}^{-1}\widehat{f_{-}}+\widehat{l_{-}}^{-1}\widehat{f_{+}}=\widehat{l_{-}}^{-1}\widehat{g_{0}}$
(71)
$\displaystyle\widehat{l_{-}}\widehat{f_{0}}+\widehat{l_{+}}^{-1}\widehat{f_{-}}+\widehat{l_{+}}^{-1}\widehat{f_{+}}=\widehat{l_{+}}^{-1}\widehat{g_{0}},$
(72)

start from the guess $\widehat{f_{+}}=0$ in Eq. (71), set

$\widehat{c_{1}}=\widehat{l_{-}}^{-1}(\widehat{g_{0}}-\widehat{f_{+}}),$ (73)

decompose $\widehat{c_{1}}=\widehat{c_{1+}}+\widehat{c_{1-}}$ with respect to
$a$, and compute the approximations

$\displaystyle\widehat{f_{0}}=\widehat{l_{+}}^{-1}\widehat{c_{1+}}$ (74)
$\displaystyle\widehat{f_{-}}=\widehat{l_{-}}\widehat{c_{1-}}$ (75)

as $+$ and $-$ functions with respect to $a$; then turn to Eq. (72), set

$\widehat{c_{2}}=\widehat{l_{+}}^{-1}(\widehat{g_{0}}-\widehat{f_{-}}),$ (76)

decompose $\widehat{c_{2}}=\widehat{c_{2+}}+\widehat{c_{2-}}$ with respect to
$b$, and compute the new approximations

$\displaystyle\widehat{f_{0}}=\widehat{l_{-}}^{-1}\widehat{c_{2-}}$ (77)
$\displaystyle\widehat{f_{+}}=\widehat{l_{+}}\widehat{c_{2+}}$ (78)

as $+$ and $-$ functions with respect to $b$; and so on until the difference
between the values of $\widehat{f_{0}}$ at an iteration and the previous falls
below a threshold. An equivalent result is obtained starting from the guess
$\widehat{f_{-}}=0$ in Eq. (72) and the computation of $\widehat{c_{2}}$.
Notice that the iterations are performed looking for a fixed point on the
variables $\widehat{f_{-}}$ and $\widehat{f_{+}}$, while $\widehat{f_{0}}$ is
a side product output at each step, but not used to compute the next step. For
details, see the function whf_gmf_filt4.m in the supplementary material.

###  4.3 Other iterative solutions

In this journal, Henery, (1977) proposed an iterative solution of the Fredholm
equation, but presented only the theory without a numerical validation. In our
tests, its literal implementation does not work. The procedure can be mapped
to ours including a missing projection and an untold detail: the $y_{n}$ found
in the residual correction scheme are corrections to the solution and thus
must be added together. However, Henery, (1977) did not express the algorithm
in terms of the Hilbert transform and thus used only the simple implementation
with the sign function, but not with a sinc function expansion.

Margetis & Choi, (2006) presented an iterative solution limited to algebraic
kernel functions. Moreover, in the example they implemented which is based on
a steady advection-diffusion problem first suggested by Choi et al., (2005),
they noted that “this choice of source function and kernel causes fortuitous
algebraic simplifications.” Therefore this method, whilst interesting as an
iterative procedure, cannot be considered to have a general validity.

###  4.4 Noble’s matrix factorisation approach

To avoid the iterations, we tried to solve the two simultaneous Wiener-Hopf
equations casted in matrix form according to the classic approach of Noble,
(1958, pp. 153–157); see also Daniele, (1984) and Daniele & Zich, (2014, Sec.
1.5.2). We write Eq. (70) multiplied once by $e^{-ia\xi}$ and once by
$e^{-ib\xi}$ as

$\begin{pmatrix}\widehat{l}&e^{i(d-a)\xi}\\\
0&e^{i(d-b)\xi}\end{pmatrix}\begin{pmatrix}\widehat{f_{0}}e^{-ia\xi}\\\
\widehat{f_{+}}e^{-id\xi}\end{pmatrix}+\begin{pmatrix}0&e^{i(c-a)\xi}\\\
\widehat{l}&e^{i(c-b)\xi}\end{pmatrix}\begin{pmatrix}\widehat{f_{0}}e^{-ib\xi}\\\
\widehat{f_{-}}e^{-ic\xi}\end{pmatrix}=\begin{pmatrix}0&1\\\
\widehat{l}&e^{i(a-b)\xi}\end{pmatrix}\begin{pmatrix}0\\\
\widehat{g_{0}}e^{-ia\xi}\end{pmatrix}$ (79)

where $a\leq c$ and $d\leq b$. Convenient choices of the parameters $c$ and
$d$ are $c=a,\ d=b$; $c=b,\ d=a$; $c=d=a$; $c=d=b$. We choose $c=a,\ d=b$ and
write for short

$\widehat{\mathbf{L}_{1}}\,\widehat{\mathbf{f}_{+}}+\widehat{\mathbf{L}_{2}}\,\widehat{\mathbf{f}_{-}}=\widehat{\mathbf{L}_{2}}\,\widehat{\mathbf{g_{+}}}.$
(80)

Multiplying from the left with $\widehat{\mathbf{L}_{2}}^{-1}$ yields a matrix
version of Eq. (56),

$\widehat{\mathbf{L}}\,\widehat{\mathbf{f}_{+}}+\widehat{\mathbf{f}_{-}}=\widehat{\mathbf{g}_{+}},$
(81)

where

$\widehat{\mathbf{L}}=\widehat{\mathbf{L}_{2}}^{-1}\,\widehat{\mathbf{L}_{1}}=\widehat{l}^{-1}\begin{pmatrix}-e^{i(a-b)\xi}&1\\\
\widehat{l}&0\end{pmatrix}\begin{pmatrix}\widehat{l}&e^{i(b-a)\xi}\\\
0&1\end{pmatrix}=\begin{pmatrix}-e^{i(a-b)\xi}&0\\\
\widehat{l}&e^{i(b-a)\xi}\end{pmatrix}$ (82)

is a triangular matrix. Swapping the elements of $\widehat{\mathbf{f}_{+}}$
and $\widehat{\mathbf{f}_{-}}$ permutes the elements of
$\widehat{\mathbf{L}}$. If we knew how to factorise
$\widehat{\mathbf{L}}=\widehat{\mathbf{L}_{-}}\,\widehat{\mathbf{L}_{+}}$,
multiplying Eq. (81) from the left with $\widehat{\mathbf{L}_{-}}^{-1}$ would
lead finally to

$\widehat{\mathbf{L}_{+}}\,\widehat{\mathbf{f}_{+}}+\widehat{\mathbf{L}_{-}}^{-1}\,\widehat{\mathbf{f}_{-}}=\widehat{\mathbf{L}_{-}}^{-1}\,\widehat{\mathbf{g}_{+}},$
(83)

which is a matrix version of Eq. (57) and can be solved in a similar fashion
decomposing its right-hand side. The same result is obtained multiplying Eq.
(80) from the left with $\widehat{\mathbf{L}_{1}}^{-1}$ or Eq. (81) from the
left with $\widehat{\mathbf{L}}^{-1}$, yielding

$\widehat{\mathbf{f}_{+}}+\widehat{\mathbf{L}}^{-1}\,\widehat{\mathbf{f}_{-}}=\widehat{\mathbf{L}}^{-1}\,\widehat{\mathbf{g}_{+}},$
(84)

where

$\widehat{\mathbf{L}}^{-1}=\widehat{\mathbf{L}_{1}}^{-1}\,\widehat{\mathbf{L}_{2}}=\widehat{l}^{-1}\begin{pmatrix}1&-e^{i(b-a)\xi}\\\
0&\widehat{l}\end{pmatrix}\begin{pmatrix}0&1\\\
\widehat{l}&e^{i(a-b)\xi}\end{pmatrix}=\begin{pmatrix}-e^{i(b-a)\xi}&0\\\
\widehat{l}&e^{i(a-b)\xi}\end{pmatrix}$ (85)

If we knew how to factorise
$\widehat{\mathbf{L}}^{-1}=\widehat{\mathbf{L}_{+}}^{-1}\,\widehat{\mathbf{L}_{-}}^{-1}$,
multiplying Eq. (84) from the left with $\widehat{\mathbf{L}_{+}}$ would lead
again to Eq. (83).

Unfortunately, though a formula to factorise triangular $2\times 2$ matrices
has been published, see Jones, (1984, Eq. (21)) and Jones, (1991, Eq. (6)), an
oscillatory element like $e^{i(b-a)\xi}$ does not fulfil the conditions to
apply it: it is assumed that $+$ or $-$ factors remain $+$ or $-$ when
inverted, but the inverse of the $+$ function $e^{i(b-a)\xi}$ is the $-$
function $e^{i(a-b)\xi}$; see also Daniele & Zich, (2014, Sec. 4.3, Example
2).

###  4.5 Voronin’s matrix factorisation approach

More recently, a different matrix form of the two simultaneous Wiener-Hopf
equations has been proposed by Voronin, (2004). We present it with slight
modifications. Start from Eq. (70) and decompose the kernel,
$\widehat{k_{0}}=\widehat{k_{-}}+\widehat{k_{+}}$ (for simplicity we drop the
0 subscript from the right-hand side), obtaining

$(\lambda-\widehat{k_{-}}-\widehat{k_{+}})\widehat{f_{0}}+\widehat{f_{-}}+\widehat{f_{+}}=\widehat{g_{0}}.$
(86)

Multiply by $e^{-ia\xi}$, take the $+$ part, thus eliminating
$\widehat{f_{-}}$, which is a $-$ function with respect to $a$, and multiply
by $e^{ia\xi}$, yielding

$(\lambda-\widehat{k_{+}})\widehat{f_{0}}-e^{ia\xi}(\widehat{k_{-}}\widehat{f_{0}}e^{-ia\xi})_{+}+\widehat{f_{+}}=\widehat{g_{0}}.$
(87)

Multiply by $e^{-ib\xi}$, take the $-$ part, thus eliminating
$\widehat{f_{+}}$, which is a $+$ function with respect to $b$, and multiply
by $e^{ib\xi}$, yielding

$\lambda\widehat{f_{0}}-e^{ia\xi}(\widehat{k_{-}}\widehat{f_{0}}e^{-ia\xi})_{+}-e^{ib\xi}(\widehat{k_{+}}\widehat{f_{0}}e^{-ib\xi})_{-}=\widehat{g_{0}}.$
(88)

Define
$\widehat{\varphi_{1}}=\big{(}\widehat{k_{-}}+\frac{1-\lambda}{2}\big{)}\widehat{f_{0}}$
and decompose it with respect to $a$,
$\widehat{\varphi_{1}}=\widehat{\varphi_{1+}}+\widehat{\varphi_{1-}}$; define
$\widehat{\varphi_{2}}=\big{(}\widehat{k_{+}}+\frac{1-\lambda}{2}\big{)}\widehat{f_{0}}$
and decompose it with respect to $b$,
$\widehat{\varphi_{2}}=\widehat{\varphi_{2+}}+\widehat{\varphi_{2-}}$; this
gives

$\widehat{f_{0}}-\widehat{\varphi_{1+}}-\widehat{\varphi_{2-}}=\widehat{g_{0}}.$
(89)

The two coupled Wiener-Hopf equations are now obtained multiplying once by
$\widehat{k_{-}}e^{-ia\xi}$ and once by $\widehat{k_{+}}e^{-ib\xi}$,

$\begin{pmatrix}1-\widehat{k_{-}}&0\\\
-\widehat{k_{+}}e^{i(a-b)\xi}&1\end{pmatrix}\begin{pmatrix}\widehat{\varphi_{1+}}e^{-ia\xi}\\\
\widehat{\varphi_{2+}}e^{-ib\xi}\end{pmatrix}+\begin{pmatrix}1&-\widehat{k_{-}}e^{i(b-a)\xi}\\\
0&1-\widehat{k_{+}}\end{pmatrix}\begin{pmatrix}\widehat{\varphi_{1-}}e^{-ia\xi}\\\
\widehat{\varphi_{2-}}e^{-ib\xi}\end{pmatrix}=\begin{pmatrix}\widehat{k_{-}}\widehat{g_{0}}e^{-ia\xi}\\\
\widehat{k_{+}}\widehat{g_{0}}e^{-ib\xi}\end{pmatrix},$ (90)

for short

$\widehat{\mathbf{M}_{\mathrm{r}-}}^{-1}\widehat{\bm{\varphi}_{+}}+\widehat{\mathbf{M}_{\mathrm{r}+}}\,\widehat{\bm{\varphi}_{-}}=\widehat{\mathbf{g}}.$
(91)

Here one can see that the parameter $\lambda$ has been inserted in the
definition of $\widehat{\bm{\varphi}_{1}}$ and $\widehat{\bm{\varphi}_{2}}$ to
avoid that it appears in place of the numbers 1 in the diagonal elements of
$\widehat{\mathbf{M}_{\mathrm{r}-}}^{-1}$ and
$\widehat{\mathbf{M}_{\mathrm{r}+}}$, which would make these matrices singular
for $\lambda=0$. Multiplying from the left by
$\widehat{\mathbf{M}_{\mathrm{r}-}}$ yields

$\widehat{\bm{\varphi}_{+}}+\widehat{\mathbf{M}}\,\widehat{\bm{\varphi}_{-}}=\widehat{\mathbf{M}_{\mathrm{r}-}}\,\widehat{\mathbf{g}},$
(92)

where

$\displaystyle\widehat{\mathbf{M}}$ $\displaystyle=$
$\displaystyle\widehat{\mathbf{M}_{\mathrm{r}-}}\widehat{\mathbf{M}_{\mathrm{r}+}}=\frac{1}{1-\widehat{k_{-}}}\begin{pmatrix}1&0\\\
\widehat{k_{+}}e^{i(a-b)\xi}&1-\widehat{k_{-}}\end{pmatrix}\begin{pmatrix}1&-\widehat{k_{-}}e^{i(b-a)\xi}\\\
0&1-\widehat{k_{+}}\end{pmatrix}$ (93) $\displaystyle=$
$\displaystyle\frac{1}{1-\widehat{k_{-}}}\begin{pmatrix}1&-\widehat{k_{-}}e^{i(b-a)\xi}\\\
\widehat{k_{+}}e^{i(a-b)\xi}&1-\widehat{k}\end{pmatrix}.$

An equivalent result is obtained multiplying Eq. (91) from the left by
$\widehat{\mathbf{M}_{\mathrm{r}+}}^{-1}$. If we knew how to factorise
$\widehat{\mathbf{M}}=\widehat{\mathbf{M}_{\mathrm{l}+}}\widehat{\mathbf{M}_{\mathrm{l}-}}$,
multiplying Eq. (92) from the left by
$\widehat{\mathbf{M}_{\mathrm{l}+}}^{-1}$ would lead finally to

$\widehat{\mathbf{M}_{\mathrm{l}+}}^{-1}\widehat{\bm{\varphi}_{+}}+\widehat{\mathbf{M}_{\mathrm{l}-}}\,\widehat{\bm{\varphi}_{-}}=\widehat{\mathbf{M}_{\mathrm{l}+}}^{-1}\widehat{\mathbf{M}_{\mathrm{r}-}}\,\widehat{\mathbf{g}},$
(94)

which, like Eq. (83), is a matrix version of Eq. (57) and can be solved in a
similar fashion decomposing its right-hand side.

Unfortunately we are stuck again: though formulas to convert left ($+-$)
factorisations of $2\times 2$ matrices into right ($-+$) ones and vice versa
have been published by Jones, (1991, Eqs. (8)–(11)), the information contained
in Eq. (93) is not sufficient to apply them, so we cannot obtain
$\widehat{\mathbf{M}_{\mathrm{l}+}}$ and $\widehat{\mathbf{M}_{\mathrm{l}-}}$
from our knowledge of $\widehat{\mathbf{M}_{\mathrm{r}-}}$ and
$\widehat{\mathbf{M}_{\mathrm{r}+}}$.

###  4.6 Iterative solution based on Voronin’s approach

An iterative solution is possible also with Voronin’s approach. Write Eq. (89)
multiplied once by $\widehat{k_{-}}/(1-\widehat{k_{-}})$ and once by
$\widehat{k_{+}}/(1-\widehat{k_{+}})$,

$\displaystyle\widehat{\varphi_{1+}}-\frac{1}{1-\widehat{k_{-}}}\widehat{\varphi_{1-}}-\frac{\widehat{k_{-}}}{1-\widehat{k_{-}}}\widehat{\varphi_{2-}}$
$\displaystyle=$
$\displaystyle\frac{\widehat{k_{-}}}{1-\widehat{k_{-}}}\widehat{g_{0}}$ (95)
$\displaystyle\widehat{\varphi_{2-}}-\frac{1}{1-\widehat{k_{+}}}\widehat{\varphi_{2+}}-\frac{\widehat{k_{+}}}{1-\widehat{k_{+}}}\widehat{\varphi_{1+}}$
$\displaystyle=$
$\displaystyle\frac{\widehat{k_{+}}}{1-\widehat{k_{+}}}\widehat{g_{0}},$ (96)

in Eq. (95) set

$\widehat{c_{1}}=\frac{\widehat{k_{-}}}{1-\widehat{k_{-}}}(\widehat{g_{0}}+\widehat{\varphi_{2-}}),$
(97)

start from the guess $\widehat{\varphi_{2-}}=0$, decompose
$\widehat{c_{1}}=\widehat{c_{1+}}+\widehat{c_{1-}}$ with respect to $a$, and
compute the approximations

$\displaystyle\widehat{\varphi_{1+}}=\widehat{c_{1+}}$ (98)
$\displaystyle\widehat{\varphi_{1-}}=(1-\widehat{k_{-}})\widehat{c_{1-}}$ (99)

as $+$ and $-$ functions with respect to $a$, as well as

$\widehat{f_{0}}=\frac{1}{\widehat{k_{-}}+\frac{1-\lambda}{2}}(\widehat{\varphi_{1+}}+\widehat{\varphi_{1-}})=\frac{1}{\widehat{k_{-}}+\frac{1-\lambda}{2}}\big{(}\widehat{c_{1+}}+(1-\widehat{k_{-}})\widehat{c_{1-}}\big{)};$
(100)

then turn to Eq. (96), set

$\widehat{c_{2}}=\frac{\widehat{k_{+}}}{1-\widehat{k_{+}}}(\widehat{g_{0}}-\widehat{\varphi_{1+}}),$
(101)

decompose $\widehat{c_{2}}=\widehat{c_{2+}}+\widehat{c_{2-}}$ with respect to
$b$, and compute the new approximations

$\displaystyle\widehat{\varphi_{2-}}=\widehat{c_{2-}}$ (102)
$\displaystyle\widehat{\varphi_{2+}}=(1-\widehat{k_{+}})\widehat{c_{2+}}$
(103)

as $+$ and $-$ functions with respect to $b$, as well as

$\widehat{f_{0}}=\frac{1}{\widehat{k_{+}}+\frac{1-\lambda}{2}}(\widehat{\varphi_{2-}}+\widehat{\varphi_{2+}})=\frac{1}{\widehat{k_{+}}+\frac{1-\lambda}{2}}\big{(}\widehat{c_{2-}}+(1-\widehat{k_{+}})\widehat{c_{2+}}\big{)};$
(104)

and so on until the difference between the values of $\widehat{f_{0}}$ at an
iteration and the previous falls below a threshold. An equivalent result is
obtained starting from the guess $\widehat{\varphi_{1+}}=0$ in Eq. (96) and
the computation of $\widehat{c_{2}}$. Similarly to Sec. 4.2, the iterations
are performed looking for a fixed point on the variables
$\widehat{\varphi_{1}}$ and $\widehat{\varphi_{2}}$, while $\widehat{f_{0}}$
is a side product output at each step, but not used to compute the next step.
For details, see the function whf_gmf_v.m in the supplementary material.

Generally, work on factorisation has concentrated on finding methods for
particular matrix classes, often corresponding to specific applications. For
example, Kisil, (2015) has developped and analysed an approximate
factorisation approach for Daniele-Khrapkov matrics but, although several
matrix classes can be reduced to this case, this is not a general solution. A
digest of constructive matrix factorisation methods was published in this
journal by Rogosin & Mishuris, (2016), who described several different matrix
classes and the most suitable methods for their factorisation. This work,
whilst important, underlines that there is no general method which can be used
for all matrices.

##  5 Test cases

As we present a general solution to the Fredholm equation, rather than one
limited to a particular application, we provide several test cases to solve
Eq. (1) for $f(x)$. Although the methods developped herein can be applied to
both Fredholm and Wiener-Hopf equations, for the numerical tests we have
chosen to concentrate on solving examples of the Fredholm equation as it is
the more challenging case and encompasses Wiener-Hopf as a special example
when $a\to-\infty$ or $b\to+\infty$.

Solutions to Eq. (1) with simple closed form expressions for $f(x)$, $g(x)$
and $k(x)$ are not readily available. However, if we limit the requirement for
simplicity to $f(x)$ and $k(x)$, then closed form solutions for $g(x)$ in Eq.
(1) can be calculated. These $g(x)$ and $k(x)$ can be used as inputs to our
numerical method, whose accuracy can be measured by comparing the result with
$f(x)$. We selected $f(x)$ and $k(x)$ to give closed form expressions for
$g(x)$ and also to have Fourier transforms which are easily calculable. We
derived three solutions, with both $f(x)$ and $k(x)$ 1\. Gaussian, 2. Cauchy,
3. Laplace (bilateral exponential). As discussed in Section 2.2 the decay of
the functions as $x,\xi\rightarrow\infty$ can influence the error performance
of Fourier-based methods. The functions were therefore selected to be
exponentially decaying in both the state space and Fourier space (Gaussian),
to be polynomially decaying in the state space and exponentially decaying in
the Fourier space (Cauchy) and to be exponentially decaying in the state space
and polynomially decaying in the Fourier space (Laplace). The derivation of
$g(x)$ is described in the following sections.

###  5.1 Gaussian

We set $f(x)=k(x)=\frac{1}{\sqrt{\pi}}e^{-x^{2}}$. The expression for $g(x)$
is then

$\displaystyle g(x)$
$\displaystyle=\frac{1}{\sqrt{\pi}}e^{-x^{2}}-\frac{1}{\pi}\int^{b}_{a}e^{-(x-y)^{2}}e^{-y^{2}}dy$
$\displaystyle=\frac{1}{\sqrt{\pi}}e^{-x^{2}}-\frac{1}{\pi}e^{-\frac{x^{2}}{2}}\int^{b}_{a}e^{-2(y-\frac{x}{2})^{2}}dy$
$\displaystyle=\frac{1}{\sqrt{\pi}}e^{-x^{2}}-\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}\big{[}\Phi(2b-x)-\Phi(2a-x)\big{]},\quad
x\in[a,b],$ (105)

where $\Phi(\cdot)$ is the standard normal cumulative distribution function.

###  5.2 Cauchy

We set $f(x)=k(x)=\frac{1}{\pi(x^{2}+1)}$. The first step in the calculation
of $g(x)$ is to solve the integral

$\displaystyle g_{\mathrm{int}}(x)$
$\displaystyle=\frac{1}{\pi^{2}}\int^{b}_{a}\frac{1}{y^{2}+1}\frac{1}{(x-y)^{2}+1}dy.$
(106)

Using partial fractions,

$\displaystyle g_{\mathrm{int}}(x)=$
$\displaystyle\frac{1}{\pi^{2}}\int^{b}_{a}\frac{1}{y^{2}+1}\frac{1}{(x-y)^{2}+1}dy$
$\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}x(x^{2}+4)}\int^{b}_{a}\left(\frac{2y}{y^{2}+1}+\frac{x}{y^{2}+1}-\frac{2(y-x)}{(y-x)^{2}+1}+\frac{x}{(y-x)^{2}+1}\right)dy$
$\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}x(x^{2}+4)}\left[\log(y^{2}+1)+x\arctan(y)-\log[(x-y)^{2}+1]+x\arctan(y-x)\right]^{b}_{a}$
$\displaystyle=$
$\displaystyle\frac{1}{\pi^{2}x(x^{2}+4)}\left\\{\log\frac{(b^{2}+1)((a-x)^{2}+1)}{(a^{2}+1)((b-x)^{2}+1)}+\right.$
$\displaystyle\left.+x\left[\arctan(b)-\arctan(a)+\arctan(b-x)-\arctan(a-x)\right]\right\\}.$
(107)

This gives $g(x)$ in closed form:

$\displaystyle g(x)=$
$\displaystyle\frac{1}{\pi(x^{2}+1)}-\frac{1}{\pi^{2}x(x^{2}+4)}\left\\{\log\frac{(b^{2}+1)((a-x)^{2}+1)}{(a^{2}+1)((b-x)^{2}+1)}+\right.$
$\displaystyle\left.+x\left[\arctan(b)-\arctan(a)+\arctan(b-x)-\arctan(a-x)\right]\right\\},\quad
x\in[a,b].$ (108)

###  5.3 Laplace or bilateral exponential

We set $f(x)=k(x)=\frac{1}{2}e^{-|x|}$. In order to make the calculation of
$g(x)$ simpler, the values of $a$ and $b$ are restricted so that
$0<a<b<\infty$. Then the formula for $g(x)$ in closed form is

$\displaystyle g(x)$
$\displaystyle=\frac{1}{2}e^{-x}-\frac{1}{4}\int^{b}_{a}e^{-|x-y|}e^{-y}dy$
$\displaystyle=\frac{1}{2}e^{-x}-\frac{1}{4}\left[\int^{b}_{x}e^{(x-y)}e^{-y}dy+\int^{x}_{a}e^{-(x-y)}e^{-y}dy\right]$
$\displaystyle=\frac{1}{2}e^{-x}-\frac{1}{4}e^{x}\left(\int^{b}_{x}e^{-2y}dy+e^{-x}\int^{x}_{a}dy\right)$
$\displaystyle=\frac{1}{2}e^{-x}+\frac{1}{8}e^{x}\left[e^{-2y}\right]^{b}_{x}-\frac{1}{4}e^{-x}\big{[}y\big{]}^{x}_{a}$
$\displaystyle=\frac{1}{2}e^{-x}+\frac{1}{8}e^{x}\left(e^{-2b}-e^{-2x}\right)-\frac{1}{4}e^{-x}(x-a)$
$\displaystyle=e^{-x}\left[\frac{3}{8}+\frac{1}{8}e^{-2(b-x)}+\frac{1}{4}(a-x)\right],\quad
x\in[a,b].$ (109)

##  6 Results

We used the following methods to recover $f(x)$ and produce the detailed
results shown in this section:

  1. 1.

4th order quadrature with preconditioner (Fusai et al.,, 2012; Press et al.,,
2007); see the Matlab functions quadrature.m and weights.m in the
supplementary material.

  2. 2.

Wiener-Hopf method using the sinc-based fast Hilbert transform with no zero
padding. In order to counteract the oscillations on the recovered function, we
used an exponential filter of order 8 on the final stage of the fixed point
algorithm. The maximum number of iterations of the fixed point algorithm is
set to 5. In fact, in most cases the final error level is achieved within 3
iterations. We discuss the use of the sinc-based fast Hilbert transform and
spectral filtering in Section 6.1.1 below.

  3. 3.

Wiener-Hopf method using the symmetric sign function in the fast Hilbert
transform, i.e. with zeros placed at both $\xi=0$ and
$\xi=\xi_{\min}=-\frac{N}{2}\Delta\xi$, similar to the method introduced by
Rino, (1970) and Henery, (1974) and tested by Fusai et al., (2016).

  4. 4.

Wiener-Hopf method with Voronin’s variant using the symmetrical sign function
for the Hilbert transform.

###  6.1 Results for the Gaussian test case

We first examine the performance of the different numerical methods with the
Gaussian test case, with particular emphasis on the method used to implement
the Hilbert transform.

####  6.1.1 Sinc-based fast Hilbert transform and spectral filtering

In the financial pricing applications described by Feng & Linetsky, (2008) and
Fusai et al., (2016) the sinc-based fast Hilbert transform has shown excellent
error convergence, especially when combined with a spectral filter as in
Phelan et al., (2019). However, when we consider its use for this application
we must take account of several ways in which the requirements differ from its
general use for finding solutions to Wiener-Hopf or Fredholm equations.

Firstly the pricing methods that were implemented using the Wiener-Hopf method
in Fusai et al., (2016), as devised by Green et al., (2010), use the analytic
continuation of $x$, i.e. they give results for values of $x$ both inside and
outside the barriers (the integration limits of Fredholm equation). This means
that for these applications there is no requirement to truncate the functions
to the integration limits of $a$ and $b$ in the state space, unlike the
problems presented as examples in this paper. The requirement to truncate the
function means that there is a jump discontinuity introduced in the state
space, meaning that, as described by Boyd, (2001), the function in the Fourier
space decays as a first order polynomial due to the Gibbs phenomenon. As
explained in Stenger, (1993) this polynomial decay means that the sinc-based
fast Hilbert transform no longer obtains an error which is exponentially
convergent with grid size but rather converges polynomially. This is in
contrast with the aforementioned finance-based papers from Green et al.,
(2010) and Fusai et al., (2016) where the Fourier domain functions subject to
the sinc-based fast Hilbert transform are exponentially decaying (or
polynomially so in the case of the VG process) and so excellent error
performance is achieved, especially in conjunction with a spectral filter to
solve the issue with the fixed point algorithm.

In contrast, here we solve the Fredholm and Wiener-Hopf equations as they were
originally formulated, i.e. the function is only defined for the range of the
integration $[a,b]$ and therefore the functions $g(x)$ and $k(x)$ must be
truncated to the ranges $[a,b]$ and $[a-b,b-a]$ respectively. This truncation
will introduce a jump in the functions which means that their Fourier
transforms now have first order polynomial decay. Therefore the truncation
error from the Hilbert transform will have a first order polynomial
convergence unless we can exploit some symmetry between the Fourier domain
functions for positive and negative $\xi$ as in Phelan et al., (2018), in
which case we may achieve second order polynomial convergence.

Moreover, there is a second important distinction to be made between the
general solution presented here and the work in the above literature. In the
finance literature the solutions to the Fredholm equation are used to
calculate the expectation of a further function, in this case the payoff
function. Therefore the exact errors in the function for individual values of
$x$ are not particularly important. Rather, the finance literature is
concerned with the average error, weighted according to the shape of the
payoff function. This also has particular importance when we are considering
the use of the sinc-based fast Hilbert transform described in Section 2.4
which was instrumental in achieving exponential error convergence with the
number of FFT grid points $N$ in Feng & Linetsky, (2008) and Fusai et al.,
(2016).

In Figures 1 and 2, we show results using the sinc-based fast Hilbert
transform with no filter for the Gaussian test case described in Section 5.1.
It is immediately obvious that, even for high values of $N$, oscillations are
visible in the numerical solution.

Figure 1: Numerical and analytical $f(x)$ using the sinc-based fast Hilbert
transform with no filter. Figure 2: Error in the numerical calculation of
$f(x)$ using the sinc-based fast Hilbert transform with no filter.

We can use a spectral filter to overcome the oscillations, however this can
have a negative effect on the accuracy of the numerical method, especially
close to the discontinuities in the state space; this is illustrated in
Figures 3–5. Figure 4 shows that the lower order filter gives a shallower
slope at the discontinuity, but has a stronger effect on the oscillations.
However, we can see from Figure 4 that, regardless of the order of the filter,
the overshoot at the discontinuity remains approximately the same. Figure 5
shows that a spectral filter removes the oscillations away from the
discontinuity and that the best results are achieved with a filter of order 8.
Although the behaviour of the numerical method using the sinc-based fast
Hilbert transform is not appropriate for a general solution to the Fredholm
equation due to the high errors at function discontinuities, it remains the
case that for applications where we are solely interested in a function value
away from any jumps this may be an appropriate method to use.

Figure 3: Numerical and analytical $f(x)$ using the sinc-based fast Hilbert
transform with an exponential filter. The parameter $p$ describes the order of
the filter. Figure 4: Numerical and analytical $f(x)$ using the sinc-based
fast Hilbert transform with an exponential filter, focussing on the
discontinuity at $x=0$. The parameter $p$ describes the order of the filter.
Figure 5: Error between the numerical and analytical calculation of $f(x)$
using the sinc-based fast Hilbert transform with an exponential filter. The
parameter $p$ describes the order of the filter. The scale has been chosen to
display the error away from the discontinuities of $f(x)$.

####  6.1.2 Sign-based fast Hilbert transform method

As an alternative to the sinc-based fast Hilbert transform, we examine the
method of Rino, (1970) and Henery, (1974), which was used also by Fusai et
al., (2016). It is based on the simple relationship between the Hilbert
transform and the Fourier transform given in Eq. (25).

The results for the Gaussian test case are shown in Figures 6–8; $f$ is the
analytic solution, $f_{\mathrm{sgn}0}$ and $f_{\mathrm{Vor}}$ are the
numerical solutions obtained with the Wiener-Hopf iterative method using the
fast Hilbert transform implemented with the symmetric sign function, the
latter in the Voronin variant. It is immediately apparent from Figure 6 that
neither implementation suffers from the overshoot that was seen using the
sinc-based fast Hilbert transform. However, looking at the discontinuity more
closely in Figure 7, we can see that we will have a peak error at a single
state-space grid point as the numerical solution increases to the final value
of $f(x)$ more slowly than the analytic function. However, unlike the sinc-
based function, where the extent of the oscillations depends not only on the
filter, but also the shape the function used, we can state here that as long
as the value of $x$ is at least one grid step away from the discontinuity, the
answer will be unaffected by the peak error. It is also interesting to note
that the error is symmetrical around the discontinuity when the iterative
Wiener-Hopf method is used, but not with the Voronin variant. This difference
is likely to account for the better performance seen in Figure 8 compared to
Figure 9. These display the error results away from the discontinuity and we
can see that, although there is some variation in the error across $x$, the
results for both methods are superior to those for the sinc-based fast Hilbert
transform.

Figure 6: Numerical and analytical $f(x)$ using the FFT based method with a
symmetrical sign function. Figure 7: Numerical and analytical $f(x)$ using the
FFT based method with a symmetrical sign function, focussing on the
discontinuity at $x=0$. Figure 8: Error of the numerical calculation of $f(x)$
for the new iterative Wiener-Hopf method using the sign-based fast Hilbert
transform. The scale has been chosen to display the error away from the
discontinuities of $f(x)$. The error is calculated by comparing the numerical
calculation to the analytic solution. Figure 9: Error of the numerical
calculation of $f(x)$ for the Voronin method using the sign-based fast Hilbert
transform. The scale has been chosen to display the error away from the
discontinuities of $f(x)$. The error is calculated by comparing the numerical
calculation to the analytic solution.

Although it is important to observe the functions which are calculated
numerically, when assessing the performance of the numerical methods, the
error convergence with CPU time and number of grid points $N$ is also
important. We measured this at 10%, 50% and 90% of the range between $a$ and
$b$; results for the Gaussian test case are shown in Figures 10–11.

The fastest converging method is the Wiener-Hopf iterative method using the
sign-based fast Hilbert transform, achieving an error of $O(1/N^{2})$. The
other methods exhibited $O(1/N)$ error convergence, with the method using the
sinc-based fast Hilbert transform with spectral filter achieving better
absolute error performance vs. $N$ but converging with CPU time almost
identically to the quadrature method. The $O(1/N)$ convergence achieved with
the sinc-based fast Hilbert transform is consistent with the error bound
described by Stenger, (1993) for a function with a first order discontinuity,
while the $O(1/N^{2})$ convergence seen for the sign-based variant is
consistent with that reported by Fusai et al., (2016).

Figure 10: Error convergence of the numerical methods vs. $N$ with the
Gaussian test case. Figure 11: Error convergence of the numerical methods vs.
CPU time with the Gaussian test case.

###  6.2 Results for Cauchy and Laplace test cases

Figures 12–15 compare the results for the test cases in Sections 5.2 and 5.3
for the iterative Wiener-Hopf method with the sinc- and sign-based fast
Hilbert transform methods; in the figures these are labelled
$f_{\mathrm{sinc}}$ and $f_{\mathrm{sgn}}$. We also compare the performance of
the iterative Voronin method with the symmetrical sign function, labelled
$f_{\mathrm{Vor}}$. An 8th order exponential filter was used with the sinc-
based fast Hilbert transform to counteract the oscillations, as described in
Section 6.1.1. As a benchmark we include results from 4th order quadrature
with preconditioner $f_{\mathrm{q}}$ (Fusai et al.,, 2012; Press et al.,,
2007), labelled $f_{\mathrm{q}}$, which was the previous state of the art. The
results for the Cauchy and Laplace test cases are consistent with those for
the Gaussian test case; the use of the sinc-based fast Hilbert transform
results in an overshoot at the function discontinuities and the sign-based
method results in a spot error at function discontinuities. We also notice
that the quadrature method has a spot error at the discontinuity, but this
effects a smaller range of $x$ than our new numerical methods. The reason for
this smaller range is that the Fourier-based methods need a truncation in
state space at $\pm 4(b-a)$ in order to avoid wrap-round effects. In contrast,
the range of $x$ for quadrature only needs truncation at the integration
limits $a$ and $b$.

Figure 12: Numerical and analytical $f(x)$ with the Cauchy test case. Figure
13: Numerical and analytical $f(x)$ with the Cauchy test case focussing on the
first discontinuity. Figure 14: Numerical and analytical $f(x)$ with the
Laplace test case. Figure 15: Numerical and analytical $f(x)$ with the Laplace
test case focussing on the first discontinuity.

We also measured the error convergence with the Cauchy and Laplace test cases
and the results are shown in Figures 16–19. These confirm the findings with
the Gaussian test case in Section 5.1 which showed that the best performing
method is the new iterative solution to the Wiener-Hopf equation with the
Hilbert transform implemented using the FFT with the symmetrical sign
function.

Figure 16: Error convergence of the numerical methods vs. $N$ with the Cauchy
test case. Figure 17: Error convergence of the numerical methods vs. CPU time
with the Cauchy test case. Figure 18: Error convergence of the numerical
methods vs. $N$ with the Laplace test case. Figure 19: Error convergence of
the numerical methods vs. CPU time with the Laplace test case.

##  7 Conclusion

We implemented four methods to solve general Fredholm equations of the second
kind and assessed their performance for three test cases with analytical
solutions. The methods are: 4th order quadrature with preconditioner Fusai et
al., (2012); Press et al., (2007); two iterative solutions that use the
Wiener-Hopf method with the sinc- or sign-based fast Hilbert transform; a
variant of the iterative method based on Voronin’s partial solution to the
Fredholm equation with the sign-based fast Hilbert transform.

Unlike an earlier application in option pricing with exponential convergence
of a weighted average error, the iterative Wiener-Hopf method with the sinc-
based fast Hilbert transform does not turn out optimal for a general solution
of the Fredholm equation, having $O(1/N)$ convergence with the number of FFT
grid points $N$ and high errors close to the function discontinuities; this
can be explained with the different requirements of the two problems, as here
the solution over the whole interval is required. Instead, the iterative
Wiener-Hopf method with the sign-based fast Hilbert transform has $O(1/N^{2})$
convergence and therefore performs better than its sinc-based sibling, the
quadrature method from the literature and the iterative method based on
Voronin’s partial solution, whose convergence is $O(1/N)$ even with the sign-
based fast Hilbert transform. So in terms of error convergence, the iterative
Wiener-Hopf method with the sign-based fast Hilbert transform reveals the new
state of the art for the numerical solution of general Fredholm equations,
achieving double the convergence speed of the known 4th order quadrature
method.

The other aspect which we must compare for the different methods is the peak
error at a discontinuity of $f(x)$ as shown in Figures 15 and 13. This error
is wider for the Wiener-Hopf method than for the quadrature method because a
wider $x$ range is required to avoid the wrap-around effects of Fourier
transform methods. Therefore, if we require an accurate answer close to a
discontinuity, the quadrature method may be best. However, the excellent CPU
time vs. error performance shown in Figures 11, 17 and 19 recommend using the
Wiener-Hopf method with the sign-based fast Hilbert transform and a larger
grid size to yield the required accuracy close to the discontinuities.

## Acknowledgments

We are thankful to Carlo Sgarra for a cue about Sec. 4.4, and to Vito Daniele
for useful discussions on the difficulties of factorising the matrices
$\mathbf{L}$ of Sec. 4.4 and $\mathbf{M}$ of Sec. 4.5, as well as for a
preprint of his book a year before it was published (Daniele & Zich,, 2014).

The support of the Economic and Social Research Council (ESRC) in funding the
Systemic Risk Centre (grant number ES/K002309/1) and of the Engineering and
Physical Sciences Research Council (EPSRC) in funding the UK Centre for
Doctoral Training in Financial Computing and Analytics (grant number 1482817)
are gratefully acknowledged.

## References

  * Bart et al., (2004) Bart, H., Gohberg, I., Kaashoek, M. A. & Ran, A. C. M. (2004) Factorization of Matrix and Operator Functions: the State Space Method. Birkhäuser, Basel. Operator Theory: Advances and Applications, Vol. 178.
  * Boyd, (2001) Boyd, J. P. (2001) Chebyshev and Fourier Spectral Methods. Springer, Heidelberg.
  * Choi et al., (2005) Choi, J., Margetis, D., Squires, T. M. & Bazant, M. Z. (2005) Steady advection-diffusion around finite absorbers in two-dimensional potential flows. Journal of Fluid Mechanics, 536, 155–184.
  * Daniele & Lombardi, (2007) Daniele, V. & Lombardi, G. (2007) Fredholm factorization of Wiener-Hopf scalar and matrix kernels. Radio Science, 42.
  * Daniele, (1984) Daniele, V. G. (1984) On the solution of two coupled Wiener-Hopf equations. SIAM Journal of Applied Mathematics, 44, 667–680.
  * Daniele & Zich, (2014) Daniele, V. G. & Zich, R. S. (2014) The Wiener-Hopf Method in Electromagnetics. SciTech Publishing, Stevenage.
  * Duffy, (2008) Duffy, D. G. (2008) Mixed Boundary Value Problems. CRC Press / Chapman & Hall (Taylor & Francis Group), Boca Raton, FL.
  * Feller, (1971) Feller, W. (1971) An Introduction to Probability Theory and Its Applications, volume 2. Wiley, New York, 3rd edition.
  * Feng & Linetsky, (2008) Feng, L. & Linetsky, V. (2008) Pricing discretely monitored barrier options and defaultable bonds in Lévy process models: a Hilbert transform approach. Mathematical Finance, 18, 337–384.
  * Feng & Linetsky, (2009) Feng, L. & Linetsky, V. (2009) Computing exponential moments of the discrete maximum of a Lévy process and lookback options. Finance and Stochastics, 13, 501–529.
  * Fredholm, (1903) Fredholm, I. (1903) Sur une classe d’équations fonctionelles. Acta Mathematica, 27, 365–390.
  * Frigo & Johnson, (2005) Frigo, M. & Johnson, S. G. (2005) The design and implementation of FFTW3. Proceedings of the IEEE, 93, 216–231.
  * Fusai et al., (2006) Fusai, G., Abrahams, I. D. & Sgarra, C. (2006) An exact analytical solution for discrete barrier options. Finance and Stochastics, 10, 1–26.
  * Fusai et al., (2016) Fusai, G., Germano, G. & Marazzina, D. (2016) Spitzer identity, Wiener-Hopf factorisation and pricing of discretely monitored exotic options. European Journal of Operational Research, 251(4), 124–134.
  * Fusai et al., (2012) Fusai, G., Marazzina, D., Marena, M. & Ng, M. (2012) $Z$-transform and preconditioning techniques for option pricing. Quantitative Finance, 12, 1381–1394.
  * Gottlieb & Shu, (1997) Gottlieb, D. & Shu, C. (1997) On the Gibbs phenomenon and its resolution. SIAM Review, 39(4), 644–668.
  * Green et al., (2010) Green, R., Fusai, G. & Abrahams, I. D. (2010) The Wiener-Hopf technique and discretely monitored path-dependent option pricing. Mathematical Finance, 20, 259–288.
  * Henery, (1974) Henery, R. J. (1974) Solution of Wiener-Hopf integral equations using the fast Fourier transform. Journal of the Institute of Mathematics and its Applications, 13, 89–96.
  * Henery, (1977) Henery, R. J. (1977) Solution of Fredholm integral equations with symmetric difference kernels. Journal of the Institute of Mathematics and its Applications, 19, 29–37.
  * Jones, (1984) Jones, D. S. (1984) Factorization of a Wiener-Hopf matrix. IMA Journal of Applied Mathematics, 32, 211–220.
  * Jones, (1991) Jones, D. S. (1991) Wiener-Hopf splitting of a $2\times 2$ matrix. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 434, 419–433.
  * King, (2009) King, F. W. (2009) Hilbert Transforms, volume 1. Cambridge University Press, Cambridge. Encyclopedia of Mathematics and its Applications, vol. 124.
  * Kisil, (2015) Kisil, A. V. (2015) Stability analysis of matrix Wiener-Hopf factorization of Daniele-Khrapkov class and reliable approximate factorization. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2178), 20150146.
  * Krein, (1962) Krein, M. G. (1962) Integral equations on a half-line with kernel depending on the difference of the arguments. In Nine Papers on Analysis, volume 22 of 2, pages 163–288. American Mathematical Society. American Mathematical Society Translations.
  * Kress & Martensen, (1970) Kress, R. & Martensen, E. (1970) Anwendung der Rechteckregel auf die reelle Hilberttransformation mit unendlichem Intervall. Zeitschrift für Angewandte Mathematik und Mechanik, 50, T61–T64.
  * Lawrie & Abrahams, (2007) Lawrie, J. B. & Abrahams, I. D. (2007) A brief historical perspective of the Wiener-Hopf technique. Journal of Engineering Mathematics, 59(4), 351–358.
  * Marazzina et al., (2012) Marazzina, D., Fusai, G. & Germano, G. (2012) Pricing credit derivatives in a Wiener-Hopf framework. In Cummins, M., Murphy, F. & Miller, J. J. H., editors, Topics in Numerical Methods for Finance, pages 139–154, New York. Springer Science+Business Media. Springer Proceedings in Mathematics & Statistics, vol. 19.
  * Margetis & Choi, (2006) Margetis, D. & Choi, J. (2006) Generalized iteration method for first-kind integral equations. Studies in Applied Mathematics, 117(1), 1–25.
  * McKechan et al., (2010) McKechan, D. J. A., Robinson, C. & Sathyaprakash, B. S. (2010) A tapering window for time-domain templates and simulated signals in the detection of gravitational waves from coalescing compact binaries. Classical and Quantum Gravity, 27(8), 084020.
  * Noble, (1958) Noble, B. (1958) Methods Based on the Wiener-Hopf Technique for the Solution of Partial Differential Equations. Pergamon Press, London. Reprinted by Chelsea, New York, 1988.
  * Pandey, (1996) Pandey, J. N. (1996) The Hilbert transform of Schwartz distributions and applications. Wiley, New York.
  * Phelan et al., (2018) Phelan, C. E., Fusai, G., Marazzina, D. & Germano, G. (2018) Fluctuation identities with continuous monitoring and their application to the pricing of barrier options. European Journal of Operational Research, 271, 210–223.
  * Phelan et al., (2019) Phelan, C. E., Marazzina, D., Fusai, G. & Germano, G. (2019) Hilbert transform, spectral filters and option pricing. Annals of Operational Research, 282(1–2), 273–298.
  * Phelan et al., (2020) Phelan, C. E., Marazzina, D. & Germano, G. (2020) Numerical pricing methods for $\alpha$-quantile and perpetual early-exercise options. Quantitative Finance, 20(6), 899–918.
  * Polyanin & Manzhirov, (1998) Polyanin, A. D. & Manzhirov, A. V. (1998) Handbook of Integral Equations. CRC Press, Boca Raton.
  * Press et al., (2007) Press, W. H., Teukolsky, S. A., Vetterling, W. T. & Flannery, B. P. (2007) Numerical Recipes in C++. Cambridge University Press, Cambridge, 3rd edition.
  * Rino, (1970) Rino, C. L. (1970) Factorization of spectra by discrete Fourier transforms. IEEE Transactions on Information Theory, 16, 484–485.
  * Rogosin & Mishuris, (2016) Rogosin, S. & Mishuris, G. (2016) Constructive methods for factorization of matrix-functions. IMA Journal of Applied Mathematics, 81(2), 365–391.
  * Ruijter et al., (2015) Ruijter, M. J., Versteegh, M. & Oosterlee, C. W. (2015) On the application of spectral filters in a Fourier option pricing technique. Journal of Computational Finance, 19(1), 75–106.
  * Spitzer, (1957) Spitzer, F. (1957) The Wiener-Hopf equation whose kernel is a probability density. Duke Mathematical Journal, 24, 327–343.
  * Stenger, (1973) Stenger, F. (1973) The approximate solution of convolution-type integral equations. SIAM Journal of Mathematical Analysis, 4, 536–555.
  * Stenger, (1993) Stenger, F. (1993) Numerical Methods Based on Sinc and Analytic Functions. Springer, Berlin.
  * Stenger, (2011) Stenger, F. (2011) Handbook of Sinc Numerical Methods. CRC Press / Chapman & Hall (Taylor & Francis Group), Boca Raton, FL.
  * Tadmor, (2007) Tadmor, E. (2007) Filters, mollifiers and the computation of the Gibbs phenomenon. Acta Numerica, 16, 305–378.
  * Tadmor & Tanner, (2005) Tadmor, E. & Tanner, J. (2005) Adaptive filters for piecewise smooth spectral data. IMA Journal of Numerical Analysis, 25(4), 635–647.
  * Vandeven, (1991) Vandeven, H. (1991) Family of spectral filters for discontinuous problems. Journal of Scientific Computing, 6(2), 159–192.
  * Vergara Caffarelli & Loreti, (1999) Vergara Caffarelli, G. & Loreti, P. (1999) Trasformata di Hilbert. Quaderni Didattici del Dipartimento di Metodi e Modelli Matematici per le Scienze Applicate, Università di Roma La Sapienza.
  * Voronin, (2004) Voronin, A. F. (2004) A complete generalization of the Wiener-Hopf method to convolution integral equations with integrable kernel on a finite interval. Differential Equations, 40(9), 1259–1267.
  * Weideman, (1995) Weideman, J. A. C. (1995) Computing the Hilbert transform on the real line. Mathematics of Computation., 64, 745–762.
  * Weisstein, (2021) Weisstein, E. W. (2021) Fourier transform — Inverse function. MathWorld, http://mathworld.wolfram.com/FourierTransformInverseFunction.html.
  * Whittaker & Watson, (1927) Whittaker, E. T. & Watson, G. N. (1927) A Course of Modern Analysis. Cambridge University Press, Cambridge, 4th edition. Reprinted 2002.
  * Wiener & Hopf, (1931) Wiener, N. & Hopf, E. (1931) Über eine Klasse singulärer Integralgleichungen. Sitzungsberichte der Preußischen Akademie der Wissenschaften, Mathematisch-Physikalische Klasse, 31, 696–706.


# Team Mia at TextVQA Challenge 2021: Vision-and-Language Representation
Learning with Pre-trained Sequence-to-Sequence Model

Yixuan Qiao1 Hao Chen1 Jun Wang1 Yihao Chen2 Xianbin Ye3 Ziliang Li4
Xianbiao Qi2 Peng Gao1 Guotong Xie1,5,6
1 SFE Deeplearning Platform  Ping An Health Technology  Beijing  China.
2 Visual Computing Group  Ping An Property & Casualty Insurance Company
Shenzhen  China.
3 Jinan University  Guangzhou  China. 4 Central University of Finance and
Economics  Beijing  China.
5 Ping An Health Cloud Company Limited  Shenzhen  China.
6 Ping An International Smart City Technology Co  Ltd  Shenzhen  China.
{qiaoyixuan528, chenhao305, wangjun916, chenyihao291,
<EMAIL_ADDRESS>
{gaopeng712<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

TextVQA requires models to read and reason about text in images to answer
questions about them. Specifically, models need to incorporate a new modality
of text present in the images and reason over it to answer TextVQA questions.
In this challenge, we use generative model T5 for TextVQA task. Based on pre-
trained checkpoint T5-3B from HuggingFace repository, two other pre-training
tasks including masked language modeling(MLM) and relative position
prediction(RPP) are designed to better align object feature and scene text. In
the stage of pre-training, encoder is dedicate to handle the fusion among
multiple modalities: question text, object text labels, scene text labels,
object visual features, scene visual features. After that decoder generates
the text sequence step-by-step, cross entropy loss is required by default. We
use a large-scale scene text dataset in pre-training and then fine-tune the
T5-3B with the TextVQA dataset only.

Figure 1: Comparison of our solution and other teams from TextVQA 2021 and
TextVQA 2020.

##  1 Pre-training Stage

###  1.1 Pre-training Dataset

Inspired by the promising results of TAP (Text-Aware Pre-training) [2],
further increasing the size of model and data are expected to have better
performance. We built a dataset contains about 9 million scene text images.
Part of it by using the similar methods as TAP, part of it by generating
synthetic text images as described in [1]. Meantime, we have corresponding
relationships between object features and scene features for both parts.

###  1.2 Pre-training Strategy

First of all, for each scene image, we have the following features.

  * •

image texts, object labels, scene texts

  * •

object visual features, scene visual features

Object labels and features are extracted from VinVL(Revisiting Visual
Representations in Vision-Language Models) [3], which can generate
representations of a richer collection of visual objects and concepts. Scene
text and scene visual features are mainly coming from public Microsoft OCR
API.

For each scene visual features, scene text pair in images, we use the former
with all other features as input, scene text as target to construct training
samples. This design comes from the fact that half of the answers in the Text-
VQA dataset came from scene text which is necessary to equip the model with
the ability to generate scene text. For MLM on the extended text input, we
randomly mask each text token with a probability of 15%, the details are the
same as BERT. The objective of RPP is to predict the relative position between
the scene visual feature and each object visual feature.

Additionally, to enable large-scale training, we adopt the adversarial
training strategy combined with KL-divergence-based regularization to promote
higher invariance in the embedding space. Specifically, we propose to perform
adversarial training in the embedding space of each modality, rather than
adding adversarial perturbations on image pixels and textual tokens. Finally,
together with the cross-entropy loss of decoder as final objective to be
optimized.

##  2 Fine-tuning stage

In fine-tuning, the model step-by-step predicts the answer using its decoding
module, and is trained with the answer classification loss and the adversarial
training strategy above in each step.

Metric | TextVQA
---|---
Number of images | 28k
Images with text | 100%
Number of questions | 45k
Number of questions with answer in OCR tokens | 39% (18k)
Number of questions with spatial words | 14%

Table 1: Details of the TextVQA dataset in Fine-tuning stage.

##  3 Post-processing Module

After we got the answer sequence through fine-tuning stage, we use a open
source fuzzy string matching method fuzzywuzzy which uses Levenshtein Distance
to calculate the differences between sequences to make the final correction.

##  4 Results and Conclusions

TextVQA requires models to read and reason about text in an image to answer
questions based on them. In order to perform well on this task, models need to
first detect and read text in the images. Models then need to reason about
this to answer the question. There are 115 and 65 submissions in validation
and test phases from the 33 teams participating in TextVQA 2021 challenge. The
final comparison between our solution and other results from TextVQA
challenges is shown in Figure 1

To summarize, the experiments from this challenge present the following key
suggestions:

  * •

Pretraining with image-text matching contrastive objectives on huge data leads
to significant gains.

  * •

Developing hard-negatives for such objectives also improve performance.

  * •

There is a scope of improvement in multi-word answers.

  * •

Using a pretrained NLP model as a base helps.

  * •

Improved object and scene-text features also lead to some gains.

## References

  * [1] Ankush Gupta, Andrea Vedaldi, and Andrew Zisserman.  Synthetic data for text localisation in natural images.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2315–2324, 2016.
  * [2] Zhengyuan Yang, Yijuan Lu, Jianfeng Wang, Xi Yin, Dinei Florencio, Lijuan Wang, Cha Zhang, Lei Zhang, and Jiebo Luo.  Tap: Text-aware pre-training for text-vqa and text-caption, 2020.
  * [3] Pengchuan Zhang, Xiujun Li, Xiaowei Hu, Jianwei Yang, Lei Zhang, Lijuan Wang, Yejin Choi, and Jianfeng Gao.  Vinvl: Revisiting visual representations in vision-language models, 2021\.