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 pacakge 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
• blur
• deriv2
• foxgood
• gravity
• heat
• phillips
• shaw
• spikes
• ursell
• wing

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 pesudo-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.

phillips

Phillips’s “famous” problem. Discretization of the “famous” Fredholm integral equation of the first kind deviced 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.