# Adaptive cross approximation for Tikhonov regularization in general form

T. Mach Institute for Mathematics, University of Potsdam, 14476 Potsdam,
Germany. E-mail<EMAIL_ADDRESS>L. Reichel Department of Mathematical
Sciences, Kent State University, Kent, OH 44242, USA. E-mail:
<EMAIL_ADDRESS>M. Van Barel Department of Computer Science, KU Leuven,
Celestijnenlaan 200A, 3001 Leuven (Heverlee), Belgium. E-mail:
<EMAIL_ADDRESS>
###### Abstract

Many problems in Science and Engineering give rise to linear integral
equations of the first kind with a smooth kernel. Discretization of the
integral operator yields a matrix, whose singular values cluster at the
origin. We describe the approximation of such matrices by adaptive cross
approximation, which avoids forming the entire matrix. The choice of the
number of steps of adaptive cross approximation is discussed. The discretized
right-hand side represents data that commonly are contaminated by measurement
error. Solution of the linear system of equations so obtained is not
meaningful because the matrix determined by adaptive cross approximation is
rank-deficient. We remedy this difficulty by using Tikhonov regularization and
discuss how a fairly general regularization matrix can be used. Computed
examples illustrate that the use of a regularization matrix different from the
identity can improve the quality of the computed approximate solutions
significantly.

###### keywords:

ill-posed problem, inverse problem, sparse discretization, regularization,
adaptive cross approximation

Dedicated to Claude Brezinski on the occasion of his 80th birthday.

##  1 Introduction

Linear integral equations of the first kind,

$\int_{\Omega_{1}}\kappa(s,t)x(t)dt=g(s),\qquad s\in\Omega_{2},$ (1)

with a smooth kernel $\kappa$ arise in many applications, including remote
sensing, computerized tomography, and image restoration. Here $\Omega_{i}$
denotes a subset of ${\mathbb{R}}^{d_{i}}$ for some positive integer $d_{i}$.
The solution of (1) is an _ill-posed problem_ , because the singular values of
the integral operator cluster at the origin; see, e.g., [13, 22] for
introductions to ill-posed problems.

Discretization of (1) yields a linear system of equations

$A\mbox{\boldmath{$x$}}=\mbox{\boldmath{$g$}},\qquad
A\in{{\mathbb{R}}}^{n\times
n},\qquad\mbox{\boldmath{$g$}}\in{{\mathbb{R}}}^{n},$ (2)

with a matrix, whose singular values coalesce at the origin. This makes the
matrix $A$ severely ill-conditioned and possibly rank-deficient; we measure
the conditioning of a matrix with its condition number, which is the ratio of
the largest and smallest singular values. Linear systems of equations with a
matrix of this kind are often referred to as _discrete ill-posed problems_ ;
see, e.g., [23]. We will for notational simplicity assume the matrix $A$ to be
square, however, the method described also can be applied, after minor
modifications, when $A$ is rectangular, in which case the linear system of
equations (2) is replaced by a least-squares problem.

In many applications, the right-hand side vector $g$ in (2) represents
measured data and is contaminated by a measurement error $e$. Due to the
severe ill-conditioning of $A$, straightforward solution of (2) typically
yields a computed solution that is severely contaminated by propagated error,
and therefore is not useful. To circumvent this difficulty, the linear system
of equations (2) commonly is replaced by a nearby problem, whose solution is
less sensitive to the error $e$ in $g$. This replacement is referred to as
_regularization_. Tikhonov regularization is possibly the most popular and
well understood regularization method. It replaces the linear system of
equations (2) by a penalized least-squares problem of the form

$\min_{\mbox{\boldmath{${\scriptstyle
x}$}}\in{\mathbb{R}}^{n}}\\{\|A\mbox{\boldmath{$x$}}-\mbox{\boldmath{$g$}}\|^{2}+\mu\|L\mbox{\boldmath{$x$}}\|^{2}\\},$
(3)

where $L\in{\mathbb{R}}^{p\times n}$ is referred to as the regularization
matrix and $\mu>0$ as the regularization parameter. The problem (3) is said to
be in _standard form_ when $L$ is the identity; otherwise (3) is in _general
form_. Throughout this paper $\|\cdot\|$ denotes the Euclidean vector norm or
the spectral matrix norm.

The choice of regularization parameter $\mu$ is important for the quality of
the computed solution: a too small value results in a computed solution that
is contaminated by needlessly much propagated error, while a too large value
yields an unnecessarily smooth solution that may lack details of interest.
Generally, a suitable value of $\mu$ is not known a priori, but has to be
computed during the solution process. This typically requires that (3) be
solved for several $\mu$-values. Methods for determining a suitable value of
$\mu$ include the L-curve criterion, generalized cross validation, and the
discrepancy principle; see, e.g., [6, 7, 15, 28, 32, 33] for discussions of
properties of these and other methods.

The matrix $L$ is assumed to be chosen so that

${\mathcal{N}}(A)\cap{\mathcal{N}}(L)=\\{\mbox{\boldmath{$0$}}\\},$ (4)

where ${\mathcal{N}}(M)$ denotes the null space of the matrix $M$. Then the
Tikhonov minimization problem (3) has the unique solution

$\mbox{\boldmath{$x$}}_{\mu}:=(A^{T}A+\mu
L^{T}L)^{-1}A^{T}\mbox{\boldmath{$g$}}$ (5)

for any $\mu>0$; see, e.g, [23] for details. Here and below the superscript T
denotes transposition. We are interested in the situation when the matrices
$A$ and $L$ are so large that it is impossible or undesirable to compute the
solution (5) by Cholesky factorization of the matrix $A^{T}A+\mu L^{T}L$. In
fact, we would like to avoid evaluating all the entries of $A$. We describe
how $A$ can be approximated by a much smaller matrix, without evaluating all
matrix entries, by applying adaptive cross approximation.

