format: text

text:
# Effects of finite-precision arithmetic on interior-point methods for nonlinear programming

## 1 Introduction

We investigate the effects of finite-precision arithmetic on the calculated steps of primal-dual interior-point (PDIP) algorithms for the nonlinear programming problem

$$\text{NLP:}\underset{z}{\mathrm{min}}\varphi (z)\text{subject to }g(z)0,$$
(1)
where $`\varphi :\text{I R }^n\text{I R }`$ and $`g:\text{I R }^n\text{I R }^m`$ are twice Lipschitz continuously differentiable functions. Optimality conditions for this problem can be derived from the Lagrangian function $`(z,\lambda )`$, which is defined as

$$(z,\lambda )=\varphi (z)+\underset{i=1}{\overset{m}{}}\lambda _ig_i(z)=\varphi (z)+\lambda ^Tg(z),$$
(2)
where $`\lambda \text{I R }^m`$ is a vector of Lagrange multipliers. When a constraint qualification (discussed below) holds at the point $`z^{}`$, first-order necessary conditions for $`z^{}`$ to be a solution of (1) are that there exists a vector of Lagrange multipliers $`\lambda ^{}\text{I R }^m`$ such that the following conditions are satisfied for $`(z,\lambda )=(z^{},\lambda ^{})`$:

$$_z(z,\lambda )=\varphi (z)+g(z)\lambda =0,g(z)0,\lambda 0,\lambda ^Tg(z)=0,$$
(3)
where

$$g(z)=[g_1(z),g_2(z),\mathrm{},g_m(z)].$$

The conditions (3) are the well-known Karush-Kuhn-Tucker (KKT) conditions. We use $`𝒮_\lambda `$ to denote the set of vectors $`\lambda ^{}`$ such that $`(z^{},\lambda ^{})`$ satisfies (3). The primal-dual solution set is defined by

$$𝒮=\{z^{}\}\times 𝒮_\lambda .$$
(4)
This paper discusses local convergence of PDIP algorithms for (1), assuming that the algorithm is implemented on a computer that performs calculations according to the standard model of floating-point arithmetic. Because of our focus on local convergence properties, we assume throughout that the current iterate $`(z,\lambda )`$ is close enough to the solution set $`𝒮`$ that superlinear convergence would occur if exact steps (uncorrupted by finite precision) were taken. In the interests of generality, we weaken an assumption that is often made in the analysis of algorithms for (1), namely, that the gradients of the active constraints are linearly independent at the solution. We replace this linear independence constraint qualification (LICQ) with the weaker Mangasarian-Fromovitz constraint qualification (MFCQ) . MFCQ allows constraint gradients to become dependent at the solution, so that the set $`𝒮_\lambda `$ of optimal Lagrange multipliers is no longer necessarily a singleton, though it remains bounded. We continue to assume that a strict complementarity (SC) condition holds, that is,

$$g_i(z^{})=0\lambda _i^{}>0,\text{for some }\lambda ^{}𝒮_\lambda .$$
(5)
In the context of rapidly convergent algorithms, the SC condition makes good sense. If SC fails to hold, superlinear convergence of Newton-like algorithms does not occur, except for specially modified algorithms such as those that identify the active constraints explicitly (see Monteiro and Wright  and El-Bakry, Tapia, and Zhang ).

The major conclusion of the paper is that the effects of roundoff errors on the rapid local convergence of the algorithm are fairly benign. When a standard second-order condition is added to the assumptions already mentioned, the steps produced under floating-point arithmetic approach $`𝒮`$ almost as effectively as do exact steps, as long as the distance to the solution set remains significantly greater than the unit roundoff $`𝐮`$. The latter condition is hardly restrictive, since the data errors made in storing the problem in a digital computer mean that the solution set is known only to within some multiple of $`𝐮`$ in any case.

The conclusions about the effectiveness of the computed steps are not obvious, because all three formulations of the linear system that must be solved to compute the step at each iteration may become highly ill conditioned near the solution. Our analysis would be significantly simpler if we were to make the LICQ assumption because, in this case, one formulation of the linear equations remains well conditioned, and stability of the three standard formulations can be proved by exploiting their relationship to this system of equations.

This work is related to earlier work of the author on finite-precision analysis of interior-point algorithms for linear complementarity problems  and linear programming . The existence of second-order effects gives the analysis here a somewhat different flavor, however. In addition, we go into more depth in checking that the computed iterates can continue to satisfy the approximate centrality conditions usually required in primal-dual algorithms, and in deriving expressions for the rate at which the computed iterates approach the solution set. Related work by Forsgren, Gill, and Shinnerl  deals with one formulation of the step equations for the nonlinear programming problem—the so-called augmented form treated here in Section 6—but makes assumptions on the pivot sequence that do not always hold in practice. M. H. Wright  recently presented an analysis of the condensed form of the step equations discussed in Section 5 under the assumption that LICQ holds, and found that the computed steps were more accurate than would be expected from a naive analysis.

For linear programming, the PDIP approach has emerged as the most powerful of the interior-point approaches. The supporting theory is strong, in terms of global and local convergence analysis and complexity theory (see the bibliography of Wright ). Implementations yield better results than pure-primal or barrier-function approaches; see Andersen et al. . Strong theory is also available for these algorithms when applied to convex programming, in which $`\varphi ()`$ and $`g_i()`$, $`i=1,\mathrm{},m`$ are all convex functions; see, for example, Wright and Ralph  and Ralph and Wright . The latter paper drops the LICQ assumption in favor of MFCQ, making the local theory stronger in one sense than the corresponding local theory for the sequential quadratic programming (SQP) algorithm. The use of MFCQ complicates the analysis considerably, however; under LICQ, the implicit function theorem can be used to prove a key technical result about the length of the step, while more complicated logic is needed to derive this same result under MFCQ.

A significant by-product of the current paper is to prove the key technical result about the length of the rapidly convergent step (Corollary 8) under MFCQ and SC, even when the problem (1) is not convex. This allows the local convergence results of Ralph and Wright  to be extended to general nonconvex nonlinear problems.

The analysis of this paper could also be applied to the recently proposed stabilized sequential quadratic programming (sSQP) algorithm (see Wright  and Hager ), in which small penalties on the change in the multiplier estimate $`\lambda `$ from one iteration to the next ensure rapid convergence even when LICQ is relaxed to MFCQ. A finite-precision analysis of the sSQP method appears in \[29, Section 3.2\], but only for the augmented form of the step equations. Analysis quite similar to that of the current paper could be applied to show that similar conclusions continue to hold when a condensed form of the step equations is used instead. We omit the details.

The remainder of this paper is structured in the following way. Section 2 contains notation, together with our basic assumptions about (1) and some relevant results from the literature. Section 3 discusses the primal-dual interior-point framework, defining the general form of each iteration and the step equations that must be solved at each iteration. Subsection 3.2 proves an important technical result about the relationship between the distance of an interior-point iterate to the solution set $`𝒮`$ and a duality measure $`\mu `$. Section 4 describes perturbed variants of the linear systems that are solved to obtain PDIP steps, and proves our key results about the effect of the perturbations on the accuracy of the steps.

Section 5 focuses on one form of the PDIP step equations: the most compact form in which most of the computational effort goes into factoring a symmetric positive definite matrix, usually by a Cholesky procedure. We trace the effect on step accuracy of errors in evaluation of the functions, formation of the system, and the factorization/solution process. Further, we show the effects of these inaccuracies on the distance that we can move along the steps before the interiority condition is violated, and on various measures of algorithmic progress. An analogous treatment of the augmented form of the step equations appears in Section 6. The conclusions of this section depend on the actual algorithm used to solve the augmented system—it is not sufficient to assume, as in Section 5, that any backward-stable procedure is used to factor the matrix. (We note that similar results hold for the full form of the step equations, but we omit the details of this case, which can be found in the technical report .) We conclude with a numerical illustration of the main results in Section 7 and summarize the paper in Section 8.

## 2 Notation, Assumptions, and Basic Results

We use $``$ to denote the set of active indices at $`z^{}`$, that is,

$$=\{i=1,2,\mathrm{},m|g_i(z^{})=0\},$$
(6)
whereas $`𝒩`$ denotes its complement

$$𝒩=\{1,2,\mathrm{},m\}\backslash .$$
(7)
The set $`_+`$ is defined as

$$_+=\{i|\lambda _i^{}>0\text{for some }\lambda ^{}\text{ satisfying (}\text{3}\text{)}\}.$$
(8)
The strict complementarity condition (5) is equivalent to

$$_+=.$$
(9)
We frequently make reference to submatrices and subvectors corresponding to the index sets $``$ and $`𝒩`$. For example, the quantities $`\lambda _{}`$ and $`g_{}(z)`$ are the vectors containing the components $`\lambda _i`$ and $`g_i(z)`$, respectively, for $`i`$, while $`g_{}(z)`$ is the matrix whose columns are $`g_i(z)`$, $`i`$.

The Mangasarian-Fromovitz constraint qualification (MFCQ) is satisfied at $`z^{}`$ if there is a vector $`\overline{y}\text{I R }^n`$ such that

$$g_{}(z^{})^T\overline{y}<0.$$
(10)
The following fundamental result about MFCQ is due to Gauvin .

###### Lemma 1.

Suppose that the first-order conditions (3) are satisfied at $`z=z^{}`$. Then $`𝒮_\lambda `$ is bounded if and only if the MFCQ condition (10) is satisfied at $`z^{}`$.

This result is crucial because it allows our (local) analysis to place a uniform bound on all $`\lambda `$ in a neighborhood of the dual solution set $`𝒮_\lambda `$.

The second-order condition used in most of the remainder of the paper is that there is a constant $`\xi >0`$ such that

$$w^T_{zz}(z^{},\lambda ^{})w\xi w^2,$$
(11)
for all $`\lambda ^{}𝒮_\lambda `$ and all $`w`$ satisfying

$$\begin{array}{cc}g_i(z^{})^Tw=0,\hfill & \text{for all }i_+,\hfill \\ g_i(z^{})^Tw0,\hfill & \text{for all }i\backslash _+.\hfill \end{array}$$
(12)
When the SC condition (5) (alternatively, (9)) is satisfied, this direction set is simply $`\text{null}g_{}(z^{})^T`$.

A simple example that satisfies MFCQ but not LICQ at the solution, and that satisfies the second-order conditions (11), (12) and the SC condition is as follows:

$$\underset{z\text{I R }^2}{\mathrm{min}}z_1\text{subject to}(z_11/3)^2+z_2^21/9,(z_12/3)^2+z_2^24/9.$$
(13)
The solution is $`z^{}=0`$, and the optimal multiplier set is

$$𝒮_\lambda =\{\lambda 0|\mathrm{\hspace{0.17em}2}\lambda _1+4\lambda _2=3\}.$$
(14)
The gradients of the two constraints are the solution are $`(2/3,0)^T`$ and $`(4/3,0)^T`$, respectively. They are linearly dependent, but the MFCQ condition (10) can be satisfied by choosing $`\overline{y}=(1,0)^T`$.

We use $`𝐮`$ to denote the unit roundoff, which we define as the smallest number such that the following property holds: When $`x`$ and $`y`$ are any two floating-point numbers, op denotes $`+`$, $``$, $`\times `$, or $`/`$, and $`\text{fl}(z)`$ denotes the floating-point approximation of a real number $`z`$, we have

$$\text{fl}(x\text{op}y)=(x\text{op}y)(1+ϵ),|ϵ|𝐮.$$
(15)
Modest multiples of $`𝐮`$ are denoted by $`\delta _𝐮`$.

We assume that the problem is scaled so that the values of $`g`$ and $`\varphi `$ and their first and second derivatives in the vicinity of the solution set $`𝒮`$, and the values $`(z,\lambda )`$ themselves, can all be bounded by moderate quantities. When multiplied by $`𝐮`$, quantities of this type are absorbed into the notation $`\delta _𝐮`$ in the analysis below.

Order notation $`O()`$ and $`\mathrm{\Theta }()`$ is used as follows: If $`v`$ (vector or scalar) and $`ϵ`$ (nonnegative scalar) are two quantities that share a dependence on other variables, we write $`v=O(ϵ)`$ if there is a moderate constant $`\beta _1`$ such that $`v\beta _1ϵ`$ for all values of $`ϵ`$ that are interesting in the given context. (The “interesting context” frequently includes cases in which $`ϵ`$ is either sufficiently small or sufficiently large, but we often use $`v=O(\mu )`$ to indicate that $`v\beta _1\mu `$ for all sufficiently small $`\mu `$ that satisfy $`\mu 𝐮`$, for some $`\beta _1`$; see later discussion.) We write $`v=\mathrm{\Theta }(ϵ)`$ if there are constants $`\beta _1`$ and $`\beta _0`$ such that $`\beta _0ϵv\beta _1ϵ`$ for all interesting values of $`ϵ`$. Similarly, we write $`v=O(1)`$ if $`v\beta _1`$, and $`v=\mathrm{\Theta }(1)`$ if $`\beta _0v\beta _1`$.

We use the notation $`\delta (z,\lambda )`$ to denote the distance from $`(z,\lambda )`$ to the primal-dual solution set, that is,

$$\delta (z,\lambda )\stackrel{\mathrm{def}}{=}\underset{(z^{},\lambda ^{})𝒮}{\mathrm{min}}(z,\lambda )(z^{},\lambda ^{}).$$
(16)
It is well known (see, for example, Theorem A.1 of Wright ) that this distance can be estimated in terms of known quantities at $`(z,\lambda )`$. We state this result formally as follows.

###### Theorem 2.

Suppose that the first-order conditions (3), the MFCQ condition (10) and the second-order conditions (11), (12) are satisfied at the solution $`z^{}`$. Then if $`\lambda 0`$, we have

$$\delta (z,\lambda )=\mathrm{\Theta }\left(\left[\begin{array}{c}_z(z,\lambda )\\ \mathrm{min}(\lambda ,g(z))\end{array}\right]\right).$$
(17)
We write the singular value decomposition (SVD) of the matrix $`g_{}(z^{})`$ of first partial derivatives as follows:

$$g_{}(z^{})=\left[\begin{array}{cc}\widehat{U}& \widehat{V}\end{array}\right]\left[\begin{array}{cc}\mathrm{\Sigma }& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}U^T\\ V^T\end{array}\right]=\widehat{U}\mathrm{\Sigma }U^T,$$
(18)
where the matrices $`\left[\begin{array}{cc}\widehat{U}& \widehat{V}\end{array}\right]`$ and $`\left[\begin{array}{cc}U& V\end{array}\right]`$ are orthogonal, and $`\mathrm{\Sigma }`$ is a diagonal matrix with positive diagonal elements.

Note that the columns of $`\widehat{U}`$ form a basis for the range space of $`g_{}(z^{})`$, while the columns of $`\widehat{V}`$ form a basis for the null space of $`g_{}(z^{})^T`$. Similarly, the columns of $`U`$ form a basis for the range space of $`g_{}(z^{})^T`$, while the columns of $`V`$ form a basis for the null space of $`g_{}(z^{})`$. These four subspaces are key to our analysis, particularly the subspace spanned by the columns of $`V`$. For the computational methods used to solve the primal-dual step equations discussed in this paper, the computed step in the $``$-components of the multipliers—that is, $`\mathrm{\Delta }\lambda _{}`$—has a larger error in the range space of $`V`$ than in the complementary subspace spanned by the columns of $`U`$. The errors in the computed primal step $`\mathrm{\Delta }z`$, the computed step of the $`𝒩`$-components of the multipliers $`\lambda _𝒩`$, and the computed step in the dual slack variables (defined later) are typically also less significant than the error in $`V^T\mathrm{\Delta }\lambda _{}`$. We show, however, that the potentially large error in $`V^T\mathrm{\Delta }\lambda _{}`$ does not affect the performance of primal-dual algorithms that use these computed steps until $`\mu `$ becomes similar to $`𝐮^{1/2}`$.

When the stronger LICQ condition holds, the matrix $`V`$ is vacuous, and the SVD (18) reduces to $`g_{}(z^{})=\widehat{U}\mathrm{\Sigma }U^T`$. Much of the analysis in this paper would simplify considerably under LICQ, in part because $`V^T\mathrm{\Delta }\lambda _{}`$—the step component with the large error—is no longer present.

We use $`\sigma _{\mathrm{min}}()`$ to denote the smallest eigenvalue, and $`\mathrm{cond}()`$ to denote the condition number, as measured by the Euclidean norm.

## 3 Primal-Dual Interior-Point Methods

### 3.1 Centrality Conditions and Step Equations

Primal-dual interior-point methods are constrained, modified Newton methods applied to a particular form of the KKT conditions (3). By introducing a vector $`s\text{I R }^m`$ of slacks for the inequality constraint, we can rewrite the nonlinear program as

$$\underset{(z,s)}{\mathrm{min}}\varphi (z)\text{subject to }g(z)+s=0\text{}s0\text{,}$$

and the KKT conditions (3) as

$$_z(z,\lambda )=0,g(z)+s=0,(\lambda ,s)0,\lambda ^Ts=0.$$
(19)
Motivated by this form of the conditions, we define the mapping $`(z,\lambda ,s)`$ by

$$(z,\lambda ,s)\stackrel{\mathrm{def}}{=}\left[\begin{array}{c}\varphi (z)+g(z)\lambda \\ g(z)+s\\ S\mathrm{\Lambda }e\end{array}\right],$$
(20)
where the diagonal matrices $`S`$ and $`\mathrm{\Lambda }`$ are defined by

$$S\stackrel{\mathrm{def}}{=}\text{diag}(s_1,s_2,\mathrm{},s_m),\mathrm{\Lambda }\stackrel{\mathrm{def}}{=}\text{diag}(\lambda _1,\lambda _2,\mathrm{},\lambda _m),$$

and $`e`$ is defined as

$$e=(1,1,\mathrm{},1)^T.$$
(21)
The KKT conditions (19) can now be stated equivalently as

$$(z,\lambda ,s)=0,(s,\lambda )0.$$
(22)
Primal-dual iterates $`(z,\lambda ,s)`$ invariably satisfy the strict bound $`(s,\lambda )>0`$, while they approach satisfaction of the condition $`()=0`$ in the limit. An important measure of progress is the duality measure $`\mu (\lambda ,s)`$, which is defined by

$$\mu (\lambda ,s)\stackrel{\mathrm{def}}{=}\lambda ^Ts/m.$$
(23)
When $`\mu `$ is used without arguments, we assume that $`\mu =\mu (\lambda ,s)`$, where $`(z,\lambda ,s)`$ is the current primal-dual iterate. We emphasize that $`\mu `$ is a function of $`(z,\lambda ,s)`$, rather than a target value explicitly chosen by the algorithm, as is the case in some of the literature.

