Analysis and Computation of Fixed Points. Proceedings of a Symposium Conducted by the Mathematics Research Center, the University of Wisconsin–Madison, May 7–8, 1979

414 Pages·1980·14.719 MB·English

by  Stephen M. Robinson

ANALYSIS
AND COMPUTARON
OF FIXED POINTS
Edited by
Stephen M. Robinson
Mathematics Research Center
University of Wisconsin—Madison
Madison, Wisconsin
Proceedings of a Symposium
Conducted by the Mathematics Research Center
The University of Wisconsin—Madison
May 7-8, 1979
ACADEMIC PRESS 1980
A Subsidiary of Harcourt Brace Jovanovich, Publishers
NEW YORK LONDON TORONTO SYDNEY SAN FRANCISCO
COPYRIGHT © 1980, BY ACADEMIC PRESS, INC.
ALL RIGHTS RESERVED.
NO PART OF THIS PUBLICATION MAY BE REPRODUCED OR
TRANSMITTED IN ANY FORM OR BY ANY MEANS, ELECTRONIC
OR MECHANICAL, INCLUDING PHOTOCOPY, RECORDING, OR ANY
INFORMATION STORAGE AND RETRIEVAL SYSTEM, WITHOUT
PERMISSION IN WRITING FROM THE PUBLISHER.
ACADEMIC PRESS, INC.
Ill Fifth Avenue, New York, New York 10003
United Kingdom Edition published by
ACADEMIC PRESS, INC. (LONDON) LTD.
24/28 Oval Road, London NW1 7DX
Library of Congress Cataloging in Publication Data
Symposium on Analysis and Computation of Fixed
Points, University of Wisconsin—Madison, 1979.
Analysis and computation of fixed points.
(Publication of the Mathematics Research Center,
University of Wisconsin—Madison ; no. 43)
"Proceedings of a symposium conducted by the
Mathematics Research Center, The University of
Wisconsin—Madison, May 7-8, 1979."
Includes index
1. Algebraic topology—Congresses. 2. Fixed
point theory—Congresses. I. Robinson, Stephen M.
II. Wisconsin. University—Madison. Mathematics
Research Center. III. Title. IV. Series:
Wisconsin. University—Madison. Mathematics Re­
search Center. Publication ; no. 43.
QA3.U45 no. 43 [QA612] 510s [514' .2] 80-20545
ISBN 0-12-590240-9
PRINTED IN THE UNITED STATES OF AMERICA
80 81 82 83 9 8 7 6 5 4 3 2 1
Contributors
Numbers in parentheses indicate the pages on which the authors' contributions begin.
Charles R. Engles (213) Booz-Allen and Hamilton, Inc., 555 California Street,
San Francisco, California 94104
Hartmut Jürgens (139) Fachbereich Mathematik, Forschungsschwerpunkt
"Dynamische Systeme," Universität Bremen, 2800 Bremen 33, West
Germany
Masakazu Kojima (93) Department of Information Sciences, Tokyo Institute
of Technology, Meguro, Tokyo 152, Japan
Tien-Yien Li (73) Department of Mathematics, Michigan State University,
East Lansing, Michigan 48824
James G. MacKinnon (197) Department of Economics, Queen's University,
Kingston, Ontario, Canada K7L 3N6
Heinz-Otto Peitgen (139) Fachbereich Mathematik, Forschungsschwerpunkt
"Dynamische Systeme," Universität Bremen, 2800 Bremen 33, West
Germany
John Piggott (183) Department of Economics, University of Western Ontario,
London, Ontario, Canada N6A 5C2
D. G. Saari (57) Department of Mathematics, Northwestern University,
Evanston, Illinois 60201
Romesh Saigal (57) Department of Industrial Engineering and Management
Sciences, Northwestern University, Evanston, Illinois 60201
Dietmar Saupe (139) Fachbereich Mathematik, Forschungsschwerpunkt
"Dynamische Systeme," Universität Bremen, 2800 Bremen 33, West
Germany
Shlomo Shamir (25) Department of Engineering—Economic Systems,
Stanford University, Stanford, California 94305
Michael J. Todd (1) School of Operations Research and Industrial Engineer­
ing, Cornell University, Ithaca, New York 14853
John Whalley (183) Department of Economics, University of Western Ontario,
London, Ontario, Canada N6A 5C2
James A. Yorke (73) Institute for Physical Science and Technology, Uni­
versity of Maryland, College Park, Maryland 20742
vii
Preface
The papers published in this book arose primarily from the Symposium on
Analysis and Computation of Fixed Points, held at the University of Wisconsin
—Madison on May 7-8, 1979. The symposium was organized and held under
the auspices of the Mathematics Research Center and was sponsored by the
National Science Foundation under Grant No. MCS-7907217 and by the United
States Army under Contract No. DAAG29-75-C-0024. Of the thirteen papers
presented there, eight are collected in this volume; the ninth paper published
here is the doctoral dissertation of Charles R. Engles, written at Stanford
University under the supervision of Professor B. Curtis Eaves. It is very much in
the spirit of the other papers, and I hope that its inclusion may help to make its
contents available to a wider group of research workers than they would
otherwise reach.
I should like to thank Mrs. Judith Siesen, for extremely capable editorial
assistance; Mrs. Gladys Moran, for excellent work in organizing the sympos­
ium; and the agencies cited above, for their support of the preparation of this
volume as well as of the symposium itself.
Stephen M. Robinson
Numerical Stability and Sparsity
in Piecewise-Linear Algorithms
Michael J. Todd
1. Introduction.
This paper demonstrates how the techniques of numerical
linear algebra that ./possess attractive stability properties
and exploit sparsity can be applied to the linear systems
that arise in algorithms that solve equations by constructing
piecewise-linear homotopies. Such algorithms are usually
called fixed-point or simplicial algorithms; however, they
are most simply described as finding zeroes rather than fixed
points, and can traverse pieces larger than simplices--hence
we shall call them piecewise-linear (PL) algorithms. General
references on these methods and their applications are Scarf
[19], Todd [20], Eaves and Scarf [6], Eaves [4] and Allgower
and Georg [1]; the latter two have comprehensive bibliogra­
phies. We shall be concerned only with the restart algorithm
of Merrill [14] but our ideas apply to most PL algorithms.
Suppose that we seek a zero of a continuous function f
from R to itself. (PL algorithms are also applicable to
the problem of approximating zeroes of upper semi-continuous
point-to-set mappings, but we will confine ourselves to the
case of functions.) We then choose a one-to-one affine func­
tion f : R -> R defined by f (x) = Gx-g and construct
the homotopy h: Rn χ [0,1] -> Rn by η(χ,θ) = 6f(x) +
(l-6)f (x). Next we make a piecewise-linear approximation £
to h; generally this is done by choosing a triangulation T
of R x [0,1], letrting £ agree with h on the vertices
ANALYSIS AND COMPUTATION OF FIXED POINTS Λ Copyright © 1980 by Academic Press, Inc.
* All rights of reproduction in any form reserved.
ISBN 0-12-590240-9
2 MICHAEL J. TODD
of T and then extending linearly on the simplices. However
it is important to note that the pieces of linearity of £
are usually much larger than simplices, even if f has no
special structure; this conglomeration of pieces can be en­
hanced if f has special structure by an appropriate choice
of triangulation T [23,24]. Then a PL path of zeroes of I
is traced from the known zero (x ,0) where x = G g.
This path either diverges to infinity or produces a zero
(x ,1) of I—x is then an approximate zero of f. The
step from x to x is one major iteration of the restart
algorithm of Merrill [14]; if desired it can then be repeated
with x replacing x and a new affine function replacing
f°.
Within each major iteration the PL algorithm generates a
sequence of pieces. In each minor iteration, some partial
information of f (or f ) is obtained, usually the vector
f(v) or f (v) for a vertex v, and a linear algebra com­
putation is performed to determine the new piece of linearity.
In the simplicial form of the algorithm the linear algebra
amounts to a linear programming pivot, but when the large
pieces are fully exploited it becomes more complicated [24].
The proof that the algorithm does not cycle relies heavily on
the uniqueness of the linear programming pivots to maintain
feasibility, and lexicographic rules are used to avoid prob­
lems of degeneracy. It is clear that a buildup of roundoff
error can easily lead to an incorrect pivot and possible cy­
cling—the situation is mere serious than in linear program­
ming since there is no objective function to "pull the iter­
ates away from a cycle" and there is no simple remedy of a
reinversion and a "crash" start.
I am unaware of any problem to date caused by cycling in
PL algorithms. (A possible case arose with the authorTs code
on one of Scarf's economic problems. The algorithm at first
failed to terminate in the limit of 7000 pivots. After ad­
justing a tolerance in the minimum ratio testing, convergence
in less than 4000 pivots was obtained. If nothing else, this
example shows the drastic effects of incorrect pivots.) How­
ever, all codes of which I have knowledge maintain and update
the basis inverse explicitly—it is known that this method in
PIECEWISE-LINEAR ALGORITHMS 3
linear programming lacks numerical stability. In this paper
we suggest the use of factorizations of the basis matrix that
do maintain numerical stability, in particular the QR and LU
decompositions.
If the only operations on the basis performed were the
straightforward column exchanges of linear programming, the
application of these factorizations would be merely an exer­
cise. There are three special properties of PL algorithms,
however, which require other operations on the basis inverse.
We show that these operations can also be performed when fac­
torizations of the basis are maintained and updated.
First, when starting with f (x) = x-c or c-x as is
usual, the initial basis has a great deal of structure en­
abling one to obtain its inverse explicitly. We show that
factorizations of this initial basis are just as straightfor-
ward. If we start with f (x) = A(x-c) and A- is known,
the initial basis inverse is again immediate; similarly we
show that if a factorization of A is known, a factorization
of the initial basis can be written down.
Second, the linear algebra computation to determine
which piece of linearity is to be entered next is simple when
the basis inverse is maintained. This computation amounts to
a column exchange together with some possible pre- and post­
processing of the basis inverse, consisting of changing signs
of rows and adding or subtracting a row from another. We
demonstrate that if the factorization is maintained and up­
dated correctly, the pre- and post-processing can be per­
formed as simply on a factorization of the basis.
Finally, when a major iteration is completed, PL algor­
ithms generate the inverse of an approximate Jacobian of f,
if the basis inverse is maintained.. This matrix can be used
for a variety of purposes: estimating the distance to the
solution and hence updating the grid size [17]; taking parall­
el chord or quasi-Newton steps [21,22]; and generating a new
affine function to be used in a restart [22]. We show that a
factorization of the final basis inverse can also be used for
all of these tasks.
Section 2 describes the factorizations we shall be con­
cerned with and how they can be obtained and updated. In
4 MICHAEL J. TODD
section 3 we present techniques for handling these factoriza­
tions to deal with the special operations required in PL al­
gorithms— this is the heart of the paper. Section 4 discusses
the application o'f the ideas of section 3 when the function
f is separable. Finally, in section 5 we consider the case
where f is sparse. We argue that a slightly different basis
matrix should be factorized in this case. This matrix is
sparse and permits the algorithm to traverse large pieces as
it finds them; the techniques here illustrate the combination
of the local strategy of [23] with the global strategy of
[24], The price paid for the sparsity of the matrix is that
two columns may have to be exchanged instead of just one, at
each minor iteration.
2. The QR and LU Factorizations
Let A be an mxm nonsingular matrix, with columns
a ,a ,...,a . We discuss here two well-known factorizations
that (i) allow easy solution of systems Ax = b and (ii) can
be readily updated when a column a is replaced by a new
column a. The use of these methods in mathematical program­
ming was proposed by Bartels, Golub and Saunders [l1].
2.1 The QR decomposition
The matrix A can be written as QR where Q is an or­
thogonal and R an upper triangular matrix of order m. The
T
matrix Q is not stored. To solve Ax = b, we form c = A b
T
and solve R y = c and Rx = y ; the number of operations
2
(multiplications and divisions) required is 2n + 0(n).
The initial factorization of A can be obtained by mul­
tiplying A on the left by elementary reflectors (or House-
T
holder matrices) I - 2uu with ||u|| =1 or by plane
rotations (Givens matrices) of the form
Γι Ί
1
cos θ sin θ
i
1
-sin θ cos θ
1
L Ί J
PIECEWISE-LINEAR ALGORITHMS 5
to create the desired zeroes. If A is sparse, the rows and
columns of A can be first ordered to make it nearly upper
triangular, so that fill-in is minimized—see [7,10,15]. In
our application we will see that the initial basis can fre­
quently be reordered so that it is upper triangular or nearly
so.
Now suppose we wish to replace column a by a new
column a and obtain the upper triangular factor R of the
new basis. We will usually already have solved the system
T
Ax = a and thus have available the vector y = Q a (i.e.,
T T
R y = A a). We can thus reorder the columns of the new basis
to write it as
-r- 1 2 p-1 p + 1 m -,
A = [ra ,a ,...,a^ ,a^ ,...,a ,a]
^ 1 2 p-1 p + 1 m _ ~
= Qr[r ,r ,...,r^ ,r^ ,...,r ,yΊ] = QτHτ
where r ,...,rm are the columns of R and H is an upper
Hessenberg matrix. In our application we will usually want
to insert the new column a in the penultimate position;
note that this also yields a Hessenberg matrix. Sometimes
the new column will be inserted in position m-2 —then we
obtain a Hessenberg matrix with an extra nonzero in position
(m, m-2). In any case we will denote the resulting matrix
H.
Plane rotations are then used to reduce H to upper tri­
angular form, yielding the required matrix R. If performed
in the natural way, this computation requires n-p square
2
roots and 2(n-p) + 0(n) multiplications, for an average of
2
n/2 square roots and 2/3 n + 0(n) multiplications if all
values of p are equally likely. It is possible to rear­
range the computations to avoid the square roots and halve
the number of muliplications, but this method may be numeri­
cally unstable—see for instance [8,25].
In the sparse case, A may have had its columns permuted
to encourage a sparse matrix R. It is possible that sparsity
of R will be enhanced by placing column a in the pth
position in A rather than in the last or penultimate posi­
tion. In this case H has a spike in position p. H is
then reduced to upper triangular form by plane rotations—in
6 MICHAEL J. TODD
the worst case, twice as many square roots and multiplica­
tions are required, but sparsity may reduce this considerably.
There are a number of related decompositions. First, we
have avoided storing Q, but if it is stored the solution of
Ax = b is slightly easier and more stable. On the other
hand, Q must be updated and this adds considerably to the
computational burden of a column exchange. Gill, Murray and
Saunders [9] have suggested decomposing A as VDR where
1/2
R is unit upper triangular, D is diagonal and VD is
T
orthogonal. V or V is then maintained in product form.
This procedure alleviates the cost of updating and if A is
sparse, the products in V may be sparse even if V is not.
Finally, the related decompositions A = LQ or A = LDV
where L is lower triangular are not used because they do
not mesh well with the pre- and post-processing we perform.
2.2 The LU Decomposition
Now suppose A is factorized as LU, with U upper
triangular and L stored implicitly as a product
E, P, E. P. , . . . E-. P., . Here each P. isa permutation matrix
k k k- 1Ί k-1 11 i ^-
and each E an elementary lower triangular matrix (E
differs from the identity in just one column). To solve
Ax = b we form c = E (P (E _ (P _ ( . . .E^^CP b) ) . . . ) and
k k k 1 k 1
then solve Ux = c. If each E differs from the identity
2
in just one entry, 1/2 n + k multiplications are required.
The initial factorization of A is obtained by Gaussian
elimination with partial pivoting. That is, A is first
premultiplied by a permutation matrix P so that its (1,1)
entry is largest in absolute value in its first column. Then
premultiplying by an elementary matrix E-. whose first col­
umn differs from the identity gives a matrix whose first
column is zero below the diagonal. Continuing for columns
2,...,m-l gives an LU representation of A with L
E ., P -....E-.P.. and U, the final transformation of A,
m-1 m-1 1 1 ' '
upper triangular. If A is sparse, its rows and columns can
be permuted to make it nearly upper triangular. In addition,
the matrices P. can be chosen at least partially on the
basis of sparsity—instead of placing the largest subdiagonal