The application of cross approximation to matrices that stem from the
discretization of Fredholm integral equations of the second kind has received
considerable attention in the literature, see, e.g., [3, 4, 17, 20, 37];
however, the use of cross approximation in the context of solving linear
discrete ill-posed problems has not been thoroughly studied.

The use of adaptive approximation to the approximate solution of (3) when $L$
is the identity is discussed in [29]. This paper extends this discussion to
general regularization matrices $L$. Our interest in this extension of the
method in [29] stems from the fact that the use of a suitably chosen
regularization matrix $L$ can deliver solutions of higher quality than $L=I$;
see, e.g., [25, 31, 35] for illustrations and discussions on various ways of
constructing regularization matrices. Roughly, $L$ should be chosen so as not
to damp known important features of the desired solution, while damping the
propagated error stemming from the error in $g$.

In the computed examples of Section 4, we use the _discrepancy principle_ to
determine $\mu>0$. Let $\widehat{\mbox{\boldmath{$g$}}}\in{\mathbb{R}}^{n}$
denote the unknown error-free vector associated with the right-hand side $g$
in (2), i.e.,
$\mbox{\boldmath{$g$}}=\widehat{\mbox{\boldmath{$g$}}}+\mbox{\boldmath{$e$}}$.
Assume that the linear system of equations with the error-free right-hand
side,

$A\mbox{\boldmath{$x$}}=\widehat{\mbox{\boldmath{$g$}}},$ (6)

is consistent and that a fairly accurate bound
$\|\mbox{\boldmath{$e$}}\|\leq\delta$ is known. The discrepancy principle
prescribes that the regularization parameter $\mu>0$ be determined so that the
Tikhonov solution (5) satisfies

$\|A\mbox{\boldmath{$x$}}_{\mu}-\mbox{\boldmath{$g$}}\|=\eta\delta,$ (7)

where $\eta>1$ is a user-specified parameter that is independent of $\delta$.
It can be shown that when $\delta$ tends to zero,
$\mbox{\boldmath{$x$}}_{\mu}$ converges to the minimal-norm solution,
$\widehat{\mbox{\boldmath{$x$}}}$, of (6); see, e.g., [13] for a proof in a
Hilbert space setting. We remark that the determination of $\mu>0$ such that
$\mbox{\boldmath{$x$}}_{\mu}$ satisfies (7) typically requires the solution of
(3) for several $\mu$-values.

The present paper is concerned with the situation when the matrix
$A\in{\mathbb{R}}^{n\times n}$ in (2) is large. Then the evaluation of all
entries of $A$ can be quite time-consuming. Cross approximation, also known as
skeleton approximation, of $A$ reduces this time by approximating $A$ by a
matrix $M_{k}\in{\mathbb{R}}^{n\times n}$ that consists of $k\ll n$ rows and
columns of $A$. We would like to choose the rows and columns of $A$ so that
$M_{k}$ approximates $A$ well and is easy to compute with.

This paper is organized as follows. Section 2 reviews the application of
adaptive cross approximation to the approximation of $A$ by a matrix of low
rank. In Section 3, we describe the application of adaptive cross
approximation to the approximation of the Tikhonov equation (3). Section 4
reports a few computed examples, and concluding remarks can be found in
Section 5.

We conclude this section with an example that leads to a large linear discrete
ill-posed problem, whose solution is difficult to compute using
straightforward discretization.

Example 1.1. Consider the solution of the Fredholm integral equation of the
first kind

$\int_{S}\frac{\sigma(\mbox{\boldmath{$y$}})}{4\pi\epsilon_{0}\mbox{$\parallel\\!\mbox{\boldmath{$x$}}-\mbox{\boldmath{$y$}}\\!\parallel$}}\,d\mbox{\boldmath{$y$}}=\phi(\mbox{\boldmath{$x$}}),\qquad\mbox{\boldmath{$x$}}\in
S,$ (8)

where $\phi$ is a given electric potential, $S$ is a surface with electrodes,
$\sigma(\mbox{\boldmath{$y$}})$ denotes the density of the charge on $S$, and
$\epsilon_{0}$ stands for the electric permittivity in vacuum. We would like
to determine $\sigma$ from $\phi$, and assume that $\phi$ is chosen so that
(8) has a solution. The computation of a solution of (8) is an ill-posed
problem. Using a weak formulation and discretization lead to a dense symmetric
matrix $K=[k_{ij}]_{i,j=1}^{n}\in{\mathbb{R}}^{n\times n}$ with entries

$k_{ij}=\int_{S}\int_{S}\frac{v_{j}(\mbox{\boldmath{$x$}})v_{i}(\mbox{\boldmath{$y$}})}{4\pi\epsilon_{0}\mbox{$\parallel\\!\mbox{\boldmath{$x$}}-\mbox{\boldmath{$y$}}\\!\parallel$}}\,d\mbox{\boldmath{$y$}}\,d\mbox{\boldmath{$x$}},\quad
i,j=1,2,\ldots,n.$

This matrix can be expensive to store and handle when the discretization is
fine. Employing a hierarchical compression with $\mathcal{H}^{2}$-matrices
reduces the required storage to $O(n)$, with a large constant hidden in the
$O(\cdot)$, and allows matrix-vector product evaluations in $O(n)$ arithmetic
floating point operations (flops); see [5].