A typical step $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda ,\mathrm{\Delta }s)`$ of the primal-dual method satisfies

$$(z,\lambda ,s)\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda \\ \mathrm{\Delta }s\end{array}\right]=(z,\lambda ,s)\left[\begin{array}{c}0\\ 0\\ t\end{array}\right],$$
(24)
where $`t`$ defines the deviation from a pure Newton step for $``$ (which is also known as a “primal-dual affine-scaling” step). The vector $`t`$ frequently contains a centering term $`\sigma \mu e`$, where $`\sigma `$ is a centering parameter in the range $`[0,1]`$. It sometimes also contains higher-order information, such as the product $`\mathrm{\Delta }\mathrm{\Lambda }_{\mathrm{aff}}\mathrm{\Delta }S_{\mathrm{aff}}e`$, where $`\mathrm{\Delta }\mathrm{\Lambda }_{\mathrm{aff}}`$ and $`\mathrm{\Delta }S_{\mathrm{aff}}`$ are the diagonal matrices constructed from the components of the pure Newton step (Mehrotra ). In any case, the vector $`t`$ usually goes to zero rapidly as the iterates converge to a solution, so that the steps generated from (24) approach pure Newton steps, which in turn ensures rapid local convergence. Throughout this paper, we assume that $`t`$ satisfies the estimate

$$t=O(\mu ^2).$$
(25)
All our major results continue to hold, with slight modification, if we replace (25) by $`t=O(\mu ^\sigma )`$, for some $`\sigma (1,2]`$. Our essential point remains unchanged; the theoretical superlinear convergence rate promised by this choice of $`t`$ is not seriously compromised by roundoff errors as long as $`\mu `$ remains significantly larger than the unit roundoff $`𝐮`$. To avoid notational clutter, however, we analyze only the case (25).

Using the definition (2), we can write the system (24) explicitly as follows:

$$\left[\begin{array}{ccc}_{zz}(z,\lambda )& g(z)& 0\\ g(z)^T& 0& I\\ 0& S& \mathrm{\Lambda }\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda \\ \mathrm{\Delta }s\end{array}\right]=\left[\begin{array}{c}_z(z,\lambda )\\ g(z)+s\\ S\mathrm{\Lambda }e+t\end{array}\right].$$
(26)
Block eliminations can be performed on this system to yield more compact formulations. By eliminating $`\mathrm{\Delta }s`$, we obtain the augmented system form, which is

$$\left[\begin{array}{cc}_{zz}(z,\lambda )& g(z)\\ g(z)^T& \mathrm{\Lambda }^1S\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda \end{array}\right]=\left[\begin{array}{c}_z(z,\lambda )\\ g(z)+\mathrm{\Lambda }^1t\end{array}\right].$$
(27)
By eliminating $`\mathrm{\Delta }\lambda `$ from this system, we obtain a system that is sometimes referred to as the condensed form (or in the case of linear programming as the normal equations form), which is

$`\left[_{zz}(z,\lambda )+g(z)\mathrm{\Lambda }S^1g(z)^T\right]\mathrm{\Delta }z`$
$`=`$ $`_z(z,\lambda )g(z)\mathrm{\Lambda }S^1[g(z)\mathrm{\Lambda }^1t].`$
We consider primal-dual methods in which each iterate $`(z,\lambda ,s)`$ satisfies the following properties:

$`r_f(z,\lambda )C\mu ,\text{where}r_f(z,\lambda )\stackrel{\mathrm{def}}{=}_z(z,\lambda ),`$ (29a)
$`r_g(z,s)C\mu ,\text{where}r_g(z,s)\stackrel{\mathrm{def}}{=}g(z)+s,`$ (29b)
$`(\lambda ,s)>0,\lambda _is_i\gamma \mu ,\text{for all }i=1,2,\mathrm{},m\text{,}`$ (29c)
for some constants $`C>0`$ and $`\gamma (0,1)`$, where $`\mu `$ is defined as in (23). (In much of the succeeding discussion, we omit the arguments from the quantities $`\mu `$, $`r_f`$, and $`r_g`$ when they are evaluated at the current iterate $`(z,\lambda ,s)`$.) These conditions ensure that the pairwise products $`s_i\lambda _i`$, $`i=1,2,\mathrm{},m`$ are not too disparate and that the first two components of $``$ in (20) can be bounded in terms of the third component. They are sometimes called the centrality conditions because they are motivated by the notion of a central path and its neighborhoods. Conditions of the type (29) are imposed in most path-following interior-point methods for linear programming (see, for example, ). For nonlinear convex programming, examples of methods that require these conditions can be found in Ralph and Wright . In nonlinear programming, we mention Gould et al.  (see Algorithm 4.1 and Figure 5.1) and Byrd, Liu, and Nocedal . In the latter paper, (29a) and (29b) are imposed explicitly, while (29c) can be guaranteed by choosing $`ϵ_\mu =(1\gamma )\mu `$. Even when the choice $`ϵ_\mu =\mu `$ is made, as in the bulk of the discussion in , their other conditions concerning positivity of $`(s,\lambda )`$ can be expected to produce iterates that satisfy (29c) in practice.

For points $`(z,\lambda ,s)`$ that satisfy (29), we can use $`\mu `$ to estimate the distance $`\delta (z,\lambda )`$ from $`(z,\lambda )`$ to the solution set $`𝒮`$ (see (16)). These results, which are proved in the following subsection, can be summarized briefly as follows. When the MFCQ condition (10) and the second-order conditions (11), (12) are satisfied, we have that $`\delta (z,\lambda )=O(\mu ^{1/2})`$. When the strict complementarity assumption (5) is added, we obtain the stronger estimate $`\delta (z,\lambda )=O(\mu )`$. We can use these estimates to obtain bounds on the elements of the diagonal matrices $`S`$, $`\mathrm{\Lambda }`$, $`S^1\mathrm{\Lambda }`$, and $`\mathrm{\Lambda }^1S`$ in the systems above; these bounds are the key to the error analysis of the remainder of the paper.

### 3.2 Using the Duality Measure to Estimate Distance to the Solution

The main result of this section, Theorem 5, shows that under certain assumptions, the distance $`\delta (z,\lambda )`$ of a primal-dual iterate $`(z,\lambda ,s)`$ to the solution set $`𝒮`$ can be estimated by the duality measure $`\mu `$. We start with a technical lemma that proves the weaker estimate $`\delta (z,\lambda )=O(\mu ^{1/2})`$. Note that this result does not assume that the SC condition (5) holds.

###### Lemma 3.

Suppose that $`z^{}`$ is a solution of (1) at which the MFCQ condition (10) and the second-order conditions (11), (12) are satisfied. Then for all $`(z,\lambda )`$ with $`\lambda 0`$ for which there is a vector $`s`$ such that $`(z,\lambda ,s)`$ satisfies (29), we have that

$$\delta (z,\lambda )=O(\mu ^{1/2}).$$
(30)
###### Proof.

We prove the result by showing that

$$\left[\begin{array}{c}_z(z,\lambda )\\ \mathrm{min}(\lambda ,g(z))\end{array}\right]=O(\mu ^{1/2})$$
(31)
and then applying Theorem 2. Since $`_z(z,\lambda )=r_f=O(\mu )`$, the first part of the vector satisfies the required estimate. For the second part, we have from (29b) that

$$g(z)=sr_g=s+O(\mu ),$$

and hence that

$$\mathrm{min}(g_i(z),\lambda _i)=\mathrm{min}(s_i,\lambda _i)+O(\mu ).$$
(32)
Because of (23) and (29c), we have that $`s_i\lambda _im\mu `$ and $`(\lambda _i,s_i)>0`$. It follows immediately that $`\mathrm{min}(\lambda _i,s_i)(m\mu )^{1/2}`$ for $`i=1,2,\mathrm{},m`$. Hence, by substitution into (32), we obtain

$$\mathrm{min}(g_i(z),\lambda _i)(m\mu )^{1/2}+O(\mu )=O(\mu ^{1/2}).$$

We conclude that the second part of the vector in (31) is of size $`O(\mu ^{1/2})`$, so the proof is complete. ∎

The following examples show the upper bound of Lemma 3 is indeed achieved and that it is not possible to obtain a lower bound on $`\delta (z,\lambda )`$ as a strictly increasing nonnegative function of $`\mu `$. To demonstrate the first claim, consider the problem

$$\mathrm{min}\frac{1}{2}z^2\text{subject to }z0\text{.}$$

The point $`(z,\lambda ,s)=(ϵ,ϵ,ϵ)`$ satisfies

$$_z(z,\lambda )=0,g(z)+s=0,s\lambda =ϵ^2,\mu =ϵ^2,$$

so that the conditions (29) are satisfied. Clearly the distance from the point $`(z,\lambda )`$ to the solution set $`𝒮=(0,0)`$ is $`\sqrt{2}ϵ=\sqrt{2}\mu ^{1/2}`$. For the second claim, consider any nonlinear program such that $`=\{1,2,\mathrm{},m\}`$ (that is, all constraints active) and strict complementarity (5) holds at some multiplier $`\lambda ^{}`$. Then for appropriate choices of $`\gamma `$ and $`C`$, the point

$$(z,\lambda ,s)=(z^{},\lambda ^{},(m\mu )/(e^T\lambda ^{})e)$$
(33)
satisfies the definition (23) and the condition (29) for any $`\mu >0`$. On the other hand, we have $`\delta (z,\lambda )=\delta (z^{},\lambda ^{})=0`$ by definition, so there are no $`\beta >0`$ and $`\sigma >0`$ that yield a lower bound estimate of the form $`\delta (z,\lambda )\beta \mu ^\sigma `$.

We now prove an extension of Lemma 5.1 of Ralph and Wright , dropping the monotonicity assumption of this earlier result.

###### Lemma 4.

Suppose that the conditions of Lemma 3 hold and in addition that the SC condition (5) is satisfied. Then for all $`(z,\lambda ,s)`$ satisfying (29), we have that

$`i`$ $``$ $`s_i=\mathrm{\Theta }(\mu ),\lambda _i=\mathrm{\Theta }(1),`$ (34a)
$`i𝒩`$ $``$ $`s_i=\mathrm{\Theta }(1),\lambda _i=\mathrm{\Theta }(\mu ).`$ (34b)
###### Proof.

By boundedness of $`𝒮`$ (Lemma 1), we have for all $`(z,\lambda ,s)`$ sufficiently close to $`𝒮`$ that

$$\lambda _i=O(1),s_i=g_i(z)+(r_g)_i=O(1).$$
(35)
Given $`(z,\lambda ,s)`$ satisfying (29), let $`P(\lambda )`$ be the projection of $`\lambda `$ onto the set $`𝒮_\lambda `$, and let $`\lambda ^{}𝒮_\lambda `$ be some strictly complementary optimal multiplier (for which (5) is satisfied). From Lemma 3 we obtain

$$zz^{}=O(\mu ^{1/2}).$$
(36)
Using this observation together with smoothness of $`\varphi ()`$ and $`g()`$, we have for the gradient of $``$ that

$`_z(z,\lambda )_z(z^{},\lambda ^{})`$
$`=`$ $`\varphi (z)\varphi (z^{})+g(z)\lambda g(z^{})\lambda ^{}`$
$`=`$ $`O(\mu ^{1/2})+g(z)[\lambda P(\lambda )]+[g(z)g(z^{})]P(\lambda )+g(z^{})[P(\lambda )\lambda ^{}].`$
Since $`P(\lambda )`$ and $`\lambda ^{}`$ are both in $`𝒮_\lambda `$, we find from (3) that the last term vanishes. From (36) and $`P(\lambda )=O(1)`$, the second-to-last term has size $`O(\mu ^{1/2})`$. For the remaining term, we have $`g(z)=O(1)`$, and $`\lambda P(\lambda )\delta (z,\lambda )=O(\mu ^{1/2})`$. By assembling all these observations, and using $`_z(z^{},\lambda ^{})=0`$, we obtain

$$_z(z,\lambda )=_z(z,\lambda )_z(z^{},\lambda ^{})=O(\mu ^{1/2}).$$
(37)
Using again that $`g(z^{})[P(\lambda )\lambda ^{}]=0`$, we have from (36) that

$`[P(\lambda )\lambda ^{}]^T[g(z)g(z^{})]`$ $`=`$ $`[P(\lambda )\lambda ^{}]^T[g(z^{})^T(zz^{})+O(zz^{}^2)]`$ (38)
$`=`$ $`O(zz^{}^2)=O(\mu ).`$
By gathering the estimates (30), (36), (37), and (38), we obtain

$`\left[\begin{array}{c}zz^{}\\ \lambda \lambda ^{}\end{array}\right]^T\left[\begin{array}{c}_z(z,\lambda )_z(z^{},\lambda ^{})\\ g(z)+g(z^{})\end{array}\right]`$ (43)
$`=`$ $`\left[\begin{array}{c}zz^{}\\ \lambda P(\lambda )\end{array}\right]^T\left[\begin{array}{c}_z(z,\lambda )_z(z^{},\lambda ^{})\\ g(z)+g(z^{})\end{array}\right]`$
$`+[P(\lambda )\lambda ^{}]^T[g(z)+g(z^{})]`$
$`=`$ $`O(\delta (z,\lambda ))O(\mu ^{1/2})+O(\mu )=O(\mu ).`$ (49)
By substituting from (29) and using (49), we obtain

$$\left[\begin{array}{c}zz^{}\\ \lambda \lambda ^{}\end{array}\right]^T\left[\begin{array}{c}r_f\\ sr_gs^{}\end{array}\right]=\left[\begin{array}{c}zz^{}\\ \lambda \lambda ^{}\end{array}\right]^T\left[\begin{array}{c}_z(z,\lambda )_z(z^{},\lambda ^{})\\ g(z)+g(z^{})\end{array}\right]=O(\mu ),$$

and therefore

$$(\lambda \lambda ^{})^T(ss^{})=(zz^{})^Tr_f+(\lambda \lambda ^{})^Tr_g+O(\mu ).$$

By using the conditions (29a), (29b), and the definition (23), we obtain

$`{\displaystyle \underset{i=1}{\overset{m}{}}}\lambda _i^{}s_i{\displaystyle \underset{i=1}{\overset{m}{}}}\lambda _is_i^{}`$
$`=`$ $`(\lambda ^{})^Ts\lambda ^Ts^{}=\lambda ^Ts+O(\mu )+O(zz^{}r_f)+O(\lambda \lambda ^{}r_g)=O(\mu ).`$
Since $`(\lambda ,s)>0`$ and $`(\lambda ^{},s^{})0`$, all terms $`\lambda _i^{}s_i`$ and $`\lambda _is_i^{}`$, $`i=1,2,\mathrm{},m`$ are nonnegative, so there is a constant $`C_1>0`$ such that

$$0\lambda _i^{}s_iC_1\mu ,0\lambda _is_i^{}C_1\mu ,\text{for all }i=1,2,\mathrm{},m\text{.}$$

For all $`i`$, we have $`\lambda _i^{}>0`$ by our strictly complementary choice of $`\lambda ^{}`$, and so

$$0<s_i\frac{C_1}{\lambda _i^{}}\mu \frac{C_1}{\mathrm{min}_i\lambda _i^{}}\mu \stackrel{\mathrm{def}}{=}C_2\mu .$$
(50)
On the other hand, we have by boundedness of $`𝒮_\lambda `$ and our assumption (29c) that

$$s_i\frac{\gamma \mu }{\lambda _i}\gamma _{\mathrm{min}}\mu ,\text{for all }i=1,2,\mathrm{},m,$$
(51)
for some constant $`\gamma _{\mathrm{min}}>0`$. By combining (50) and (51), we have that

$$s_i=\mathrm{\Theta }(\mu ),\text{for all }i.$$

For the $`\lambda _{}`$ component, we have that

$$s_i\lambda _i\gamma \mu \lambda _i\frac{\gamma \mu }{s_i}\frac{\gamma }{C_2},\text{for all }i.$$

Hence, by combining with (35), we obtain that

$$\lambda _i=\mathrm{\Theta }(1),\text{for all }i.$$

This completes the proof of (34a). We omit the proof of (34b), which is similar. ∎

Next, we show that when the strict complementarity assumption is added to the assumptions of Lemma 3, the upper bound on the distance to the solution set in (30) can actually be improved to $`O(\mu )`$.

###### Theorem 5.

Suppose that $`z^{}`$ is a solution of (1) at which the MFCQ condition (10), the second-order conditions (11), (12), and the SC condition (5) are satisfied. Then for all $`(z,\lambda )`$ with $`\lambda 0`$ for which there is a vector $`s`$ such that $`(z,\lambda ,s)`$ satisfies (29), we have that

$$\delta (z,\lambda )=O(\mu ).$$
(52)
###### Proof.

From (29a), we have directly that $`r_f=O(\mu )`$. We have from (29) and (34a) that

$$g_i(z)=s_i+(r_g)_i=O(\mu ),\lambda _i=\mathrm{\Theta }(1),\lambda _i>0\text{for all }i,$$

so that

$$\mathrm{min}(g_i(z),\lambda _i)=g_i(z)=O(\mu ),\text{for all }i\text{,}$$
(53)
whenever $`\mu `$ is sufficiently small. For the remaining components, we have

$$\mathrm{min}(g_i(z),\lambda _i)=\lambda _i=O(\mu ),\text{for all }i𝒩\text{.}$$
(54)
By substituting (29a), (53), and (54) into (17), we obtain the result. ∎

Similar conclusions to Lemma 4 and Theorem 5 can be reached in the case of linear programming algorithms. The second-order conditions (11), (12) are not relevant for this class of problems, and the SC assumption (5) holds for every linear program that has a solution.

## 4 Accuracy of PDIP Steps: General Results

By partitioning the constraint index set $`\{1,2,\mathrm{},m\}`$ into active indices $``$ and inactive indices $`𝒩`$, we can express the system (27) without loss of generality as follows:

$$\left[\begin{array}{ccc}_{zz}(z,\lambda )& g_{}(z)& g_𝒩(z)\\ g_{}(z)^T& D_{}& 0\\ g_𝒩(z)^T& 0& D_𝒩\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda _{}\\ \mathrm{\Delta }\lambda _𝒩\end{array}\right]=\left[\begin{array}{c}_z(z,\lambda )\\ g_{}(z)+\mathrm{\Lambda }_{}^1t_{}\\ g_𝒩(z)+\mathrm{\Lambda }_𝒩^1t_𝒩\end{array}\right],$$
(55)
where $`D_{}`$ and $`D_𝒩`$ are positive diagonal matrices defined by

$$D_{}=\mathrm{\Lambda }_{}^1S_{},D_𝒩=\mathrm{\Lambda }_𝒩^1S_𝒩.$$
(56)
When the SC condition (5) is satisfied, we have from Lemma 4 that the diagonals of $`D_{}`$ have size $`\mathrm{\Theta }(\mu )`$ while those of $`D_𝒩`$ have size $`\mathrm{\Theta }(\mu ^1)`$. By eliminating $`\mathrm{\Delta }\lambda _𝒩`$ from (55), we obtain the following intermediate stage between (27) and (3.1):

$`\left[\begin{array}{cc}H(z,\lambda )& g_{}(z)\\ g_{}(z)^T& D_{}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda _{}\end{array}\right]`$ (61)
$`=`$ $`\left[\begin{array}{c}_z(z,\lambda )g_𝒩(z)D_𝒩^1[g_𝒩(z)\mathrm{\Lambda }_𝒩^1t_𝒩]\\ g_{}(z)+\mathrm{\Lambda }_{}^1t_{}\end{array}\right],`$ (64)
where we have defined

$$H(z,\lambda )\stackrel{\mathrm{def}}{=}_{zz}(z,\lambda )+g_𝒩(z)D_𝒩^1g_𝒩(z)^T.$$
(65)
In this section, we start by proving a key result about the solutions of perturbed forms of the system (61). Subsequently, we use this result as the foundation for proving results about the three alternative formulations (26), (27), and (3.1) of the PDIP step equations. The principal reason for our focus on (61) is that the proof of the main result can be derived from fairly standard linear algebra arguments. Gould \[13, Section 6\] obtains a system similar to (61) for the Newton equations for the primal log-barrier function, and notes that the matrix approaches a nonsingular limit when certain optimality conditions, including LICQ, are satisfied. Because we replace LICQ by MFCQ, the matrix in (61) may approach a singular limit in our case.

We note that the form (61) is also relevant to the stabilized sequential quadratic programming (sSQP) method of Wright  and Hager ; that is, slight modifications to the results of this paper can be used to show that the condensed and augmented formulations of the step equations for the sSQP algorithm yield good steps even in the presence of roundoff errors and cancellation. We omit further details in this paper.

Errors in the step equations arise from cancellation and roundoff errors in evaluating both the matrix and right-hand side and from roundoff errors that arise in the factorization/solution process. We discuss these sources of error further and quantify them in the next section. In this section, we consider the following perturbed version of (61):

$$\left[\begin{array}{cc}H(z,\lambda )+\stackrel{~}{E}_{11}& g_{}(z)+\stackrel{~}{E}_{12}\\ g_{}(z)^T+\stackrel{~}{E}_{21}& D_{}+\stackrel{~}{E}_{22}\end{array}\right]\left[\begin{array}{c}w\\ y\end{array}\right]=\left[\begin{array}{c}r_1\\ g_{}(z^{})^Tr_3+r_4\end{array}\right].$$
(66)
Here, $`\stackrel{~}{E}`$ is the perturbation matrix (appropriately partitioned and not assumed to be symmetric) and $`r_1`$, $`r_3`$, and $`r_4`$ represent components of a general right-hand side. Note the partitioning of the second right-hand side component into a component $`g_{}(z^{})^Tr_3`$ in the range space of $`g_{}(z^{})^T`$ and a remainder term $`r_4`$. When LICQ is satisfied, the range space of $`g_{}(z^{})^T`$ spans the full space, so we can choose $`r_4`$ to be zero. Under MFCQ, however, we have in general that $`r_4`$ must be nonzero. The main interest of the results below is in isolating the component of the solution of (66) that is sensitive to $`r_4`$.

To make the results applicable to a wider class of linear systems, we do not impose the assumptions that were needed in the preceding section to ensure that the matrices $`D_{}`$ and $`D_𝒩`$ defined by (56) have diagonals of the appropriate size. Instead, we assume that the diagonals have the given size, and derive the application to the linear systems of interest (those that arise in primal-dual interior-point methods) as a special case.

Our results in this and later sections make extensive use of the SVD (18) of $`g_{}(z^{})`$. They also make assumptions about the size of the smallest singular value of this matrix, and about the size of the smallest eigenvalue of $`\widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V}`$, the two-sided projection of the Lagrangian Hessian onto the active constraint manifold.

###### Theorem 6.

Let $`(z,\lambda )`$ be an approximate primal-dual solution of (1) with $`\delta (z,\lambda )=O(\mu )`$, and suppose the diagonal matrices $`D_{}`$ and $`D_𝒩^1`$ defined by (56) have all their diagonal elements of size $`\mathrm{\Theta }(\mu )`$. Suppose that the perturbation submatrices in (66) satisfy

$$\stackrel{~}{E}_{11}=\delta _𝐮/\mu +O(\mu ),\stackrel{~}{E}_{21},\stackrel{~}{E}_{12},\stackrel{~}{E}_{22}=\delta _𝐮,$$
(67)
and that the following conditions hold for some $`\beta >0`$:

$`𝐮/\mu 1,𝐮1,`$ (68a)
$`\sigma _{\mathrm{min}}(\mathrm{\Sigma })\beta \mathrm{max}(\mu ^{1/3},𝐮/\mu ),`$ (68b)
$`\sigma _{\mathrm{min}}(\widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V})\beta \mathrm{max}(\mu ^{1/3},𝐮/\mu ),\text{for all }\lambda ^{}𝒮_\lambda .`$ (68c)
Then if $`\beta `$ is sufficiently large (in a sense to be specified in the proof), the step $`(w,y)`$ computed from (66) satisfies

$`(U^Ty,\widehat{V}^Tw,\widehat{U}^Tw)`$ $`=`$ $`O(r_1+r_3+r_4),`$
$`V^Ty`$ $`=`$ $`O(r_1+r_3+r_4/\mu ).`$
###### Proof.

If we define

$$y_U=U^Ty,y_V=V^Ty,w_{\widehat{U}}=\widehat{U}^Tw,w_{\widehat{V}}=\widehat{V}^Tw,$$

we have

$$y=Uy_U+Vy_V,w=\widehat{U}w_{\widehat{U}}+\widehat{V}w_{\widehat{V}}.$$

Using this notation, we can rewrite (66) as

$`\left[\begin{array}{cccc}\widehat{U}^TM_{11}\widehat{U}& \widehat{U}^TM_{11}\widehat{V}& \widehat{U}^TM_{12}U& \widehat{U}^TM_{12}V\\ \widehat{V}^TM_{11}\widehat{U}& \widehat{V}^TM_{11}\widehat{V}& \widehat{V}^TM_{12}U& \widehat{V}^TM_{12}V\\ U^TM_{21}\widehat{U}& U^TM_{21}\widehat{V}& U^TM_{22}U& U^TM_{22}V\\ V^TM_{21}\widehat{U}& V^TM_{21}\widehat{V}& V^TM_{22}U& V^TM_{22}V\end{array}\right]\left[\begin{array}{c}w_{\widehat{U}}\\ w_{\widehat{V}}\\ y_U\\ y_V\end{array}\right]`$ (77)
$`=\left[\begin{array}{c}\widehat{U}^Tr_1\\ \widehat{V}^Tr_1\\ U^Tg_{}(z^{})^Tr_3+U^Tr_4\\ V^Tg_{}(z^{})^Tr_3+V^Tr_4\end{array}\right],`$ (82)
where we have defined

$`M_{11}=H(z,\lambda )+\stackrel{~}{E}_{11},M_{12}=g_{}(z)+\stackrel{~}{E}_{12},`$ (83)
$`M_{21}=g_{}(z)^T+\stackrel{~}{E}_{21},M_{22}=D_{}+\stackrel{~}{E}_{22},`$
and $`H(,)`$ is defined in (65). From (18), we have

$$V^Tg_{}(z^{})^T=0,U^Tg_{}(z^{})^T=\mathrm{\Sigma }\widehat{U}^T.$$

The fact that $`V^T`$ annihilates $`g_{}(z^{})^T`$ is crucial, because it causes the term with $`r_3`$ to disappear from the last component of the right-hand side of (77), which becomes

$$\left[\begin{array}{c}\widehat{U}^Tr_1\\ \widehat{V}^Tr_1\\ \mathrm{\Sigma }\widehat{U}^Tr_3+U^Tr_4\\ V^Tr_4\end{array}\right].$$
(84)
From the definitions (83) and (65), the perturbation bound (67), our assumptions that $`D_𝒩^1=O(\mu )`$ and $`\delta (z,\lambda )=O(\mu )`$, compactness of $`𝒮`$, and the fact that $`_{zz}`$ is Lipschitz continuous, we have that

$$M_{11}=_{zz}(z^{},\lambda ^{})+\delta _𝐮/\mu +O(\mu ),$$
(85)
for some $`\lambda ^{}𝒮_\lambda `$. Using these same facts, we have likewise that

$$M_{21}=g_{}(z^{})^T+\delta _𝐮+O(\mu ),$$

so by substituting from (18), we have that

$`U^TM_{21}\widehat{U}=\mathrm{\Sigma }+\delta _𝐮+O(\mu ),U^TM_{21}\widehat{V}=\delta _𝐮+O(\mu ),`$ (86a)
$`V^TM_{21}\widehat{U}=\delta _𝐮+O(\mu ),V^TM_{21}\widehat{V}=\delta _𝐮+O(\mu ).`$ (86b)
Similarly, from the definition of $`M_{12}`$, we have

$`\widehat{U}^TM_{12}U=\mathrm{\Sigma }+\delta _𝐮+O(\mu ),\widehat{U}^TM_{12}V=\delta _𝐮+O(\mu ),`$ (87a)
$`\widehat{V}^TM_{12}U=\delta _𝐮+O(\mu ),\widehat{V}^TM_{12}V=\delta _𝐮+O(\mu ).`$ (87b)
For the $`M_{22}`$ block, we have from (83) and (67) that

$`U^TM_{22}U=U^TD_{}U+\delta _𝐮=O(\mu )+\delta _𝐮,`$ (88a)
$`U^TM_{22}V=O(\mu )+\delta _𝐮,V^TM_{22}U=O(\mu )+\delta _𝐮,`$ (88b)
$`V^TM_{22}V=V^TD_{}V+\delta _𝐮=\stackrel{~}{M}_{VV}+\delta _𝐮,`$ (88c)
where $`\stackrel{~}{M}_{VV}\stackrel{\mathrm{def}}{=}V^TD_{}V`$ has all its singular values of size $`\mathrm{\Theta }(\mu )`$, so that

$$\stackrel{~}{M}_{VV}^1=\mathrm{\Theta }(\mu ^1).$$
(89)
Using these estimates together with (84), we can rewrite (77) as

$$\left\{\left[\begin{array}{cc}Q& 0\\ 0& \stackrel{~}{M}_{VV}\end{array}\right]+\left[\begin{array}{cc}\widehat{E}_{11}& \widehat{E}_{12}\\ \widehat{E}_{21}& \widehat{E}_{22}\end{array}\right]\right\}\left[\begin{array}{c}w_{\widehat{U}}\\ w_{\widehat{V}}\\ y_U\\ \\ y_V\end{array}\right]=\left[\begin{array}{c}\widehat{U}^Tr_1\\ \widehat{V}^Tr_1\\ \mathrm{\Sigma }\widehat{U}^Tr_3+U^Tr_4\\ \\ V^Tr_4\end{array}\right],$$
(90)
where

$`Q`$ $`=`$ $`\left[\begin{array}{ccc}\widehat{U}^T_{zz}(z^{},\lambda ^{})\widehat{U}& \widehat{U}^T_{zz}(z^{},\lambda ^{})\widehat{V}& \mathrm{\Sigma }\\ \widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{U}& \widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V}& 0\\ \mathrm{\Sigma }& 0& 0\end{array}\right]`$ (98)
$`+\left[\begin{array}{ccc}\delta _𝐮/\mu +O(\mu )& \delta _𝐮/\mu +O(\mu )& \delta _𝐮+O(\mu )\\ \delta _𝐮/\mu +O(\mu )& \delta _𝐮/\mu +O(\mu )& 0\\ \delta _𝐮+O(\mu )& 0& 0\end{array}\right]`$
$`\stackrel{\mathrm{def}}{=}`$ $`\left[\begin{array}{ccc}N_{UU}& N_{UV}& \overline{\mathrm{\Sigma }}_1\\ N_{VU}& N_{VV}& 0\\ \overline{\mathrm{\Sigma }}_2& 0& 0\end{array}\right],`$ (102)
while

$$\widehat{E}_{11}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& \delta _𝐮+O(\mu )\\ 0& \delta _𝐮+O(\mu )& \delta _𝐮+O(\mu )\end{array}\right],$$
(103)
and

$$\widehat{E}_{12},\widehat{E}_{21}=\delta _𝐮+O(\mu )=O(\mu ),\widehat{E}_{22}=\delta _𝐮.$$
(104)
For purposes of specifying the required conditions on $`\beta `$ in (68b) and (68c), we define $`\kappa `$ to be a constant such that expressions of size $`\delta _𝐮`$ and $`O(\mu )`$ that arise in the perturbation terms in the coefficient matrix in (90) can be bounded by $`\kappa 𝐮`$ and $`\kappa \mu `$, respectively. For example, we suppose that the perturbations in $`\overline{\mathrm{\Sigma }}_1`$, $`\overline{\mathrm{\Sigma }}_2`$, and $`N_{VV}`$ can bounded as follows:

$`\overline{\mathrm{\Sigma }}_1\mathrm{\Sigma }\kappa (\mu +𝐮),\overline{\mathrm{\Sigma }}_2\mathrm{\Sigma }\kappa (\mu +𝐮),`$ (105a)
$`N_{VV}\widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V}\kappa (𝐮/\mu +\mu ),`$ (105b)
and that

$$\widehat{E}_{11}\kappa (𝐮+\mu ),\widehat{E}_{12}\kappa (𝐮+\mu ),\widehat{E}_{21}\kappa (𝐮+\mu ),\widehat{E}_{22}\kappa 𝐮.$$
(106)
From (105a) and (68b), we have that

$$\overline{\mathrm{\Sigma }}_1\mathrm{\Sigma }\kappa \mathrm{max}(\mu ^{1/3},𝐮/\mu )(\kappa /\beta )\sigma _{\mathrm{min}}(\mathrm{\Sigma })(\kappa /\beta )\mathrm{\Sigma }.$$

It is therefore easy to show that if $`\beta `$ can be chosen large enough that $`\beta >2\kappa `$ (while still satisfying (68b) and (68c)), then

$$\overline{\mathrm{\Sigma }}_12\mathrm{\Sigma },\overline{\mathrm{\Sigma }}_1^12\mathrm{\Sigma }^1.$$
(107)
Similarly, we can show that

$$\overline{\mathrm{\Sigma }}_22\mathrm{\Sigma },\overline{\mathrm{\Sigma }}_2^12\mathrm{\Sigma }^1,$$
(108)
$$N_{VV}2\widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V},N_{VV}^12(\widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V})^1.$$
(109)
Note, too, that because of Lipschitz continuity of $`_{zz}`$ and compactness of $`𝒮`$, and the bounds (68a), the norms of $`N_{UU}`$, $`N_{UV}`$, $`N_{VU}`$, $`N_{VV}`$, and $`\mathrm{\Sigma }`$ are all $`O(1)`$. Hence the matrix $`Q`$ is itself invertible, and we have

$$Q^1=\left[\begin{array}{ccc}0& 0& \overline{\mathrm{\Sigma }}_2^1\\ 0& N_{VV}^1& N_{VV}^1N_{VU}\overline{\mathrm{\Sigma }}_2^1\\ \overline{\mathrm{\Sigma }}_1^1& \overline{\mathrm{\Sigma }}_1^1N_{UV}N_{VV}^1& \overline{\mathrm{\Sigma }}_1^1(N_{UU}N_{UV}N_{VV}^1N_{VU})\overline{\mathrm{\Sigma }}_2^1\end{array}\right].$$
(110)
Noting that

$$(Q+\widehat{E}_{11})^1=(I+Q^1\widehat{E}_{11})^1Q^1,$$
(111)
we examine the size of $`Q^1\widehat{E}_{11}`$. Note first from (68b) and (68c) together with (107), (108), and (109) that

$`\overline{\mathrm{\Sigma }}_1^1{\displaystyle \frac{2}{\beta }}(𝐮/\mu )^1,\overline{\mathrm{\Sigma }}_2^1{\displaystyle \frac{2}{\beta }}(𝐮/\mu )^1,N_{VV}^1{\displaystyle \frac{2}{\beta }}(𝐮/\mu )^1,`$ (112a)
$`\overline{\mathrm{\Sigma }}_1^1{\displaystyle \frac{2}{\beta }}\mu ^{1/3},\overline{\mathrm{\Sigma }}_2^1{\displaystyle \frac{2}{\beta }}\mu ^{1/3},N_{VV}^1{\displaystyle \frac{2}{\beta }}\mu ^{1/3}.`$ (112b)
By forming the product of (110) with (103) and using the bounds in (112), we can show that the norm of $`Q^1\widehat{E}_{11}`$ can be made less than $`1/2`$ provided that $`\beta `$ in (68b), (68c) is sufficiently large. The $`(3,3)`$ block of $`Q^1\widehat{E}_{11}`$, for instance, has the form

$$\overline{\mathrm{\Sigma }}_1^1N_{UV}N_{VV}^1(\delta _𝐮+O(\mu ))+\overline{\mathrm{\Sigma }}_1^1(N_{UU}N_{UV}N_{VV}^1N_{VU})\overline{\mathrm{\Sigma }}_2^1(\delta _𝐮+O(\mu )).$$

Because of (106), its norm can be bounded by a quantity of the form

$$C\kappa \left(\overline{\mathrm{\Sigma }}_1^1N_{VV}^1+\overline{\mathrm{\Sigma }}_1^1\overline{\mathrm{\Sigma }}_2^1N_{VV}^1+\overline{\mathrm{\Sigma }}_1^1\overline{\mathrm{\Sigma }}_2^1\right)\left((𝐮/\mu )\mu +\mu \right),$$

(for some $`C`$ that depends on $`_{zz}(z^{},\lambda ^{})`$), which in turn because of (112) is bounded by the following quantity:

$$8C\kappa \left(\frac{1}{\beta ^2}\mu ^{2/3}+\frac{1}{\beta ^3}\mu ^{1/3}\right)+8C\kappa \left(\frac{1}{\beta ^2}\mu ^{1/3}+\frac{1}{\beta ^3}\right).$$

Provided that $`\beta `$ is large enough that this and the other blocks of $`Q^1\widehat{E}_{11}`$ can be bounded appropriately, we have that $`Q^1\widehat{E}_{11}1/2`$, and therefore from (111) we have

$$(Q+\widehat{E}_{11})^1=2Q^1.$$

Our conclusion is that for $`\beta `$ satisfying the conditions outlined in this paragraph, the inverse of the $`(1,1)`$ block of the matrix in (90) can be bounded in terms of $`Q^1`$, which because of (107), (108), (109), and (110) can in turn be bounded by a finite quantity that depends only on the problem data and not on $`\mu `$.

Returning to (90), and using (104), we have that

$`\left[\begin{array}{c}w_{\widehat{U}}\\ w_{\widehat{V}}\\ y_U\end{array}\right]`$ $`=`$ $`(Q+\widehat{E}_{11})^1\widehat{E}_{12}y_V+(Q+\widehat{E}_{11})^1\left[\begin{array}{c}\widehat{U}^Tr_1\\ \widehat{V}^Tr_1\\ \mathrm{\Sigma }\widehat{U}^Tr_3+U^Tr_4\end{array}\right]`$
$`=`$ $`O(\widehat{E}_{12}y_V)+O(r_1+r_3+r_4)`$
$`=`$ $`O(\mu )y_V+O(r_1+r_3+r_4).`$ (120)
Meanwhile, for the second block row of (90), we obtain

$$y_V=(\stackrel{~}{M}_{VV}+\widehat{E}_{22})^1\widehat{E}_{21}\left[\begin{array}{c}w_{\widehat{U}}\\ w_{\widehat{V}}\\ y_U\end{array}\right]+(\stackrel{~}{M}_{VV}+\widehat{E}_{22})^1V^Tr_4.$$
(121)
Since from (89), (104), and (68a), we have

$$(\stackrel{~}{M}_{VV}+\widehat{E}_{22})^1=(I+\stackrel{~}{M}_{VV}^1\widehat{E}_{22})^1\stackrel{~}{M}_{VV}^1=(I+\delta _𝐮/\mu )\stackrel{~}{M}_{VV}^1=O(\mu ^1),$$

it follows from (121) and (104) that

$$y_V=O(\mu ^1)O(\mu )\left[\begin{array}{c}w_{\widehat{U}}\\ w_{\widehat{V}}\\ y_U\end{array}\right]+O(\mu ^1)O(r_4).$$

By substituting from (120), we obtain

$$y_V=O(\mu )y_V+O(r_1+r_3+r_4)+O(r_4/\mu ),$$

and therefore

$$y_V=O(r_1+r_3+r_4/\mu ),$$

as claimed. The estimate for $`(w_{\widehat{U}},w_{\widehat{V}},y_U)`$ is obtained by substituting into (120). ∎

The conditions (68) need a little explanation. For the typical value $`𝐮=10^{16}`$, the minimum value of the quantity $`\mathrm{max}(\mu ^{1/3},𝐮/\mu )`$ is $`10^4`$, achieved at $`\mu ^{12}`$. Moreover, we have $`\mathrm{max}(\mu ^{1/3},𝐮/\mu )10^2`$ only for $`\mu `$ in the range $`[10^{14},10^6]`$. It would seem, then, that the problem would need to be quite well conditioned for (68b) and (68c) to hold and that $`\mu `$ may have to become quite small before the results apply. We note, however, that the purpose of the bounds (68b) and (68c) is to ensure that the inverse of $`Q+\widehat{E}_{11}`$ can be bounded independently of $`\mu `$, and that for this purpose they are quite conservative. That is, we would expect to find that $`(Q+\widehat{E}_{11})^1`$ is not too much larger than the norm of the inverse of the corresponding exact matrix (the first term on the right-hand side of (98)) for $`\mu `$ not much less than the smallest eigenvalues of $`\mathrm{\Sigma }`$ and $`\widehat{V}^T_{zz}(z^{},\lambda ^{})\widehat{V}`$.

The requirement that $`𝐮/\mu `$ and $`\mu `$ both be small in (68) may not seem to sit well with expressions such as $`O(\mu )`$ and $`O(\mu ^2)`$, which crop up repeatedly in the analysis and which assert that certain bounds hold “for all sufficiently small $`\mu `$.” As noted in the preceding paragraph, this requirement implies that the analysis holds for $`\mu `$ in a certain range, or “window,” of values. Similar windows are used in the analysis of S. Wright , and M. H. Wright , and numerical experience indicates that such a window does indeed exist in most practical cases. We expect the same to be true of the problem and algorithms discussed in this paper.

At this point, we assemble the assumptions that are made in the remainder of the paper into a single catch-all assumption.

###### Assumption 4.1.

* $`z^{}`$ is a solution of (1), so that the condition (3) holds. The MFCQ condition (10), the second-order conditions (11), (12), and the SC condition (5) are satisfied at this solution. The current iterate $`(z,\lambda ,s)`$ of the PDIP algorithm satisfies the conditions (29), and the right-hand side modification $`t`$ satisfies (25).
* The quantities $`\mu `$, $`𝐮`$ (15), $`_{zz}(z^{},\lambda ^{})`$, $`\mathrm{\Sigma }`$, and $`\widehat{V}`$ (18) satisfy the conditions (68).

From our observations following (56), we have under this assumption that

$$D_{}=O(\mu ),D_{}^1=O(\mu ^1),D_𝒩=O(\mu ^1),D_𝒩^1=O(\mu ).$$
(122)
Our next result considers a perturbed form of the system (55), with a general right-hand side. By eliminating one component to obtain the form (61), we can apply Theorem 6 to obtain estimates of the dependence of the solution on the right-hand side components.

###### Theorem 7.

Suppose that Assumption 4.1 holds. Consider the linear system

$`\left[\begin{array}{ccc}_{zz}(z,\lambda )+E_{11}& g_{}(z)+E_{12}& g_𝒩(z)+E_{13}\\ g_{}(z)^T+E_{21}& D_{}+E_{22}& E_{23}\\ g_𝒩(z)^T+E_{31}& E_{32}& D_𝒩+E_{33}\end{array}\right]\left[\begin{array}{c}w\\ y\\ q\end{array}\right]`$ (129)
$`=\left[\begin{array}{c}r_5\\ g_{}(z^{})^Tr_6+r_7\\ r_8\end{array}\right],`$ (133)
where

$`E_{11}=\delta _𝐮/\mu ,E_{33}=\delta _𝐮/\mu ^2,`$ (134a)
$`E_{12},E_{21},E_{22}=\delta _𝐮,E_{13},E_{31},E_{23},E_{32}=\delta _𝐮/\mu .`$ (134b)
Then the step $`(w,y,q)`$ satisfies the following estimates:

$`(U^Ty,w)`$ $`=`$ $`O(r_5+r_6+r_7+\mu r_8),`$
$`V^Ty`$ $`=`$ $`O(r_5+r_6+r_7/\mu +(\delta _𝐮/\mu +O(\mu ))r_8),`$
$`q`$ $`=`$ $`O(\mu )\left[r_5+r_6+r_8\right]+(\delta _𝐮/\mu +O(\mu ))r_7.`$
###### Proof.

From (122) and the assumed bound (134a) on the size of $`E_{33}`$, we have that

$`(D_𝒩+E_{33})^1`$ (135)
$`=`$ $`(ID_𝒩^1E_{33})^1D_𝒩^1=(I+O(\mu )\delta _𝐮/\mu ^2)O(\mu )=O(\mu ).`$
By eliminating $`q`$ from (133), we obtain the reduced system

$$\left[\begin{array}{cc}H(z,\lambda )+\stackrel{~}{E}_{11}& g_{}(z)+\stackrel{~}{E}_{12}\\ g_{}(z)^T+\stackrel{~}{E}_{21}& D_{}+\stackrel{~}{E}_{22}\end{array}\right]\left[\begin{array}{c}w\\ y\end{array}\right]=\left[\begin{array}{c}r_5+O(\mu )r_8\\ g_{}(z^{})^Tr_6+r_7+\delta _𝐮r_8\end{array}\right],$$

where from (68) and (65), we obtain

$`\stackrel{~}{E}_{11}`$ $`=`$ $`E_{11}(g_𝒩(z)+E_{13})(D_𝒩+E_{33})^1(g_𝒩(z)^T+E_{31})g_𝒩(z)D_𝒩^1g_𝒩(z)^T`$
$`=`$ $`\delta _𝐮/\mu +O(\mu ),`$
$`\stackrel{~}{E}_{12}`$ $`=`$ $`E_{12}(g_𝒩(z)+E_{13})(D_𝒩+E_{33})^1E_{32}=\delta _𝐮+O(1)O(\mu )\delta _𝐮/\mu =\delta _𝐮,`$
$`\stackrel{~}{E}_{21}`$ $`=`$ $`E_{21}E_{23}(D_𝒩+E_{33})^1(g_𝒩(z)^T+E_{31})=\delta _𝐮+(\delta _𝐮/\mu )O(\mu )O(1)=\delta _𝐮,`$
$`\stackrel{~}{E}_{22}`$ $`=`$ $`E_{22}E_{23}(D_𝒩+E_{33})^1E_{32}=\delta _𝐮+(\delta _𝐮/\mu )^2O(\mu )=\delta _𝐮.`$
These perturbation matrices satisfy the assumptions of Theorem 6, which can be applied to give

$`(U^Ty,\widehat{V}^Tw,\widehat{U}^Tw)`$ $`=`$ $`O(r_5+r_6+r_7+\mu r_8),`$ (136a)
$`V^Ty`$ $`=`$ $`O(r_5+r_6+r_7/\mu )+(\delta _𝐮/\mu +O(\mu ))r_8.`$ (136b)
From the last block row of (133), and using (68), (135), (136), we obtain

$`q`$ $`=`$ $`(D_𝒩+E_{33})^1\left[r_8(g_𝒩(z)^T+E_{31})wE_{32}y\right]`$
$`=`$ $`O(\mu )\left[r_8+w+(\delta _𝐮/\mu )y\right]`$
$`=`$ $`O(\mu )\left[r_5+r_6+r_7+r_8\right]+`$
$`\delta _𝐮\left[r_5+r_6+r_7/\mu +(\delta _𝐮/\mu +O(\mu ))r_8\right]`$
$`=`$ $`O(\mu )\left[r_5+r_6+r_8\right]+(\delta _𝐮/\mu +O(\mu ))r_7.`$
An estimate for the solution of the exact system (26) follows almost immediately from this result. This is the key technical result used by Ralph and Wright  to prove superlinear convergence of PDIP algorithms for convex programming problems. The result below, however, does not require a convexity assumption.

###### Corollary 8.

Suppose that Assumption 4.1(a) holds. Then the (exact) solution $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda ,\mathrm{\Delta }s)`$ of the system (26) satisfies

