Group Title: Department of Computer and Information Science and Engineering Technical Reports
Title: A sparse multilevel implementation of the LP Dual Active Set Algorithm
CITATION PDF VIEWER THUMBNAILS PAGE IMAGE
Full Citation
STANDARD VIEW MARC VIEW
Permanent Link: http://ufdc.ufl.edu/UF00095598/00001
 Material Information
Title: A sparse multilevel implementation of the LP Dual Active Set Algorithm
Series Title: Department of Computer and Information Science and Engineering Technical Reports
Physical Description: Book
Language: English
Creator: Davis, Timothy A.
Hager, William W.
Publisher: Department of Computer and Information Science and Engineering, University of Florida
Place of Publication: Gainesville, Fla.
Publication Date: July 26, 2003
Copyright Date: 2003
 Record Information
Bibliographic ID: UF00095598
Volume ID: VID00001
Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.

Downloads

This item has the following downloads:

2003345 ( PDF )


Full Text









A Sparse Multilevel Implementation of the

LP Dual Active Set Algorithm *


Timothy A. Davist William W. Hagerl

July 26, 2003



Abstract
We present a factorization-based 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 32611-6120. Phone (352) 392-1481. Fax (352) 392-1220.
Shager@math.ufl.edu, http://www.math.ufl.edu/-hager, PO Box 118105, Depart-
ment of Mathematics, University of Florida, Gainesville, FL 32611-8105. Phone (352)
392-0281. Fax (352) 392-8357.










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 factorization-based 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 factorization-based 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 factorization-based 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 factorization-based 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 first-order 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


1-1
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 first-order 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 first-order 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
first-order 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 i-th 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 row-dropping 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 i-th 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 i-th 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 i-th 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 NP-hard 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 i-th 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 real-valued 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 (AX-A. .
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 = rP-l(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(A-R.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 = 7q-1(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 = rP-l(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 7-l(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. i-i, [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) : A-R = (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 7r-1(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
row-oriented 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
floating-point operations per second, than the supernodal-based 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 row-by-row, 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 matrix-matrix 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 well-structured
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 7-node tree (Figure 9); the permuted
matrix is seen in Figure 10.
In Tables 1-3 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 built-in 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 out-of-


































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 1-3, 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 well-suited 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 single-level 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 single-level implementation. In Table 4, we compare the solu-
tion times of single-level and multi-level 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 fill-in 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: \ 1il-iii ;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: \ nil-iii ;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


S||II = 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\ 3-112
2
(A.,x) +2lx-y|

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)+ ||x-y||- = (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. 221-245.

[2] M. BERGOUNIOUX AND K. KUNISCH, Primal-dual strategy for state-
constrained optimal control problems, Computational Optimization and
Applications, 22 (2002), pp. 193-224.

[3] R. .E. BIxBY, Progress in linear progrn, ,,, ;i,; ORSA J. on Computing,
6 (1994), pp. 15-22.

[4] I. Bi: A\'.::.! AN AND S. TOLEDO, Nested-dissection 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. 247-262.

[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 TR-00-005,
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. 606-627.

[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. 997-1013.

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) 1-17.

12] I. EKELAND AND R. TEMAM, Convex Analysis and Variational Prob-
lems, North-Holland, 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. 35-62.

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. 99-138.

15] C. D. HA, Generalization of the proximal point u1,li,,,;11i, SIAM J. Con-
trol Optim., 28 (1990) 503-512.

16] W. IT C:-.'ii SCH, .1,///Ii,, ;. methods and applications, Springer-Verlag,
Berlin, I''-,

17] W. W. HAGER, Convex control and dual approximations, part I, Con-
trol and Cybernetics, 8 (1979), 5-22.









[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. 189-202.

[19] W. W. HAGER, U1ii,,;,.i the inverse of a matrix, SIAM Review, 31
(lI' '), pp. 221-239.

[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. 137-142.

[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. 243-254.

[22] W. W. HAGER, Iterative methods for nearly .;,spii, linear systems,
SIAM J. Sci. Comput., 22 (2000), pp. 747-766.

[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. 263-275.

[24] W. W. HAGER, The dual active set iJii ;11i,, and the iterative solution
of linear pro'tr-am -. in Novel Approaches to Hard Discrete Optimization,
P. M. Pardalos and H. Wolkowicz, Eds., Fields Institute Communica-
tions, Vol. 37 (2003), 95-107.

[25] W. W. HAGER AND D. HEARN, The dual active set u1rilimi,,, applied
to quadratic networks, Computational Optimization and Applications,
1 (1993), pp. 349-373.

[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, CHUN-LIANG 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. 235-253.









[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' '),
359-392.

[32] G. KARYPIS AND V. KUMAR, Multilevel k-way pari.,.,,, scheme for
,, I,,l,, ,, i,. J. Parallel Distrib. Comput., 48 (1999), 96-129.

[33] G. KARYPIS AND V. KUMAR, Parallel multilevel k-way ilii,,,i;,,,,
scheme for ;,,i ,,, ri1ph. SIAM Review, 41 (1999), 278-300.

[34] J. W. H. LIu, A ,'I.1 ,1 envelope method for sparse factorization by
rows, AC'.I Transactions on Mathematical Software, 17 (1991) 112-129.

[35] F. J. LUQUE, Asymptotic (i,, I ,,i ,i analysis of the proximal point
(,l,,;i, ;i,, SIAM J. Control Optim., 22 (1984), pp. 277-293.

[36] B. MARTINET, !i ,',,,,;s l.,., d'inequations variationelles par ap-
proximations successives, Revue Francaise Informatique et Recherche
Operationnelle, 4 (1970) 154-159.

[37] B. MARTINET, Determination approachee d'un point fixe d'une appli-
cation pseudo-contractante, Comptes Rendus des Seances de l'Academie
des Sciences, Paris, 274 (1972) 163-165.

[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), 97-116.









[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. 26-2.

[43] C. ZHU, Asymptotic convergence analysis of some inexact proximal point
(,l,,7,i,/,,,s for minimization, SIAM J. Optim., 6 (1996), pp. 626-637.




University of Florida Home Page
© 2004 - 2010 University of Florida George A. Smathers Libraries.
All rights reserved.

Acceptable Use, Copyright, and Disclaimer Statement
Last updated October 10, 2010 - - mvs