Stabilization of Mehrotra's primal–dual algorithm and its implementation

Abstract

In this paper we apply a stabilization procedure proposed by Kovačević-Vujčić and Ašić to the Mehrotra's primal–dual interior-point algorithm for linear programming. Transformations of the dual problem corresponding to the stabilization procedure are considered. The stabilization procedure and Mehrotra's algorithm are implemented in the package MATHEMATICA. A number of highly degenerate test examples are used to compare the modified Mehrotra's method with respect to the original one and the codes PCx, LIPSOL and MOSEK. We also provide numerical results on some examples from the Netlib test set.

Introduction

Primal–dual algorithms are designed for linear programming problems in the standard form:

$$\text{minimize}\text{c}{}^{\mathrm{<mtext>T</mtext>}}\text{x}\text{subject}\text{to}\text{Ax=b,}\text{x⩾0,}$$

where $\text{c,x∈}\text{R}{}^{\mathrm{n}}$, $\text{b∈}\text{R}{}^{\mathrm{m}}$, A is an m×n real matrix and c^T is transpose of the vector c. The dual problem for (1.1) is

$$\text{maximize}\text{b}{}^{\mathrm{<mtext>T</mtext>}}\text{λ}\text{subject}\text{to}\text{A}{}^{\mathrm{<mtext>T</mtext>}}\text{λ+s=c,}\text{s⩾0,}$$

where $\text{λ∈}\text{R}{}^{\mathrm{m}}$ and $\text{s∈}\text{R}{}^{\mathrm{n}}$ and b^T, A^T denote transpose of the vector b and the matrix A, respectively. It is well-known that the vector $\text{x}{}^{\mathrm{\ast }}\text{∈}\text{R}{}^{\mathrm{n}}$ is a solution of (1.1) if and only if there exist vectors $\text{s}{}^{\mathrm{\ast }}\text{∈}\text{R}{}^{\mathrm{n}}$ and $\text{λ}{}^{\mathrm{\ast }}\text{∈}\text{R}{}^{\mathrm{m}}$ such that the following conditions hold:

$$\text{A}{}^{\mathrm{<mtext>T</mtext>}}\text{λ}{}^{\mathrm{\ast }}\text{+s}{}^{\mathrm{\ast }}\text{=c,}$$

$$\text{Ax}{}^{\mathrm{\ast }}\text{=b,}$$

$$\text{x}{}_{\mathrm{i}}{}^{\mathrm{\ast }}\text{s}{}_{\mathrm{i}}{}^{\mathrm{\ast }}\text{=0,}\text{i=1,…,n,}$$

$$\text{(x}{}^{\mathrm{\ast }}\text{,s}{}^{\mathrm{\ast }}\text{)⩾0.}$$

The feasible sets of , will be denoted by X and Y, respectively. Further, let X⁰={x|Ax=b,x>0},Y⁰={(λ,s)|A^Tλ+s=c,s>0}. As it is usual in the context of interior-point methods we shall assume that the sets X⁰ and Y⁰ are nonempty. This assumption implies that both sets $\text{X}{}^{\mathrm{\ast }}$ and $\text{Y}{}^{\mathrm{\ast }}$ are nonempty and bounded, where $\text{X}{}^{\mathrm{\ast }}$ and $\text{Y}{}^{\mathrm{\ast }}$ are the sets of optimal solutions to problems , , respectively.

All primal–dual methods generate iterates (x^t,λ^t,s^t) that satisfy the bounds (1.3d) strictly, that is, x^t>0 and s^t>0, and instead the condition (1.3c) deal with the condition x_i^ts_i^t=τ,i=1,…,n, where τ→0. Primal–dual strictly feasible set will be denoted by $\text{F}{}^0$, i.e. $\text{F}{}^0\text{={(x,λ,s)|Ax=b,A}{}^{\mathrm{<mtext>T</mtext>}}\text{λ+s=c,(x,s)>0}}$. It is well-known that there exist at least one primal solution $\text{x}{}^{\mathrm{\ast }}\text{∈X}{}^{\mathrm{\ast }}$ and one dual solution $\text{(λ}{}^{\mathrm{\ast }}\text{,s}{}^{\mathrm{\ast }}\text{)∈Y}{}^{\mathrm{\ast }}$ such that $\text{x}{}^{\mathrm{\ast }}\text{+s}{}^{\mathrm{\ast }}\text{>0}$ and such solutions are called strictly complementary solutions. Note that if (x,λ,s) and $\text{(x}{}^{\mathrm{\ast }}\text{,λ}{}^{\mathrm{\ast }}\text{,s}{}^{\mathrm{\ast }}\text{)}$ are two strictly complementary solution pairs, then x_i>0 if and only if $\text{x}{}_{\mathrm{i}}{}^{\mathrm{\ast }}\text{>0}$, i=1,…,n, and similarly, s_i>0 if and only if $\text{s}{}_{\mathrm{i}}{}^{\mathrm{\ast }}\text{>0}$, i=1,…,n, i.e. the sets of indices of positive coordinates are the same for all strictly complementary optimal pairs. Let $\text{B}\text{={j|x}{}_{\mathrm{j}}{}^{\mathrm{\ast }}\text{>0},}\text{N}\text{={j|s}{}_{\mathrm{j}}{}^{\mathrm{\ast }}\text{>0}}$. It is easy to see that $\text{B}\text{∩}\text{N}\text{=∅}$. The partition of the set {1,…,n} into $\text{B}$ and $\text{N}$ induces the partition of A into $\text{A}{}_{\mathrm{<mtext>B</mtext>}}$ and $\text{A}{}_{\mathrm{<mtext>N</mtext>}}$, where $\text{A}{}_{\mathrm{<mtext>B</mtext>}}\text{(A}{}_{\mathrm{<mtext>N</mtext>}}\text{)}$ is the submatrix of A with columns whose indices are in $\text{B}\text{(}\text{N}\text{)}$. Similar notations will be used for partitions of x, s and c. The sets $\text{X}{}^{\mathrm{\ast }}$ and $\text{Y}{}^{\mathrm{\ast }}$ can now be represented as

