Menke W. Geophysical Data Analysis: Discrete Inverse Theory
Подождите немного. Документ загружается.
I28
7
Applications
of
Vector Spaces
Here
y
can be partitioned into two parts
yE
and
y,
(possibly requiring
reordering of the constraints) that satisfy
F,m
-
hE
=
0
(7.50)
Y=[’~>’]
Ys
=
0
and
F,m
-
h,
>
0
The first group of equality-inequality constraints are satisfied in the
equality sense (thus the subscript
E
for “equality”). The rest are
satisfied more loosely in the inequality sense (thus the subscript
S
for
“slack”).
The theorem states that any feasible solution
m
is the minimum
solution only
if
the direction in which one would have to perturb
m
to
decrease the total error
E
causes the solution to cross some constraint
hyperplane and become infeasible. The direction of decreasing error is
-
VE
=
GT[d
-
Gm].
The constraint hyperplanes have normals
+
V[Fm]
=
F’
which point into the feasible side. Since
F,m
-
h,
=
0,
the solution lies exactly on the bounding hyperplanes of the
F,
con-
straints but within the feasible volume of the
F,
constraints. An
infinitesimal perturbation
6m
of
the solution can, therefore, only
violate the
F,
constraints.
If
it
is not to violate these constraints, the
perturbation must be made in the direction of feasibility,
so
that it
must be expressible as a nonnegative combination of hyperplane
normals
6m
-
V[Fm]
2
0.
On the other hand,
if
it is to decrease the
total prediction error it must satisfy
6m
*
VE
5
0.
For solutions that
satisfy the Kuhn
-
Tucker theorem these two conditions are incompat-
ible and
6m
*
VE
=
dm
-
V[Fm]
*
y
2
0
since both
6m
-
V[Fm]
and
y
are positive. These solutions are indeed minimum solutions to the
constrained problem (Fig.
7.5).
To
demonstrate how the Kuhn-Tucker theorem can be used, we
consider the simplified problem
Minimize
Ild
-
Gmll,
subject to
m
2
0
(7.51)
We find the solution using an iterative scheme of several steps [Ref.
141.
Sfep
I.
Start with an initial guess for
m.
Since
F
=
I
each model
parameter
is
associated with exactly one constraint. These model
parameters can be separated into a set
mE
that satisfies the constraints
in the equality sense and a set
m,
that satisfies the constraints in the
inequality sense. The particular initial guess
m
=
0
is clearly feasible
and has all its elements in
mE.
7.9
Inequality Constraints
I
129
Fig.
7.5.
Error
E(m)
(contoured) has a single minimum (dot). The linear inequality
constraint (straight line) divides
S(m)
into a feasible and infeasible halfspace (arrows
point into feasible halfspace). Solution (triangle) lies on the boundary of the halfspaces
and therefore satisfies the constraint in the equality sense. At this point the normal ofthe
constraint hyperplace is antiparallel
to
-
OE.
Step
2.
Any model parameter
m,
in
m,
that has associated with it a
negative gradient
[VE],
can be changed both to decrease the error and
to remain feasible. Therefore, if there is no such model parameter in
m,,
the Kuhn
-
Tucker theorem indicates that this
m
is
the solution to
the problem.
Step
3.
If
some model parameter
mi
in
mE
has a corresponding
negative gradient, then the solution can be changed to decrease the
prediction error.
To
change the solution, we select the model parame-
ter corresponding to the most negative gradient and move it to the set
ms
.
All
the model parameters in
m,
are now recomputed by solving the
system
G,m&
=
d,
in the least squares sense. The subscript
S
on the
matrix indicates that only the columns multiplying the model parame-
ters in
m,
have been included in the calculation.
All
the
mE)s
are still
zero.
If
the new model parameters are all feasible, then we set
m
=
m’
and return to Step
2.
Step
4.
If
some of the
mL’s
are infeasible, however, we cannot use this
vector as a new guess for the solution. Instead, we compute the change
in the solution
6m=
mL-m,
and add
as
much of this vector as
130
7
Applications
of
Vector Spaces
possible to the solution
m,
without causing the solution to become
infeasible. We therefore replace
m,
with the new guess
m,
+
adm,
where
a
=
min,(m,,/[ms,
-
mk,])
is the largest choice that can be
made without some
m,
becoming infeasible. At least one of the
ms,’s
has its constraint satisfied in the equality sense and must be moved
back to
mE.
The process then returns to Step
3.
This algorithm contains two loops, one nested within the other. The
outer loop successively moves model parameters from the group that is
constrained to the group that minimizes the prediction error. The
inner loop ensures that the addition of a variable to this latter group
has not caused any of the constraints to be violated. Discussion of the
convergence properties of this algorithm can be found in Ref. 14.
This algorithm can also be used to solve the problem [Ref. 141
Minimize
Ilmll,
subject to
Fm
L
h
(7.52)
and, by virtue of the transformation described in Section
7.9.1,
the
completely general problem. The method consists of forming the
(7.53)
and finding the
m’
that minimizes
Ild’
-
G’m’ll,
subject to
m’
2
0
by
the above algorithm. If the prediction error
e‘=
d‘- G’m’
is
identically zero, then the constraints
Fm
2
h
are inconsistent.
Otherwise, the solution is
mi
=
el/e$+
I
.
We shall show that this method does indeed solve the indicated
problem [adapted from Ref. 141. We first note that the gradient of the
error is
VE’
=
GtT[d’
-
G’m’]
=
-GITe‘,
and that because of the
Kuhn
-
Tucker theorem,
m’
and
VE’
satisfy
mf:
=
0
mk
>
0
[VE’IE
<
0
[VE’IS
=
0
(7.54)
The length of the error is therefore
If the error is not identically zero,
e&+
I
is greater than zero. We can use
this result to
show
that the constraints
Fm
2
h
are consistent and
m
is
7.9
Inequality Constraints
131
feasible since
0
I
-VE‘
=
GITe’
=
[F, h][m’,
-
lITeh+,
=
[Fm
-
h]ea+,
(7.56)
Since
eh+,
>
0,
we have
[Fm
-
h]
2
0.
The solution minimizes
Ilmll,
because the Kuhn
-
Tucker condition that the gradient
V
II
mll
,
=
m
be
represented as a nonnegative combination
of
the rows
of
F
is
satisfied:
Vllmll,
=
m
=
[e;,
.
.
.
,
ehIT/eh+,
=
FTm’/eh+,
(7.57)
Here
m’
and
eh+,
are nonnegative. Finally, assume that the error is
identically zero but that a feasible solution exists. Then
0
=
elTe’/eh+,
=
[mT,
-
l][d’
-
G’m’]
=
-
1
-
[Fm
-
hITm’
(7.58)
Since
m’
2
0,
the relationship
Fm
<
h
is implied. This contradicts the
constraint equations
Fm
2
h,
so
that an identically zero error implies
that no feasible solution exists and that the constraints are inconsis-
tent.
This page intentionally left blank
LINEAR INVERSE
GAUSSIAN DISTRIBUTIONS
PROBLEMS
AND
NON-
8.1
L,
Norms
and Exponential
Distributions
In Chapter
5
we showed that the method
of
least squares and the
more general use of
L,
norms could be rationalized through the
assumption that the data and a priori model parameters followed
Gaussian statistics. This assumption
is
not always appropriate, how-
ever; some data sets follow other distributions. The
exponential
distrz-
bution
is one simple alternative. When Gaussian and exponential
distributions of the same mean
(d)
and variance
o2
are compared, the
exponential distribution is found to be much longer-tailed (Fig.
8.1
and Table
8.1):
Gaussian
ExDonential
(8.1)
I33
134
8
Non
-Gaussian Distributions
TABLE
8.1
The area beneath
a
portion
of
the Gaussian and
exponential distributions, centered on the mean and
having the given half-width‘
Gaussian Exponential
Half-width Area
(%)
Half-width Area
(%)
la
68.2
la
16
20 95.4 2a 94
3a 99.7 30 98.6
4a 99.999+ 4a 99.1
~~~
~
Note that the exponential distribution is longer-
tailed.
Note that the probability of realizing data far from
(d)
is much higher
for the exponential distribution than for the Gaussian distribution.
A
few data several standard deviations from the mean are reasonably
probable in a data set of, say,
1000
samples drawn from an exponential
distribution but very improbable for data drawn from a Gaussian
distribution. We therefore expect that methods based on the exponen-
tial distribution will be able to handle a data set with a few “bad” data
(outliers) better than Gaussian methods. Methods that can tolerate a
few outliers are said to be
robust
[Ref.
41.
0.75r
0.50
-
h
Y
‘0
n
0.25
-
d
Fig.
8.1.
Gaussian (curve A) and exponential (curve
B)
distribution with zero mean
and unit variance. The exponential distribution has the longer tail.
8.2
Maximum Likelihood Estimate
of
Mean
135
8.2
Maximum Likelihood Estimate
of
the
Mean
of
an Exponential Distribution
Exponential distributions bear the same relationship to
L,
norms as
Gaussian distributions bear to
L,
norms.
To
illustrate this relation-
ship, we consider the joint distribution for
N
independent data, each
with the same mean
(d)
and variance
02.
Since the data are indepen-
dent, the joint distribution is just the product of
N
univanate distribu-
tions:
To maximize the likelihood of
P(d)
we must maximize the argument
of the exponential, which involves minimizing the sum of absolute
residuals as
N
Minimize
E
=
2
Id,
-
(d)
I
(8.3)
i-
I
This
is
the
L,
norm of the prediction error in a linear inverse problem
of the form
Gm
=
d,
where
M=
1,
G
=
[I,
1,
.
. .
,
1IT,
and
m
=
[(
d)].
Applying the principle of maximum likelihood, we obtain
Maximize L=logP=--log2 -Nloga-T
xIdl--
(d)l
N 2112
N
2
1-
I
(8.4)
Setting the derivatives to zero yields
where the sign function sign@) equals
+
1
ifx
3
0,
-
1
ifx
<
0
and
0
if
x
=
0.
The first equation yields the implicit expression for
(d)cst
for
which
Xi
sign(d,
-
(d))
=
0.
The second equation can then be solved
for an estimate of the variance as
136
8
Non-Gaussian
Distributions
The equation for
(
d)est
is exactly the sample median; one finds a
(d)
such that half the
d:s
are less than
(d)
and half are greater than
(d).
There is then an equal number
of
negative and positive signs and the
sum
of
the signs is zero. The median is a robust property of a set of
data. Adding one outlier can at worst move the median from one
central datum to another nearby central datum. While the maximum
likelihood estimate of a Gaussian distribution's true mean is the
sample arithmetic mean, the maximum likelihood estimate of an
exponential distribution's true mean is the sample median.
The estimate
of
the variance also differs between the two distribu-
tions: in the Gaussian case it is the square
of
the sample standard
deviation but in the exponential case it is not. If there is an odd number
of samples, then
(
d)est
equals the middle
d;:
if there is an even number
of
samples, any
(d)est
between the two middle data maximizes the
likelihood. In the odd case the error
E
attains a minimum only at the
middle sample, but in the even case it is flat between the two middle
samples (Fig.
8.2).
We see, therefore, that
L,
problems of minimizing
the prediction error of
Gm
=
d
can possess nonunique solutions that
are distinct from the type
of
nonuniqueness encountered in the
L2
problems. The
L,
problems can still possess nonuniqueness owing to
the existence of null solutions since the null solutions cannot change
the prediction error under any norm. That kind of nonuniqueness
leads to a completely unbounded range
of
estimates. The new type of
nonuniqueness,
on
the other hand, permits the solution to take on any
values between
finite
bounds.
E
____r__)
<dest)
datum d
range
of(dDTt)/+
datum
d
Fig.
8.2.
The
L,
error
in determining the mean
of
(a) an odd number
of
data and
(b)
an
even number
of
data.
8.3
The General Linear Problem
137
We also note that regardless of whether
N
is even or odd, we can
choose
(
d)est
so
that one of the equations
Gm
=
d
is satisfied exactly
(in this case,
(
d)est
=
dmid,
where mid denotes “middle”). This can be
shown to be a general property
of
L,
norm problems. Given
N
equations and Munknowns related by
Gm
=
d,
it is possible to choose
m
so that the
L,
prediction error is minimized and
so
that
M
of the
equations are satisfied exactly.
8.3
The
General Linear Problem
Consider the linear inverse problem
Gm
=
d
in which the data and a
priori model parameters are uncorrelated with known means
dobs
and
(m)
and known variances
03
and
02,
respectively. The joint distribu-
tion is then
i=
I
1=l
where the prediction error is given by
e
=
d
-
Gm
and the solution
length by
t
=
m
-
(m).
The maximum likelihood estimate of the
model parameters occurs when the exponential is a minimum, that is,
when the sum of the weighted
L,
prediction error and the weighted
L,
solution length is minimized:
N
M
Minimize
E
+
L
=
2
lel(c;ll
+
2
ltlla;:
(8.8)
1-1
1-
I
In this case the weighting factors are the reciprocals of the standard
deviations- this in contrast to the Gaussian case, in which they are
the reciprocals of the variances. Note that linear combinations of
exponentially distributed random variables are not themselves expo-
nential (unlike Gaussian variables, which give rise to Gaussian combi-
nations). The covariance matrix of the estimated model parameters is,
therefore, difficult both to calculate and to interpret since the manner
in which it
is
related to confidence intervals varies from distribution to
distribution. We shall focus our discussion on estimating only the
model parameters themselves.