Test Problems for Regularization Methods

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

\[\int_{0}^1 K(s,t) f(t) dt = g(s), \quad 0 \leq s \leq 1,\]

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.