$$(\mathrm{\Delta }z,\mathrm{\Delta }\lambda ,\mathrm{\Delta }s)=O(\mu ).$$
(137)
###### Proof.

Note first that Assumption 4.1(b) holds trivially in this case for $`\mu `$ sufficiently small, because our assumption of exact computations is equivalent to setting $`𝐮=0`$. We prove the result by identifying the system (55) with (133) and then applying Theorem 7.

For the right-hand side, we note first that, because of smoothness of $`g`$, Taylor’s theorem, the definition (6) of $``$, and Theorem 5,

$`g_{}(z)`$ $`=`$ $`g_{}(z^{})+g_{}(z^{})^T(zz^{})+O(zz^{}^2)`$ (138)
$`=`$ $`g_{}(z^{})^T(zz^{})+O(\mu ^2).`$
We now identify the right-hand side of (55) with (133) by setting

$`r_5`$ $`=`$ $`_z(z,\lambda ),`$
$`r_6`$ $`=`$ $`(zz^{}),`$
$`r_7`$ $`=`$ $`g_{}(z^{})^T(zz^{})g_{}(z)+\mathrm{\Lambda }_{}^1t_{},`$
$`r_8`$ $`=`$ $`g_𝒩(z)+\mathrm{\Lambda }_𝒩^1t_𝒩.`$
The sizes of these vectors can be estimated by using (29), Lemma 4, (138), Theorem 5, and the assumption (25) on the size of $`t`$ to obtain