$$\text{X}{}^{\mathrm{\ast }}\text{={x∈}\text{R}{}^{\mathrm{n}}\text{|A}{}_{\mathrm{<mtext>B</mtext>}}\text{x}{}_{\mathrm{<mtext>B</mtext>}}\text{=b,x}{}_{\mathrm{<mtext>B</mtext>}}\text{⩾0,x}{}_{\mathrm{<mtext>N</mtext>}}\text{=0},}$$

$$\text{Y}{}^{\mathrm{\ast }}\text{={(λ,s)∈}\text{R}{}^{\mathrm{m}}\text{×}\text{R}{}^{\mathrm{n}}\text{|s}{}_{\mathrm{<mtext>B</mtext>}}\text{=c}{}_{\mathrm{<mtext>B</mtext>}}\text{−A}{}_{\mathrm{<mtext>B</mtext>}}{}^{\mathrm{<mtext>T</mtext>}}\text{λ=0,s}{}_{\mathrm{<mtext>N</mtext>}}\text{=c}{}_{\mathrm{<mtext>N</mtext>}}\text{−A}{}_{\mathrm{<mtext>N</mtext>}}{}^{\mathrm{<mtext>T</mtext>}}\text{λ⩾0}.}$$

According to [3], [9] the optimal face $\text{X}{}^{\mathrm{\ast }}$ is degenerate if $\text{rank}\text{(A}{}_{\mathrm{<mtext>B</mtext>}}\text{)<m}$.

Let us mention that interior-point methods require to solve at tth iteration a system of linear equations with the coefficient matrix of the form AD²_tA^T. One of the most serious sources of numerical difficulties is the fact that in the case when the optimal face $\text{X}{}^{\mathrm{\ast }}$ is degenerate the matrices involved tend to become extremely ill-conditioned as the optimal points are approached. There are several papers which address properties of the related system

$$\text{(AD}{}_{\mathrm{t}}{}^2\text{A}{}^{\mathrm{<mtext>T</mtext>}}\text{)u=AD}{}_{\mathrm{t}}\text{d.}$$

Its solution u=(AD_t²A^T)⁻¹AD_td has an interesting property that it is uniformly bounded in spite of possible ill-conditionedness of the matrix AD_t²A^T. The first result of that type was obtained by Dikin [6], who showed that

$$\text{∥u∥⩽}\text{max}\text{J∈Ω(A)}\text{∥(A}{}_{\mathrm{J}}{}^{\mathrm{<mtext>T</mtext>}}\text{)}{}^{\mathrm{-1}}\text{d}{}_{\mathrm{J}}\text{∥,}$$

where $\text{Ω(A)}$ is the set of column indices of all nonsingular m×m submatrices of A and A_J is a submatrix of A with columns whose indices are in J. This result has been rediscovered by Stewart [11] and Ben-Tal and Teboulle [4]. Stewart moreover proves that

$$\text{∥(AD}{}_{\mathrm{t}}{}^2\text{A}{}^{\mathrm{<mtext>T</mtext>}}\text{)}{}^{\mathrm{-1}}\text{AD}{}_{\mathrm{t}}\text{∥}$$

is uniformly bounded with respect to t and gives an explicit bound in terms of singular values of A. Nevertheless, the fact that solution is bounded does not imply that it can be computed in a stable way. Computational techniques most commonly used are Cholesky decomposition and preconditioned conjugate gradients. In presence of degeneracy both approaches have numerical problems. In the case of conjugate gradients the problem is how to choose a good preconditioner. Cholesky decomposition is usually stabilized by a heuristic rule that pivots smaller than a threshold value ε should be skipped. Under suitable assumptions Wright [16] shows that this heuristic approach has theoretical grounds, i.e. the distance between the exact and the approximate solution is bounded by $\text{ε}$. Nevertheless, in the papers [3], [9] it is shown that in the case of degeneracy, the interior-point methods can be seriously affected by ill-conditioning, even in their most robust implementations, such as PCx [5], [10], [15] and HOPDM [7]. Also, in [3], [9] a stabilization technique based on Gaussian elimination is presented which reformulates the original problem in an equivalent way such that ill-conditionedness is avoided. Our goal is to develop an implementation of this stabilization algorithm. The developed software solves a class of ill-conditioned test examples from [3], [9] which could not be solved by the codes PCx and HOPDM.

The paper is organized as follows. In Section 2 we restate the Mehrotra's primal–dual interior-point method [10] and the corresponding algorithm for the stabilization of the Mehrotra's method, established in [3], [9]. We also investigate transformations of the dual problem corresponding to the stabilization procedure.

In Section 3 we illustrate the effects of the stabilization procedure on a set of highly degenerate test problems. We also report results of the developed software on some of the Netlib test problems.