A Sparse Multilevel Implementation of the
LP Dual Active Set Algorithm *
Timothy A. Davist William W. Hagerl
July 26, 2003
Abstract
We present a factorizationbased implementation of the LP Dual
Active Set Algorithm. This implementation exploits a proximal ap
proximation to the dual function, a multilevel solution strategy, and
recently developed algorithms for updating a sparse Cholesky factor
ization after a small rank change. We compare solution times to those
of simplex and barrier algorithms, as implemented in CPLEX, using
the Netlib LP Test Problems.
Key words. Dual active set algorithm, linear programming, numerical
experimentation, simplex method, barrier method, interior point method
AMS(MOS) subject classifications. 90C05, 90C06, 65Y20
*This material is based upon work supported by the National Science Foundation under
Grant No. CC'r[,i'.i'70. Any opinions, Ifl,,11,_ and conclusions or recommendations
expressed in this material are those of the author and do not necessarily reflect the views
of the National Science Foundation.
tdavis@cise.ufl. edu, http://www.cise.ufl.edu/~davis, PO Box 116120, Department
of Computer and Information Science and EHIli._ I i.i, University of Florida, Gainesville,
FL 326116120. Phone (352) 3921481. Fax (352) 3921220.
Shager@math.ufl.edu, http://www.math.ufl.edu/hager, PO Box 118105, Depart
ment of Mathematics, University of Florida, Gainesville, FL 326118105. Phone (352)
3920281. Fax (352) 3928357.
1 Introduction
Suppose that : A x X 1* R where A C " and X C P' are closed and
convex. In this paper, we develop techniques for solving the problem
sup inf (A,x). (1)
AEA xEX
We refer to in (1) as the "Lagrangian," while the "dual function," defined
on A, is given by
(A)= inf (A,x). (2)
xEX
Hence, the maximin problem (1) is equivalent to the maximization of the
possibly nonsmooth dual function. We focus in particular on the case where
the dual function is obtained from a proximal regularization (e. g. [15, 35,
36, 37, 39, 40, 41, 43]) of the linear programming problem
min cTx subject to Ax = b, x > 0, (3)
where A is m by n, and all vectors are of compatible size. In this case,
X= {x PE :x> 0}, A= (4)
and the associated regularized Lagrangian is
(A, x) = c'x + AX(b Ax) IXA pi (5)
2
where 5 > 0 is a scalar and pI E '" is the proximall shift."
Even though the regularized Lagrangian (5) is strongly concave in A, the
associated dual function is nonsmooth, and computing a maximizer may not
be easy. Consequently, we also consider a regularization in x given by
,(A;y) = min (A,x) + x (6)
xex 2
where e > 0 is a scalar and y E X. We solve (1) by solving a sequence of
smooth optimization problems corresponding to the iteration:
j = arg max L,(A;yj), (7)
AEA
Yj+l = arg min (Aj,x) + y,'. (8)
xEX 2
The dual active set algorithm (DASA) [20, 21, 23, 24, 25, 27] is used to im
plement the iteration (7)(8). In the context of linear programming (LP),
we show in [24] that the resulting proximal/DASA scheme closely approx
imates, to within O(5), the search directions of the LP DASA [23] gotten
by orthogonal projection. Based on preliminary numerical experiments in
[27], it appeared that a direct implementation of the projective LP DASA
could not compete with the simplex method implemented with an LU fac
torization. With the proximal approach, the orthogonal projections of the
LP DASA are approximated by the solution of a symmetric, positive definite
linear system.
The numerical results presented here use factorizationbased sparse ma
trix techniques [8, 9, 10] to solve the linear systems that arise in each step of
the DASA. In contrast [24] gives numerical results based on an SSOR iter
ative solver [22]. This iterative implementation is attractive in cases where
the constraint matrix A is sparse, while the Clii1,,ky factorization of AA'
is relatively dense. In this situation, the product between A and a vector
can be evaluated much faster than a factorization or an update of a factor
ization. Moreover, the storage requirements of this iterative approach are
much less than that of factorizationbased approaches. In particular, we saw
in [24] that this iterative approach could solve large LPs connected with
lower bounds to the quadratic assignment problem which were intractable
for factorizationbased simplex, barrier, or dual active set methods.
This iterative approach, however, is not as efficient as a factorization
based approach when both A and the Cl,,li.ky factorization of AA' are
sparse. In this paper we compare the efficiency of our factorizationbased LP
DASA to simplex and barrier algorithms as implemented in the commercial
code CPLEX [3]. Another variation of the LP DASA appears in [38], where
the bound constrained least squares problem appearing in the first step of
the algorithm [21, 23] is replaced by an unconstrained least squares problem.
As a result, a nondegeneracy assumption is needed to ensure convergence,
while we obtain in [23] convergence without any assumptions on the problem.
Previous applications of the dual active set algorithm have been to classes
of bound constrained control problems [2, 26, 42] and to quadratic network
optimization [25].
Our paper is organized as follows: In Section 2 we prove convergence of
the iteration (7)(8). In Section 3 we review the dual active set approach for
solving the proximal subproblems, focusing in particular on the case of linear
programming. Section 4 explains how we can delete entire rows from a linear
programming problem, using our dual active set approach. These dropped
rows are identified dynamically, as the iterations progress, as opposed to
resolve techniques that identify rows that can be dropped before starting the
solution process. Section 5 develops the multilevel DASA strategy. In this
approach, we subdivide the original optimization problem into a sequence
of smaller subproblems; the solution of the subproblems at any level are
used as starting guess for the solution at the next higher level. Multilevel
techniques have obvious benefits in a parallel computing environment since
the subproblems can be solved on separate processors. We show, however,
that even on a single processor, multilevel techniques can yield faster solution
times. \Niil iii ;I comparisons with simplex and barrier codes in CPLEX are
given in Section 6. For fixed e, the convergence of (7)(8) is extremely slow
in some situations. One way to deal with these cases is to let e tend to zero.
A completely different approach, developed in the Appendix, is based on an
active set strategy applied to the proximal shift vector. This update of the
proximal shift can be extremely efficient in a neighborhood of a solution.
2 Proximal approximation
Since the extremand in (6) is strongly convex in x, the minimizing argument,
denoted x(A), is unique. As a consequence of [6, Thm. 2.1], if the derivative
VB(A, x) exists and is continuous with respect A and x, then e( y ;y) is
differentiable and
VA'C(A; y) = VA'(A, x(A)). (9)
Moreover, if the mixed derivative VxVAf(A, x) exists and is continuous, then
by [17, Lemma 2.2] or [18, Theorem 4.1], x(A) is Lipschitz continuous in A,
and by (9), ,(A; y) is Lipschitz continuously differentiable in A. Hence, by
introducing the c term in (6), we are able to transform a (possibly) nonsmooth
extremand in (1) into a smooth function e.
In the smooth proximal iteration (7)(8), Aj is uniquely determined when
the associated extremand is uniformly, strongly concave in A; that is, there
exists 5 > 0, independent of x, such that
(A, ) < (1, x) + VA(, x)(A ) I (10)
2
for each A, /I E A. Similarly, yj+i is uniquely determined since the associated
extremand is strongly convex in x. Since the iteration (7)(8) does not seem
to conform to any of the proximal schemes appearing in the literature, we
provide a convergence proof:
Theorem 1. Suppose the L.,iI.Ii;.,,, is continuously differentiable,
( ,x) is uniformly, strongly concave for each fixed x E X, (A, ) is
convex for each fixed A E A, and the dual function is not identically equal to
oo. Then the Aj ',1 rated by (7)(8) (,,. I ,, to the maximizer A* of (1),
and we have
(Axj+l; yj+i) < (Aj; Y)  2 Ij+i
2
(11)
A11' 1. yj+i
2
Proof. We focus on the iterates associated with j = 1 and j = 2. Putting
x = Y2 in (6) when A = A2 and y = Y2, and exploiting the uniform, strong
concavity of ( x), gives:
f(A2; y2) < (A1, Y2) + VA(Al, Y2)(2 ) A2
2
A1 '. (12)
Also, by the firstorder optimality conditions at the maximizer in (7), it
follows that
VA(Ajy )(A Aj) <0
for all A E A. By the definition of y2,
(A,; yI) = (A, y2)+
2
(13)
(14)
y ll'.
Combining (12), (13) with j = 1, and (14) gives
6(A2; 2) < (A1; y) A2
2
SI
Analogous to (15), the iterates (Aij+, yj+i) and (Aj, yj) satisfy the rela
tion
(Aj+1; yj+i) < (A; Yj) 2 2 Aj+l
It follows that
11
2i1
1=1
21 .j+i
I,+Xii 1 I+ + (i.,+iA Y Hi) .
y l .
(15)
(16)
f,(x,; Yj) <_ r'(A1 YI)
For any A E A, we have
(Axj; yj) > 6(A; yj) > o(A; yj) = (A). (17)
Since the dual function is not identically oo, it follows from (16) and (17)
that both the dual iterates Aj and the primal iterates yj form a Cauchy
sequence.
Let A* be a limit of a subsequence of the Aj. Further prune this subse
quence as necessary so that the associated yj+i approach a limit y*. From
the continuity of the derivative, (13) implies that
VA(A*, y*)(A A*)< 0 (18)
for all A E A. The firstorder optimality conditions for the minimizer in (8)
yield
(Vx(AX, yj,+) + E(yj+1 yj)') (x yj+i) > 0
for all x E X. Again, (16) implies that lyj+i yjl tends to 0, from which
we deduce that
Vx(A*,y*)(x y*) > 0 (19)
for all x E X. Since the firstorder optimality conditions are sufficient for
optimality in our convex setting, we conclude from (19) that y* minimizes
(A*,x) over x E X:
(A*,y*) = min (A*,x)
xEX
Again, by [6, Thm. 2.1], the derivative of the dual function is given by
V(A*) = V (A*, y*). (20)
Strong concavity of the Lagrangian (10) implies strong concavity of the as
sociated dual function; that is, take 0 E [0, 1], and add 0 times the inequality
(10) with A = AX and yI = 0AI + (1 O)A2, to 1 0 times the inequality (10)
with A = A2 and / = OAX + (1 0)A2, to obtain
(0AX + (1 0)A2, x) > (Al,x) + (1 8)(A2, x) + 0(1 0) A1 A
> 0L(Al) + (1 0)C(A2) + 60(1 0) AI AIl.
Minimizing over x E X gives
(0AI + (1 0)A2) > 061(Al) + (1 0)C(X2) + 60(1 )AX A,11. (21)
Combining (18), (20), the concavity of the dual function, and the fact that the
firstorder optimality conditions are sufficient for optimality in the concave
setting, X* maximizes the dual function. By (21), the maximizer of the dual
function is unique and the entire Xj sequence convergences to A*. O
Remark. Although e is independent of j in the statement of Theorem 1,
it could depend on j without any changes in the analysis. When e depends
on j, (11) becomes:
,+1(Aj+1, yj+l) _< 6(Aj, y) j+1 2j ( I ,j+l Yj
3 The dual active set algorithm
As mentioned in the introduction, we will implement (7) using the dual active
set algorithm. Hence, we consider the maximin problem (1) in the case that
the Lagrangian is strongly convex in x and the primal constraint set X is
the nonnegative cone (4). If B C {1, 2,..., n}, let XB be the subvector of
x consisting of those components xi associated with i E B, and define the
functions
CB+(A) = min (A,x), (22)
XB>O
and
BO(A) = min (A,x).
XB=O
In carrying out the minimizations above, the components of x corresponding
to indices in the complement of B are unconstrained. By the strong convexity
assumption, there exists a unique minimizer x(A, B) in (22) for each choice
of A and B. Let x(A) denote x(A, B) in the case B = {1,2,..., n}.
In presenting the DASA iteration below, A denotes the current DASA
iterate while A+ is the next iterate.
Dual Active Set Algorithm
Convergence test: If A is a solution of (1), then stop.
Initialization: Set k = 0, vo = A, Bo = {j [1,n] : xj(A ) = 0}.
Subiteration:
k = arg max Bo(A) and vk+= arg max B+(A).
A I 11,. ,0. ] k
Here Ivk, Wk] denotes the line segment connecting vk and wk.
Constraint deletion:
Bk+ l j E Bk : yj = 0, (Vx(Vk+l, y)) > 0}
where y = x(vk+l, Bk).
Stopping criterion: If vk+1 = Wk, then set A+ = vk+l, and proceed
to the next iteration, starting with the convergence test. Otherwise,
increment k and continue the subiteration.
By [23, Thm. 1], the DASA solves (7) in a finite number of iterations. The
proof is based on the fact that the subiteration strictly increases the value of
the dual function; as a result, the final set Bk generated in the subiteration
cannot repeat. Since the sets Bk are chosen from a finite set, convergence
occurs in a finite number of steps. Note that the indices contained in the final
Bk of the subiteration will be included in the initial Bo of the next iteration.
Hence, the dual initialization step adds indices to the current "bound set"
B, while the constraint deletion step frees bound indices.
The computationally intensive operation in the DASA is the maximiza
tion of gBo. In the case of linear programming, with proximal regularization
of primal and dual variables, the maximization of C,.. in the subiteration
amounts to solving a problem of the form:
max min c'x + A (b Ax) + x  6  (23)
A XB=o 2 2
Given any set B C {1, 2,..., n}, let F denote its complement. The set F
corresponds to "free" indices (components not at a bound) of the primal
approximation. The maximin problem (23) is equivalent to the following
linear systems:
(AFA + oI)A = arp + AFc + e(b AFYF), and (24)
1
XF = YF (CF AA) (25)
where a = c5 and AF denotes the submatrix of A associated with columns
corresponding to indices in F. We solve (24) for A and substitute in (25) to
obtain the associated XF that attains the minimum in (23).
4 Row dropping in linear programming
In this section, we focus on linear programming, and a particular simplifi
cation that we call "row dropping." With the proximal version of the LP
dual active set algorithm, the subiteration requires the solution of (24), for
which the matrix AFA' is typically much sparser than the full matrix AA'.
Additional sparsity in linear programming is achieved by dropping not only
binding columns, but "inactive rows." We say that row i is inactive if there
exists a column j of A that is completely nonzero except for an element aij,
and there exists an optimal solution x of the LP with xj > 0. In this case,
it follows from complementary slackness that at a corresponding solution to
the dual problem, Aiaij = cj. Since A, = cjlaij is known, we simply hold
it fixed during the solution of the dual problem. In essence, we eliminate
the ith equation in (24) and solve a system with one less row and one less
column.
Although presolvers [1] may identify some inactive equations, other inac
tive equations may not be apparent to even the best presolvers. With the
LP DASA, there is a natural way to dynamically delete inactive equations
or add newly activated equations. FilI some notation: Let C denote the
set of column indices j corresponding to column singletons. That is, for each
j E C, there is precisely one i such that aij : 0. When smoothing the dual
functions as in (8), we do not smooth the components of x corresponding to
column singletons:
B+(A) = min (A, x) A + ( and
XB>o 2 2 ( and
jECc
BO(A) = min (A, x) A  + ( (26)
XB=o 2 2 (26)
jECc
Let J7 denote the set of column singletons in row i. In other words, J7
is the set of column indices with the property that aj : 0 and aj is the
only nonzero in column j. If for some i, there exist p and q E J7 with
cp/aip = cq/aiq, then columns p and q can be combined together into a single
column (see [1]). To simplify the discussion, we assume that cp/aip # cq/aiq
for all p, q E J7. At subiteration k, we define a projection Pk as follows:
af if A, > a,
(PkX)i = k if Ai < k,
Ai otherwise,
where
ai = min {cj/aij j E J Cj/aij > (ik)i}
g = max {cj/ai j : J cj/aij < (Vk)i 
If the argument of the min above is empty, then set a = oo; if the argument
of the max is empty, then set /3 = oo. With this notation, we now give
the rowdropping version of the LP DASA. The dropped rows correspond to
indices in the set Dk below.
DASA for LP (with row dropping)
Convergence test: If A is a solution of (1), then stop.
Initialization: Set k = 0, vo = A, Bo = {j E Cc : xj(A) = 0}, and
Do is the set of indices i with the property that for some j E J, we
have
A = jl/aij and (b 5(, 1) 6 ,( p(i))ia < 0.
1j
Subiteration:
Wk = arg max{Bo(A) : AD = (Vk)Dk},
Vk+i = arg max{lB+(A) A E Pk [k, Wk]}
Constraint deletion:
Bk+1 = {j Bk : yj = 0, (Vxf(Vk+1,y))j > 0},
Dk+ = {i" (Vk+1)i = v or },
where y = x(vk+l, Bk).
Stopping criterion: If Vk+1 = k, then set A+ = wk, and proceed
to the next iteration, starting with the convergence test. Otherwise,
increment k and continue the subiteration.
Observe that in the constraint deletion step, we delete those rows for
which the associated component of the dual multiplier has reached one of the
bounds ao or P/3. In the dual initialization step, we not only bind columns (as
we did before), but also add rows which may have been previously dropped.
Both of these operations ensure that when we are at a nonoptimal point, the
dual subiteration yields strict ascent of the dual function. The dual function
is nonsmooth since we only smooth with the set CC in (26). At points of
differentiability, the derivative of the LP dual function is
V(A) = b Ax(A) 5(X p).
In other words, compute the partial derivative of the Lagrangian with respect
to A and evaluate at x(A).
Now suppose Ai = cj/aij where aij is a column singleton in row i. Observe
that
= 0 when Aiaij < cj, 7
S( A ) =I (27)
xo when Aiaij > cj.
It follows, due to the assumption cp/aip : cq/aiq for all p, q E J7, that at
a point of nondifferentiability where A, = cjlaij, the subdifferential of the
dual function with respect to its ith component is one of the intervals below:
(ooy] when aij > 0 where = (A ,) ~ il X).
[7, +oo) when aij < 0, j'
In the dual initialization step, we include in Do those rows that were pre
viously dropped, and where zero does not lie in the subdifferential interval.
When the subdifferential interval corresponding to the ith component of A
is strictly positive, a small increase in A, will increase the dual function; like
wise, when the subdifferential interval corresponding to the ith component
of A is strictly negative, a small decrease in A, will increase the dual function.
By freeing these components of A for which zero lies outside the subdifferen
tial interval, we ensure strict increase in the dual function when the current
dual iterate is not optimal.
The proof of finite convergence is the same given previously for the DASA
(e. g. [20, 23, 25]), but with the following adjustment: Again, finite con
vergence is due to the fact that certain sets cannot repeat. With the row
dropping algorithm, the entity that cannot repeat is the pair (Bk, (Wk)Dk)
associated with the final dual subiteration. The vector (wk)Dk comes from
a finite set since all components have the form cj/aij for some i E Dk and
j E JG.
5 A multilevel decomposition
Multilevel techniques have proved to be very effective for solving elliptic par
tial differential equations [5, 16], and for solving NPhard graph partitioning
problems [29, 30, 31, 32, 33]. We now describe a multilevel version of the
dual active set algorithm. The decomposition that we use is related to the ar
rowhead partitioning presented in [13], however, instead of working towards
an arrowhead form, we utilize a more general recursive partitioning, closer in
spirit to the nested dissection ordering [28, 30] used for symmetric matrices.
Our discussion is in the context of the general dual active set algorithm of
Section 3, not the LP version of Section 4. The multilevel scheme is described
using a tree T with r nodes, where r is the root. First, some terminology:
The parent of any node i, denoted 7r(i), is the node above i in T. The
children of i are the inverse of the parent map:
r1(i= {j: (j) = i}.
These are the nodes beneath i in T. We denote the path from node i to the
root of the tree as the sequence P,, where
P,, = (i, 7r(i), 7r(7r(i)), 73i), r) .
The ancestors of i, denoted as the set Ai, are the nodes in the path P7 but
excluding i itself. If i is an ancestor of j, then j is a descendant of i.
In a multilevel decomposition of the dual function, we relate each node i
of the tree to components of the dual variables associated with a set Ri C
{1, 2,..., m}. The ith subproblem will be expressed in terms of Aj for j E
Ri. These sets should have the following properties:
M1. Rr = {1, 2, ... m}.
M2. For each node i of T, Rj c Ri when i = 7(j).
l:i For each node i of T, Rj n 7k = 0 whenever i = 7r(j) = ir(k) and
j#k.
In other words, as we move up the tree, an increasing number of dual vari
ables appear in the subproblems (property M2), and when we reach the root
of the tree, the subproblem is the entire problem, involving all the dual vari
ables (property Ml). Finally, at any level of the tree, the variables in the
various subproblems are all disjoint, so that the subproblems can be solved
independently (property '.1:3).
The final multilevel property explains how the dual functions decompose.
Again, some notation: The separator indices S, for node i are those elements
of Ri not belonging to any of the children:
S, = ,Ri \U Z
where \ denotes the usual set difference. Thus, we can also write
Ri s,u S U ) Z
Let Si denote the union of the separators for all ancestors of i,
s = U sj.
jcAi
Note that ,i is not included in Si since i is not an ancestor of i. The set Sr
is empty since the root of the tree has no ancestors. The fourth property of
a multilevel decomposition is the following:
..14 There exist realvalued functions 1, 2..., r with
,(A) = (A)= inf (A, x).
xEX
Moreover, the function i associated with any node i can be expressed
as a sum of the functions for its children:
A(AiAXS,)= SC (AXA. .
i='(c)
We now explore some properties of a multilevel decomposition.
Lemma 1. Let i and j be nodes in a tree T associated with a multilevel
decomposition. If i is neither an ancestor nor descendant of j, then RinRj =
0.
Proof. Consider the sequence of ancestors 'z"i), p = 1, 2,... and 7q(j),
q = 1,2,.... Let k = '"(i) = 7q(j) be the first common ancestor of i
and j. Since k is the first common ancestor, the children c = rPl(i) and
c2 = 1 (j) are distinct. By (r.1:), 7z, n R,, = 0. By (M2) 7, C R,, and
Rj C R2,. Hence, Z, n Rj = 0. o
The usual definition of the subtree rooted at node i is the subtree of T
induced by node i and all of its descendants. We use another kind of subtree
in our algorithm, which we refer to as a special subtree. Given a subset AF
of the nodes of T, the induced special subtree is obtained by pruning from
T all the descendants of the nodes AF. We denote OT as the leaves of the
tree T (that is, the childless nodes). The root nodes of T and T7 are always
identical, and neither tree can be empty.
Lemma 2. Given a multilevel decomposition associated with a tree T,
and a special subtree Tr, we have
(A)= Y &(AR".A81). (28)
IEa'Tr
Proof. Suppose that (28) is violated for some special subtree TAo of
T. Let 10 E 07To, and let po = 7r(lo) be its parent. Let TAg be the special
subtree of T obtained by pruning from TA, all the children r (Po). That is,
AI = No U {Po}. Either (28) holds for TA, and we stop the pruning process,
or it is violated for Tr,, in which case we choose l1 E 07A,, pi = (lI),
and TA2 is the special subtree of T induced by Ai2 = AFi U {pi}. Since
the decomposition (28) holds for the special subtree {r}, which is contained
in any special subtree, this pruning process must terminate with a special
subtree TAk for which (28) is violated, while (28) holds for 7Ak+,. Since 7Tk+i
is obtained by pruning the children of some node p from TAk, it follows from
(.114) that
Cr ,As,)= r c(AR.AS). (29)
P=5(c)
Since (28) holds for TAg+I, we have
Since p & 7A+I, and Tr, only differs from T7A+, by the children of p, it
follows from (29) that (28) holds for TAg, which is a contradiction. Hence,
(28) holds for all special subtrees of T. o
Lemma 3. If TA is a special subtree of T and i and j E T7r, then
7i n zj = 0 and Sin Rz = 0.
Proof. Since leaves of a tree cannot be ancestors of each other, the
condition zi, n Rj = 0 follows from Lemma 1. Let p and q be chosen so
that k = .''(i) = z'(j) is the first common ancestor of i and j. If I E R,,
where c2 = 7q1(j) is the child of k on the path down to node j, then I S,,
since the separator set St for any node t on the path between k and the root
excludes the indices in the sets R, for the children s E 7r(t). Moreover, if
I E R'2, then by (.13:),
I Rc1, (30)
where cl = rPl(i). Since the separator sets St for all nodes t on the path
between cl and i are all contained in Rc, it follows from (30) that I ( St.
To summarize, I St for all nodes t on the path from cl to i and I St
for all nodes t on the path from k to the root. It follows that I S,. Since
R C ',,2 we conclude that ~j n s, = 0. O
By Lemmas 1, 2, and 3, it follows that if TA is a special subtree of
T, then maximizing the dual function in (28) over AX for each I E &Tg,
while holding the other components of A fixed, decomposes into independent
maximizations:
max ( ,,,)
AX71
Before stating the multilevel DASA, we illustrate a decomposition using
the Lagrangian associated with a linear program with constraint matrix
110000
001000
A= 1 1 1 1 00 (31)
0 0 0 0 1 1
101010
For the tree shown in Figure 1, properties (M1)(:. 13) are satisfied since a set
7Ri at any level includes the union of the sets R, for its children c E 7l(i),
and at any level in the tree, the sets ,R are disjoint. With regard to ('14),
the complete dual function is
(A) = b'A + inf (c A'A)x.
x>O
Z5= {1,2,3,4,5}
R4 = {4}
n2 = {2}
Z3 = {1,2,3}
R1 = {1}
Figure 1: Multilevel decomposition associated with (31)
The product ATAx involves terms of the form Aiaijxj. Due to the zeros
in columns 5 and 6 of A, the terms involving products between {A1, A2, A3}
and {xs, x6} vanish. Hence the dual function decomposes into two separate
terms:
(A) = 3(A, A2, A3, A5) + 4(A4, A5)
The last arguments of 3 and 4 correspond to the separator {5}. The first
arguments correspond to the sets Z3 and 74, respectively. The decomposed
dual functions are:
C3(A1, A2, A3, A5)
4(A4, A5)
cl > AI + A3 + A5
bIA + b2A2 + b3A3 if > A + A3 + A5
oo otherwise
b4A4 + b55 if C5 A4 + A5
Sc6 > A4
oo otherwise
By the structure of 3, it can be decomposed as follows:
3(A) = f1(A A3, A5) + 2(A2, A3, A5),
Figure 2: Nonzero structure of AAT corresponding to multilevel partitioning
tree in Figure 1
where
AC(A,, A3, A5,
and
2(A2 ,3, A5)
SblA if cl > A1,+A3 +A5
= c2 > AI + A3
oo otherwise
b2A2 + bA if ca > A2 + A3 + A5
oo othe 4 > A3
oo otherwise
The smoothed dual functions L as well as the functions CB+ and Bo de
compose exactly as the dual function.
The multilevel structure indicated in Figure 1 is connected with the
nonzero structure of the symmetric matrix AAT. The right side of Figure 2
shows the nonzero pattern of AAT corresponding to the tree of Figure 1.
The regions labeled 0 correspond to (i, j) pairs with the property that o I,,
is zero for all k. Figure 3 shows the nonzero pattern of A with this partition.
A permutation of the rows of A to achieve a multilevel structure can be
accomplished using graph partitioning codes such as 1. ii, [31, 32, 33] and
Q,3 S1
S
11 02\5
IIC^
I I
II &
II ^T1
I I I
2i 0 0 0 0t
I0 0 > R5
24 0 0 0 0
Figure 3: Nonzero structure of A corresponding to multilevel partitioning
tree in Figure 1
Cli1n, [29, 30]. Given the nonzero pattern of a symmetric matrix M, these
codes can be used to partition the row indices into two sets Q1 and Q2,
chosen to minimize the number of nonzeros mij with i E Q1 and j E Q2
The separator set is gotten by extracting elements from Qi or Q2 in order to
obtain a set S for which mij = 0 whenever i E Q1\S and j E Q2\S. To obtain
an ordering like that shown in Figure 2, we can first apply one of the graph
partitioning codes to the matrix M = AAT. After permuting M, putting
the rows corresponding to Q1 \ S first, followed by the rows corresponding
to Q2 \ 8, followed by the rows corresponding to S, the nonzero structure
of AAT is that shown in the left side of Figure 2. Again, applying the
partitioning code to the upper left diagonal block on the left side of Figure 2,
we obtain the final partitioning shown on the right side of Figure 2, after a
suitable permutation of rows and columns
The size of the separator set measures, in some sense, the degree to which
an LP can be uncoupled into two independent problems. The extreme cases
are (a) when the the separator is empty and the problem uncouples into two
independent problems, or (b) when the separator includes all the rows, so
that the decomposed problem is the original one.
The multilevel version of DASA, based on a multilevel decomposition of
the dual function, and its modifications, is the following:
Multilevel Dual Active Set Algorithm
Convergence test: If A maximizes the dual function, then stop.
Dual initialization: Set No = 0, k = 0, vo = A B0 = {j : xj(A)
0}.
Dual subiteration: For each each I E O7rk,
wz = arg max{(,. (A) : AR = (vk)R},
v1 = arg max { B+(A) : A [vk,W']}
Set (vk+i)R (v)R, for each I e OT7N and (vk+i)i = (v)i if i
Constraint deletion, tree modification:
Bk+1 = fi E Bk : Yi = 0, (Vxf (Vk+I, Y))i > 0} ,
Ak+1 = kA/ U {(E O rk (vk+l)>R = (wk)R for all c E 7r1(r(l))}
where y = x(vk+,, Bk).
Stopping criterion: If TA = {r} and Vk+l = Wr, then set A+ = Vk+l,
and proceed to the next iteration, starting with the convergence test.
Otherwise, increment k and continue the dual subiteration.
As seen after Lemmas 1, 2, and 3, the maximizations over the leaves of Trk
in the dual subiteration decomposes into independent maximizations of local
dual function, which could be done in parallel. These local maximizations
all increase the value of the dual function. The proof of finite convergence
is the same given previously for the DASA (e. g. [20, 23, 25]). That is, in
the multilevel setting, we work up to the root of the tree T increasing the
value of the dual function, and by (Ml), the set R, contains all the indices
of A. Hence, in the multilevel setting, the final dual subiterate again takes
us to a maximizer of Cj.c for some set Bk, which can never repeat without
contradicting the strict ascent of the dual iterates.
6 Numerical comparisons
In this section, we compare the performance of the LP DASA to that of
simplex and barrier methods as implemented in the commercial code ILOG
CPLEX version 7.0 [3], focusing on problems from the Netlib LP test set with
sparse factorizations. In particular, we consider all the problems except (i)
the LPs associated with lower bounds for the quadratic assignment problem,
which we studied in [24], and (ii) the 5 problems which utilize the "ranges"
specification in MPS. Since the CPLEX barrier code automatically applies
a presolver to the input problem, we use the CPLEX resolved problem as
input to our code.
In our implementation of LP DASA, we apply the proximal iteration
Ik+1 = arg max inf cTx + AT(b Ax) A pi 11, (32)
A x>o 2
which is implemented using (7)(8). By [39, Thm. 1], the iteration (32)
converges to a solution of the dual to the linear program (3), assuming a so
lution exists. The factorization of the matrix AFAT in (24) is done using a
roworiented sparse Clii, 1,lky factorization algorithm [34]. We do not exploit
supernodes in this code, and thus it runs about 10 times slower, in terms of
floatingpoint operations per second, than the supernodalbased code used
in the CPLEX barrier method. Fortunately, we refactor the matrix infre
quently, and after a column is freed or a row is dropped, we use the sparse
modification techniques developed in [8, 9, 10] to modify the factorization.
Since the CIli1.ky factorization is computed rowbyrow, it can be done in
steps, in the multilevel setting. The submatrix associated with a subtree
(nodes 1, 2 and 3 in Figure 1, for example) can be factored without having
to compute the entire factorization. New rows corresponding to separators
can be added to the factorization without touching the previously computed
partial factorization. The codes for modifying a factorization after column
or row updates, take advantage of speedup associated with multiple rank
updates. As a result, a rank 16 update of a sparse matrix achieves flop rates
comparable to those of a BLAS dot product, and about 1/3 the speed of a
BLAS matrixmatrix multiply [11], in the experiments given in [9].
In coding LP DASA, we use COLAMD [7] to reorder the rows of A to
minimize the fill generated during a Ci, 1blky factorization of each diagonal
block of AAT corresponding to a leaf of the partitioning tree. This ordering
does not take into account the fill generated outside the diagonal blocks, so
the ordering is not optimal with respect to the fill generated in the separator
rows. A constrained COLAMD (such as that developed in [4] for a sparse
LU factorization) would likely result in better orderings.
Our convergence test is based on the LP optimality conditions. Com
plementary slackness is always satisfied, and the primal approximation x(A)
Figure 4: Nonzero structure of AA' for KEN_07
always satisfies the bound constraints x > 0. Hence, our convergence cri
terion is expressed in terms of the residuals in the primal and dual linear
systems:
1 AxII + <,A I 108_
I A IX II 
Here . denotes the maximum absolute component of a vector. Since
the test problems are given in MPS format, for which data often has 5 digit
accuracy, our 8 digit error tolerance for the equation is, in a sense, much more
accurate than the input data. With this convergence tolerance, the primal
approximations often yield a smaller cost than the published optimal values.
This is due to the fact that we are not satisfying the equation Ax = b with
16 digit accuracy, but only with 8 digit accuracy. With this small reduction
in the error tolerance in the equation, which seems reasonable when the input
data has only 5 digit accuracy, we are able to reduce the cost, sometimes in
the 5th significant digit.
The test problems were solved on a single processor (200 Mhz) in an II 1
SP cluster. We solved all the test problems with a single level DASA, while 13
problems, which could be partitioned into nearly independent subproblems,
were solved by the multilevel approach. The partitioning was done in the
following way: The problems KEN_07, KEN_11, KEN_13, KEN_18, PDS_02,
PDS_06, PDS_10, PDS_20 have the property that the nonzeros of AAT fall in
a series of diagonal blocks and along the outside border. For example, Fig
ure 4 displays the sparsity pattern for KEN_07, which has 7 x 7 = 49 diagonal
Figure 5: Nonzero structure of AAT for PDS_02
blocks. The PDS problems are similar except there are 11 diagonal blocks
(Figure 5 shows the sparsity for PDS_02). For this set of 8 wellstructured
problems, we specified that there were two levels, with the leaf nodes corre
sponding to the diagonal blocks of AAT. The remaining 6 problems, TRUSS,
D2Q06C, GREENBEA, GREENBEB, CYCLE, and FFFFF800, were partitioned
using the graph partitioning technique of Section 5, and the code Clin; 1I to ap
proximate a graph bisection. After each bisection of the adjacency graph for
AA', the separator was evaluated. If the number of nodes in the separator
was at most 10% of the size of the graph being bisected, we accepted the bi
section. Otherwise, the partition was rejected. The separator was evaluated
using a MATLAB code of John Gilbert based on maximum matching (MAT
LAB's dmperm), and the entire partitioning process was carried out through a
MATLAB interface which generated the partitioning tree and an associated
reordering of the rows. As an example, the problem D2Q06C (Figure 6) was
partitioned into a 32 node tree depicted in Figure 7. A more regular matrix
TRUSS (Figure 8) was partitioned into a 7node tree (Figure 9); the permuted
matrix is seen in Figure 10.
In Tables 13 we compare the execution times in seconds for the test
problems. The execution time is shown in bold if it is the fastest for that
problem, and the tables are sorted according to the fastest run time. Since
the builtin timer used for these comparisons has an error on the order of
hundredths of seconds, small time measurements contain significant error.
For problems OSA_60 and KEN_18, the barrier code generated an outof
Figure 6: Nonzero structure of AAT for D2Q06C
Figure 7: Partitioning tree associated with D2Q06C
Figure 8: Nonzero structure of AAT for TRUSS
23
:. '' I~ ~ :r ~ i "r r(g
F
I
.!. "`:
.i '
o
I
'";, Ji~
o
1;; :
O
'
;i : ; i
T~:
or.. ....
oa~L  k B'C :
ic is
r
Rf+lLI
D~'~r~ `
k" "
I i
Figure 9: Partitioning tree associated with TRUSS
Figure 10: Nonzero pattern of TRUSS after partitioning
memory message, and hence, no times are given. In Tables 13, we see that
in 24 out of the 106 test problems, LP DASA achieved the fastest measured
time; similarly, in 80 out of 106 problems, the simplex code (by default,
dual simplex) achieved the fastest time, while in 17 out of 106 problems, the
barrier code achieved the fastest time.
In the current version of the LP DASA code, the best results, relative to
the other codes, were obtained for the OSA and FIT problems for which many
rows could be dropped. The worst relative results were obtained with the
PILOT problems. For these problems, there was significant growth in error
in the CIi, 1bky factorization during downdates. Hence, further refinements
in the downdate formulas to reduce error growth should help with these
problems.
As noted earlier, multilevel LP DASA is wellsuited for a parallel com
puting environment since the subproblems can be assigned to separate pro
cessors. Nonetheless, even with a single processor, we obtained significant
speedup, relative to singlelevel LP DASA, for problems that nearly decom
pose into independent subproblems; the small subproblems in the multilevel
approach can be solved quickly and the overall flop count is much smaller
than that of a singlelevel implementation. In Table 4, we compare the solu
tion times of singlelevel and multilevel LP DASA for the 14 problems cited
earlier.
Enhancements to the LP DASA code that should further increase its
speed include the development of a constrained COLAMD similar to [4] that
takes into account fillin in the separators as well as in the diagonal blocks,
a factorization code that exploits supernodes, and the implementation of the
proximal shift technique developed in the appendix.
Problem
sc50A
sc50B
KB2
sc105
STOCFOR1
I;[ I:(:' IPE
VTP_BASE
BOI: :.i[ )
SCORPION
BEACONFD
STANDATA
STANDGUB
ADLITTLE
BLEND
SCAGR7
IlI AIl: 2 n
LOTFI
STANDMPS
sc205
i: A I \ 1 ~1
CAPRI
SCSD1
AGG
TUFF
AGG2
AGG3
SCTAP1
BRANDY
I i 'I [..
iIiO['04s
liiP04L
SAN C; IlS
 Ii: i. i.S
FINNIS
LP DASA
0.00
0.01
0.00
0.02
0.00
0.01
0.00
0.01
0.00
0.01
0.00
0.01
0.00
0.02
0.01
0.02
0.02
0.03
0.04
0.02
0.07
0.17
0.07
0.03
0.08
0.04
0.05
0.05
0.08
0.06
0.07
0.12
0.13
0.10
0.10
Table 1: \N ill ii ;il Comparisons, part 1
Simplex
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.00
0.01
0.01
0.01
0.01
0.01
0.01
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.02
0.03
0.03
0.03
0.03
0.03
0.03
0.04
0.04
Barrier
0.00
0.00
0.00
0.02
0.02
0.02
0.02
0.01
0.04
0.02
0.01
0.05
0.00
0.01
0.02
0.02
0.03
0.04
0.13
0.04
0.06
0.10
0.05
0.07
0.17
0.15
0.15
0.08
0.09
0.15
0.11
0.15
0.21
0.14
0.13
Problem
GFRD_PNC
SCTAP2
Mt ([);': [,, 1
SHIPO8S
PDS_02
GROW7
SCFXM1
FFFFF800
SHIP12S
SCAGR25
BANDM
. 11: I:: A
SCSD6
SCTAP3
CYCLE
STAIR
FIT1D
SHIP08L
KEN_07
FIT1P
ETAMACRO
DEGEN2
WOOD1P
SHIP12L
SCSD8
CZPROB
SCFXM2
GROW15
BNL1
STOCFOR2
GROW22
SCFX:.I!.
WOODW
LP DASA
0.33
0.19
0.10
0.23
0.08
0.33
0.06
0.17
0.12
0.17
0.17
0.07
0.13
0.23
0.15
0.25
2.78
0.53
0.10
0.23
0.16
0.11
0.23
0.30
1.70
0.37
0.35
0.17
0.37
0.19
2.83
0.73
0.32
0.58
0.48
Table 2: \ 1iliii ;il Comparisons, part 2
Simplex
0.05
0.05
0.05
0.05
0.05
0.05
0.07
0.06
0.06
0.06
0.06
0.07
0.07
0.07
0.08
0.08
0.09
0.10
0.15
0.10
0.11
0.41
0.13
0.15
0.15
0.16
0.56
0.20
0.18
0.54
0.30
0.30
1.58
0.33
0.37
Barrier
0.11
0.16
0.34
0.27
0.15
0.92
0.09
0.12
0.10
0.33
0.19
0.08
0.09
0.40
0.10
0.47
0.83
0.25
0.22
0.27
0.34
0.33
0.23
0.27
1.14
0.41
0.17
0.39
0.23
0.26
0.36
0.44
0.41
0.36
0.72
Problem
MAROS
BNL2
PILOT
25.1:47
PEROLD
D6CUBE
CRE_C
CREA
OSA_07
FIT2P
PILOT_WE
PILOTNOV
TRUSS
KEN_11
PDS_06
GREENBEB
FIT2D
GREENBEA
PILOT_JA
80r[) .In
OSA_14
D2Q06C
STOC(1 I.I1
PDS_10
MAROSR7
KEN_13
OSA_30
PILOT
CRE_D
CRE_B
OSA_60
PILOT87
PDS_20
KEN_18
DFL001
LP DASA
0.83
5.71
2.73
1.22
&, ,Sr
19.57
1.93
2.41
0.81
0.82
10.74
9.26
6.08
1.33
9.71
2.98
1.55
9.03
35.19
4.97
3.79
2.67
8.54
16.99
31.68
6.25
8.42
6.66
41.64
50.04
67.36
18.79
132.41
307.29
62.05
517.18
Table 3: \ niliii ;il Comparisons, part 3
28
Simplex
0.39
0.44
0.47
1.57
0.90
0.77
0.79
0.81
1.81
21.58
3.37
1.26
49.08
1.50
1.48
5.15
3.23
3.76
1.79
1.98
2.09
11.47
10.47
10.20
6.06
12.22
6.45
45.50
13.80
7.76
13.17
228.52
97.14
42.38
44.98
89.71
Barrier
0.42
1.51
0.57
0.65
0.77
2.02
1.30
1.51
2.03
2.30
0.99
1.10
1.18
2.94
9.67
1.52
2.84
1.73
1.81
2.23
2.70
5.04
3.35
5.76
32.04
8.41
9.71
11.58
6.79
13.69
16.98
*.**
23.91
157.16
*.**
121.92
Problem
PDS_02
PDS_06
PDS_10
PDS_20
i. .: .: _07
KEN_ 11
KEN_13
KEN_18
TRUSS
GREENBEA
GREENBEB
D2Q06C
CYCLE
FFFFF800
One Level
0.26
14.48
56.78
653.05
0.25
1.89
16.00
155.70
9.14
11.30
2.76
15.86
3.07
0.30
Multilevel
0.33
9.71
31.68
307.29
0.16
1.33
8.42
62.05
6.08
9.03
2.98
8.54
2.78
0.17
Table 4: Single level versus multilevel LP DASA
A Appendix: Optimizing the proximal shift
Although we show in Theorem 1 that the iteration (7)(8) converges, we
observe in some of the LP test problems, that convergence can be slow when
e is fixed. One way to accelerate convergence is to push e towards zero (by
the remark after Theorem 1, we can use a different e in each step). This
is the approach used for the numerical experiments reported in Section 6.
Here we develop a better approach for dealing with slow convergence near a
solution. NIi.ly, the proximal shift y is updated through an optimization
approach, rather than using (8). Although we observed in initial implemen
tations of the LP DASA that this proximal shift technique was more efficient
than letting e tend to zero near an optimum, time constraints soon made it
infeasible to both refine the primary code implementing (7)(8), and make
the corresponding adjustments to the proximal shift code. In this appendix,
we explain the proximal shift approach.
Let Q be defined by
Q(y) = sup inf (A,x) + I
AEA xEX 2
SII = sup e(A; y).
AEA
(33)
Proposition 1. Suppose ( ,x) is uniformly, strongly concave for each
fixed x E X, (A, ) is convex for each fixed X E A, A is closed and
convex, and X is compact and convex. Tl,, there exists a minimizer x* of
Q over X, there exists a unique maximizer X* of the dual function (2), and
Q(x*) = (A*). If y* is any minimizer of Q over X, then X* is the unique
maximizer of L ( ;y*) over A.
Proof. For any fixed x and y E X, the function
(AX, X) +( I\ 3112
2
(A.,x) +2lxy
is strongly concave in A. It follows from the analysis given at the end of the
proof of Theorem 1 that e( ; y) is strongly concave. We now show that for
any fixed A, f(A; ) is convex. FilI observe that the function
(A,x, y) = (X,x) + , '
2
for fixed A, is convex with respect to the pair (x, y) since the Lagrangian
(A, ) is convex by assumption while the term Ix y' is convex with
respect to the pair (x, y). Given any xi, yi E X, i = 1, 2, and 0 E [0, 1], we
have
,(A;0y + (1 )y2) = min (AX,x, 0y + (1 )y2)
xEX
< (A, Ox + (1 )x2, 0y + (1 )y2)
< 0(A,x1,y) +(1 0)(AX,2,y2).
Minimizing over xl and x2 E X gives
4(A; 0yi + (1 0)y2) < 8C(A; yi) + (1 0)4(; y2).
Since L( ; y) is strongly concave, L(A; ) is convex, and X is convex
and compact, it follows from [12, Prop. 2.4, p. 176] that L possesses a saddle
point (A*, x*):
,(A; x*) < L(A*; x*) < L(A*; y) for all A E A, y E X (34)
Since Q is given by a maximum over A in (33), we have, for any A E A and
y X,
Q(y) > inf (A, x) + x y' > (A). (35)
xEX 2
By the saddle point condition (34) and the definition of Q, we have
Q(x*) = sup L,(A; x*) < inf (A*;y)= inf (A*,x)+ xy = (A*).
AEA yEX x,yEX 2
(36)
Combining (35) and (36) gives Q(x*) = (A*) where x* minimizes Q over X
and A* maximizes over A.
If y* is an arbitrary minimizer of Q, then we have
Q(y*) = max ,(A;y*) > L,(A*;y*)
AEA
=min (A*,x) + Ix y I
xEX 2
> (XA*) = Q(x*) = Q(y*)
Consequently, A* maximizes ( ;y*). O
Based on Proposition 1, we can solve (1) by minimizing Q over X and
then maximizing C(A; x*) over A E A. To summarize, even though the
iteration (7)(8) is convergent in theory by Theorem 1, we find in some test
problems that e must be small for reasonable convergence. But taking e
small can lead to slower convergence of the DASA iterations used to solve
(7). As an alternative to taking e small, we can hold e fixed and minimize
the function Q over X.
We now give an active set strategy for minimizing Q, under the hypotheses
of Proposition 1, in the case that X is the nonnegative cone (4). For any
B C {1,2,...,n}, define
B(Y) = max mmin (A,x) + I ', (37)
AEA XB=O 2
and let x(y, B) denote the x that attains the minimum in (37). Let x(y) de
note x(y, B) in the case B = {1, 2,..., n}. In the statement of the algorithm
below, y denotes the current primal iterate while y+ is the next iterate.
Primal/Dual Active Set Algorithm
Convergence test: If y minimizes Q, then stop.
Initialization: Set k = 0, xo = y, Bo = {j E [1, n] : xj(y) = 0}.
* Subiteration: Set
zk = arg min QB,(y)
y>O
Xk+1 is the first y E [Xk, k], the line segment connecting xk and zk,
with the property that x(y, Bk)j = 0 for some j E B'.
Constraint addition: Bk+ = Bk U {j}.
Stopping criterion: If xk+1 zk, then set y+ = Xk+l, and proceed
to the next iteration, starting with the convergence test. Otherwise,
increment k and continue the subiteration.
With the primal/dual active set scheme, we achieve convergence in a finite
number of steps, unlike the iteration (7)(8) for which convergence takes
place in the limit as j tends to infinity. On the other hand, to implement
the primal/dual scheme, we first need a code to carry out the maximization
in (7) since the vector x(y) appearing in the initialization step is obtained
by solving a problem of the form (7). Hence, the code to implement (7) is
critical for the implementation of the primal/dual algorithm.
References
[1] E. D. ANDERSEN AND K. D. ANDERSEN, Pi. ,. I;,I in linear pro
grnuiil,,, i, Mathematical Programming, 71 (1995), pp. 221245.
[2] M. BERGOUNIOUX AND K. KUNISCH, Primaldual strategy for state
constrained optimal control problems, Computational Optimization and
Applications, 22 (2002), pp. 193224.
[3] R. .E. BIxBY, Progress in linear progrn, ,,, ;i,; ORSA J. on Computing,
6 (1994), pp. 1522.
[4] I. Bi: A\'.::.! AN AND S. TOLEDO, Nesteddissection orderings for sparse
LU with partial I,. .;;,, SIAM J. Matrix Anal. Appl., 23 (2002)
pp. ;.' 1012.
[5] J. H. BRAMBLE, 11,ill;,, i;! methods, John Wiley, New York, NY, 1993.
[6] F. H. CLARKE, Generalized gradients and applications, Transactions
of the American Mathematical Society, 205 (1975), pp. 247262.
[7] T. A. DAVIS, J. R. Gi [.i:i:i', S. I. LARIMORE, AND E. G. NG, A
column approximate minimum Ji,, oi ,, i, ,;, uli ,,,,;i, AC'.1 Transac
tions on Mathematical Software, (under submission, see also TR00005,
Univ. of Florida, CISE Dept., www.cise.ufl.edu).
[8] T. A. DAVIS AND W. W. HAGER, l,,;f,;,,,l/ a sparse C'li Ii/ fac
torization, SIAM J. Matrix Anal. Appl., 20 (1999), pp. 606627.
[9] T. A. DAVIS AND W. W. HAGER, Multiple rank modifications of a
sparse C(',i, ,/ j factorization, SIAM J. Matrix Anal. Appl., 22 (2001),
pp. 9971013.
10] T. A. DAVIS AND W. W. HAGER, The effect of row additions and
deletions on a sparse (C'li / j factorization, in preparation.
11] J. DONGARRA, J. DU CROZ, S. II\:.I:.I IAI :1NG AND I. DUFF, A Set
of Level 3 Basic Linear Algebra S;dprrilran AC(.I Transactions on
Mathematical Software, 16 (1990) 117.
12] I. EKELAND AND R. TEMAM, Convex Analysis and Variational Prob
lems, NorthHolland, Amsterdam, 1976.
13] M. C. Fl:1:i:is AND J. D. HORN, Par;i;ll ;,,,. mathematical programs
for parallel solution, Math. Prog., 80 (1'1'), pp. 3562.
14] P. E. GILL AND W. MURRAY, A numerically stable form of the simplex
(l,,, Ii/1,,, Linear Algebra and Its Applications, 7 (1973), pp. 99138.
15] C. D. HA, Generalization of the proximal point u1,li,,,;11i, SIAM J. Con
trol Optim., 28 (1990) 503512.
16] W. IT C:.'ii SCH, .1,///Ii,, ;. methods and applications, SpringerVerlag,
Berlin, I'',
17] W. W. HAGER, Convex control and dual approximations, part I, Con
trol and Cybernetics, 8 (1979), 522.
[18] W. W. HAGER, Inequalities and approximation, Constructive Ap
proaches to Mathematical Models, C. V. Coffman and G. J. Fix, eds.,
Academic Press, New York, (1979), pp. 189202.
[19] W. W. HAGER, U1ii,,;,.i the inverse of a matrix, SIAM Review, 31
(lI' '), pp. 221239.
[20] W. W. HAGER, The dual active set lpi.J/i,,ii, Advances in Optimiza
tion and Parallel Computing, P.M. Pardalos, ed., North Holland, Ams
terdam, (1992), pp. 137142.
[21] W. W. HAGER, The LP dual active set Jl, i;,,, High Performance
Algorithms and Software in Nonlinear Optimization, R. De Leone,
A. Murli, P. M. Pardalos, and G. Toraldo, eds., Kluwer, Dordrecht,
l'I', pp. 243254.
[22] W. W. HAGER, Iterative methods for nearly .;,spii, linear systems,
SIAM J. Sci. Comput., 22 (2000), pp. 747766.
[23] W. W. HAGER, The dual active set lri ,I,;11i, and its application to
linear progn i,, ,, in,; Computational Optimization and Applications, 21
(2002), pp. 263275.
[24] W. W. HAGER, The dual active set iJii ;11i,, and the iterative solution
of linear pro'tram . in Novel Approaches to Hard Discrete Optimization,
P. M. Pardalos and H. Wolkowicz, Eds., Fields Institute Communica
tions, Vol. 37 (2003), 95107.
[25] W. W. HAGER AND D. HEARN, The dual active set u1rilimi,,, applied
to quadratic networks, Computational Optimization and Applications,
1 (1993), pp. 349373.
[26] W. W. HAGER AND G. IANCU[.[:ECU, Dual approximations in optimal
control, SIAM J. Control Optim., 22 (1984), pp. 423 4t11
[27] W. W. HAGER, CHUNLIANG SHIH, AND ERIK O. LUNDIN, Active set
strategies and the LP dual active set ulI .,;li,,, Department of Mathe
matics, University of Florida, Gainesville, Florida, 1996.
[28] M. T. HEATH AND P. RAGHAVAN, A Cartesian parallel nested dissec
tion u',I, ;in,,,. SIAM J. Matrix Anal. Appl., 16 (1995), pp. 235253.
[29] B. HENDRICKSON AND R. LELAND, A multilevel (,,,,,;lii,, for parti
i;,r,,,, irplh Proc. Supercomputing '95, A(C'.1 1995.
[30] B. HENDiC:'I.ON AND E. ROTHBERG, Impro, ,,.I the runtime and
quality of nested dissection orderi, SIAM J. Sci. Comput., 20 (1999),
[31] G. KARYPIS AND V. KUMAR, A fast and high quality multilevel scheme
for paril;I,' ,,, i, , 1.,, i/,ii/, SIAM J. Sci. Stat. Comput., 20 ('1' '),
359392.
[32] G. KARYPIS AND V. KUMAR, Multilevel kway pari.,.,,, scheme for
,, I,,l,, ,, i,. J. Parallel Distrib. Comput., 48 (1999), 96129.
[33] G. KARYPIS AND V. KUMAR, Parallel multilevel kway ilii,,,i;,,,,
scheme for ;,,i ,,, ri1ph. SIAM Review, 41 (1999), 278300.
[34] J. W. H. LIu, A ,'I.1 ,1 envelope method for sparse factorization by
rows, AC'.I Transactions on Mathematical Software, 17 (1991) 112129.
[35] F. J. LUQUE, Asymptotic (i,, I ,,i ,i analysis of the proximal point
(,l,,;i, ;i,, SIAM J. Control Optim., 22 (1984), pp. 277293.
[36] B. MARTINET, !i ,',,,,;s l.,., d'inequations variationelles par ap
proximations successives, Revue Francaise Informatique et Recherche
Operationnelle, 4 (1970) 154159.
[37] B. MARTINET, Determination approachee d'un point fixe d'une appli
cation pseudocontractante, Comptes Rendus des Seances de l'Academie
des Sciences, Paris, 274 (1972) 163165.
[38] P.Q. PAN, I ,,;fi ,',;in of the dual projective pivot ,lJ,'i I,,, ;ii for linear
progruiiiiio,,, Department of Applied Mathematics, Southeast Univer
sity, China, 2003.
[39] R.. T. ROCKAFELLAR, Monotone operators and the proximal point al
,i,,II,,, SIAM J. Control Optim., 14 (1976), pp. 877 '',
[40] R. T. ROCKAFELLAR, Ai,,,,, rl ,1 ,"1,'"';,n's and applications of the
proximal point ul1,i, ; ,, in convex progrnii,,ii,;, Mathematics of Oper
ations Research, 2 (1976), 97116.
[41] J. E. S['I.:(:ARN, Submonotone mappings and the proximal point (,lI,
rithm, Niliii ;il Functional Analysis and Optimization, 4 (1982) 123
150.
[42] S. Volkwein, /11n,,n i SQP techniques for the control constrained opti
mal boundary control for Burger's equation, to appear, Computational
Optimization and Applications, Vol. 262.
[43] C. ZHU, Asymptotic convergence analysis of some inexact proximal point
(,l,,7,i,/,,,s for minimization, SIAM J. Optim., 6 (1996), pp. 626637.