$$r_5=O(\mu ),r_6=O(\mu ),r_7=O(\mu ^2),r_8=O(1).$$
(139)
(By choosing $`r_6`$ and $`r_7`$ in this way, we ensure that the terms involving $`r_7/\mu `$ in the estimates of the solution components in Theorem 7 are not grossly larger than the other terms in these expressions.) We complete the identification of (55) with (133) by setting all the perturbation matrices $`E_{11},E_{12},\mathrm{},E_{33}`$ to zero and by identifying the solution vector components $`\mathrm{\Delta }z`$, $`\mathrm{\Delta }\lambda _{}`$, and $`\mathrm{\Delta }\lambda _𝒩`$ with $`w`$, $`y`$, and $`q`$, respectively. By directly applying Theorem 7, substituting the estimates (139), and setting $`\delta _𝐮=0`$ (since we are assuming exact computations), we have that

$$(U^T\mathrm{\Delta }\lambda _{},\mathrm{\Delta }z)=O(\mu ),V^T\mathrm{\Delta }\lambda _{}=O(\mu ),\mathrm{\Delta }\lambda _𝒩=O(\mu ).$$

To show that the remaining solution component $`\mathrm{\Delta }s`$ of (26) is also of size $`O(\mu )`$, we write the second block row in (26) as

$$\mathrm{\Delta }s=(g(z)+s)g(z)^T\mathrm{\Delta }z,$$

from which the desired estimate follows immediately by substituting from (29b) and $`\mathrm{\Delta }z=O(\mu )`$. ∎

The next result uses Theorem 7 to compare perturbed and exact solutions of the system of the system (55).

###### Corollary 9.

Suppose that Assumption 4.1 holds. Let $`(w,y,q)`$ be obtained from the following perturbed version of (27):

$`\left[\begin{array}{ccc}_{zz}(z,\lambda )+E_{11}& g_{}(z)+E_{12}& g_𝒩(z)+E_{13}\\ g_{}(z)^T+E_{21}& D_{}+E_{22}& E_{23}\\ g_𝒩(z)^T+E_{31}& E_{32}& D_𝒩+E_{33}\end{array}\right]\left[\begin{array}{c}w\\ y\\ q\end{array}\right]`$ (146)
$`=\left[\begin{array}{c}_z(z,\lambda )+f_1\\ g_{}(z)+\mathrm{\Lambda }_{}^1t_{}+f_2\\ g_𝒩(z)+\mathrm{\Lambda }_𝒩^1t_𝒩+f_3\end{array}\right],`$ (150)
where $`E_{ij}`$, $`i,j=1,2,3`$, satisfy the conditions (134) and $`f_1`$, $`f_2`$, and $`f_3`$ are all of size $`\delta _𝐮`$. Then if $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda ,\mathrm{\Delta }s)`$ is the (exact) solution of the system (26), we have

$`(\mathrm{\Delta }zw,U^T(\mathrm{\Delta }\lambda _{}y))`$ $`=`$ $`\delta _𝐮,`$
$`V^T(\mathrm{\Delta }\lambda _{}y)`$ $`=`$ $`\delta _𝐮/\mu ,`$
$`\mathrm{\Delta }\lambda _𝒩q`$ $`=`$ $`\delta _𝐮.`$
###### Proof.

By combining (150) with (55), we obtain

$`\left[\begin{array}{ccc}_{zz}(z,\lambda )+E_{11}& g_{}(z)+E_{12}& g_𝒩(z)+E_{13}\\ g_{}(z)^T+E_{21}& D_{}+E_{22}& E_{23}\\ g_𝒩(z)^T+E_{31}& E_{32}& D_𝒩+E_{33}\end{array}\right]\left[\begin{array}{c}w\mathrm{\Delta }z\\ y\mathrm{\Delta }\lambda _{}\\ q\mathrm{\Delta }\lambda _𝒩\end{array}\right]`$ (157)
$`=\left[\begin{array}{c}f_1\\ f_2\\ f_3\end{array}\right]\left[\begin{array}{ccc}E_{11}& E_{12}& E_{13}\\ E_{21}& E_{22}& E_{23}\\ E_{31}& E_{32}& E_{33}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda _{}\\ \mathrm{\Delta }\lambda _𝒩\end{array}\right].`$ (167)
From the bounds on the perturbations $`E`$ in (134) and the result of Corollary 8, we have for the right-hand side of this expression that

$`\left[\begin{array}{c}r_5\\ r_7\\ r_8\end{array}\right]`$ $`\stackrel{\mathrm{def}}{=}`$ $`\left[\begin{array}{c}f_1\\ f_2\\ f_3\end{array}\right]\left[\begin{array}{ccc}E_{11}& E_{12}& E_{13}\\ E_{21}& E_{22}& E_{23}\\ E_{31}& E_{32}& E_{33}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda _{}\\ \mathrm{\Delta }\lambda _𝒩\end{array}\right]`$ (180)
$`=`$ $`\left[\begin{array}{c}\delta _𝐮+(\delta _𝐮/\mu )\mu +\delta _𝐮\mu +(\delta _𝐮/\mu )\mu \\ \delta _𝐮+\delta _𝐮\mu +\delta _𝐮\mu +(\delta _𝐮/\mu )\mu \\ \delta _𝐮+(\delta _𝐮/\mu )\mu +(\delta _𝐮/\mu )\mu +(\delta _𝐮/\mu ^2)\mu \end{array}\right]=\left[\begin{array}{c}\delta _𝐮\\ \delta _𝐮\\ \delta _𝐮/\mu \end{array}\right].`$ (187)
Using these estimates, we can simply apply Theorem 7 to (167) (with $`r_6=0`$) to obtain the result. ∎

For later reference, we show how the estimates of Corollary 9 can be modified when the perturbations have a special form. Suppose that

$$E_{23}=0,E_{33}=\delta _𝐮/\mu ,f_2=Uf_2^U+O(\mu ^2),\text{where}f_2^U=\delta _𝐮,$$
(188)
where $`U`$ is the matrix from (18). Instead of setting $`r_6=0`$ as in the proof above, we set

$$r_6=\widehat{U}\mathrm{\Sigma }f_2^U=\delta _𝐮$$

(using (18) to obtain an $`r_6`$ for which $`g_{}(z^{})^Tr_6=Uf_2^U`$). By modifying (187) to account for the remaining perturbations, we can identify (167) with (133) by setting

$`\left[\begin{array}{c}r_5\\ r_7\\ r_8\end{array}\right]`$ $`\stackrel{\mathrm{def}}{=}`$ $`\left[\begin{array}{c}f_1\\ f_2Uf_2^U\\ f_3\end{array}\right]\left[\begin{array}{ccc}E_{11}& E_{12}& E_{13}\\ E_{21}& E_{22}& E_{23}\\ E_{31}& E_{32}& E_{33}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }z\\ \mathrm{\Delta }\lambda _{}\\ \mathrm{\Delta }\lambda _𝒩\end{array}\right]`$ (201)
$`=`$ $`\left[\begin{array}{c}\delta _𝐮+(\delta _𝐮/\mu )\mu +\delta _𝐮\mu +(\delta _𝐮/\mu )\mu \\ O(\mu ^2)+\delta _𝐮\mu +\delta _𝐮\mu \\ \delta _𝐮+(\delta _𝐮/\mu )\mu +(\delta _𝐮/\mu )\mu +(\delta _𝐮/\mu ^2)\mu \end{array}\right]=\left[\begin{array}{c}\delta _𝐮\\ O(\mu ^2)\\ \delta _𝐮/\mu \end{array}\right].`$ (208)
Using these modified right-hand side estimates, we can apply Theorem 7 to obtain the following improved bound on one of the components:

$$V^T(\mathrm{\Delta }\lambda _{}y)=O(\mu ).$$
(209)
The bounds on the other components remain unchanged.

We emphasize that the conditions (29), and in particular (29c), are critical to the results of this and all the remaining sections of the paper. These conditions enable Lemma 4, which in turn enable us to assert that the diagonals of $`D_{}`$ all have size $`\mathrm{\Theta }(\mu )`$ while those of $`D_𝒩`$ all have size $`\mathrm{\Theta }(\mu ^1)`$ (see (122)). This neat classification of the diagonals of $`D`$ into two categories drives all the subsequent analysis. The motivation for conditions like (29) in the literature for path-following methods (with exact steps) is not unrelated: It allows us to obtain bounds on the steps and to show that we can move a significant distance along this direction while ensuring that (29) continues to be satisfied at the new iterate. (See, for example, \[26, Chapters 5 and 6\] and its bibliography for the case of linear programming and  for the case of nonlinear convex programming.) In the analysis above, we obtain bounds on the errors (rather than the steps themselves) when perturbation terms of a certain structure appear in the matrix and right-hand side.

Many practical implementations of path-following methods for linear programming do not explicitly check that the condition (29c) is satisfied by the calculated iterates (see, for example,  and ). However, the heuristics for “stepping back” from the boundary of the nonnegative orthant by a small but significant quantity are motivated by this condition, and it is observed to hold in practice on all but the most recalcitrant problems.

## 5 The Condensed System

Here we consider an algorithm in which the condensed linear system (3.1) is formed and solved to obtain $`\mathrm{\Delta }z`$, and the remaining step components $`\mathrm{\Delta }\lambda `$ and $`\mathrm{\Delta }s`$ are recovered from (26). We obtain expressions for the errors in the calculated step $`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ and discuss the effects of these errors on certain measures of step quality. We also derive conditions under which the Cholesky factorization applied to (3.1) is guaranteed to run to completion.

