# Test Problems for Regularization Methods¶

A Fredholm integral equation of the first kind (in 1-dimensional) can be written as

where \(g\) and \(K\) (called kernel) are known functions and \(f\) is the unknown solution. This is a classical example of a linear ill-posed problem, i.e., an arbitrary small perturbation of the data can cause an arbitrarily large perturbation of the solution. For example, in computerized tomography, \(K\) is an X-ray source, \(f\) is the object being scanned, and \(g\) is the measured damping of the X-rays. The goal here is to reconstruct the scanned object from information about the locations of the X-ray sources and measurements of their damping.

After discretizations (by the quadrature method or the Galerkin
method), we obtain a linear system of equations \(Ax=b\).
All the regularization test problems are derived from
discretizations of a Fredholm integral equation of the first kind.
Each generated test problem has type `RegProb`

, which is
defined as:

```
immutable RegProb{T}
A::AbstractMatrix{T} # matrix of interest
b::AbstractVector{T} # right-hand side
x::AbstractVector{T} # the solution to Ax = b
end
```

Here is an example:

```
julia> matrixdepot("deriv2")
Computation of the Second Derivative:
A classical test problem for regularization algorithms.
Input options:
1. [type,] n, [matrixonly]: the dimension of the matrix is n.
If matrixonly = false, the linear system A, b, x will be generated.
(matrixonly = true by default.)
Reference: P.C. Hansen, Regularization tools: A MATLAB package for
analysis and solution of discrete ill-posed problems.
Numerical Algorithms, 6(1994), pp.1-35
julia> A = matrixdepot("deriv2", 4) # generate the test matrix
4x4 Array{Float64,2}:
-0.0169271 -0.0195313 -0.0117188 -0.00390625
-0.0195313 -0.0481771 -0.0351563 -0.0117188
-0.0117188 -0.0351563 -0.0481771 -0.0195313
-0.00390625 -0.0117188 -0.0195313 -0.0169271
julia> A = matrixdepot("deriv2", Float32, 4)
4x4 Array{Float32,2}:
-0.0169271 -0.0195313 -0.0117188 -0.00390625
-0.0195313 -0.0481771 -0.0351563 -0.0117188
-0.0117188 -0.0351563 -0.0481771 -0.0195313
-0.00390625 -0.0117188 -0.0195313 -0.0169271
julia> r = matrixdepot("deriv2", 3, false) # generate the test problem
Test problems for Regularization Method
A:
3x3 Array{Float64,2}:
-0.0277778 -0.0277778 -0.00925926
-0.0277778 -0.0648148 -0.0277778
-0.00925926 -0.0277778 -0.0277778
b:
[-0.01514653483985129,-0.03474793286789414,-0.022274315940957783]
x:
[0.09622504486493762,0.28867513459481287,0.48112522432468807]
julia> r.A # generate the matrix A
3x3 Array{Float64,2}:
-0.0277778 -0.0277778 -0.00925926
-0.0277778 -0.0648148 -0.0277778
-0.00925926 -0.0277778 -0.0277778
```

Here is a list of test problems in the collection:

- baart
Discretization of an artificial Fredholm integral equation of the first kind [baart82]. The kernel \(K\) is given by

\[K(s,t) = \exp(s \cos (t)).\]The right-hand side \(g\) and the solution \(f\) are given by

\[g(s)=2\frac{\sin (s)}{s}, \quad f(t) = \sin(t).\][baart82] M.L. Baart, The use of auto-correlation for pseudo-rank determination in noisy ill-conditioned linear least-squares problems, IMA, J. Numer. Anal. 2 (1982), 241-247. - blur
Image deblurring test problem. It arises in connection with the degradation of digital images by atmospheric turbulence blur, modelled by a Gaussian point-spread function

\[h(x,y) = \frac{1}{2\pi\sigma^2}\exp(-\frac{x^2+y^2}{2\sigma^2}).\]The matrix A is a symmetric \(n^2\times n^2\) doubly block Toeplitz matrix, stored in sparse format.

- deriv2
Computation of the second derivative. The kernel \(K\) is Green’s function for the second derivative

\[\begin{split}K(s,t) = \begin{cases} s(t - 1), \quad s < t, \\ t(s - 1), \quad s \geq t, \\ \end{cases}\end{split}\]and both integration intervals are \([0,1]\). The function \(g\) and \(f\) are given by

\[g(s) = (s^3 - s)/6, \quad f(t) = t.\]The symmetric matrix \(A\) and vectors \(x\) and \(b\) are computed from \(K,f\) and \(g\) using the Galerkin method.

- foxgood
A severely ill-posed problem suggested by Fox & Goodwin. This is a model problem which does not satisfy the discrete Picard condition for the small singular values [baker77].