A fine discretization with $n=262,146$ nodes results in a large, $n\times n$,
dense matrix. We used the H2Lib library [27] for the computations and base
this example on one of the standard examples provided in this library. Without
compression, 512 GB of memory are needed to store the matrix. On a laptop
computer with an Intel Core i710710U CPU and 16 GB of RAM it took 1103 s to
assemble the matrix $K$ in the compressed $\mathcal{H}^{2}$-matrix format. The
matrix required 15.45 GB of storage, thus almost all the available RAM.
Carrying out one matrix-vector product evaluation required 1596 s, that is 44%
more time than for assembling the matrix. The reason for this is that the
$O(n)$ flops require a significant amount of communication between faster and
slower storage. This example illustrates that there are linear discrete ill-
posed problems of interest that are difficult to solve on a laptop computer,
even if a significant amount of memory is available. It therefore is important
to develop methods that are able to determine approximations of dense matrices
that requires less computer storage and less CPU time for the evaluation of
matrix-vector products. $\Box$

##  2 Adaptive cross approximation

Cross approximation of a large matrix $A\in{\mathbb{R}}^{n\times n}$
determines an approximation $M_{k}\in{\mathbb{R}}^{n\times n}$ of rank at most
$k\ll n$. All entries of $M_{k}$ can be evaluated much faster than all entries
of $A$, because $M_{k}$ is constructed from only $k$ rows and columns of $A$.
We would like $M_{k}$ to be an accurate approximation of $A$. This is achieved
by a careful choice of the $k$ rows and columns of $A$ that define the matrix
$M_{k}$. A cross approximation method is said to be _adaptive_ when the rows
and columns of $A$ that determine $M_{k}$ (and $k$) are chosen depending on
properties of $A$ revealed during the computations; see [2, 4].

We outline the adaptive cross approximation method for a general square
nonsymmetric matrix described in [29]. This method is an adaptation of the
scheme in [17] to the approximation of the matrix of linear discrete ill-posed
problems. When $A$ is symmetric, the matrix $M_{k}$ can be chosen to be
symmetric. This roughly halves the storage requirement for $M_{k}$. Both the
situations when $A$ is symmetric positive definite or symmetric indefinite are
discussed in [29]. We therefore will not dwell on these special cases in the
present paper.

Let the matrix $A\in{\mathbb{R}}^{n\times n}$ be nonsymmetric and choose $k$
rows of $A$ with indices in ${\mathbb{N}}=\\{1,2,\dots,n\\}$. We let
$\mathbf{i},\mathbf{j}\in{\mathbb{N}}^{k}$ denote index vectors with $k$
entries in ${\mathbb{N}}$. The submatrices $A_{(\mathbf{i},:)}$ and
$A_{(:,\mathbf{j})}$ of $A$ are made up of the $k$ rows with indices
$\mathbf{i}$ and the $k$ columns with indices $\mathbf{j}$, respectively.
Moreover, the _core matrix_ $A_{(\mathbf{i},\mathbf{j})}$ is made up of $k$
rows and columns of $A$. Assume that this matrix is nonsingular. Then the rows
and columns of the matrix

$M_{k}=A_{(:,\mathbf{j})}A_{(\mathbf{i},\mathbf{j})}^{-1}A_{(\mathbf{i},:)}$

are equal to the corresponding rows and columns of $A$; when $A$ is of rank
$k$, we have $M_{k}=A$.

Goreinov et al. [21] show that it is beneficial to choose the index vectors
$\mathbf{i}$ and $\mathbf{j}$ so that $A_{(\mathbf{i},\mathbf{j})}$ is a
submatrix of $A$ of maximal volume, i.e., so that the modulus of the
determinant of $A_{(\mathbf{i},\mathbf{j})}$ is maximal. However, it is
difficult to determine such index vectors $\mathbf{i}$ and $\mathbf{j}$. We
therefore seek to determine a low-rank matrix $M_{k}$ that is a sufficiently
accurate approximation of $A$ by a greedy algorithm. Suppose that we already
have computed an approximation

$M_{k-1}=\sum_{\ell=1}^{k-1}\mbox{\boldmath{$w$}}^{(c)}_{\ell}(\mbox{\boldmath{$w$}}^{(r)}_{\ell})^{T},\qquad\mbox{\boldmath{$w$}}^{(c)}_{\ell},\mbox{\boldmath{$w$}}^{(r)}_{\ell}\in{\mathbb{R}}^{n},$

of rank at most $k-1$ of $A$. To compute the next approximation, $M_{k}$, of
$A$ of rank at most $k$, we determine a row index $i^{*}$ and a column index
$j^{*}$ by looking for the index of the maximum element in magnitude in the
previously computed vectors $\mbox{\boldmath{$w$}}^{(c)}_{k-1}$ (for index
$i^{*}$) and $\mbox{\boldmath{$w$}}^{(r)}_{k-1}$ (for index $j^{*}$). The
vector $(\mbox{\boldmath{$w$}}^{(r)}_{1})^{T}$ can be chosen as an arbitrary
row of $A$. We will let $(\mbox{\boldmath{$w$}}^{(r)}_{1})^{T}$ be the first
row of $A$ in the computed examples of Section 4. The vector
$\mbox{\boldmath{$w$}}^{(c)}_{1}$ can be chosen in a similar way.

In the simplest form of cross approximation, the determination of the vectors
$\mbox{\boldmath{$w$}}^{(c)}_{k}$ and $\mbox{\boldmath{$w$}}^{(r)}_{k}$ only
requires the entries in row $i^{*}$ and column $j^{*}$ of $A$ and the elements
of already computed vectors $\mbox{\boldmath{$w$}}^{(c)}_{\ell}$ and
$(\mbox{\boldmath{$w$}}^{(r)}_{\ell})^{T}$, $\ell=1,2,\ldots,k-1$:

$\displaystyle(\mbox{\boldmath{$w$}}^{(r)}_{k})_{j}$ $\displaystyle=$
$\displaystyle
A_{i_{k}^{*},j}-\sum_{\ell=1}^{k-1}(\mbox{\boldmath{$w$}}^{(c)}_{\ell})_{i_{k}^{*}}(\mbox{\boldmath{$w$}}^{(r)}_{\ell})_{j},\hskip
14.22636pt\delta=(\mbox{\boldmath{$w$}}^{(r)}_{k})_{j^{*}_{k}},$
$\displaystyle(\mbox{\boldmath{$w$}}^{(c)}_{k})_{i}$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{\delta}\Big{(}A_{i,j_{k}^{*}}-\sum_{\ell=1}^{k-1}(\mbox{\boldmath{$w$}}^{(c)}_{\ell})_{i}(\mbox{\boldmath{$w$}}^{(r)}_{\ell})_{j_{k}^{*}}\Big{)}.$

A new skeleton is obtained from the remainder,

$R_{k}=A-\sum_{\ell=1}^{k}\mbox{\boldmath{$w$}}^{(c)}_{\ell}(\mbox{\boldmath{$w$}}^{(r)}_{\ell})^{T},$

without explicitly computing all entries of the matrix $R_{k}$.

The required number of rank-one matrices, $k$, that make up $M_{k}$ is
generally not known a priori. We would like the difference $A-M_{k}$ to be of
small norm. However, we cannot evaluate this difference, because most entries
of $A$ are not known. Following [17], we include $t$ randomly chosen matrix
entries $A_{i_{\ell},j_{\ell}}$, for $\ell=1,2,\dots,t$, with
$i_{\ell},j_{\ell}\in{\mathbb{N}}$. Define for future reference the set

$\displaystyle\Pi=\\{(i_{\ell},j_{\ell}),~{}\forall~{}\ell=1,2,\dots,t\\}$ (9)

When a new skeleton is determined, the values of these entries are updated by
Subtraction from the available skeletons,

$(R_{k})_{i_{\ell},j_{\ell}}=(R_{k-1})_{i_{\ell},j_{\ell}}-(\mbox{\boldmath{$w$}}^{(c)}_{k})_{i_{\ell}}(\mbox{\boldmath{$w$}}^{(r)}_{k})_{j_{\ell}},$
(10)

with $(R_{0})_{i_{\ell},j_{\ell}}=A_{i_{\ell},j_{\ell}}$. The values
$(R_{k})_{i_{\ell},j_{\ell}}$ are used in subsection 3.3 as part of the
stopping criterion to determine the final value $k^{*}$ for $k$. The value of
$t$ is a percentage of the total number of entries. The choice of $t$ should
depend on properties of the matrix $A$; see [17, 29] for further details. An
algorithm is presented in [29].

##  3 Tikhonov regularization in general form

This section discusses how to combine adaptive cross approximation with
Tikhonov regularization in general form.

###  3.1 Using adaptive cross approximation

The matrix $M_{k}$, whose computation was outlined in the previous section, is
of the form

$M_{k}=W_{k}^{(c)}(W_{k}^{(r)})^{T},$ (11)

where

$W_{k}^{(c)}=[\mbox{\boldmath{$w$}}_{1}^{(c)},\mbox{\boldmath{$w$}}_{2}^{(c)},\ldots,\mbox{\boldmath{$w$}}_{k}^{(c)}]\in{\mathbb{R}}^{n\times
k},\qquad
W_{k}^{(r)}=[\mbox{\boldmath{$w$}}_{1}^{(r)},\mbox{\boldmath{$w$}}_{2}^{(r)},\ldots,\mbox{\boldmath{$w$}}_{k}^{(r)}]\in{\mathbb{R}}^{n\times
k}.$

Compute the skinny QR factorizations

$W_{k}^{(c)}=Q_{k}^{(c)}R_{k}^{(c)},\qquad
W_{k}^{(r)}=Q_{k}^{(r)}R_{k}^{(r)},$ (12)

where the matrices $Q_{k}^{(c)},Q_{k}^{(r)}\in{\mathbb{R}}^{n\times k}$ have
orthonormal columns and the matrices
$R_{k}^{(c)},R_{k}^{(r)}\in{\mathbb{R}}^{k\times k}$ are upper triangular. The
factorizations (12) can be computed by the Householder-QR method or by
factorization methods that are designed to perform efficiently on modern
computers, such as the methods described in [9, 14, 38].

Combining (11) and (12) yields

$M_{k}=Q_{k}^{(c)}R_{k}^{(c)}(R_{k}^{(r)})^{T}(Q_{k}^{(r)})^{T}.$ (13)

Replacing $A$ by $M_{k}$ in (3) gives the minimization problem

$\min_{\mbox{\boldmath{${\scriptstyle
x}$}}\in{\mathbb{R}}^{n}}\\{\|M_{k}\mbox{\boldmath{$x$}}-\mbox{\boldmath{$g$}}\|^{2}+\mu\|L\mbox{\boldmath{$x$}}\|^{2}\\},$
(14)

which can be solved in several ways. If the matrix $L$ has a special
structure, such as being banded with small bandwidth, then it may be
attractive to transform (14) to standard form by a technique described by
Eldén [12]. Regularization matrices $L$ with a small bandwidth arise, e.g.,
when $L$ represents a finite difference approximation of a differential
operator in one space-dimension. It also is easy to transform (14) to standard
form when $L$ is an orthogonal projector; see [30].

In the remainder of this section, we discuss the situation when $L$ is such
that transformation of (14) to standard form as described in [12] is too
expensive to be attractive. This is the case, for instance, when $L$
represents a finite difference approximation of a differential operator in two
or more space-dimensions. This kind of matrices $L$ will be used in computed
examples of Section 4.

We describe an approaches to compute a solution of (14) and start with the
simplest one. The matrix $M_{k}$ has a null space of large dimension (at least
$n-k$). Therefore the Tikhonov minimization problem (14) is not guaranteed to
have a unique solution. To remedy this difficulty, we require the solution of
(14) to live in a subspace of fairly low dimension. A simple solution method
is obtained when using the solution subspace ${\mathcal{R}}(Q_{k}^{(r)})$.
Then we obtain the minimization problem