Formally, the complete procedure can be described as follows:

| procedure condensed |
| --- |
| given the current iterate $`(z,\lambda ,s)`$ |
|  | form the coefficient matrix and right-hand side for (3.1); |
|  | solve (3.1) using a backward stable algorithm to obtain $`\mathrm{\Delta }z`$; |
|  | set $`\mathrm{\Delta }\lambda =D^1[g(z)\mathrm{\Lambda }^1t+g(z)^T\mathrm{\Delta }z]`$; |
|  | set $`\mathrm{\Delta }s=(g(z)+s)g(z)^T\mathrm{\Delta }z`$. |

We have used the definition (56) of the matrix $`D`$. For convenience, we restate the system (3.1) here as follows:

$$\left[_{zz}(z,\lambda )+g(z)D^1g(z)^T\right]\mathrm{\Delta }z=_z(z,\lambda )g(z)D^1[g(z)\mathrm{\Lambda }^1t].$$
(210)
Note that this procedure requires evaluation of $`D^1=S^1\mathrm{\Lambda }`$, rather than $`D`$ itself.

### 5.1 Quantifying the Errors

When implemented in finite-precision arithmetic, solution of (210) gives rise to errors of three types:

* cancellation in evaluation of the matrix and right-hand side;
* roundoff errors in evaluation of the matrix and right-hand side;
* roundoff errors that accumulate during the process of factoring the matrix and using triangular substitutions to obtain the solution.

Cancellation may be an issue in the evaluation of the nonlinear functions $`_{zz}(z,\lambda )`$, $`_z(z,\lambda )`$, $`g(z)`$, and $`g(z)`$, because intermediate terms computed during the additive evaluation of these quantities may exceed the size of the final result (see Golub and Van Loan \[12, p. 61\]). The intermediate terms generally contain rounding error (which occurs whenever real numbers are represented in finite precision). Cancellation becomes a significant phenomenon whenever we take a difference of two nearly equal quantities, since the error in the computed result due to roundoff in the two arguments may be large relative to the size of the result. If, as we can reasonably assume, intermediate quantities in the calculations of our right-hand sides remain bounded, the absolute size of the errors in the result is $`\delta _𝐮`$. In the case of $`_z(z,\lambda )`$ and $`g_{}(z)`$, the final result in exact arithmetic has size $`O(\mu )`$, so that the error $`\delta _𝐮`$ takes on a large relative significance for small values of $`\mu `$. This fact causes the error bound in some components of the solution to be larger than in others, as we see in (225c) below. In summary, the computed versions of the quantities discussed above differ from their exact values in the following way:

$`\text{computed }_{zz}(z,\lambda )_{zz}(z,\lambda )+\overline{F},`$ (211a)
$`\text{computed }_z(z,\lambda )_z(z,\lambda )+\overline{f},`$ (211b)
$`\text{computed }g(z)g(z)+F=\left[\begin{array}{c}g_{}(z)\\ g_𝒩(z)\end{array}\right]+\left[\begin{array}{c}F_{}\\ F_𝒩\end{array}\right],`$ (211g)
$`\text{computed }g(z)g(z)+f=\left[\begin{array}{c}g_{}(z)\\ g_𝒩(z)\end{array}\right]+\left[\begin{array}{c}f_{}\\ f_𝒩\end{array}\right],`$ (211l)
where $`\overline{F}`$, $`\overline{f}`$, $`F`$, and $`f`$ are all of size $`\delta _𝐮`$. Earlier discussion of cancellation in similar contexts can be found in the papers of S. Wright  and M. H. Wright .

The second source of error is evaluation of the matrix $`D^1`$. From the model (15) of floating-point arithmetic and the estimates (34) of Lemma 4, we have that

$`\text{computed }D_{}^1`$ $``$ $`(D_{}+G_{})^1,G_{}=\mu \delta _𝐮,`$ (212a)
$`\text{computed }D_𝒩^1`$ $``$ $`(D_𝒩+G_𝒩)^1,G_𝒩=\delta _𝐮/\mu ,`$ (212b)
where $`G_{}`$ and $`G_𝒩`$ are both diagonal matrices that can be composed into a single diagonal matrix $`G`$.

Third, we account for the error in forming the matrix and right-hand side of (210) from the computed quantities described in the last two paragraphs. Since we are now dealing with floating-point numbers, the model (15) applies; that is, any additional errors that arise during the combination of these floating-point quantities have size $`𝐮`$ relative to the size of the result of the calculation. Since the norm of the coefficient matrix is of size $`O(\mu ^1)`$ while the right-hand side has size $`O(1)`$ (see (29)), we represent these errors by a matrix $`\widehat{F}`$ of size $`\delta _𝐮/\mu `$ and a vector $`\widehat{f}`$ of size $`\delta _𝐮`$.

Finally, we account for the error that arises in the application of a backward-stable method to solve (210). Specifically, we assume that the method yields a computed solution that is the exact solution of a nearby problem whose data contains relative perturbations of size $`𝐮`$. The absolute sizes of these terms would therefore be $`\delta _𝐮/\mu `$ in the case of the matrix and $`\delta _𝐮`$ in the case of the right-hand side. Since these errors are the same size as those discussed in the preceding paragraph, we incorporate them into the matrix $`\widehat{F}`$ and the vector $`\widehat{f}`$.

Summarizing, we find that the computed solution $`\widehat{\mathrm{\Delta }z}`$ of (210) satisfies the following system:

$`\left[_{zz}(z,\lambda )+\overline{F}+(g(z)+F)(D+G)^1(g(z)+F)^T+\widehat{F}\right]\widehat{\mathrm{\Delta }z}`$
$`=`$ $`_z(z,\lambda )\overline{f}(g(z)+F)(D+G)^1[g(z)+f\mathrm{\Lambda }^1t]+\widehat{f},`$
where the perturbation terms $`\overline{F}`$, $`F`$, $`\widehat{F}`$, $`G`$, $`\overline{f}`$, $`\widehat{f}`$, and $`f`$ are described in the paragraphs above. By “unfolding” this system and using the partitioning of $`F`$, $`G`$, and $`f`$ defined in (211) and (212), we find that $`\widehat{\mathrm{\Delta }z}`$ also satisfies the following system, for some vectors $`y`$ and $`q`$:

$`\left[\begin{array}{ccc}_{zz}(z,\lambda )+\overline{F}+\widehat{F}& g_{}(z)+F_{}& g_𝒩(z)+F_𝒩\\ g_{}(z)^T+F_{}^T& D_{}G_{}& 0\\ g_𝒩(z)^T+F_𝒩^T& 0& D_𝒩G_𝒩\end{array}\right]\left[\begin{array}{c}\widehat{\mathrm{\Delta }z}\\ y\\ q\end{array}\right]`$ (220)
$`=\left[\begin{array}{c}_z(z,\lambda )\overline{f}+\widehat{f}\\ g_{}(z)+\mathrm{\Lambda }_{}^1t_{}f_{}\\ g_𝒩(z)+\mathrm{\Lambda }_𝒩^1t_𝒩f_𝒩\end{array}\right].`$ (224)
This system has precisely the form of (150) (in particular, the perturbation matrices satisfy the appropriate bounds). Hence, by a direct application of Corollary 9, we conclude that

$`\mathrm{\Delta }z\widehat{\mathrm{\Delta }z}`$ $`=`$ $`\delta _𝐮,`$ (225a)
$`U^T(\mathrm{\Delta }\lambda _{}y)`$ $`=`$ $`\delta _𝐮,`$ (225b)
$`V^T(\mathrm{\Delta }\lambda _{}y)`$ $`=`$ $`\delta _𝐮/\mu .`$ (225c)
We return now to the recovery of the remaining solution components $`\widehat{\mathrm{\Delta }\lambda }`$ and $`\widehat{\mathrm{\Delta }s}`$ in the procedure condensed. We have from Assumption 4.1 together with (29b), Lemma 4, (137), (225a), (212a), (68), and (122) that

$`g_{}(z)=r_g(z,s)_{}s_{}=O(\mu ),\mathrm{\Lambda }_{}^1=\mathrm{\Theta }(1),\widehat{\mathrm{\Delta }z}=\mathrm{\Delta }z+\delta _𝐮=O(\mu ),`$ (226a)
$`(D_{}+G_{})^1=(I+D_{}^1G_{})D_{}^1=(I+\delta _𝐮)^1O(\mu ^1)=O(\mu ^1).`$ (226b)
Since $`t=O(\mu ^2)`$, we have from our model (15) that the floating-point version of the calculation of $`\widehat{\mathrm{\Delta }\lambda }_{}`$ in the procedure condensed satisfies the following:

$$\widehat{\mathrm{\Delta }\lambda }_{}=(D_{}+G_{})^1\left[g_{}(z)+f_{}\mathrm{\Lambda }_{}^1t_{}+(g_{}(z)+F_{})^T\widehat{\mathrm{\Delta }z}+\mu \delta _𝐮\right]+\delta _𝐮.$$

(The final term $`\delta _𝐮`$ arises from (15) because our best estimate of the quantity in the brackets at this point of the analysis is $`O(\mu )`$, so from (226b) the result has size $`O(1)`$.) Meanwhile, we have from the second block row of (220) that

$$y=(D_{}+G_{})^1\left[g_{}(z)+f_{}\mathrm{\Lambda }_{}^1t_{}+(g_{}(z)+F_{})^T\widehat{\mathrm{\Delta }z}\right].$$

By a direct comparison of these two expressions, and using $`(D_{}+G_{})^1=O(\mu )`$, we find that

$$\widehat{\mathrm{\Delta }\lambda }_{}y=\delta _𝐮.$$
(227)
By combining (227) with (225b) and (225c), we find that

$$U^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{})=\delta _𝐮,V^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{})=\delta _𝐮/\mu .$$
(228)
For the “nonbasic” part $`\widehat{\mathrm{\Delta }\lambda }_𝒩`$, we have from (29b), Lemma 4, (137), (225a), (212b), (68), and (122) that

$`g_𝒩(z)=O(1),\mathrm{\Lambda }_𝒩^1=O(\mu ^1),\widehat{\mathrm{\Delta }z}=O(\mu ),`$ (229a)
$`(D_𝒩+G_𝒩)^1=(I+D_𝒩^1G_𝒩)^1D_𝒩^1=D_𝒩^1+\mu \delta _𝐮=O(\mu ).`$ (229b)
By using $`t_𝒩=O(\mu ^2)`$ and applying the model (15) to the appropriate step in the procedure condensed, we obtain

$$\widehat{\mathrm{\Delta }\lambda }_𝒩=(D_𝒩+G_𝒩)^1\left[g_𝒩(z)+f_𝒩\mathrm{\Lambda }_𝒩^1t_𝒩+(g_𝒩(z)+F_𝒩)^T\widehat{\mathrm{\Delta }z}+\delta _𝐮\right]+\mu \delta _𝐮.$$

By comparing this expression with the corresponding exact formula, which is

$$\mathrm{\Delta }\lambda _𝒩=D_𝒩^1\left[g_𝒩(z)\mathrm{\Lambda }_𝒩^1t_𝒩+g_𝒩(z)^T\mathrm{\Delta }z\right],$$

and by using the bounds (229) and the fact that $`f_𝒩`$ and $`F_𝒩`$ have size $`\delta _𝐮`$, we obtain

$`\widehat{\mathrm{\Delta }\lambda }_𝒩\mathrm{\Delta }\lambda _𝒩`$ $`=`$ $`\mu \delta _𝐮+(D_𝒩+G_𝒩)^1[f_𝒩+\delta _𝐮]+`$ (230)
$`[(D_𝒩+G_𝒩)^1D_𝒩^1][g_𝒩(z)\mathrm{\Lambda }_𝒩^1t_𝒩]+`$
$`(D_𝒩+G_𝒩)^1(g_𝒩(z)+F_𝒩)^T\widehat{\mathrm{\Delta }z}D_𝒩^1g_𝒩(z)^T\mathrm{\Delta }z`$
$`=`$ $`\mu \delta _𝐮+(D_𝒩+G_𝒩)^1[g_𝒩(z)^T(\widehat{\mathrm{\Delta }z}\mathrm{\Delta }z)+F_𝒩^T\widehat{\mathrm{\Delta }z}]`$
$`[(D_𝒩+G_𝒩)^1D_𝒩^1]g_𝒩(z)^T\mathrm{\Delta }z`$
$`=`$ $`\mu \delta _𝐮.`$
Finally, for the recovered step $`\widehat{\mathrm{\Delta }s}`$, we have from the last step of procedure condensed, together with (29b), (211l), (226b), and (15) that

$$\widehat{\mathrm{\Delta }s}=(g(z)+f+s)(g(z)+F)^T\widehat{\mathrm{\Delta }z}+\delta _𝐮,$$

where the final term accounts for the rounding error (15) that arises from accumulating the terms in the sum, which are all bounded. By substituting the expression for the exact $`\mathrm{\Delta }s`$ together with the estimates (211l) and (226b) on the sizes of the perturbation terms, we obtain

$`\widehat{\mathrm{\Delta }s}`$ $`=`$ $`(g(z)+s)g(z)^T\mathrm{\Delta }zfg(z)^T(\widehat{\mathrm{\Delta }z}\mathrm{\Delta }z)F^T\widehat{\mathrm{\Delta }z}+\mu \delta _𝐮`$ (231)
$`=`$ $`\mathrm{\Delta }s+\delta _𝐮.`$
We summarize the results obtained so far in the following theorem.

###### Theorem 10.

Suppose that Assumption 4.1 holds. Then when the step
$`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ is calculated in a finite-precision environment by using the procedure condensed (and where, in particular, a backward stable method is used to solve the linear system for the $`\widehat{\mathrm{\Delta }z}`$ component), we have that

$`(\mathrm{\Delta }z\widehat{\mathrm{\Delta }z},U^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{}),\mathrm{\Delta }s\widehat{\mathrm{\Delta }s})`$ $`=`$ $`\delta _𝐮,`$ (232a)
$`V^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{})`$ $`=`$ $`\delta _𝐮/\mu ,`$ (232b)
$`\mathrm{\Delta }\lambda _𝒩\widehat{\mathrm{\Delta }\lambda }_𝒩`$ $`=`$ $`\mu \delta _𝐮.`$ (232c)
This theorem extends the result of M. H. Wright  for accuracy of the computed solution of the condensed system by relaxing the LICQ assumption to MFCQ. When LICQ holds, the matrix $`V`$ is vacuous, so the absolute error in all components is of size at most $`\delta _𝐮`$. The higher accuracy (232c) of the components $`\widehat{\mathrm{\Delta }\lambda }_𝒩`$ (also noted in ) does not contribute significantly to the progress that can be made along the inexact direction $`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$, in the sense of Section 5.3.

We return briefly to the case discussed immediately after Corollary 9, in which the perturbations have the special form (188), using these results to show that the bound (232b) can be strengthened when $`f_{}`$ satisfies

$$V^Tf_{}=O(\mu ^2).$$
(233)
This case is of interest when the cancellation errors in computing $`g_{}(z)`$ are smaller than the estimate we made following (211l), possibly because of use of higher-precision arithmetic or the fact that the computation did not require differencing of quantities whose size is large relative to the final result. When (233) holds, we see by comparing (150) with (220) that

$$E_{23}=0,E_{33}=G_𝒩=\delta _𝐮/\mu ,f_2=UU^Tf_{}+O(\mu ^2),\text{where }f_{}=\delta _𝐮.$$

Therefore, we deduce from (209) that (225c) can be replaced by

$$V^T(\mathrm{\Delta }\lambda _{}y)=O(\mu ).$$

Using (227) and $`\mu \delta _𝐮`$, we can therefore replace (232b) in this case by

$$V^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{})=O(\mu ).$$
(234)
### 5.2 Termination of the Cholesky Algorithm

In deriving the estimate (225), we have assumed that a backward stable algorithm is used to solve (210). Because of (11), (12), and the SC condition, and the estimates of the sizes of the diagonals of $`D`$ (from (56) and Lemma 4), it is easy to show that the matrix in (210) is positive definite for all sufficiently small $`\mu `$. The Cholesky algorithm is therefore an obvious candidate for solving this system. However, the condition number of the matrix in (210) usually approaches $`\mathrm{}`$ as $`\mu 0`$, raising the possibility that the Cholesky algorithm may break down when $`\mu `$ is small. A simple argument, which we now sketch, suffices to show that successful completion of the Cholesky algorithm can be expected under the assumptions we have used in our analysis so far.

We state first the following technical result. Since it is similar to one proved by Debreu \[6, Theorem 3\], its proof is omitted.

###### Lemma 11.

Suppose that $`M`$ and $`A`$ are two matrices with the properties that $`M`$ is symmetric and

$$A^Tx=0x^TMx\alpha Mx^2,$$

for some constant $`\alpha >0`$. Then for all $`\mu `$ such that