[baker77] C.T.H Baker, The Numerical Treatment of Integral Equations, Clarendon Press, Oxford, 1977, p. 665. - gravity
One-dimensional gravity surveying model problem. Discretization of a 1-D model problem in gravity surveying, in which a mass distribution f(t) is located at depth d, while the vertical component of the gravity field g(s) is measured at the surface. The resulting problem is a first-kind Fredholm integral equation with kernel

\[K(s,t) = d(d^2 + (s-t)^2)^{-3/2}.\]- heat
Inverse heat equation [carasso82]. It is a Volterra integral equation of the first kind with integration interval \([0,1]\). The kernel \(K\) is given by

\[K(s,t) = k(s-t),\]where

\[k(t) = \frac{t^{-3/2}}{2\kappa \sqrt{\pi}}\exp\big(-\frac{1}{4\kappa^2t}\big).\]\(\kappa\) controls the ill-conditioning of the matrix \(A\). \(\kappa = 1\) (default) gives an ill-conditioned matrix and \(\kappa = 5\) gives a well-conditioned matrix.

[carasso82] A.S. Carasso, Determining surface temperatures from interior observations, SIAM J. Appl. Math. 42 (1982), 558-574. - parallax
Stellar parallax problem with 26 fixed, real observations. The underlying problem is a Fredholm integral equation of the first kind with kernel

\[K(s,t) = \frac{1}{\sigma \sqrt{2\pi}}\exp\Big(-\frac{1}{2}\big(\frac{s-t}{\sigma}\big)^2\Big),\]with \(\sigma = 0.014234\) and it is discretized by means of a Galerkin method with n orthonormal basis functions. The right-hand side b consists of a measured distribution function of stellar parallaxes, and its length is fixed at 26; i.e., the matrix \(A\) is \(26\times n\). The exact solution, which represents the true distribution of stellar parallaxes, is unknown.

- phillips
Phillips’s “famous” problem. Discretization of the “famous” Fredholm integral equation of the first kind devised by D.L. Phillips [phillips62]. The kernel \(K\) and solution \(f\) are given by

\[K(s,t) = \theta(s-t), \quad f(t) = \theta(t),\]where

\[\begin{split}\theta(x) = \begin{cases} 1+\cos(\frac{\pi x}{3}), & |x| < 3, \\ 0, & |x| \geq 3. \\ \end{cases}\end{split}\]The right-hand side \(g\) is given by

\[g(s) = (6 - |s|)\Big( 1 + \frac{1}{2}\cos\big(\frac{\pi s}{3}\big)\Big) + \frac{9}{2 \pi}\sin\Big(\frac{\pi |s|}{3}\Big).\]Both integration intervals are \([-6,6]\).

[phillips62] D.L. Phillips, A technique for the numerical solution of certain integral equations of the first kind, J. ACM 9 (1962), 84-97. - shaw
One-dimensional image restoration model. This test problem uses a first-kind Fredholm integral equation to model a one-dimensional image restoration situation. The kernel \(K\) is given by

\[K(s,t) = (\cos(s)+\cos(t))^2\big(\frac{\sin(u)}{u}\big)^2,\]where

\[u = \pi(\sin(s) + \sin(t)).\]Both integration intervals are \([-\pi/2, \pi/2]\). The solution \(f\) is given by

\[f(t) = a_1 \exp(-c_1(t-t_1)^2) + a_2 \exp(-c_2(t-t_2)^2).\]\(K\) and \(f\) are discretized by simple quadrature to produce the matrix \(A\) and the solution vector \(x\). The right-hand \(b\) is computed by \(b=Ax\).

- spikes
- Artificially generated discrete ill-posed problem.
- ursell
Discretization of a Fredholm integral equation of the first kind with kernel K and right-hand side g given by

\[K(s,t) = \frac{1}{s+t+1}, \quad g(s) = 1,\]where both integration intervals are \([0,1]\) [ursell].

[ursell] F. Ursell, Introduction to the theory of linear integral equations, Chapter 1 in L. M. Delves and J. Walsh, Numerical Solution of Integral Equations, Clarendon Press, Oxford, 1974. - wing
A problem with a discontinuous solution. The kernel \(K\) is given by

\[K(s,t) = t \exp(-st^2),\]with both integration intervals are \([0,1]\). The functions \(f\) and \(g\) are given as

\[\begin{split}f(t) = \begin{cases} 1, \quad t_1 < t < t_2, \\ 0, \quad \mbox{otherwise},\\ \end{cases} \quad g(s) = \frac{\exp(-st_1^2) - \exp(-st_2^2)}{2s}.\end{split}\]Here \(0 < t_1 < t_2 < 1\). The matrix \(A\) and two vectors \(x\) and \(b\) are obtained by Galerkin discretization with orthonormal basis functions defined on a uniform mesh.