$\min_{\mbox{\boldmath{${\scriptstyle
y}$}}\in{\mathbb{R}}^{k}}\\{\|M_{k}Q_{k}^{(r)}\mbox{\boldmath{$y$}}-\mbox{\boldmath{$g$}}\|^{2}+\mu\|LQ_{k}^{(r)}\mbox{\boldmath{$y$}}\|^{2}\\},$
(15)

which has a unique solution if and only if

${\mathcal{N}}(M_{k}Q_{k}^{(r)})\cap{\mathcal{N}}(LQ_{k}^{(r)})=\\{\mbox{\boldmath{$0$}}\\}.$
(16)

This holds, in particular, when the triangular matrices $R_{k}^{(c)}$ and
$R_{k}^{(r)}$ in (13) are nonsingular. We found (16) to hold in all computed
examples that we solved.

Introduce the QR factorization

$LQ_{k}^{(r)}=Q_{k}^{(L)}R_{k}^{(L)},$ (17)

where the matrix $Q_{k}^{(L)}\in{\mathbb{R}}^{n\times k}$ has orthonormal
columns and $R_{k}^{(L)}\in{\mathbb{R}}^{k\times k}$ is upper triangular. We
note that since the matrix $L$ typically is very sparse and $k$ is not large,
the left-hand side of (17) generally can be evaluated quite quickly also when
$n$ is large. The minimization problem (15) yields the small problem

$\min_{\mbox{\boldmath{${\scriptstyle
y}$}}\in{\mathbb{R}}^{k}}\\{\|R_{k}^{(c)}(R_{k}^{(r)})^{T}\mbox{\boldmath{$y$}}-(Q_{k}^{(c)})^{T}\mbox{\boldmath{$g$}}\|^{2}+\mu\|R_{k}^{(L)}\mbox{\boldmath{$y$}}\|^{2}\\}.$
(18)

This problem can be solved in several ways: We may compute a generalized SVD
(GSVD) of the matrix pair $\\{R_{k}^{(c)}(R_{k}^{(r)})^{T},R_{k}^{(L)}\\}$
(see, e.g., [10, 23]), or apply a cheaper reduction of the matrix pair that
can be used when the generalized singular values of the matrix pair are not
explicitly required; see [11].

When the matrix $R_{k}^{(L)}$ is nonsingular and not very ill-conditioned,
which is the case in many applications, one may consider transforming the
minimization problem (18) to standard form by the substitution
$\mbox{\boldmath{$z$}}=R_{k}^{(L)}\mbox{\boldmath{$y$}}$. This yields the
problem

$\min_{\mbox{\boldmath{${\scriptstyle
z}$}}\in{\mathbb{R}}^{k}}\\{\|R_{k}^{(c)}(R_{k}^{(r)})^{T}(R_{k}^{(L)})^{-1}\mbox{\boldmath{$z$}}-(Q_{k}^{(c)})^{T}\mbox{\boldmath{$g$}}\|^{2}+\mu\|\mbox{\boldmath{$z$}}\|^{2}\\},$
(19)

which easily can be solved, e.g., by computing the singular value
decomposition of the matrix $R_{k}^{(c)}(R_{k}^{(r)})^{T}(R_{k}^{(L)})^{-1}$.
The solution $\mbox{\boldmath{$z$}}_{\mu}$ of (19) yields the solution
$\mbox{\boldmath{$y$}}_{\mu}=(R_{k}^{(L)})^{-1}\mbox{\boldmath{$z$}}_{\mu}$ of
(15), from which we determine the approximate solution
$\widetilde{\mbox{\boldmath{$x$}}}_{\mu}=Q_{k}^{(r)}\mbox{\boldmath{$y$}}_{\mu}$
of (3). The solution of (19) is cheaper than the solution of (18) with the aid
of the GSVD; see [19] for counts of the arithmetic floating point operations
necessary to compute the GSVD of a pair of $k\times k$ matrices, and the SVD
of a $k\times k$ matrix.

###  3.2 The discrepancy principle

We turn to the computation of the regularization parameter $\mu>0$ by the
discrepancy principle. Assume for the moment that the matrix $A$ is available.
Then we can solve equation (7) for $\mu>0$ by using a zero-finder such as
Newton’s method or one of the zero-finders described in [8, 34]. The
theoretical justification of the discrepancy principle requires that the
unavailable error-free vector $\widehat{\mbox{\boldmath{$g$}}}$ associated
with the available error-contaminated vector $g$ satisfies
$\widehat{\mbox{\boldmath{$g$}}}\in{\mathcal{R}}(A)$.

Now consider the application of the discrepancy principle to the determination
of the regularization parameter in (14). Generally,
$\widehat{\mbox{\boldmath{$g$}}}\not\in{\mathcal{R}}(M_{k})$ and, therefore,
the discrepancy principle cannot be applied when solving (14) without
modification. In the computed examples, we determine $\mu$ so that the
computed solution ${\mbox{\boldmath{$x$}}}_{\mu}$ of (14) satisfies

$\|M_{k}\mbox{\boldmath{$x$}}_{\mu}-Q_{k}^{(c)}(Q_{k}^{(c)})^{T}\mbox{\boldmath{$g$}}\|=\eta\delta;$
(20)

cf. (7). If the matrices $R_{k}^{(c)}$ and $R_{k}^{(r)}$ in (13) are
nonsingular, which generally is the case, then $Q_{k}^{(c)}(Q_{k}^{(c)})^{T}$
is an orthogonal projector onto ${\mathcal{R}}(M_{k})$, and
$Q_{k}^{(c)}(Q_{k}^{(c)})^{T}\widehat{\mbox{\boldmath{$g$}}}$ lives in
${\mathcal{R}}(M_{k})$. Equation (20) can be solved for $\mu\geq 0$ by using a
zero-finder.

