Menke W. Geophysical Data Analysis: Discrete Inverse Theory
Подождите немного. Документ загружается.
248
13
Numerical
Algorithms
Program
PLS,
continued.
11
12
13
14
do
17
j
=
ntr, mcur
/*
triangularize new columns
alfa
=
0.0
/*
build transformation
k
=
igp(j)
/*
location of j-th column of triangle
do
11
i=j, n
alfa
=
alfa
t
gp(i,k)**2
continue
alfa
=
sqrt(
alfa
)
if( gp(j,k) .It.
0.0
)
then
/*
choose sign
alfa
=
-alfa
end if
u(k)
=
-alfa
t
=
gp(j,k)
+
alfa
if( abs(t) .It. test
)
then
/*
div by
O?
ierr=l
return
end if
gp(j,k)
=
t
beta
=
-u(k)
*
gp(j,k)
do
14
jj
=
j+l, mcur
/*
apply transformation to columns
kk
=
igpcjj)
gama
=
0.0
do 12 i=j, n
gama
=
gama
+
gp(i,kk)*gp(i,k)
continue
gama
=
gama
/
beta
do
13
i=j, n
gp(i,kk)
=
gp(i,kk)
-
gama*gp(i,k)
continue
continue
gama
=
0.0
/*
apply transformation to data vector
do 15 i
=
j, n
gama
=
gama
+
dp(i)*gp(i,k)
15 continue
gama
=
gama
/
beta
do
16
i
=
j, n
dp(i)
=
dp(i)
-
gama*gp(i,k)
16
continue
17
continue
18
19
xm(iv)
=
0.0
do
19
i=mcur,l,-1
/*
back-solve for unknowns
t
=
dp(i)
do
18
j=itl,rncur
k
=
igp(j)
/*
location of j-th column of triangle
t
=
t
-
gp(i,k)
*
m(k)
continue
xm(igp(i))
=
t
/
u(igp(i))
continue
return
end
13.4
L,
Problems with Inequality Constraints
249
FORTRAN
Code (MLIC)
for
the Minimum-Length Solution with Inequality
Constraints
subroutine
mlic
(f,xm,h,n,m,nst,mst, fp,test,ierr)
real
f(nst,mst),
xm(mst),
h(nst), test, fp((m+l)*n)
integer n,
m,
nst,
mst,
npst, mpst,
ierr
/*
subroutine
MLIC
finds the minimum
L2
solution length
/*
subject to linear inequality constraints f*xm .ge. h
/*
protocol:
/*
/*
h: (sent, altered) the n by
m
constraint matrix. This matrix
/*
must be dimensioned at least (m+l)*n in
size
/*
/*
xm:
(returned)
a
vector
of
length
m
of
unknowns. Must be
/*
dimensioned
at
least
n in length
/*
/*
h: (sent, altered)
a
vector of length n of constraints. must be
/*
at
least
m+l
in length
/*
/*
n: (sent) the number of
data
(rows of g, currently limited to be
/*
less
than
or
equal to
99)
/*
/*
m:
(sent) the number
of
unknowns (columns of g)
/*
/*
fp: (scratch)
a
scratch vector
at
least
n*(m+l) in length
/*
/*
test: (sent) division by numbers
smaller
than
test
generates
/*
an error condition, solutions within
test
of zero are considered
/*
feasible
/*
/*
ierr: (returned) an error flag, zero on no error,
-1
on
/*
no feasible solution, positive on other errors.
/*
subroutines needed:
ADD,
SUB,
PLS
real
hp(100), e(100)
/*
must be redimensioned
if
n
gt
99
do
2
i=l,
n
/*
build
(m+l)
by n matrix with f-transpose at top
/*
and h-transpose at bottom
j
=
(m+l)*(i-l)
do
1
k=l,
m
fp(j+k)
=
f(i,k)
1
continue
2
continue
fp(j+m+l)
=
h(i)
do
3
i=l,
m
/*
build
m+l
vector hp
hp(i)
=
0.0
3
continue
hp(m+l)
=
1.0
/*
solve positive
least
squares problem fp*xm=hp
call
pls (fp,xm, hp,m+l, n,m+l, n, f, h, e,m+l,n,
test,
ierr)
if
(
ierr .ne.
0
)
then
return
end if
xnorm
=
0.0
250
13
Numerical
Algorithms
Program
MLIC,
continued.
do
5
i=l, m+l
/*
compute residuals and their
L-2
norm
e(i)
=
-hp(i)
do
4
j=1,
n
e(i)
=
e(i)
+
fp((j-l)*(m+l)+i)*xm(j)
4
continue
5
continue
xnorm
=
xnom
+
e
(i)
**2
xnorm
=
sqrt( xnorm
/
float(m+l)
)
6
if( xnorm .It. test
)
then
else
ierr=-1
/*
constraints incompatible
do
6
i=l, m
xm(i)
=
-e(i)/e(m+l)
continue
ierr=O
end if
return
end
13.5
Finding Eigenvalues
and
Eigenvectors
25
1
13.5
Finding the Eigenvalues and
Eigenvectors
of
a Real
Symmetric Matrix
If
the real symmetric matrix
A
is diagonal, then the process of
finding its eigenvalues and eigenvectors is trivial: the eigenvalues are
just the diagonal elements, and the eigenvectors are just unit vectors in
the coordinate directions. If it were possible to find a transformation
that diagonalized a matrix while preserving lengths, then the eigen-
value problem would essentially be solved. The eigenvalues in both
coordinate systems would be the same, and the eigenvectors could
easily be rotated back into the original coordinate system. (In fact,
since they are unit vectors in the transformed coordinate system, they
simply
are
the elements of the inverse transform.)
Unfortunately, there is no known method of finding a finite number
of unitary transformations that accomplish this goal. It is possible,
however, to reduce an arbitrary matrix to either bi-diagonal or tri-diag-
onal form by pre- and postmultiplying by a succession of Householder
transformations. The zeros are placed alternately in the rows and
columns, starting two elements from the main diagonal
(Fig.
13.4).
The tri-diagonal form is symmetric, and therefore simpler to work
with. We shall assume that this triangularization has been performed
on the matrix
A.
At this point there are two possible ways to proceed. One method is
to try to further reduce the tri-diagonal matrix into a diagonal one.
Although this cannot be done with a finite number of transformations,
it can be accomplished with an infinite succession
of
Givens rotations.
In practice, only a finite number of rotations need be applied, since at
some point the off-diagonal elements become sufficiently small to be
ignored. We shall discuss this method in the section on singular-value
decomposition (Section
13.6).
post
L
Fig.
13.4. Tri-diagonalization of a square matrix
by
alternate pre- and post-multiplica-
tion
by
Householder transformations.
252
13
Numerical Algorithms
The second method
is
to find the eigenvalues of the tri-diagonal
matrix, which are the roots of its characteristic polynomial
det(A
-
AI)
=
0.
When
A
is
tri-diagonal, these roots can be found by a recursive
algorithm.
To
derive this algorithm, we consider the
4
X
4 tri-diagonal matrix
’
1
(13.6)
where
u,,
b,,
and
c,
are the three nonzero diagonals. When the determi-
nant
is
computed by expanding in minors about the fourth row, there
are only two nonzero terms:
det(A
-
A,)=
(u4-
A)det[
;
u2;A
h”,
]
-
c4
det
[
(‘2
U2
-
A
0
a,
-A
b2
0
c2
a2
-A
b,
0
c3
a3
-
A
b4
0
0
c,
a,
-
A
uI
-A
h2
0I-A
h2
0
uj
-
i
0
(‘3
h4
(I
3.7)
Because all but one element in the last column of the second matrix are
zero, this expression can be reduced to
det(A
-
U)
=
(a,
-
A)
det
[
u2,
A
h”,
1
-
b4h
det
[
[
A-AAI=
U,
-2
b2
UI
-A
h2
(‘2
u2
-
I
U)
-
I
(1
3.8)
The determinant
of
the
4
X
4
matrix
A
is related to the determinants
of
the
3
X
3
and
2
X
2
submatrices in the upper left-hand comer
of
A.
This relationship can easily be generalized to a matrix
of
any size as
PN
=
(‘N
-
A)PN-
I
-
bNcNpN-2
(1
3.9)
where
pN
is the determinant
of
the
N
X
N
submatrix. By defining
p-l
=
0
and
p,,
=
1,
this equation provides a recursive method for
computing the value of the characteristic polynomial for any specific
value
of
A.
This sequence
of
polynomials is a
Sturm
sequence,
which
13.5 Finding Eigenvalues and Eigenvectors 253
has two properties that facilitate the locating of the roots:
(I)
the roots
of
PNare separated by the roots
of
PN-I ; (2) the
number
of roots of
PN(A) less than A
o
is equal to the
number
of
sign changes in the
sequence
[Po(A
o),
PI(A
o),
. . . , PN(A
o)].
The first property is proved in
most elementary numerical analysis texts.
The
second follows from
the fact that for large positive
A,PNandPN-I have opposite signs, since
PN= (aN- A)PN-l (Fig. 13.5).
Bounds on the locations
of
the roots
of
this polynomial can be found
by computing the norm
of
the eigenvalue equation
AX
=
Ax:
IAlllxll
2
= IIAxll
2
~
IIAII
211xll
2
(13.10)
The
roots must lie on the interval
[-IIAI1
2
,
IIAII
2
] . If the Sturm se-
quence is evaluated for some
Ao
on this interval, then the number of
roots in the subintervals
[-IIAII
2
,
A
O
]
and
[A
o
,
IIA11
2
]
can be determined
by counting the sign changes in the Sturm sequence. By successive
subdivision
of
intervals the position
of
the roots (the eigenvalues) can
be determined to any desired accuracy. This method
of
solving the
eigenvalue problem is known as the
Givens-
Householder algorithm.
Once the eigenvalues have been determined, the eigenvectors x
can be computed by solving the equation (A -
AiI)x
i
= 0 for each
eigenvalue. These equations are underdetermined, since only the
directions and not the overall magnitudes
of
the eigenvectors are
determined. If
Gauss-Jordan
reduction isused to solve the equations,
then some pivoting algorithm must be employed to avoid division
by near-zero numbers. Another possibility is to solve the equation
[A -
(Ai+
<5)I]x{
= x{-l iteratively, where the initial guess
of
the
eigenvector x? is any random vector and
<5
-e: Aisome small number,
PI
(A)
I
+
• A
•
I
P
2(A)
I
+ +
• A
•
I
P3
(A)
I
+
+
• A
•
P4
(A)
J
+ + +
• A
I
I
-IIAII
2
+IIAII
2
Fig. 13.5. Roots (dots) of the Sturm ploynomials p.().). Note that the roots of
p._
1
separate (interleave with) the roots of pn :
254
13
Numerical Algorithms
to prevent the equations from becoming completely singular. This
iteration typically converges to the eigenvector in only three iterations.
Once the eigenvectors have been determined, they are transformed
back into the original coordinate system by inverting the sequence of
Householder transformations.
13.6
The Singular-Value Decomposition
of
a
Matrix
The first step in this decomposition process is to apply Householder
transformations to reduce the matrix
G
to bi-diagonal form (in the
4
X
4
case) (see Section 13.5):
pl
"2
0
0
1
(13.1
1)
LO
0
0
44-I
where
4,
and
0,
are the nonzero diagonals. The problem is essentially
solved, since the squares ofthe singular values of
G
are the eigenvalues
of
GTG,
which is tri-diagonal with diagonal elements
a,
=
q'
and
off-diagonal elements
b,
=
c,
=
e,q,-,
.
We could therefore proceed by
the symmetric matrix algorithm of Section 13.5. However, since this
algorithm would
be
computing the squares of the singular values,
it
would not
be
very precise. Instead, it is appropriate to deal with the
matrix
G
directly, using the
Golith-
Reinch
mcrhod.
We shall examine the computation of the eigenvalues
of
a
tri-diago-
nal matrix by transformation into
a
diagonal matrix first. Suppose that
the symmetric, tri-diagonal
N
X
N
matrix
A
is transformed by a
unitary transformation into another symmetric tri-diagonal matrix
A'
=
TATT.
It is possible
to
choose this transformation
so
that the
h;
=
ch
elements become smaller in absolute magnitude. Repeated
application of this procedure eventually reduces these elements to
zero. leaving
a
matrix consisting of an
N
-
I
X
N
-
1
tri-diagonal
block and
a
1
X
I
block. The latter is its own eigenvalue. The whole
process can then be repeated on the
N
-
1
X
N
-
I
matrix to yield an
N
-
2
X
N
-
2
block and another eigenvalue. The transformation
A'
=
?'ATT
is accomplished by a coordinate shift. For any
sh$ppuram-
c"cr
a,
A'
=
TATT
=
T[A
-
aI]TT
+
aI.
It can be shown that
if
a
is
13.6
The Singular-Value Decomposition
of
a Matrix
255
chosen to be the eigenvalue of the lower-right
2
X
2
submatrix of
A
that is closest in value to
A,,
and
if
T
is chosen to make
T[A
-
a11
upper triangular, then
A’
is triangular and the
ch
=
6;
elements of
A’
are smaller than the corresponding elements of
A.
The singular values of
G
are found by using the Householder-
Reinch algorithm to diagonalize
A
=
GTG.
All of the operations are
performed on the matrix
G,
and
A
is never explicitly computed. The
shift parameter is first computed by finding the eigenvalues of the
2
X
2
matrix
and selecting the eigenvalue closest to
(4;
+
pi).
A
transformation
must then be found that makes
T(A
-
crT)
upper-triangular. It can be
shown [Ref. 141 that if
A
and
TATT
are symmetric, tri-diagonal
matrices and
if
all but the first element of the first column of
T(A
-
01)
are zero, then
T(A
-
01)
is upper triangular. The effect
T
has on
G
can
be controlled
by
introducing another unitary transformation
R
and
noting that if
G’
=
RTGTT,
then
A’
=
GT’G’
=
TGTRRTGTT
=
TGTGTT
is not a function of
R.
The transformations
T
and
R
are constructed out of a series
of
Givens transformations. These transformations are a special case of
Householder transformations designed to change only two elements of
a vector, one of which is annihilated. Suppose that the elements
u,
and
uJ
are to be altered,
uJ
being annihilated. Then the transformation
TI,
=
q,
=
c,
7;,
=
-
=
s
(all other elements being zero) will
transform
v:
=
(of
+
$)I/*,
u(
=
0
if
c
=
~,/(uf
+
uf!)II2
and
s
=
vJ/(vf
+
u;)’/~.
The first Givens transformation in
T
is chosen to
annihilate the subdiagonal element of the first column
ofT(A
-
GI),
[4:
-
(T.
4,c,,
0,
0,
.
. .
.
OIT.
Now
GT:
is no longer bi-diagonal,
so
a
series of pre- and postmultiplied Givens transformations are applied
which bi-diagonalize it. These transformations define the rest of
’I’
and
R.
Their application is sometimes called
cliusing
since. after each
application
of
a transformation, a nonzero element appears at a new
(off-bi-diagonal) position in
G.
This step assures that
TATT
is tri-diag-
onal. The whole process is iterated until the last subdiagonal element
of
TATT
becomes sufficiently small to be ignored. At this stage
ch
must also be
zero.
so
one singular value
4,
has been found. The whole
process is then repeated on a
G
the size of which has been reduced
by
one row and one column.
256
13
Numerical
Algorithms
13.7
The Simplex Method and the Linear
Programming Problem
In
the linear programming problem, we must determine a vector
x
that maximizes the
ohjecfive
function
z
=
cTx
with the equality and
inequality constraints
Ax[<,
=,
z]b
and with the nonnegativity con-
straint
x
2
0.
As
discussed in Section
7.9,
the constraints (presuming
they are consistent) define a convex polyhedron in solution space,
within which the solution must lie. Since the objective function is
linear, it can have no maxima or minima. The solution must therefore
lie on the boundary of the polyhedron, at one of its vertices (or
“extreme points”).
In
many instances there is a unique solution, but
there may
be
no solution ifthe constraints are inconsistent, an infinite
solution
if
the polyhedron is not closed in the direction
of
increasing
z,
or an infinite number of finite solutions
if
the solution lies on an edge
or face of the polyhedron that is normal to
Dz
(Fig.
13.6).
The Simplex method is typically used to solve the linear program-
ming problem.
It
has two main steps. The first step (which, for
historical reasons, is usually called the “second phase” consists of loca-
ting one extreme point on the polyhedron defined by
Ax
[s,
=,
2j
b
that is also feasible in the sense that it satisfies
x
2
0.
The second step
(the “first phase”) consists of jumping from this extreme point to
another neighboring one that is also feasible and that has a larger value
of the objective function. When no such extreme point can be found,
the solution of the linear programming problem has been found.
There is no guarantee that this search procedure will eventually reach
m2
kq
mkt,
&,
m/b
Fig.
13.6.
Depending
on
the objective function (dashed contours) and the shape
of
the
feasible volume defined by the inequality constraints (solid lines with arrows pointing
into feasible halfspace), linear programming problems can have (a) no solution,
(b)
one
solution (triangle), (c) an infinite number of finite solutions (on line segment bounded
by triangles),
or (d) an infinite solution.
u
m,
ml
m,
13.7
The Simplex Method and the Linear Programming Problem
257
the “optimum” extreme point (closed loops are possible), but it
usually does.
The inequality constraints
Ax
[s,
=,
21
b
can be converted to
equality constraints by changing all the
2
constraints to
5
constraints
by multiplication by
-
1,
and then changing all the
5
constraints to
equality constraints by increasing the number of unknowns. For
example, if the problem originally has two unknowns and a constraint
x,
+
2x2
I
4,
it can be converted into the equality constraint
x,
+
2x2
+
x3
=
4
by adding the unknown
x3
and requiring that
x3
2
0.
The
transformed linear programming problem therefore has, say,
M
un-
knowns (original ones plus the new “slack” variables)
x,
objective
function
z
=
cTx
(where
c
has been extended with zeros), and, say,
N
linear equality constraints
Ax
=
b.
Since one slack variable was added
for each inequality constraint, these equations will usually be under-
determined
(M
>
N).
If
the ith column
of
A
is denoted by the vector
vi,
the equations
Ax
=
b
can be rewritten as
(13.13)
The vector
b
is of length
N.
It can be represented by a linear combina-
tion of only
N
<
M
of the
vi.
We can find one solution to the constraint
equation by choosing
N
linearly independent columns
of
A,
expand-
ing
b
on this basis, and noting that
N
components of the solution
vector are the expansion coefficients (the rest being zero):
Any solution that can be represented by a sum of precisely Mcolumns
of
A
is
said to be a basic
solution,
and the
N
vectors
vi
are said to form
the basis of the problem. If all the
x,’s
are also nonnegative, the
problem is then said to be basicTeasible.
The key observations used by the Simplex algorithm are that all
basic-feasible solutions are extreme points of the polyhedron and that
all extreme points of the polyhedron are basic-feasible solutions. The
first assertion is proved by contradiction. Assume that
x
is both
basic-feasible and
not
an extreme point. If
x
is not an extreme point,
then it lies along a line segment bounded by two distinct vectors
xI
and
x2
that are feasible points (that is, are within the polyhedron):
x=Ax,
+(I
-A)x2
where
O<A<
1
(1
3.15)
Since
x
is basic, all but Nof its elements are zero. Let the unknowns
be