$$0<\mu <\overline{\mu }\stackrel{\mathrm{def}}{=}\mathrm{min}(\frac{\alpha A^2}{4M},\frac{A}{\alpha M}),$$

we have that

$$x^T(M+\mu ^1AA^T)x\frac{\alpha }{2}x^2,\text{for all }x\text{.}$$

We apply this result to (210) by setting

$`M`$ $`=`$ $`_{zz}(z,\lambda )+g_𝒩(z)D_𝒩^1g_𝒩(z)^T=_{zz}(z,\lambda )+O(\mu ),`$
$`A`$ $`=`$ $`\mu ^{1/2}g_{}(z)D_{}^{1/2}`$
(where again we use (56) and Lemma 4 to derive the order estimates). The conditions (11), (12), and strict complementarity ensure that this choice of $`M`$ and $`A`$ satisfies the assumptions of Lemma 11. The result then implies that the smallest singular value of the matrix in (210) is positive and of size $`\mathrm{\Theta }(1)`$ for all values of $`\mu `$ below a threshold that is also of size $`\mathrm{\Theta }(1)`$. Since $`D=O(\mu ^1)`$, the largest eigenvalue of this matrix is of size $`O(\mu ^1)`$, so we have

$$\mathrm{cond}(_{zz}(z,\lambda )+g(z)D^1g(z)^T)=O(\mu ^1).$$
(235)
(An estimate similar to this is derived by M. H. Wright \[23, Theorem 3.2\], under the LICQ assumption.) It is known from a result of Wilkinson (cited by Golub and Van Loan \[12, p. 147\]) that the Cholesky algorithm runs to completion if $`q_n\delta _𝐮\mathrm{cond}()1`$, where $`q_n`$ is a modest quantity that depends polynomially on the dimension $`n`$ of the matrix. By combining this result with (235), we conclude that for the matrix in (210), we can expect completion of the Cholesky algorithm whenever $`\mu \delta _𝐮`$. That is, no new assumptions need to be added to those made in deriving the results of earlier sections.

We note that this situation differs a little from the case of linear programming where, because second-order conditions are not applicable, it is usually necessary to modify the Cholesky procedure to ensure that it runs to completion (see ).

### 5.3 Local Convergence with Computed Steps

We begin this section by showing how the quantities $`r_f`$, $`r_g`$, and $`\mu `$ change along the computed step $`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ obtained from the finite-precision implementation of the procedure condensed. We compare these with the changes that can be expected along the exact direction $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda ,\mathrm{\Delta }s)`$. We then consider the effects of these perturbations on an algorithm of the type in which the iterates are expected to satisfy the conditions (29). Rapidly convergent variants of these algorithms for linear programming problems usually allow the values of $`C`$ and $`\gamma `$ in these conditions to be relaxed, so that a near-unit step can be taken. We address the following question: If similar relaxations are allowed in an algorithm for nonlinear programming, are near-unit steps still possible when the steps contain perturbations of the type considered above?

We show in particular that for the computed search direction, the maximum step length that can be taken without violating the nonnegativity conditions on $`\lambda `$ and $`s`$ satisfies

$$1\widehat{\alpha }_{\mathrm{max}}=\delta _𝐮/\mu +O(\mu ),$$
(236)
while the reductions in pairwise products, $`\mu `$, $`r_f`$, and $`r_g`$ satisfy

$`(\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i)(s_i+\alpha \widehat{\mathrm{\Delta }s}_i)`$ $`=`$ $`(1\alpha )\lambda _is_i+\delta _𝐮+O(\mu ^2),i=1,2,\mathrm{},m,`$ (237a)
$`\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s})`$ $`=`$ $`(1\alpha )\mu +\delta _𝐮+O(\mu ^2),`$ (237b)
$`r_f(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })`$ $`=`$ $`(1\alpha )r_f+\delta _𝐮+O(\mu ^2),`$ (237c)
$`r_g(z+\alpha \widehat{\mathrm{\Delta }z},s+\alpha \widehat{\mathrm{\Delta }s})`$ $`=`$ $`(1\alpha )r_g+\delta _𝐮+O(\mu ^2).`$ (237d)
The corresponding maximum steplength for the exact direction satisfies

$$1\alpha _{\mathrm{max}}=O(\mu ),$$
(238)
while the reductions in $`r_f`$, $`r_g`$, and $`\mu `$ satisfy

$`(\lambda _i+\alpha \mathrm{\Delta }\lambda _i)(s_i+\alpha \mathrm{\Delta }s_i)`$ $`=`$ $`(1\alpha )\lambda _is_i+O(\mu ^2),i=1,2,\mathrm{},m,`$ (239a)
$`\mu (\lambda +\alpha \mathrm{\Delta }\lambda ,s+\alpha \mathrm{\Delta }s)`$ $`=`$ $`(1\alpha )\mu +O(\mu ^2),`$ (239b)
$`r_f(z+\alpha \mathrm{\Delta }z,\lambda +\alpha \mathrm{\Delta }\lambda )`$ $`=`$ $`(1\alpha )r_f+O(\mu ^2),`$ (239c)
$`r_g(z+\alpha \mathrm{\Delta }z,s+\alpha \mathrm{\Delta }s)`$ $`=`$ $`(1\alpha )r_g+O(\mu ^2).`$ (239d)
Our proof of the estimates (236) and (237) is tedious but not completely straightforward, and we have included it in the Appendix.

It is clear from (236) and (237) that the direction $`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ makes good progress toward the solution set $`𝒮`$. If the actual steplength $`\alpha `$ is close to its maximum value $`\widehat{\alpha }_{\mathrm{max}}`$, in the sense that

$$\widehat{\alpha }_{\mathrm{max}}\alpha =\delta _𝐮/\mu +O(\mu ),$$
(240)
we have by direct substitution in (236) and (237) that

$`\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s})`$ $`=`$ $`\delta _𝐮+O(\mu ^2),`$
$`r_f(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })`$ $`=`$ $`\delta _𝐮+O(\mu ^2),`$
$`r_g(z+\alpha \widehat{\mathrm{\Delta }z},s+\alpha \widehat{\mathrm{\Delta }s})`$ $`=`$ $`\delta _𝐮+O(\mu ^2).`$
These formulae suggest that finite precision does not have an observable effect on the quadratic convergence rate of the underlying algorithm until $`\mu `$ drops below about $`\sqrt{𝐮}`$. Stopping criteria for interior-point methods usually include a condition such as $`\mu 10^4𝐮`$ or $`\mu \sqrt{𝐮}`$ (see, for example, ), so that $`\mu `$ is not allowed to become so small that the assumption $`\mu 𝐮`$ made in (68) is violated.

In making this back-of-the-envelope assessment, however, we have not taken into account the approximate centrality conditions (29), which must continue to hold (possibly in a relaxed form) at the new iterate. These conditions play a central role both in the analysis above and in the convergence analysis of the underlying “exact” algorithms, and also appear to be important in practice. Typically (see, for example, Ralph and Wright ), the conditions (29) are relaxed by allowing a modest increase in $`C`$ and a modest decrease in $`\gamma `$ on the rapidly convergent steps. We show in the next result that enforcement of these relaxed conditions is not inconsistent with taking a step length $`\alpha `$ that is close to $`\widehat{\alpha }_{\mathrm{max}}`$, so that rapid convergence can still be observed even in the presence of finite-precision effects.

###### Theorem 12.

Suppose that Assumption 4.1 holds. Then when the step
$`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ is calculated in a finite-precision environment by using the procedure condensed, there is a constant $`\widehat{C}`$ such that for all $`\tau [0,1/2]`$ and all $`\alpha `$ satisfying

$$\alpha [0,1\widehat{C}\tau ^1(𝐮/\mu +\mu )],$$
(241)
the following relaxed form of the approximate centrality conditions holds:

$`r_f(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })`$ $``$ $`C(1+\tau )\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s}),`$ (242a)
$`r_g(z+\alpha \widehat{\mathrm{\Delta }z},s+\alpha \widehat{\mathrm{\Delta }s})`$ $``$ $`C(1+\tau )\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s}),`$ (242b)
$`(\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i)(s_i+\alpha \widehat{\mathrm{\Delta }s}_i)`$ $``$ $`\gamma (1\tau )\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s}),`$
for all $`i=1,2,\mathrm{},m`$,
where $`C`$ is the constant from conditions (29).
Moreover, when we set $`\alpha `$ to its upper bound in (241), we find that

$$\delta (z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })\tau ^1(\delta _𝐮+O(\mu ^2)).$$
(243)
###### Proof.

From (29) and (237), we have that

$`r_f(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })`$
$`=`$ $`(1\alpha )r_f+\delta _𝐮+O(\mu ^2)`$
$``$ $`C(1\alpha )\mu +\delta _𝐮+O(\mu ^2)`$
$`=`$ $`C(1+\tau )(1\alpha )\mu C\tau (1\alpha )\mu +\delta _𝐮+O(\mu ^2)`$
$`=`$ $`C(1+\tau )\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s})C\tau (1\alpha )\mu +\delta _𝐮+O(\mu ^2).`$
We deduce that the required condition (242a) will hold provided that

$$\delta _𝐮+O(\mu ^2)C\tau (1\alpha )\mu .$$

Since by definition we have that $`\delta _𝐮+O(\mu ^2)\overline{C}(𝐮+\mu ^2)`$ for some positive constant $`\overline{C}`$, we find that a sufficient condition for the required inequality is that

$$(1\alpha )(\overline{C}/C)\tau ^1(𝐮/\mu +\mu ),$$

which is equivalent to (241) for an obvious definition of $`\widehat{C}`$. Identical logic can be applied to $`r_g`$ to yield a similar condition on $`\alpha `$.

For the condition (242), we have from (29) and (237) that

$`(\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i)(s_i+\alpha \widehat{\mathrm{\Delta }s}_i)`$
$`=`$ $`(1\alpha )\lambda _is_i+\delta _𝐮+O(\mu ^2)`$
$``$ $`(1\alpha )\gamma \mu +\delta _𝐮+O(\mu ^2)`$
$`=`$ $`\gamma (1\tau )(1\alpha )\mu +\gamma \tau (1\alpha )\mu +\delta _𝐮+O(\mu ^2)`$
$`=`$ $`\gamma (1\tau )\mu (\lambda +\alpha \widehat{\mathrm{\Delta }\lambda },s+\alpha \widehat{\mathrm{\Delta }s})+\gamma \tau (1\alpha )\mu +\delta _𝐮+O(\mu ^2).`$
Hence, the condition (242) holds provided that

$$\gamma \tau (1\alpha )\mu +\delta _𝐮+O(\mu ^2)0.$$

Similar logic can be applied to this inequality to derive a bound of the type (241), after a possible adjustment of $`\widehat{C}`$.

Finally, we obtain (243) by substituting $`\alpha =1\widehat{C}\tau ^1(𝐮/\mu +\mu )`$ into (237) and applying Theorem 5. (Note that, despite the relaxation of the centrality conditions (242), the result of Theorem 5 still holds; we simply modify the proof to replace $`C`$ by $`(3/2)C`$ in (29a) and (29b), and $`\gamma `$ by $`\gamma /2`$ in (29c).) ∎

## 6 The Augmented System

In this section, we consider the case in which the augmented system (27) (equivalently, (55)) is solved to obtain $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda )`$, while the remaining step component $`\mathrm{\Delta }s`$ is recovered from (26). The formal specification for this procedure is as follows:

| procedure augmented |
| --- |
| given the current iterate $`(z,\lambda ,s)`$ |
|  | form the coefficient matrix and right-hand side for (55); |
|  | solve (55) to obtain $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda )`$; |
|  | set $`\mathrm{\Delta }s=(g(z)+s)g(z)^T\mathrm{\Delta }z`$. |

Much of our work in analyzing the augmented system form (55) has already been performed in Section 4; the main error result is simply Corollary 9. However, we can apply this result only if the floating-point errors made in evaluating and solving this system satisfy the assumptions of this corollary. In particular, we need to show that the perturbation matrices $`E_{ij}`$, $`i,j=1,2,3`$ in (150) satisfy the estimates (134).

This task is not completely straightforward. Unlike the condensed and full-system cases, it is not simply a matter of assuming that a backward-stable algorithm has been used to solve the system (55). The reason is that the largest terms in the coefficient matrix in (27)—the diagonal elements in the matrix $`D_𝒩`$—have size $`O(\mu ^1)`$. The usual analysis of backward-stable algorithms represents the floating-point errors as a perturbation of the entire coefficient matrix whose size is bounded by $`\delta _𝐮`$ times the matrix norm—in this case, $`\delta _𝐮/\mu `$. However, Corollary 9 requires some elements of the perturbation matrix to be smaller than this estimate; in particular, the submatrices $`E_{12}`$, $`E_{21}`$, and $`E_{22}`$ need to be of size $`\delta _𝐮`$. Therefore, we need to look closely at the particular algorithms used to solve (55) to see whether they satisfy the following condition.

###### Condition 6.1.

The solution obtained by applying the algorithm in question to the system (55) in floating-point arithmetic is the exact solution of a perturbed system in which the perturbations of the coefficient matrix satisfy the estimates (134), while the right-hand side is unperturbed.

We focus on diagonal pivoting methods, which take a symmetric matrix $`T`$ and produce a factorization of the form

$$PTP^T=LYL^T,$$
(244)
where $`P`$ is a permutation matrix, $`L`$ is unit lower triangular, and $`Y`$ is block diagonal, with a combination of $`1\times 1`$ and symmetric $`2\times 2`$ blocks. The best-known methods of this class are due to Bunch and Parlett  and Bunch and Kaufman , while Duff et al.  and Fourer and Mehrotra  have described sparse variants. These algorithms differ in their selection criteria for the $`1\times 1`$ and $`2\times 2`$ pivot blocks. In our case, the presence of the diagonal elements of size $`\mathrm{\Theta }(\mu ^1)`$ (from the submatrix $`D_𝒩=\mathrm{\Lambda }_𝒩^1S_𝒩`$) and their place in these pivot blocks are crucial to the result.

We start by stating a general result of Higham  concerning backward stability that applies to all diagonal pivoting schemes. We then examine the Bunch-Kaufman scheme, showing that the large diagonals appear only as $`1\times 1`$ pivots and that this algorithm satisfies Condition 6.1. (In \[17, Theorem 4.2\], Higham actually proves that the Bunch-Kaufman scheme is backward stable in the normwise sense, but this result is not applicable to our context, for the reasons mentioned above.)

Next, we briefly examine the Bunch-Parlett method, showing that it starts out by selecting all the large diagonal elements in turn as $`1\times 1`$ pivots, before going on to factor the remaining matrix, whose elements are all $`O(1)`$ in size. This method also satisfies Condition 6.1. We then examine the sparse diagonal pivoting approaches of Duff et al.  and Fourer and Mehrotra , which may not satisfy Condition 6.1, because of the possible presence of $`2\times 2`$ pivots in which one of the diagonals has size $`\mathrm{\Theta }(\mu ^1)`$. These algorithms can be modified in simple ways to overcome this difficulty, possibly at the expense of higher density in the $`L`$ factor. We then mention Gaussian elimination with pivoting and refer to previous results in the literature to show that this approach satisfies Condition 6.1. Finally, we state a result like Theorem 12 about convergence of a finite-precision implementation of an algorithm based on the augmented system form.

Higham \[17, Theorem 4.1\] proves the following result.

###### Theorem 13.

Let $`T`$ be an $`\overline{n}\times \overline{n}`$ symmetric matrix, and let $`\widehat{x}`$ be the computed solution to the linear system $`Tx=b`$ produced by a method that yields a factorization of the form (244), with any diagonal pivoting strategy. Assume that, during recovery of the solution, the subsystems that involve the $`2\times 2`$ diagonal blocks are solved via Gaussian elimination with partial pivoting. Then we have that

$$(T+\mathrm{\Delta }T)\widehat{x}=b,|\mathrm{\Delta }T|\delta _𝐮(|T|+P^T|\widehat{L}||\widehat{Y}||\widehat{L}^T|P)+\delta _𝐮^2,$$
(245)
where $`\widehat{L}`$ and $`\widehat{Y}`$ are the computed factors, and $`|A|`$ denotes the matrix formed from $`A`$ by replacing all its elements by their absolute values.

In Higham’s result, the coefficient of $`𝐮`$ in the $`\delta _𝐮`$ term is actually a linear polynomial in the dimension of the system. The partial pivoting strategy for the $`2\times 2`$ systems can actually be replaced by any method for which the computed solution of $`Ry=d`$ satisfies $`(R+\mathrm{\Delta }R)\widehat{y}=d`$, where $`R`$ is the $`2\times 2`$ matrix in question and $`|\mathrm{\Delta }R|\delta _𝐮|R|`$. This property was also key in an earlier paper of S. Wright , who derived a result similar to Theorem 13 in the context of the augmented systems that arise from interior-point methods for linear programming.

All the procedures below have the property that the growth in the maximum element size in the remaining submatrix is bounded by a modest quantity at each individual step of the factorization. (In the case of Bunch-Kaufman and Bunch-Parlett, this bound averages about 2.6 per elimination step; see Golub and Van Loan \[12, Section 4.4.4\].) As with Gaussian elimination with partial pivoting, exponential element growth is possible, so that $`L`$ and $`Y`$ in (244) contain much larger elements than the original matrix $`T`$. Such behavior is, however, quite rare and is confined to pathological cases and certain special problem classes. In our analysis below, we make the safe assumption that catastrophic growth of this kind does not occur.

### 6.1 The Bunch-Kaufman Procedure

At each iteration, the Bunch-
Kaufman procedure chooses either a $`1\times 1`$ or $`2\times 2`$ pivot by examining at most two columns of the remaining matrix, that is, the part of the matrix that remains to be factored at this stage of the process. It makes use of quantities $`\chi _i`$ defined by

$$\chi _i=\underset{j|ji}{\mathrm{max}}|T_{ij}|,$$

where in this case $`T`$ denotes the remaining matrix. We define the pivot selection strategy for the first step of the factorization process. The entire algorithm is obtained by applying this procedure recursively to the remaining submatrix.

| set $`\nu =(1+\sqrt{17})/8`$; |
| --- |
| calculate $`\chi _1`$, and store the index $`r`$ for which $`\chi _1=|T_{r1}|`$; |
| if $`|T_{11}|\nu \chi _1`$ |
|  | choose $`T_{11}`$ as a $`1\times 1`$ pivot; |
| else |
|  | calculate $`\chi _r`$; |
|  | if $`\chi _r|T_{11}|\nu \chi _1^2`$ |
|  |  | choose $`T_{11}`$ as a $`1\times 1`$ pivot; |
|  | else if $`|T_{rr}|\nu \chi _r`$ |
|  |  | choose $`T_{rr}`$ as a $`1\times 1`$ pivot; |
|  | else |
|  |  | choose a $`2\times 2`$ pivot with diagonals $`T_{11}`$ and $`T_{rr}`$; |
|  | end if |
| end if. |