###  3.3 Stopping criterion for the adaptive cross approximation algorithm

In view of (7), we would like to determine a value of the regularization
parameter $\mu>0$ such that

$\|A\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|=\eta\delta.$
(21)

Even though the matrix $A$ is not available, we can determine an approximate
upper bound for the left-hand side of (21) as follows:

$\displaystyle\|A\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|$
$\displaystyle=$
$\displaystyle\|(A-M_{k})\mbox{\boldmath{$x$}}_{\mu}(M_{k})+M_{k}\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|$
(22) $\displaystyle\leq$
$\displaystyle\|A-M_{k}\|\|\mbox{\boldmath{$x$}}_{\mu}(M_{k})\|+\|M_{k}\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|$
$\displaystyle\lessapprox$ $\displaystyle
S_{k}\|\mbox{\boldmath{$x$}}_{\mu}(M_{k})\|+\|M_{k}\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|,$

where $S_{k}$ is an approximation of $\|A-M_{k}\|$. Based on the values of
$(R_{k})_{i_{\ell},j_{\ell}}$ from (10) such an approximation $S_{k}$ can be
computed as

$\displaystyle\|A-M_{k}\|$ $\displaystyle\leq$ $\displaystyle\|A-M_{k}\|_{F}$
(23) $\displaystyle\approx$
$\displaystyle\sqrt{\sum_{(i_{\ell},j_{\ell})\in\Pi}|(A-M_{k})_{i_{\ell},j_{\ell}}|^{2}}(mn)/|\Pi|$
$\displaystyle=$
$\displaystyle\sqrt{\sum_{(i_{\ell},j_{\ell})\in\Pi}|(R_{k})_{i_{\ell},j_{\ell}}|^{2}}(mn)/|\Pi|$
$\displaystyle=:$ $\displaystyle S_{k},$

where $\|\cdot\|_{F}$ denotes the matrix Frobenius norm, and the set $\Pi$ is
defined by (9) with $|\Pi|=t$ elements.

The number of step $k=k^{*}$ of the adaptive cross approximation algorithm is
chosen to be as small as possible such that there exists a $\mu>0$ such that

$S_{k}\|\mbox{\boldmath{$x$}}_{\mu}(M_{k})\|=\eta_{1}\delta\mbox{ and
}\|M_{k}\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|\leq\eta_{2}\delta.$

In this case we have, based on (22)

$\displaystyle\|A\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|$
$\displaystyle\lessapprox$ $\displaystyle
S_{k}\|\mbox{\boldmath{$x$}}_{\mu}(M_{k})\|+\|M_{k}\mbox{\boldmath{$x$}}_{\mu}(M_{k})-\mbox{\boldmath{$g$}}\|$
$\displaystyle\lessapprox$ $\displaystyle(\eta_{1}+\eta_{2})\delta.$

##  4 Numerical experiments

This section describes a few computed examples with the adaptive cross
approximation method. For problems in one space-dimension, we will use the
regularization matrices $L_{0}=I$, where $I\in{\mathbb{R}}^{n\times n}$
denotes the identity matrix, and

$L_{1}=\frac{1}{2}\left[\begin{array}[]{cccccc}1&-1&&&&\mbox{\Large 0}\\\
&\phantom{-}1&-1&&&\\\ &&\phantom{-}1&-1&&\\\ &&&\ddots&\ddots&\\\
\mbox{\Large
0}&&&&\phantom{-}1&-1\end{array}\right]\in{{\mathbb{R}}}^{(n-1)\times n}$ (24)

or

$L_{2}=\frac{1}{4}\left[\begin{array}[]{cccccc}-1&\phantom{-}2&-1&&&\mbox{\Large
0}\\\ &-1&\phantom{-}2&-1&&\\\ &&\ddots&\ddots&\ddots&\\\ \mbox{\Large
0}&&&-1&\phantom{-}2&-1\end{array}\right]\in{{\mathbb{R}}}^{(n-2)\times n},$
(25)

which approximate a multiple of the first and second order derivative
operators, respectively, assuming that $\mbox{\boldmath{$x$}}_{\mu}$ is the
discretization of a function $x_{\mu}$ at equidistant points on a bounded
interval $\Omega_{1}$. For problems in two space-dimensions, we use the
regularization matrices

$L_{1,\otimes}=\left[\begin{array}[]{ccc}I&\otimes&L_{1}\\\
L_{1}&\otimes&I\end{array}\right]$ (26)

or

$L_{2,\otimes}=\left[\begin{array}[]{ccc}I&\otimes&L_{2}\\\
L_{2}&\otimes&I\end{array}\right],$ (27)

where $\otimes$ stands for the Kronecker product. These choices of
regularization matrices are fairly common; see, e.g., [23, 26, 32] for
illustrations.

The examples are taken from the Regularization Tools [24] and from IR Tools
[18]. In all examples, the number of elements in the set $\Pi$ is $50n$. The
values of the parameters $\eta,\eta_{1},\eta_{2}$ are chosen as
$\eta=\eta_{1}=\eta_{2}=1.0$.

Experiment 1: We consider the problem “gravity” of size $n=1024$ from [24].
The aim of this experiment is to illustrate that the quantities $S_{k}$
defined by (23) provide quite accurate approximations of $\|A-M_{k}\|$. This
is displayed by Figure 1.

Fig. 1: For the “gravity” problem the values of $S_{k}$, $\|A-M_{k}\|$, and
$\|A-M_{k}\|_{F}$ are plotted as a function of the number of steps $k$ of the
adaptive cross approximation algorithm.

