The Backus-Gilbert (BG) method, Backus & Gilbert (1968), also
known as the method of optimally localized averages, has been used
successfully in seismology, Parker (1977), helioseismology,
Christensen-Dalsgaard et al. (1990) and other areas that
involve inverse problems, Craig & Brown (1986). The BG method
belongs to a class of methods known as mollifier methods. The
underlying idea in these methods is that if we are given m
measurements 
which represent (noisy) values of m functionals
on an unknown function f, i.e.,
where 
represent the noise,
then we estimate f in a given point 
by a linear combination of the measurements:
Here, we have introduced the two m-vectors
If the noise is unbiased and has covariance matrix 
, then the
variance 
of the error in 
due to the noise
in 
is given by
The BG method is a particular technique for
computing the coefficients 
explicitly. If we compare this
approach with Tikhonov regularization of a discrete ill-posed problem
, whose solution is formally given by
for some 
, then we
notice that 
in a mollifier method conceptually corresponds to
a row in the matrix 
.
If 
denotes the exact solution, then for any regularization
method it is easy to see from Eq. (1 (click here)) that 
is given by
The function
is called the averaging kernel at 
, and it shows how the noise-free
part of the regularized estimate 
is obtained as an average of the
exact solution 
. Thus, 
determines the
resolving power of the particular regularization method at 
,
and for 
to be meaningful the function 
should be peaked around 
and be normalized such that
Hence, the averaging kernel
often plays a role which is at least as important
as the computed estimate 
.
In the BG method the peakedness of 
is
controlled by specifying a criteria function 
and requiring that the quantity
be sufficiently small, subject to the constraint in Eq.
(5 (click here)). Often, a regularization term is also added in
Eq. (6 (click here)). A commonly used criteria function is
. Various other criteria functions are
discussed in the original paper by Backus &\
Gilbert (1968), while the additional regularization term is not
used there. See Hansen (1994) for an SVD-analysis and
computational details.
The BG method is computationally very expensive: If n is the number of quadrature points used to evaluate the integrals then 
operations are required for each 
. This is prohibitively
expensive in large-scale applications. In, e.g., 2D rotational
inversion in helioseismology m is of the order 
and n is of the order 
.
Therefore, there is a need for computationally more efficient
methods. One such method is subtractive optimally localized
averages (SOLA), which has been discovered independently by several
researchers; see, e.g., Louis & Maass (1990) and Pijpers &\
Thompson (1992). The idea in the SOLA method is, instead of
minimizing (6 (click here)), to minimize the squared difference
between the averaging kernel and some chosen target function 
,
and it is common practice to retain
the constraint in Eq. (5 (click here)), see , Pijpers &\
Thompson (1992, 1994).
The peakedness of 
is controlled by the choice
of 
, and therefore to compute 
one would choose a
target function peaked around 
.
In practice it is often necessary to add a regularizing term, i.e., to
replace Eq. (7 (click here)) with
Here 
is the error covariance matrix of the measured
data points 
and 
is a trade-off or regularization
parameter whose value determines the relative importance of the two
terms in Eq. (8 (click here)). Thus, 
controls the propagation
of errors from the measurements to 
.
The great advantage of the SOLA method compared to the BG method is
that the 
coefficient matrix involved does not depend on
and therefore its factorization, costing 
operations, is performed only once. The additional cost to compute
is 
for each 
.
In this paper we discuss two different approaches to implementing the
SOLA method: The standard approach using a Lagrange multiplier, and a
new approach that uses an orthogonal transformation to incorporate the
constraint in Eq. (5 (click here)) and bidiagonalization
algorithms to solve the resulting least squares problem. We
demonstrate that the latter approach is more efficient when many
values of the regularization parameter 
must be tried in
order to find the optimal one. This is important, since the optimal
is very rarely known in advance - instead it must be
computed by means of parameter-choice methods, such as
generalized cross validation or the L-curve criterion, that
involve the optimization of a function that depends on 
;
see Hanke & Hansen (1993).
Our paper is organized as follows. In Sect. 2 (click here) we introduce our notation and present the algorithms. These algorithms are based on standard linear algebra routines which are available in most scientific software libraries, and in Sect. 3 (click here) we point to available Fortran software. The computational load in the algorithms, and the efficiency of our new approach, is discussed in Sect. 4 (click here). Finally, in Sect. 5 (click here) we illustrate our algorithms with an example from helioseismology.