For each choice of pivot, the permutation matrix $`P_1`$ is chosen so that the desired $`1\times 1`$ or $`2\times 2`$ pivot is in the upper left of the matrix $`P_1TP_1^T`$. If one writes

$$P_1TP_1^T=\left[\begin{array}{cc}R& C^T\\ C& \widehat{T}\end{array}\right],$$

where $`R`$ is the chosen pivot, the first step of the factorization yields

$$P_1TP_1^T=\left[\begin{array}{cc}I& \\ CR^1& I\end{array}\right]\left[\begin{array}{cc}R& \\ & \overline{T}\end{array}\right]\left[\begin{array}{cc}I& R^1C^T\\ & I\end{array}\right],$$
(246)
where $`\overline{T}=\widehat{T}CR^1C^T`$ is the matrix remaining after this stage of the factorization.

At the first step of the factorization, the quantities $`\chi _1`$ and $`\chi _r`$ (if calculated) both have size $`O(1)`$, since the large elements of this matrix occur only on the diagonal. Since a $`2\times 2`$ pivot is chosen only if

$$|T_{11}|<\nu \chi _1\text{and}|T_{rr}|<\nu \chi _r,$$

it follows immediately that both diagonals in a $`2\times 2`$ pivot must be $`O(1)`$. Hence, the pivot chosen by this procedure is one of three types:

$`1\times 1`$ pivot of size $`O(1)`$; (247a)
$`2\times 2`$ pivot in which both diagonals have size $`O(1)`$; (247b)
$`1\times 1`$ pivot of size $`\mathrm{\Theta }(\mu ^1)`$. (247c)
In fact, the pivots are one of the types (247) at all stages of the factorization, not just the first stage. The reason is that the updated matrix $`\overline{T}`$ in (246) has the same essential form as the original matrix $`T`$—its elements are all of size $`O(1)`$ except for some large diagonal elements of size $`\mathrm{\Theta }(\mu ^1)`$. We demonstrate this claim by showing that the update $`CR^1C^T`$ that is applied to the remaining matrix in (246) is a matrix whose elements are of size at most $`O(1)`$, regardless of the type of pivot, so that it does not disturb the essential structure of the remaining matrix. When the pivots are of type (247a) and (247b), the standard argument of Bunch and Kaufman  can be applied to show that the norm of $`CR^1C^T`$ is at most a modest multiple of $`C`$. We know that $`C=O(1)`$, since $`C`$ consists only of off-diagonal elements, so we conclude that $`CR^1C^T=O(1)`$ in this case as claimed. For the other pivot type (247c), we have $`R=\mathrm{\Theta }(\mu ^1)`$ and $`C=O(1)`$, so the elements of $`CR^1C^T`$ have size $`O(\mu )`$, and the claim holds in this case too.

In the rest of this subsection, we show by using Theorem 13 that Condition 6.1 holds for the Bunch-Kaufman algorithm. In fact, we prove a stronger result: When $`T`$ in Theorem 13 is the matrix (55), the perturbation matrix $`\mathrm{\Delta }T`$ contains elements of size $`\delta _𝐮`$, except in those diagonal locations corresponding to the elements of $`D_𝒩`$, where they may be as large as $`\delta _𝐮/\mu `$. Given the bound on $`|\mathrm{\Delta }T|`$ in (245), we need only to show that $`P^T|\widehat{L}||\widehat{Y}||\widehat{L}|^TP`$ has the desired structure. In fact, it suffices to show that the exact factor product $`P^T|L||Y||L|^TP`$ has the structure in question, since the difference between these two products is just $`\delta _𝐮`$ in size.

We demonstrate this claim inductively, using a refined version of the arguments from Higham \[17, Section 4.3\] for some key points, and omitting some details. For simplicity, and without loss of generality, we assume that $`P=I`$.

When $`\overline{n}=1`$ (that is, $`T`$ is $`1\times 1`$), we have that $`L=1`$ and $`Y=T`$, and the result holds trivially. When $`\overline{n}=2`$, there are three cases to consider. If the matrix contains no elements of size $`\mathrm{\Theta }(\mu ^1)`$, then the analysis for general matrices can be used to show that $`|L||Y||L|^T=O(1)`$, as required. If either or both diagonal elements have size $`\mathrm{\Theta }(\mu ^1)`$, then both pivots are $`1\times 1`$, and the factors have the form

$$L=\left[\begin{array}{cc}1& 0\\ T_{21}/T_{11}& 1\end{array}\right],Y=\left[\begin{array}{cc}T_{11}& 0\\ 0& T_{22}T_{21}^2/T_{11}\end{array}\right].$$
(248)
Two cases arise.

* A diagonal of size $`O(1)`$ is accepted as the first pivot and moved (if necessary) to the $`(1,1)`$ position. We then have

$$|T_{11}|\nu \chi _1=\nu \chi _2=\nu |T_{21}|,$$

and therefore $`|T_{21}/T_{11}|1/\nu `$ and hence $`|T_{21}^2/T_{11}||T_{21}|/\nu =O(1)`$. If the $`(2,2)`$ diagonal is also $`O(1)`$, we have that $`L=O(1)`$ and $`Y=O(1)`$, and we are done. Otherwise, $`T_{22}=\mathrm{\Theta }(\mu ^1)`$, and so the $`(2,2)`$ element of $`Y`$ satisfies this same estimate. We conclude from (248) that $`|L||Y||L|^T`$ also has an $`\mathrm{\Theta }(\mu ^1)`$ element in the $`(2,2)`$ position and $`O(1)`$ elements elsewhere.
* A diagonal of size $`\mathrm{\Theta }(\mu ^1)`$ is accepted as the first pivot and moved (if necessary) to the $`(1,1)`$ position. We then have

$$|T_{21}/T_{11}|=O(\mu ),|T_{21}^2/T_{11}|=O(\mu ).$$

It follows from (248) that

$$|L||Y||L|^T=\left[\begin{array}{cc}|T_{11}|& |T_{21}|\\ |T_{21}|& |T_{22}|+O(\mu )\end{array}\right],$$

which obviously has the desired structure.

We now assume that our claim holds for some dimension $`\overline{n}2`$, and we prove that it continues to hold for dimension $`\overline{n}+1`$. Using the notation of (246) (assuming that $`P_1=I`$), and denoting the factorization of the Schur complement $`\overline{T}`$ in (246) by $`\overline{T}=\overline{L}\overline{Y}\overline{L}^T`$, we have that

$$T=LYL^T=\left[\begin{array}{cc}I& \\ CR^1& \overline{L}\end{array}\right]\left[\begin{array}{cc}R& \\ & \overline{Y}\end{array}\right]\left[\begin{array}{cc}I& R^1C^T\\ & \overline{L}^T\end{array}\right].$$
(249)
It follows that

$$|L||Y||L|^T=\left[\begin{array}{cc}|R|& |R||R^1C^T|\\ |CR^1||R|& |CR^1||R||R^1C^T|+|\overline{L}||\overline{Y}||\overline{L}|^T\end{array}\right].$$
(250)
Since, as we mentioned above, the norm of $`CR^1C^T`$ is at most $`O(1)`$, the Schur complement $`\overline{T}=\widehat{T}CR^1C^T`$ has size $`O(1)`$ except for large $`\mathrm{\Theta }(\mu ^1)`$ elements in the same locations as in the original matrix. Hence, by our inductive hypothesis, $`|\overline{L}||\overline{Y}||\overline{L}|^T`$ has a similar structure, and we need to show only that the effects of the first step of the factorization (246) do not disturb the desired structure.

For the case in which $`R`$ is a pivot of type either (247a) and (247b), Higham \[17, Section 4.3\] shows all elements of both $`|CR^1||R|`$ and $`|CR^1||R||R^1C^T|`$ are bounded by a modest multiple of either $`\chi _1`$ (if $`T_{11}`$ was selected as the pivot because it passed the test $`|T_{11}|\nu \chi _1`$) or $`(\chi _1+\chi _r)`$, where $`r`$ is the “other” column considered during the selection process. In our case, this observation implies that both $`|CR^1||R|`$ and $`|CR^1||R||R^1C^T|`$ have size $`O(1)`$. By combining these observations with those of the preceding paragraph, we conclude that for pivots of types (247a) and (247b), “large” elements of the matrix in (250) occur only in the diagonal locations originally occupied by $`D_𝒩`$.

For the remaining case—pivots of type (247c)—we have that $`C`$ has size $`O(1)`$ while $`R^1`$ has size $`O(\mu )`$. Therefore, $`|CR^1||R|`$ has size $`O(1)`$ and $`|CR^1||R||R^1C^T|`$ has size $`O(\mu )`$, while $`|R|`$, which occupies the $`(1,1)`$ position in the matrix (250), just as it did in the original matrix $`T`$, has size $`\mathrm{\Theta }(\mu ^1)`$. We conclude that the desired structure holds in this case as well.

We conclude from this discussion that Condition 6.1 holds for the Bunch-Kaufman procedure. We show later that the perturbations arising from other sources, namely, roundoff and cancellation in the evaluation of the matrix and right-hand side, also satisfy the conditions of Corollary 9, so this result can be used to bound the error in the computed steps.

Finally, we note that it is quite possible for pivots of types (247a) and (247b) to be chosen while diagonal elements of size $`\mathrm{\Theta }(\mu ^1)`$ still remain in the submatrix. Therefore, a key assumption of the analysis of Forsgren, Gill, and Shinnerl \[9, Theorem 4.4\]—namely, that all the diagonals of size $`\mathrm{\Theta }(\mu ^1)`$ are chosen as $`1\times 1`$ pivots before any of the other diagonals are chosen—may not be satisfied by the Bunch-Kaufman procedure.

### 6.2 The Bunch-Parlett Procedure

The Bunch-Parlett procedure is conceptually simpler but more expensive to implement than Bunch-Kaufman, since it requires $`O(n^2)`$ (rather than $`O(n)`$) comparisons at each step of the factorization. The pivot selection strategy is as follows.

| set $`\nu =(1+\sqrt{17})/8`$; |
| --- |
| calculate $`\chi _{\mathrm{off}}=|T_{rs}|=\mathrm{max}_{ij}|T_{ij}|`$, $`\chi _{\mathrm{diag}}=|T_{pp}|=\mathrm{max}_i|T_{ii}|`$; |
| if $`\chi _{\mathrm{diag}}\nu \chi _{\mathrm{off}}`$ |
|  | choose $`T_{pp}`$ as the $`1\times 1`$ pivot; |
| else |
|  | choose the $`2\times 2`$ pivot whose off-diagonal element is $`T_{rs}`$; |
| end if. |

The elimination procedure then follows as in (246).

It is easy to show that the Bunch-Parlett procedure starts by selecting all the diagonals of size $`\mathrm{\Theta }(\mu ^1)`$ in turn as $`1\times 1`$ pivots. (Because of this property, it satisfies the key assumption of  mentioned at the end of the preceding section.) The update $`CR^1C^T`$ generated by each of these pivot steps has size only $`O(\mu )`$, so the matrix that remains after this phase of the factorization contains only $`O(1)`$ elements. The remaining pivots are then a combination of types (247a) and (247b).

By using the arguments of the preceding subsection in a slightly simplified form, we can show that Condition 6.1 holds for this procedure as well.

### 6.3 Sparse Diagonal Pivoting

For large instances of (1), the Bunch-
Kaufman and Bunch-Parlett procedures are usually inefficient because they do not try to maintain sparsity in the lower triangular factor $`L`$. Sparse variants of these algorithms, such as those of Duff et al.  and Fourer and Mehrotra , use pivot selection strategies that combine stability considerations with Markowitz-like estimates of the amount of fill-in that a candidate pivot will cause in the remaining matrix.

At each stage of the factorization, both algorithms examine a roster of possible $`1\times 1`$ and $`2\times 2`$ pivots, starting with those that would create the least fill-in. As soon as a pivot is found that meets the stability criteria described below, it is accepted. Both algorithms prefer to use $`1\times 1`$ pivots where possible.

For candidate $`1\times 1`$ pivots, Duff et al. \[7, p. 190\] use the following stability criterion:

$$|R^1|C_{\mathrm{}}\rho ,$$
(251)
where the notation $`R`$ and $`C`$ is from (246) and $`\rho [2,\mathrm{})`$ is some user-selected parameter that represents the tolerable growth factor at each stage of the factorization. For a $`2\times 2`$ pivot, the criterion is

$$|R^1|\left[\begin{array}{c}C_{,1}_{\mathrm{}}\\ C_{,2}_{\mathrm{}}\end{array}\right]\left[\begin{array}{c}\rho \\ \rho \end{array}\right],$$
(252)
where $`C_{,1}`$ and $`C_{,2}`$ are the two columns of $`C`$. The stability criteria of Fourer and Mehrotra  are similar.

As they stand, the stability tests (251) and (252) do not necessarily restrict the choice of pivots to the three types (247). If a $`1\times 1`$ pivot of size $`\mathrm{\Theta }(\mu ^1)`$ is ever considered for structural reasons, it will pass the test (251) (the left-hand side of this expression will have size $`O(\mu )`$) and therefore will be accepted as a pivot. However, it is possible that $`2\times 2`$ pivots in which one or both diagonals have size $`\mathrm{\Theta }(\mu ^1)`$ may pass the test (252) and may therefore be accepted. Although the test (252) ensures that the size of the update $`CR^1C^T`$ is modest (so that the update $`\overline{T}=\widehat{T}CR^1C^T`$ does not disturb the large-diagonal structure of $`\widehat{T}`$), there is no obvious assurance that the matrix $`|L||Y||L|^T`$ in (250) mirrors the structure of $`|T|`$, in terms of having the large diagonal elements in the same locations. The terms $`|CR^1||R|`$ and $`|CR^1||R||R^1C^T|`$ in (250) may not have size $`O(1)`$, as they do for pivots of the three types (247) arising from the Bunch-Kaufman and Bunch-Parlett selection procedures.

The Fourer-Mehrotra algorithm does, however, rule out the possibility of a $`2\times 2`$ pivot in which both diagonals are of size $`\mathrm{\Theta }(\mu ^1)`$. It considers a $`2\times 2`$ candidate only if one of its diagonal elements has previously been considered as a $`1\times 1`$ pivot but failed the stability test. However, if either of the diagonals had been subjected to the test (251), they would have been accepted, as noted in the preceding paragraph, so this situation cannot occur.

If the sparse algorithms are modified to ensure that all pivots have one of the three types (247), and all continue to satisfy the stability tests (251) or (252), then simple arguments (simpler than those of Section 6.1!) can be applied to show that Condition 6.1 is satisfied. One possible modification that achieves the desired efffect is to require that a $`2\times 2`$ pivot be allowed only if both its diagonals have previously been considered as $`1\times 1`$ pivots but failed the stability test (251).

### 6.4 Gaussian Elimination

Another possibility for solving the system (55) is to ignore its symmetry and apply a Gaussian elimination algorithm, with row and/or column pivoting to preserve sparsity and prevent excessive element growth. Such a strategy satisfies Condition 6.1. In , the author uses a result of Higham  to show that the effects of the large diagonal elements are essentially confined to the columns in which they appear. Assuming that the pivot sequence is chosen to prevent excessive element growth in the remaining matrix, and using the notation of (133) and (134), we can account for the effects of roundoff error in Gaussian elimination with perturbations in the coefficient matrix that satisfy the following estimates:

$$E_{11},E_{21},E_{31},E_{12},E_{22},E_{32}=\delta _𝐮,E_{13},E_{23},E_{33}=\delta _𝐮/\mu .$$

These certainly satisfy the conditions (134), so Condition 6.1 holds.

### 6.5 Local Convergence with the Computed Steps

We can now state a formal result to show that when the evaluation errors are taken into account as well as the roundoff errors from the factorization/solve procedure discussed above, the accuracies of the computed steps obtained from the procedure augmented, implemented in finite precision, satisfy the same estimates as for the corresponding steps obtained from the procedure condensed. The result is analogous to Theorem 10.

###### Theorem 14.