Experiment 2: We again consider the example “gravity” from [24] of size
$n=1024$. Let $\delta=10^{-2}$. This example illustrates that for certain
problems only a fairly small number of steps of the adaptive cross
approximation algorithm suffices to yield a satisfactory result. The example
also shows that it may be very beneficial to use a regularization matrix
different from the identity matrix. The maximum number of adaptive cross
approximation steps is $30$. Results are shown for
$L\in\\{L_{0},L_{1},L_{2}\\}$. The quality of the computed solution
$\mbox{\boldmath{$x$}}_{\rm computed}$ is measured by the relative error
$\|\mbox{\boldmath{$x$}}_{\rm
computed}-\hat{\mbox{\boldmath{$x$}}}\|/\|\hat{\mbox{\boldmath{$x$}}}\|$. The
horizontal axis of Figure 2 shows the number of steps of adaptive cross
approximation; the vertical line indicates that for each one of the choices of
$L$, the stopping criterion for the method is satisfied at step $k^{*}=20$.

Fig. 2: For the “gravity” problem the relative error between the exact
solution and the computed approximate solution is plotted as a function of the
iteration step $k$ and the chosen regularization matrix $L_{0},L_{1},L_{2}$.
The size of the problem is $n=1024$ with $\delta=1.0e-2$ and
$\eta=\eta_{1}=\eta_{2}=1.0$. The vertical line indicates where the stopping
criterion is satisfied for adaptive rank approximation.

Experiment 3: This experiment is similar to Experiment 2, but for problem
“baart” from [24]. Results are displayed in Figure 3. Also in this example it
is beneficial to use a regularization matrix different from the identity.

Fig. 3: For the “baart” problem the relative error between the exact solution
and the computed approximate solution is plotted as a function of the
iteration step $k$ and the chosen regularization matrix $L_{0},L_{1},L_{2}$.
The size of the problem is $n=1024$ with $\delta=1.0e-2$ and
$\eta=\eta_{1}=\eta_{2}=1.0$. The vertical line indicates where the stopping
criterion is satisfied for adaptive rank approximation.

Experiment 4: This expperiment is similar to Experiment 2, but for problem
“phillips” from [24]. The result is shown in Figure 4. For this example all
three regularization matrices used perform about equally well. The singular
values of the matrix $A$ decay to zero slower for this example than for the
previous examples. Therefore more steps with the adaptive cross approximation
algorithm have to be carried out.

Fig. 4: For the “phillips” problem the relative error between the exact
solution and the computed approximate solution is plotted as a function of the
iteration step $k$ and the chosen regularization matrix $L_{0},L_{1},L_{2}$.
The size of the problem is $n=1024$ with $\delta=1.0e-2$ and
$\eta=\eta_{1}=\eta_{2}=1.0$. The vertical line indicates where the stopping
criterion is satisfied for adaptive rank approximation.

Experiment 5: We consider the example EXdiffusion_rrgmres from the IR Toolbox
[18]. The size of the problem is $4096$ and
$\delta=5.0e-3\|\widehat{\mbox{\boldmath{$g$}}}\|$. The other parameters are
as in Experiment 2. The “best” solution determined by the example script in
the IR Toolbox has relative error (when compared with the exact solution)
$0.1875$, while for $L=L_{2}$ our algorithm in iteration step $k=76$ reaches a
relative error of $0.1910$.

Fig. 5: For the “diffusion” problem the relative error between the exact
solution and the computed approximate solution is plotted as a function of the
iteration step $k$ and the chosen regularization matrix $L$. The size of the
problem is $n=4096$ with $\delta=5.0e-3\|\widehat{\mbox{\boldmath{$g$}}}\|$
and $\eta=\eta_{1}=\eta_{2}=1.0$. The vertical line indicates where the
stopping criterion is satisfied for adaptive rank approximation.

Experiment 6: The same as in Experiment 2 where we consider $A=B\otimes B$
with $B$ the matrix of the “baart” regularization problem. The matrix $B$ is
of order $40$. We use (27) for $L_{0}$, $L_{1}$ and $L_{2}$ as regularization
matrices. The relative errors of the computed solutions are displayed in
Figure 6.

Fig. 6: For the 2D “baart” problem the relative error between the exact
solution and the computed approximate solution is plotted as a function of the
iteration step $k$ and the chosen 2D regularization matrix based on
$L_{0},L_{1},L_{2}$. The size of the problem is $n=40^{2}$ with
$\delta=1.0e-3$ and $\eta=\eta_{1}=\eta_{2}=1.0$. The vertical line indicates
where the stopping criterion is satisfied for adaptive rank approximation.

##  5 Conclusion

This paper discusses the application of adaptive cross approximation to
Tikhonov regularization problems in general form. The computed examples
illustrate that often only quite few cross approximation steps are required to
yield useful approximate solutions. Particular attention is given to the
stopping criterion for adaptive cross approximation.

## Acknowledgment

This research was partially supported by the Fund for Scientific
Research–Flanders (Belgium), Structured Low-Rank Matrix/Tensor Approximation:
Numerical Optimization-Based Algorithms and Applications: SeLMA, EOS 30468160,
the KU Leuven Research Fund, Numerical Linear Algebra and Polynomial
Computations, OT C14/17/073.

