Table Of ContentSolution of Linear Programming
and Non-Linear Regression Problems
Using Linear M-Estimation Methods
Ove Edlund
Department of Mathematics,
Lule(cid:23)a University of Technology,
S-971 87 Lule(cid:23)a, Sweden
July 1999
1
2
Abstract
This Ph.D. thesis is devoted to algorithms for two optimization problems, and
their implementation. The algorithmsarebasedonsolvinglinear M-estimation
problems.
First, an algorithm for the non-linear M-estimation problem is considered.
The main idea of the algorithm is to linearize the residual function in each
iteration and thus calculate the iteration step by solving a linear M-estimation
problem. A 2-norm bound on the variables restricts the step size, to guarantee
convergence.
The other algorithm solves the linear programming problem. If a variable
in the primal problem has both a lower and an upper bound, it gives rise to
an edge in the dual objective function. This edge is \smoothed" by replacing
it and its neighbourhood with a quadratic function, thus making it possible to
solvetheapproximatedproblemwithNewton’smethod. Forvariableswithonly
lower or only upper bounds, a quadratic penalty function is used on the dual
problem. In this way also variables with one-sided bounds can be handled. A
crucialpropertyofthealgorithmisthatoncetherightactivecolumnsetforthe
optimal solution is identi(cid:12)ed, the optimal solution is found in one step.
The implementation uses sparse matrix techniques. Since it is an active set
method, it is possible to reuse the old factor when calculating the new step.
This is accomplished by up- and downdating the old factor, thus saving much
computation time. It is only occasionally, when the downdating fails, that the
factor instead has to be found with a sparse multifrontal LQ-factorization.
3
Acknowledgements
First of all, I would like to thank my supervisor, associate professor H(cid:23)akan
Ekblom. One way or another, he managed to get me a position as a Ph.D.
student in numerical analysis at a time when the university did not recognize
this as a research area in its own right. Indeed, professor Lars-Erik Persson,
GerdBrandellandthedepartmentofmathematicsalsowereveryhelpfulinthis
respect. Things have improved since then, as the number of people involved in
this research area has grown. Also, the name of the research area has changed
from numerical analysis to scienti(cid:12)c computing.
During my studies, H(cid:23)akan has always been very generous with his time,
listening to whatever I had in mind, reading my manuscripts, and giving much
help in return.
I would also like to express my gratitude to professor Kaj Madsen and as-
sociate professor Hans Bruun Nielsen at the department of mathematical mod-
elling, DTU, Denmark. They have had a substantial input in all the articles
in the thesis and I have learned a lot by the discussions we have had. Nordisk
Forskerutdanningsakademi(NorFA)has(cid:12)nanciallysupportedmyvisitstoDen-
mark, making this collaborationpossible.
Finally, I would like to thank my wife, Susanne. The work on the thesis
has taken more time than I anticipated, but still Susanne has put up with my
late nights and done more than her share of taking careof our childrenand the
household. Perhaps there will be a time when I can pay her back.
4
Contents
1 Introduction 7
2 Linear M-Estimation 8
3 Non-Linear M-Estimation 11
4 Using a Modi(cid:12)ed Huber-criterion for Solving Linear Program-
ming Problems 12
5 The Concept of Sparse Matrices 14
Article I: Algorithms for Non-Linear M-Estimation 17
Article II: Linear M-Estimation with Bounded Variables 31
Article III: A Piecewise Quadratic Approach for Solving Sparse
Linear Programming Problems 51
Article IV: A Software Package for Sparse Orthogonal Factoriza-
tion and Updating 79
5
6
1 Introduction
ThisPh.D.thesisinscienti(cid:12)ccomputingisdevotedtothesolutionoftwodi(cid:11)er-
ent optimization problems. The non-linear M-estimation problem comes from
the(cid:12)eldofrobustregression. Theobjectiveistoidentifyparametersinamath-
ematicalmodelsothatthemodelandreallifeobservationsmatch. Themeaning
of \robust" is that a few erroneous observations should not alter the solution
in a signi(cid:12)cant way. The other problem considered is the linear programming
problem,thatwasformulatedandsolvedbyDanziginthelate1940:s,andsince
then has found its use in many applications, including such di(cid:11)erent areas as
integrated circuit design, chemistry and economics.
Thoughtheseproblemsareverydi(cid:11)erent,thetwoalgorithmssharetheprop-
erty that the solution is found by solving a sequence of linear M-estimation
problems.
Linear programming problems are frequently formulated with sparse ma-
trices, and much computational work can be saved by taking the sparsity into
account. Thelinearprogrammingalgorithmhasbeenimplementedusingsparse
matrixtechniques. The sparsematrixpartofthe implementationisthe subject
of one of the articles in this thesis.
The thesis consists of four articles. The (cid:12)rst two deal with solving the
non-linearM-estimationproblem,andthefollowingtwodescribethelinearpro-
gramming algorithm and implementation:
Article I: Edlund, O., Ekblom,H. & Madsen,K.(1997),‘Algorithms for non-
linear M-estimation’, Computational Statistics 12, 373{383.
Article II: Edlund, O. (1997), ‘Linear M-estimation with bounded variables’,
BIT 37(1), 13{23.
Article III: Edlund, O., Madsen K. & Nielsen H. B. (1999), ‘A Piecewise
Quadratic Approach for Solving Sparse Linear Programming Problems’.
To be submitted.
Article IV: Edlund, O. (1999), ‘A Software Package for Sparse Orthogonal
Factorization and Updating’. To be submitted.
Beside these articles, the author has published two papers at conferences.
One is a short note on solving linear Huber problems using a vector processor
(Edlund 1994). The other paper is on algorithms for robust error-in-variables
problems (Ekblom & Edlund 1998). This second paper describes a work in
progress. Unfortunately it did not evolve far enough to make it to the thesis.
The articles are preceded by a general introduction where the basic ideas
of the algorithms are introduced. Section 2 describes briefly the concept of
linear M-estimation, the core in both of the algorithms. The case of nonlinear
M-estimation is the subject of section 3. This serves as an introduction to
Article I{II. In section 4 the linear programming algorithm of Article III is
introduced. Section 5 gives a brief introduction to the area of sparse matrices,
as a prelude to Article IV.
7
0.6
0.5
(t, b)
3 3
0.4
(t, b)
4 4 (t, b) (t, b)
(t, b) 8 8 9 9
2 2
(t, b)
b0.3 (t5, b5) 6 6
0.2
0.1
(t, b)
7 7
0 (t1, b1)
0 1 2 3 4 5 6 7 8
t
Figure 1: The observations, and some randomly chosen model functions.
2 Linear M-Estimation
TheM-estimatesariseinrobuststatistics. Theirpurposeistomakeastatistical
modelthatisnotsensitivetogrosserrorsinthedata. Theterm\M-estimate"is
tobeinterpretedas\maximumlikelihoodtypeestimate",andisjusti(cid:12)edbythe
fact that the de(cid:12)nition of an M-estimate is somewhat similar to the maximum
likelihood problem (Huber 1981).
ToexplaintheprinciplesoflinearM-estimation,wewilllookatanexample.
Suppose that we have some observations (t ;b ) and there is a mathematical
i i
model that describes the relation between t and b with the equation
b=x t+x te−t (1)
1 2
where x and x are unknown. In Figure 1 the observations are plotted, as
1 2
well as b as a function of t for some random [x x ]. The objective is to (cid:12)nd a
1 2
vector [x x ] such that the equation is as close to the observations as possible.
1 2
As seen in the (cid:12)gure, the observation (t ;b ) is very di(cid:11)erent compared to the
7 7
others and should be regardedas an erroneous observation.
We introduce a measurement of the distance from each observation to the
function, called the residual as shown in Figure 2. The residual is a vector
containing as many elements as the number of observations. For our example
an element in the residual is de(cid:12)ned as
ri =bi−x1ti−x2tie−ti;
8
0.6
r
0.5 3
0.4 r r
2 4
0.3
b r
5
r r r
6 8 9
0.2
0.1
r
7
0 r
1
0 1 2 3 4 5 6 7
t
Figure 2: The elements of the residual.
thus the residual vector can be calculated as
2 3 2 3 2 3
66 rr12 77 66 bb12 77 66 tt12 tt12ee−−tt12 77(cid:20) x (cid:21)
r =6 . 7=6 . 7−6 . . 7 1 ;
4 . 5 4 . 5 4 . . 5 x
. . . . 2
r9 b9 | t9 {t9ze−t9 }
=A
i.e. r =b−Ax where A is a constantmatrix and b is a constantvector,and x
is the unknown vector that is sought. The reason why we could formulate the
residualwiththehelpofamatrixisthattheelementsofxarelinearin(1). The
mathematical model is often disregarded in the linear case, since it is su(cid:14)cient
toknowthe matrixAandthevectorbtosolvethe problem. Themoredi(cid:14)cult
non-linear case is introduced in section 3 and is described in more detail in
Article I.
Looking at Figure 2, it is not di(cid:14)cult to imagine that the model equation
is close to the observations if the elements of the residual are small, so our
objective is to minimize the residual. This is done by constructing an objective
function that sums up the residual entries in a special way:
Xm
G(x)= (cid:26)(r (x)):
i
i=1
Here(cid:26)(t)isapositivefunctionthatisdecreasingast<0andincreasingast>0.
By(cid:12)nding the solutionofmin G(x), the residualis minimized. This minimum
x
9
Table 1: (cid:26)-functions for some M-estimates.
‘ estimate
1
(cid:26)(t)=jtj
Huber estimate (k>0)
t2=2; jtj(cid:20)k
(cid:26)(t)=
kjtj−k2=2; jtj>k
(cid:26)
Fair estimate (k>0)
(cid:26)(t)=k2(jtj=k−ln(1+jtj=k))
Tukey estimate (k>0)
(k2=6)[1−f1−(t=k)2g3]; jtj(cid:20)k
(cid:26)(t)=
k2=6; jtj>k
(cid:26)
Welsh estimate (k>0)
(cid:26)(t)=k2=2 1−e−(t=k)2
(cid:16) (cid:17)
varies though with the choice of (cid:26). The most common choice is (cid:26)(t)=t2 which
gives the least squares solution. The least squares solution is easy to (cid:12)nd by
solving the normal equations, though other methods have better stability from
anumericalpointofview. Theleastsquaressolutionisthemaximumlikelihood
solutionifthe errorsarenormallydistributed, but it is notverygoodathandle
erroneous observations, since there is a high penalty on large residual entries
duetothesquareinthesum. Somerobust choicesof(cid:26),thatarelesssensitiveto
erroneous observations, are shown in Table 1. Solutions for di(cid:11)erent choices of
(cid:26)arecalledM-estimates. Manyofthe(cid:26)-functionshaveaconstantparameterk,
so the solution found does not only vary with the choice of(cid:26), but also with the
choice of k.
Using some di(cid:11)erent (cid:26)-functions to solve our example give the solutions in
Figure 3. From the (cid:12)gure it is obvious that the Huber and Welsh estimates are
less disturbed by the dissentient observation than the least squares estimate.
The choice of k should reflect a threshold between good residual elements
and bad ones corresponding to erroneous observations. To be able to do this
choice in a consistent way, the residual entries are often rescaled by a factor (cid:27),
giving the following problem
Xm
minimize G(x)= (cid:26)(r (x)=(cid:27)): (2)
i
i=1
Finding the solution of (2) then comes in two flavours. Either the scale (cid:27)
is known, or it is not. In the second case (cid:27) can be found in the process of
solving (2).
ThefastestalgorithmsforsolvinglinearM-estimationproblemsarebasedon
Newton’s method with linesearch, though it demands that (cid:26)(t) is continuously
di(cid:11)erentiable. This rules out Newton’s method for (cid:12)nding ‘ estimates. Other
1
peculiar things happen if the (cid:26)-functions are non-convex, as the case is for the
10