Suppose that Assumption 4.1 holds. Then when the step
$`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ is calculated in a finite-precision environment by using the procedure augmented, where the algorithm used to solve (55) satisfies Condition 6.1, we have

$`(\mathrm{\Delta }z\widehat{\mathrm{\Delta }z},U^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{}),\mathrm{\Delta }s\widehat{\mathrm{\Delta }s})`$ $`=`$ $`\delta _𝐮,`$ (253a)
$`V^T(\mathrm{\Delta }\lambda _{}\widehat{\mathrm{\Delta }\lambda }_{})`$ $`=`$ $`\delta _𝐮/\mu ,`$ (253b)
$`\mathrm{\Delta }\lambda _𝒩\widehat{\mathrm{\Delta }\lambda }_𝒩`$ $`=`$ $`\delta _𝐮.`$ (253c)
###### Proof.

The proof follows from Corollary 9 when we show that the perturbations to (55) from all sources—evaluation of the matrix and right-hand side as well as the factorization/solution procedure—satisfy the bounds required by this earlier result.

Because of Condition 6.1, perturbations arising from the factorization/solution procedure satisfy the bounds (134). The expressions (211) show that the errors arising from evaluation of $`_{zz}(z,\lambda )`$, $`_z(z,\lambda )`$, $`g(z)`$, and $`g(z)`$ are all of size $`\delta _𝐮`$, and hence they too satisfy the required bounds. Similarly to (212), evaluation of $`D_{}`$ and $`D_𝒩`$ yields errors of relative size $`\delta _𝐮`$, that is,

$`\text{computed }D_{}`$ $``$ $`D_{}+G_{},G_{}=\mu \delta _𝐮,`$ (254a)
$`\text{computed }D_𝒩`$ $``$ $`D_𝒩+G_𝒩,G_𝒩=\delta _𝐮/\mu ,`$ (254b)
where $`G_{}`$ and $`G_𝒩`$ are diagonal matrices.

We now obtain all the estimates in (253) by a direct application of Corollary 9, with the exception of the estimate for $`(\mathrm{\Delta }s\widehat{\mathrm{\Delta }s})`$. Since the expressions for recovering $`\mathrm{\Delta }s`$ are identical in procedures condensed and augmented, we can apply expression (231) from Section 5.1 to deduce that the desired estimate holds for this component as well. ∎

The only difference between the error estimates of Theorem 10 for the condensed system and those obtained above for the augmented system is that the $`\widehat{\mathrm{\Delta }\lambda }_𝒩`$ components are slightly less accurate in the augmented case. If we work through the analysis of Section 5.3 with the estimate (253c) replacing (232c), we find that the main results are unaffected. Therefore, we conclude this section by stating without proof a result similar to Theorem 12.

###### Theorem 15.

Suppose that all the assumptions of Theorem 12 hold, except that the step $`(\widehat{\mathrm{\Delta }z},\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})`$ is calculated by using the procedure augmented with a factorization/solution algorithm that satisfies Condition 6.1. Then the conclusions of Theorem 12 hold.

## 7 Numerical Illustration

We illustrate the results of Sections 5 and 6 using the two-variable example (13). Consider a simple algorithm that takes steps satisfying (26) with $`t`$ set rather arbitrarily to $`t=\mu ^2e`$. (The search directions thus used are like those generated in the later stages of a practical primal-dual algorithm such as Mehrotra’s algorithm .) We start this algorithm from the point

$$z_0=(1/30,1/9)^T,\lambda _0=(1,1/5)^T,s_0=(1/10,1/2).$$

(It is easy to check that the conditions (29) are satisfied at this point for a modest value of $`C`$.) At each step we calculated $`\widehat{\alpha }_{\mathrm{max}}`$, defined in Section 5.3, and took an actual step of $`.99\widehat{\alpha }_{\mathrm{max}}`$.

We programmed the method in Matlab, using double-precision arithmetic. In our first experiment, we solved the formulation (55) of the linear equations using Matlab’s standard Gaussian elimination solver for general systems of linear equations, which was analyzed in Section 6.4. From Theorem 14, the estimates (253) apply to this case.

Results are tabulated in Table 1. Note first the size of the component $`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$, which grows as $`\mu `$ decreases below $`𝐮^{1/2}`$, in accordance with (253b). (We cannot tabulate the difference $`V^T(\widehat{\mathrm{\Delta }\lambda }_{}\mathrm{\Delta }\lambda _{})`$ because of course we do not know the true step $`(\mathrm{\Delta }z,\mathrm{\Delta }\lambda ,\mathrm{\Delta }s)`$, but since the true step has size $`O(\mu )`$ (Corollary 8), the error is dominated by the term $`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$ in any case.) As predicted by (236), the maximum step $`\widehat{\alpha }_{\mathrm{max}}`$ becomes significantly smaller than $`1`$ as $`\mu `$ is decreased below $`𝐮^{1/2}`$. As indicated by (237), however, good progress still can be made along this direction (in the sense of reducing $`\mu `$ and the norms of the residuals $`r_f`$ and $`r_g`$) almost until $`\mu `$ reaches the level of $`𝐮`$. In fact, between iterations 5 and 8 we see the reduction factor of $`100`$ that we would expect by moving a distance of $`.99`$ along a direction that is close to the pure Newton direction. The component with the large error—$`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$—does not interfere significantly with rapid convergence, but only causes the $`\lambda `$ iterates to move tangentially to $`𝒮_\lambda `$. This effect may be noted in the final iterate where the value of $`\lambda `$ changes significantly. In some cases, however, when the current $`\lambda `$ is near the edge of the set $`𝒮_\lambda `$, this error may result in a severe curtailment of the step length.

Next, we performed the same experiment using the condensed formulation (3.1) of the linear system, as described in Section 5. Results are shown in Table 2. The main difference with Table 1 is that there is no increase in the value $`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$ as $`\mu `$ approaches unit roundoff; this component appears to decrease at the same rate as the other step components. This observation can be explained by our analysis of the case in which the cancellation error term $`f_{}`$ incurred in the evaluation of $`g_{}(z)`$ satisfies (233). We calculated the product $`V^T(g_{}(z)+f_{})`$ (the product of $`V`$ with our computed version of $`g_{}(z)`$) and found it to be exactly zero on iterations 7, 8, and 9. Therefore, using Taylor’s theorem, (18), and Theorem 5, we have

$$V^Tf_{}=V^Tg_{}(z)=V^Tg_{}(z^{})(zz^{})+O(zz^{}^2)=O(\mu ^2).$$

Hence, (234) together with Corollary 8 shows that $`V^T\widehat{\mathrm{\Delta }\lambda }_{}=O(\mu )`$, which is consistent with the results in Table 2. Note too that because of the higher accuracy in the $`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$ component, the maximum step length stays very close to $`1`$ during the last few iterations. By comparing Tables 1 and 2, however, we can verify that the convergence of $`\mu `$ to zero, and of the iterates to the solution set, is not materially affected by the presence or absence of the large error in $`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$.

To show that the lack of cancellation effects in Table 2 cannot be assumed in general, we modified problem (13) slightly, changing the second constraint to

$$g_2(z)\stackrel{\mathrm{def}}{=}\frac{2}{3\sqrt{5}}(z_1\sqrt{5})^2+z_2^2\frac{2\sqrt{5}}{3}0.$$
(255)
The primal and dual solutions remain unchanged, and we ran the condensed-equations version of the algorithm from the same starting point as above. Results are shown in Table 3. We observed that $`g_{}(z)`$ did not escape cancellation errors in this instance and, as in Table 1, we observe significant errors in $`V^T\widehat{\mathrm{\Delta }\lambda }_{}`$ that do not materially affect the convergence of the algorithm to the solution set.

## 8 Summary and Conclusions

In this paper, we have analyzed the finite-precision implementation of a primal-dual interior point method whose convergence rate is theoretically superlinear. We have made the standard assumptions that appear in most analyses of local convergence of nonlinear programming algorithms and path-following algorithms, with one significant exception: The assumption of linearly independent active constraint gradients is replaced by the weaker Mangasarian-Fromovitz constraint qualification, which is equivalent to boundedness of the set of optimal Lagrange multipliers. Because of this assumption, it is possible that all reasonable formulations of the step equations—the linear system that needs to be solved to obtain the search direction—are ill conditioned, so it is not obvious that the numerical errors that occur when this system is solved in finite precision do not eventually render the computed search direction useless. We show that although the error in the computed step may indeed become large as $`\mu `$ decreases, most of the error is restricted to a subspace that does not matter, namely, the null space of the matrix $`g_{}(z^{})`$ of first derivatives of the active constraints. Although this error causes the computed iterates to “slip” in a tangential direction to the optimal Lagrange multiplier set, it does not interfere with rapid convergence of the iterates to the primal-dual solution set.

We found that the centrality conditions (29), which are usually applied in path-following methods, played a crucial role in the analysis, since they enabled us to establish the estimates (34) in Lemma 4 concerning the sizes of the basic and nonbasic components of $`s`$ and $`\lambda `$ near the solution set. The analysis of Section 4, culminating in Corollary 9, finds bounds on the errors induced in step components by certain structured perturbations of the step equations. We show in the same section that the exact step is $`O(\mu )`$, allowing the local convergence analysis of Ralph and Wright  to be extended from convex programs to nonlinear programs.

In Sections 5 and 6 we apply the general results of Section 4 to the two most obvious ways of formulating and solving the step equations; namely, as a “condensed” system involving just the primal variables $`z`$, or as an “augmented” system involving both $`z`$ and the Lagrange multipliers $`\lambda `$. In each case, the errors introduced in finite-precision implementation have the structure of the perturbations analyzed in Section 4, so the error bounds obtained in Corollary 9 apply. In Section 5.3 (whose analysis also applies to the computed solutions analyzed in Section 6), we show that the potentially large error component discussed above does not interfere appreciably with the near-linear decrease of the quantities $`\mu `$, $`r_f`$, and $`r_g`$ to zero along the computed steps, indicating that until $`\mu `$ becomes quite close to $`𝐮`$, the convergence behavior predicted by the analysis of the “exact” algorithm will be observed in the finite-precision implementation. We conclude in Section 7 with a numerical illustration of our major observations on a simple problem with two variables and two constraints, first introduced in Section 2.

## Acknowledgments

Many thanks are due to an anonymous referee for close and careful readings of various versions of the paper and for many helpful suggestions.

##

Justification of the Estimates (236) and (237).

To prove (236), we use analysis similar to that of S. Wright . From the definition (23) of $`\mu `$, and the centrality condition (29c), we have that

$$\lambda _is_i=\mathrm{\Theta }(\mu ),\text{for all }i=1,2,\mathrm{},m\text{.}$$

Hence, from the third block row of (26) and the assumption (25) on the size of $`t`$, we have that

$$\frac{\mathrm{\Delta }\lambda _i}{\lambda _i}+\frac{\mathrm{\Delta }s_i}{s_i}=1\frac{t_i}{s_i\lambda _i}=1+O(\mu ),\text{for all }i=1,2,\mathrm{},m\text{.}$$
(256)
We have from Lemma 4 and (137) that $`\mathrm{\Delta }\lambda _i/\lambda _i=O(\mu )`$ for all $`i`$. Hence, by using (34a) from (4) together with (256), we obtain

$$\mathrm{\Delta }s_i=s_i+O(\mu ^2),\text{for all }i\text{.}$$
(257)
For the computed step components $`\widehat{\mathrm{\Delta }s}_{}`$, we have by combining (232a) with (257) that

$$\widehat{\mathrm{\Delta }s}_i=s_i+\delta _𝐮+O(\mu ^2),\text{for all }i\text{.}$$
(258)
Therefore, if $`s_i+\alpha \widehat{\mathrm{\Delta }s}_i=0`$ for some $`i`$ and some $`\alpha [0,1]`$, we have by using (34a) again that

$`s_i+\alpha (s_i+\delta _𝐮+O(\mu ^2))=0`$ (259)
$``$ $`(1\alpha )s_i=\delta _𝐮+O(\mu ^2)`$
$``$ $`(1\alpha )=\delta _𝐮/\mu +O(\mu ),\text{for any }i\text{.}`$
Meanwhile, for $`i𝒩`$, we have from Lemma 4, (137), and (232a) that

$$s_i+\alpha \widehat{\mathrm{\Delta }s}_i>0,\text{for all }i𝒩\text{ and all }\alpha [0,1]\text{,}$$
(260)
so the components $`\widehat{\mathrm{\Delta }s}_𝒩`$ do not place a limit on the step length bound $`\widehat{\alpha }_{\mathrm{max}}`$. For the components $`\widehat{\mathrm{\Delta }\lambda }_𝒩`$, we have by using Lemma 4, (137), (232c), and (256) that

$$\widehat{\mathrm{\Delta }\lambda }_i=\lambda _i+\mu \delta _𝐮+O(\mu ^2),\text{for all }i𝒩\text{.}$$

Therefore, if $`\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i=0`$ for some $`i𝒩`$ and some $`\alpha [0,1]`$, we have by arguing as in (259) that

$$1\alpha =\delta _𝐮+O(\mu ).$$
(261)
Finally, for $`i`$, we have from Lemma 4 that $`\lambda _i=\mathrm{\Theta }(1)`$, while from (137), (232a), and (232b), we have that

$$\mathrm{\Delta }\lambda _i=O(\mu ),\widehat{\mathrm{\Delta }\lambda }_i=O(\mu )+\delta _𝐮/\mu ,\text{for all }i\text{.}$$
(262)
Therefore, we have for $`\mu 𝐮`$ that

$$\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i>0,\text{for all }i\text{ and all }\alpha [0,1]\text{.}$$
(263)
By combining the observations (259), (260), (261), and (263), we conclude that there is a value $`\widehat{\alpha }_{\mathrm{max}}`$ satisfying

$$\widehat{\alpha }_{\mathrm{max}}[0,1],1\widehat{\alpha }_{\mathrm{max}}=\delta _𝐮/\mu +O(\mu )$$

such that

$$(\lambda ,s)+\alpha (\widehat{\mathrm{\Delta }\lambda },\widehat{\mathrm{\Delta }s})>0,\text{for all }\alpha [0,\widehat{\alpha }_{\mathrm{max}}]\text{,}$$

proving the claim (236). By making various simplifications to the analysis above, it is easy to show that (238) holds as well.

We now prove the claims (237) concerning the changes in the feasibility and duality measures along the computed step.

From (2), (29a), and the first block row of (26), we have

$`r_f(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })`$ (264)
$`=`$ $`_z(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })`$
$`=`$ $`_z(z,\lambda )+\alpha _{zz}(z,\lambda )\widehat{\mathrm{\Delta }z}+\alpha g(z)\widehat{\mathrm{\Delta }\lambda }+O(\alpha ^2\widehat{\mathrm{\Delta }z}^2)`$
$`=`$ $`(1\alpha )_z(z,\lambda )+\alpha _{zz}(z,\lambda )(\widehat{\mathrm{\Delta }z}\mathrm{\Delta }z)+\alpha g_{}(z)(\widehat{\mathrm{\Delta }\lambda }_{}\widehat{\mathrm{\Delta }\lambda }_{})`$
$`+\alpha g_𝒩(z)(\widehat{\mathrm{\Delta }\lambda }_𝒩\mathrm{\Delta }\lambda _𝒩)+O(\alpha ^2\widehat{\mathrm{\Delta }z}^2).`$
From (137) and (232a), we have $`\widehat{\mathrm{\Delta }z}=\delta _𝐮+O(\mu )`$, so for $`\mu 𝐮`$ and $`\alpha [0,1]`$, we have

$$\alpha ^2\widehat{\mathrm{\Delta }z}^2=O(\mu ^2).$$
(265)
From the definition (18) of the SVD of $`g_{}(z^{})`$, Theorem 5, and (232a), we have that

$`g_{}(z)(\widehat{\mathrm{\Delta }\lambda }_{}\mathrm{\Delta }\lambda _{})`$ $`=`$ $`g_{}(z^{})(\widehat{\mathrm{\Delta }\lambda }_{}\mathrm{\Delta }\lambda _{})+O(zz^{}\widehat{\mathrm{\Delta }\lambda }_{}\mathrm{\Delta }\lambda _{})`$ (266)
$`=`$ $`\widehat{U}\mathrm{\Sigma }U^T(\widehat{\mathrm{\Delta }\lambda }_{}\mathrm{\Delta }\lambda _{})+O(\mu )\delta _𝐮/\mu `$
$`=`$ $`\delta _𝐮.`$
Note that the larger error (232b) in the component $`V^T(\widehat{\mathrm{\Delta }\lambda }_{}\mathrm{\Delta }\lambda _{})`$, which is present when MFCQ is satisfied but not when LICQ is satisfied, does not enter into the estimate (266). By substituting this estimate into (264) together with estimates for $`\widehat{\mathrm{\Delta }z}\mathrm{\Delta }z`$ and $`\widehat{\mathrm{\Delta }\lambda }_𝒩\mathrm{\Delta }\lambda _𝒩`$ from (232), we obtain that

$$r_f(z+\alpha \widehat{\mathrm{\Delta }z},\lambda +\alpha \widehat{\mathrm{\Delta }\lambda })=(1\alpha )r_f+\delta _𝐮+O(\mu ^2),$$

verifying our claim (237c). The potentially large error (232b) does not affect rapid decrease of the $`r_f`$ component along the computed search direction.

For the second feasibility measure $`r_g`$, we have from (29b), the second block row of (26), and the estimates (232a) and (265) that

$`r_g(z+\alpha \widehat{\mathrm{\Delta }z},s+\alpha \widehat{\mathrm{\Delta }s})`$
$`=`$ $`g(z+\alpha \widehat{\mathrm{\Delta }z})+s+\alpha \widehat{\mathrm{\Delta }s}`$
$`=`$ $`g(z)+\alpha g(z)^T\widehat{\mathrm{\Delta }z}+s+\alpha \widehat{\mathrm{\Delta }s}+O(\alpha ^2\widehat{\mathrm{\Delta }z}^2)`$
$`=`$ $`(1\alpha )(g(z)+s)+\alpha g(z)^T(\widehat{\mathrm{\Delta }z}\mathrm{\Delta }z)+\alpha (\widehat{\mathrm{\Delta }s}\mathrm{\Delta }s)+O(\mu ^2)`$
$`=`$ $`(1\alpha )r_g+\delta _𝐮+O(\mu ^2),`$
verifying (237d).

To examine the change in $`\mu `$, we look at the change in each pairwise product $`\lambda _is_i`$, $`i=1,2,\mathrm{},m`$. We have

$`(\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i)(s_i+\alpha \widehat{\mathrm{\Delta }s}_i)`$
$`=`$ $`\lambda _is_i+\alpha (s_i\widehat{\mathrm{\Delta }\lambda }_i+\lambda _i\widehat{\mathrm{\Delta }s}_i)+\alpha ^2\widehat{\mathrm{\Delta }s}_i\widehat{\mathrm{\Delta }\lambda }_i`$
$`=`$ $`\lambda _is_i+\alpha (s_i\mathrm{\Delta }\lambda _i+\lambda _i\mathrm{\Delta }s_i)+\alpha s_i(\widehat{\mathrm{\Delta }\lambda }_i\mathrm{\Delta }\lambda _i)+\alpha \lambda _i(\widehat{\mathrm{\Delta }s}_i\mathrm{\Delta }s_i)`$
$`+\alpha ^2\widehat{\mathrm{\Delta }\lambda }_i\widehat{\mathrm{\Delta }s}_i.`$
From the last block row in (26), the estimate $`t=O(\mu ^2)`$ (25), and the estimate (137) of the exact step, we have

$$\lambda _is_i+\alpha (s_i\mathrm{\Delta }\lambda _i+\lambda _i\mathrm{\Delta }s_i)=(1\alpha )\lambda _is_i+O(\mu ^2).$$
(268)
From (137) and (232), we have

$$\widehat{\mathrm{\Delta }\lambda }_i\widehat{\mathrm{\Delta }s}_i=(\delta _𝐮/\mu +O(\mu ))(O(\mu )+\delta _𝐮)=\delta _𝐮+O(\mu ^2),$$
(269)
since $`\mu 𝐮`$. For $`i`$, we have from Lemma 4, (232a), and (232b) that

$$s_i(\widehat{\mathrm{\Delta }\lambda }_i\mathrm{\Delta }\lambda _i)=O(\mu )\delta _𝐮/\mu =\delta _𝐮,\text{for all }i.$$
(270)
For $`i𝒩`$, we have from Lemma 4 and (232c) that

$$s_i(\widehat{\mathrm{\Delta }\lambda }_i\mathrm{\Delta }\lambda _i)=\mu \delta _𝐮,\text{for all }i𝒩.$$
(271)
For the remaining term $`\lambda _i(\widehat{\mathrm{\Delta }s}_i\mathrm{\Delta }s_i)`$, we have from Lemma 4 and (232a) that

$$\lambda _i(\widehat{\mathrm{\Delta }s}_i\mathrm{\Delta }s_i)=\delta _𝐮,\text{for all }i=1,2,\mathrm{},m\text{.}$$
(272)
By substituting (268)–(272) into (Effects of finite-precision arithmetic on interior-point methods for nonlinear programming), we obtain

$$(\lambda _i+\alpha \widehat{\mathrm{\Delta }\lambda }_i)(s_i+\alpha \widehat{\mathrm{\Delta }s}_i)=(1\alpha )\lambda _is_i+\delta _𝐮+O(\mu ^2),\text{all }i=1,2,\mathrm{},m\text{.}$$
(273)
Therefore, by summing over $`i`$ and using (23), we obtain (237b).