## References

  * [1] M. L. Baart, The use of auto-correlation for pseudo-rank determination in noisy ill-conditioned least-squares problems, IMA J. Numer. Anal., 2 (1982), pp. 241–247.
  * [2] M. Bebendorf, Approximation of boundary element matrices, Numer. Math., 86 (2000), pp. 565–589.
  * [3] M. Bebendorf and R. Grzibovski, Accelerating Galerkin BEM for linear elasticity using adaptive cross approximation, Math. Methods Appl. Sci., 29 (2006), pp. 1721–1747.
  * [4] M. Bebendorf and S. Rjasanow, Adaptive low-rank approximation of collocation matrices, Computing, 70 (2003), pp. 1–24.
  * [5] S. Börm, Efficient Numerical Methods for Non-Local Operators: $\mathcal{H}_{2}$-Matrix Compression, Algorithms and Analysis, European Math. Society, Zürich, 2010.
  * [6] C. Brezinski, G. Rodriguez, and S. Seatzu, Error estimates for linear systems with applications to regularization, Numer. Algorithms, 49 (2008), pp. 85–104.
  * [7] C. Brezinski, M. Redivo-Zaglia, G. Rodriguez, and S. Seatzu, Multi-parameter regularization techniques for ill-conditioned linear systems, Numer. Math., 94 (2003), pp. 203–224.
  * [8] A. Buccini, M. Pasha, and L. Reichel, Generalized singular value decomposition with iterated Tikhonov regularization, J. Comput. Appl. Math., 373 (2020), Art. 112276
  * [9] D. Calvetti, J. Petersen, and L. Reichel, A parallel implementation of the GMRES algorithm, in Numerical Linear Algebra, eds. L. Reichel, A. Ruttan, and R. S. Varga, de Gruyter, Berlin, 1993, pp. 31–46.
  * [10] L. Dykes, S. Noschese, L. Reichel, Rescaling the GSVD with application to ill-posed problems, Numer. Algorithms, 68 (2015), pp. 531–545.
  * [11] L. Dykes, L. Reichel, Simplified GSVD computations for the solution of linear discrete ill-posed problems, J. Comput. Appl. Math., 255 (2013), pp. 15–27.
  * [12] L. Eldén, A weighted pseudoinverse, generalized singular values, and constrained least squares problems, BIT Numer. Math., 22 (1982), pp. 487–501.
  * [13] H. W. Engl, M. Hanke, and A. Neubauer, Regularization of Inverse Problems, Kluwer, Dordrecht, 1996.
  * [14] J. Erhel, A parallel GMRES version for general sparse matrices, Electron. Trans. Numer. Anal., 3 (1995), pp. 160–176.
  * [15] C. Fenu, L. Reichel, G. Rodriguez, and H. Sadok, GCV for Tikhonov regularization by partial SVD, BIT Numer. Math., 57 (2017), pp. 1019–1039.
  * [16] L. Fox and E. T. Goodwin, The numerical solution of non-singular linear integral equations, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 245:902 (1953), pp. 501–534.
  * [17] K. Frederix and M. Van Barel, Solving a large dense linear system by adaptive cross approximation, J. Comput. Appl. Math., 234 (2010), pp. 3181–3195.
  * [18] S. Gazzola, P. C. Hansen, and J. G. Nagy, IR Tools: A MATLAB package of iterative regularization methods and large-scale test problems, Numer. Algorithms, 81 (2019), pp. 773–811.
  * [19] G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, Baltimore, 2013.
  * [20] S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, A theory of pseudo-skeleton approximation, Linear Algebra Appl., 261 (1997), pp. 1–21.
  * [21] S. A. Goreinov, N. L. Zamarashkin, and E. E. Tyrtyshnikov, Pseudo-skeleton approximations by matrices of maximal volume, Math. Notes, 62 (1997), pp. 515–519.
  * [22] C. W. Groetsch, The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind, Pitman, Boston, 1984.
  * [23] P. C. Hansen, Rank-Deficient and Discrete Ill-Posed Problems, SIAM, Philadelphia, 1998\.
  * [24] P. C. Hansen, Regularization Tools version 4.0 for Matlab 7.3, Numer. Algorithms, 46 (2007), pp. 189–194. Software is available in Netlib at http://www.netlib.org.
  * [25] G. Huang, S. Noschese, and L. Reichel, Regularization matrices determined by matrix nearness problems, Linear Algebra Appl., 502 (2016), pp. 41–57.
  * [26] G. Huang, L. Reichel, and F. Yin, On the choice of subspace for large-scale Tikhonov regularization problems in general form, Numer. Algorithms, 81 (2019), pp. 33–55.
  * [27] H2Lib, http://www.h2lib.org/, 2015–2020.
  * [28] S. Kindermann and K. Raik, A simplified L-curve method as error estimator, Electron. Trans. Numer. Anal., 53 (2020), pp. 217–238.
  * [29] T. Mach, L. Reichel, M. Van Barel, and R. Vandebril, Adaptive cross approximation for ill-posed problems, J. Comput. Appl. Math., 303 (2016), pp. 206–217.
  * [30] S. Morigi, L. Reichel, and F. Sgallari, Orthogonal projection regularization operators, Numer. Algorithms, 44 (2007), pp. 99–114.
  * [31] S. Noschese and L. Reichel, Inverse problems for regularization matrices, Numer. Algorithms, 60 (2012), pp. 531–544.
  * [32] Y. Park, L. Reichel, G. Rodriguez, and X. Yu, Parameter determination for Tikhonov regularization problems in general form, J. Comput. Appl. Math., 343 (2018), pp. 12–25.
  * [33] L. Reichel and G. Rodriguez, Old and new parameter choice rules for discrete ill-posed problems, Numer. Algorithms, 63 (2013), pp. 65–87.
  * [34] L. Reichel and A. Shyshkov, A new zero-finder for Tikhonov regularization, BIT Numer. Math., 48 (2008), pp. 627–643.
  * [35] L. Reichel and Q. Ye, Simple square smoothing regularization operators, Electron. Trans. Numer. Anal., 33 (2009), pp. 63–83.
  * [36] C. B. Shaw, Jr., Improvements of the resolution of an instrument by numerical solution of an integral equation, J. Math. Anal. Appl., 37 (1972), pp. 83–112.
  * [37] E. E. Tyrtyshnikov, Incomplete cross approximation in the mosaic-skeleton method, Computing, 64 (2000), pp. 367–380.
  * [38] Y. Yamamoto, Y. Nakatsukasa, Y. Yanagisawa, and T. Fukaya, Roundoff error analysis of the CholeskyQR2 algorithm, Electron. Trans. Numer. Anal., 44 (2015), pp. 306–326.