Accurate Solution of Weighted Least Squares by Iterative Methods.
We consider the weighted least-squares (WLS) problem with a very ill-conditioned weight matrix. WLS problems arise in many applications including linear programming, electrical networks, boundary value problems, and structures.  Because of roundoff errors, standard iterative methods for solving a WLS problem with ill-conditioned weights may not give the correct answer.  Indeed, the difference between the true and computed solution (forward error) may be large.  We propose an iterative algorithm, called MINRES-L, for solving WLS problems.  The MINRES-L method is the application of MINRES, a Krylov-space method due to Paige and Saunders [SIAM J. Numer. Anal., 12 (1975), pp. 617--629], to a certain layered linear system.  Using a simplified model of the effects of roundoff error, we prove that MINRES-L ultimately yields answers with small forward error. We present computational experiments for some applications.
Introduction
Consider the weighted least-squares (WLS) problem
In this formula and for
the remainder of this article, k \Delta k indicates the 2-norm. We make the following
assumptions: D is a diagonal positive definite matrix and rank A = n. These
assumptions imply that (1) is a nonsingular linear system with a unique
solution. The normal equations for (1) have the form
A T
Weighted least-squares problems arise in several application domains including
linear programming, electrical power networks, elliptic boundary
value problems and structural analysis, as observed by Strang [21]. This
article focuses on the case when matrix D is severely ill-conditioned. This
happens in certain classes of electrical power networks. In this case, A is
a node-arc adjacency matrix, D is matrix of load conductivities, b is the
vector of voltage sources, and x is the vector of voltages of the nodes. Ill-conditioning
occurs when resistors are out of scale, for instance, when modeling
leakage of current through insulators.
Ill-conditioning also occurs in linear programming when an interior-point
method is used. To compute the Newton step for an interior-point method,
we need to solve a weighted least-squares equation of the form (2). Since
some of the slack variables become zero at the solution, matrix D always becomes
ill-conditioned as the iterations approach the boundary of the feasible
region. In Section 9, we cover this application in more detail. Ill-conditioning
also occurs in finite element methods for certain classes of boundary value
problems, for example, in the heat equilibrium equation r \Delta
thermal conductivity field c varies widely in scale.
An important property of problem (1) or (2) is the norm bound on the
solution, which was obtained independently by Stewart [20], Todd [22] and
several other authors. See [6] for a more complete bibliography. Here we
state this result as in the paper by Stewart.
Theorem 1 Let D denote the set of all positive definite m \Theta m real diagonal
matrices. Let A be an m \Theta n real matrix of rank n. Then there exist constants
-A and -
-A such that for any D 2 D
Note that the matrix appearing in (3) is the solution operator for the normal
equations (2), in other words, (2) can be rewritten as
Since the bounds (3), (4) exist, we can hope that there exist algorithms
for (2) that possess the same property, namely, the forward error bound does
not depend on D. We will call these algorithms stable, where stability, as
defined by Vavasis [23], means that forward error in the computed solution
x satisfies
where ffl is machine precision and f(A) is some function of A not depending
on D. Note that the underlying rationale for this kind of bound is that the
conditioning problems in (1) stem from an ill-conditioned D rather than an
ill-conditioned A.
This stability property is not possessed by standard direct methods such
as QR factorization, Cholesky factorization, symmetric indefinite factoriza-
tion, range-space and null-space methods, nor by standard iterative methods
such as conjugate gradient applied to (2). The only two algorithms in literature
that are proved to have this property are the NSH algorithm by Vavasis
[23] and the complete orthogonal decomposition (COD) algorithm by Hough
and Vavasis [12], both of them direct. See Bj-orck [1] for more information
about algorithms for least-squares problems.
We would like to have stable iterative methods for this problem because
iterative methods can be much more efficient than direct methods for large
sparse problems, which is the common setting in applications.
This article presents an iterative algorithm for WLS problems called
MINRES-L. MINRES-L consists of applying the MINRES algorithm of Paige
and Saunders [14] to a certain layered linear system. We prove that MINRES-
satisfies (5). This proof of the forward error bound for MINRES-L is based
on a simplified model of how roundoff error affects Krylov space methods.
This analysis is then confirmed with computational experiments in Section 8.
(The simplified model itself is described in Section 5.) An analysis of round-off
in MINRES-L starting from first principles is not presented here because
the effect of roundoff on the MINRES iteration is still not fully understood.
MINRES-L imposes the additional assumption on the WLS problem instance
that D is "layered." This assumption is made without loss of generality
(i.e., every weighted least-squares problem can be rewritten in layered
form), but the MINRES-L algorithm is inefficient for problems with many
layers.
This article is organized as follows. In Section 2 we state the layering
assumption, and also the layered least-squares (LLS) problem. In Section 3
we consider previous work. In Section 4 we describe the MINRES-L method
for two-layered WLS problems. In Section 5 we analyze the convergence in
the two-layered case using the simplifying assumptions about roundoff error.
In Section 6 and Section 7 we extend the algorithm and analysis to the case
of p layers. In Section 8 we present some computational experiments in
support of our claims. In Section 9 we consider application of MINRES-L to
interior-point methods for linear programming.
2 The Layering Assumption
Recall that we have already assumed that the weight matrix D appearing
in (1) is diagonal, positive definite and ill-conditioned. For the rest of this
article we impose an additional "layering" assumption: we assume, after
a suitable permutation of the rows of (A; b) and corresponding symmetric
permutation of D, that D has the structure
where each D k is well-conditioned and scaled so that its smallest diagonal
entry is 1, and where denote the maximum
diagonal entry among D . The layering assumption is that - is not
much larger than 1.
Note that this assumption is made without any loss of generality (and we
could assume since we could place each diagonal entry of D in its own
layer. Unfortunately, the complexity of our algorithm grows quadratically
with p. Furthermore, our upper bound on the forward error degrades as p
increases (see (39) below). Thus, a tacit assumption is that the number of
layers p is not too large.
From now on, we write A in partitioned form as
A pC C A
to correspond with the partitioning of D. We partition
similarly.
Under this assumption, we say that (1) is a "layered WLS" problem. In
the context of electrical networks, this assumption means that there are several
distinct classes of wires in the circuit, where the resistance of wires in
class l is of order 1=ffi l . For instance, one class of wires might be transmission
lines, whereas the other class might consist of broken wires (open lines)
where the resistance is much higher. In the context of the heat equilibrium
equation, the layering assumption means that the object under consideration
is composed of a small number of different materials. Within each material
the conductivity ffi l is constant, but the different materials have very different
conductivities. In linear programming, taking means that the some
of the slack variables at the current interior-point iterate are "small" while
others are "large."
A limiting case of layered WLS occurs when the gaps between the ffi l 's
tend to infinity, that is, ffi 1 is infinitely larger than ffi 2 and so on. As the
weight gaps tend to infinity, the solution to (1) tends to the solution of the
following problem, which we refer to as layered least squares (LLS). Construct
a sequence of nested affine subspaces L 0 oe L 1 oe \Delta \Delta \Delta oe L p of R n . These spaces
are defined recursively: L
fminimizers of kD 1=2
l
Finally, x, the solution to the LLS problem, is the unique element in L p . The
layered least-squares problem was first introduced by Vavasis and Ye [25] as
a technique for accelerating the convergence of interior-point methods. They
also established the result mentioned above in this paragraph: the solution
to the WLS problem in the limit as ffi l+1 =ffi l ! 0 for all l converges to the
solution of the LLS problem.
Combining this result with Theorem 1 yields the following corollary, also
proved by Vavasis and Ye.
Corollary 1 Let x be the solution to the LLS problem posed with matrix A
and right-hand side vector b. Then kxk -Akbk and kAxk -
-Akbk for
any choice of diagonal positive definite weight matrices D
3 Previous Work
The standard iterative method for least-squares problems, including WLS
problems, is conjugate gradient (see Golub and Van Loan [7] or Saad [18])
applied to the normal equations (2). This algorithm is commonly referred to
as CGNR, which is how we will denote it here. There are several variants of
CGNR in the literature; see, e.g., Bj-orck, Elfving, and Strako-s [2]. Note that
in most variants one does not form the triple product A T DA when applying
CG to (2); instead, one forms matrix-vector products involving matrices A T ,
D and A. This trick can result in a substantial savings in the running time
since A T DA could be much denser than A alone. The same trick is applicable
to our MINRES-L method and was used in our computational experiments.
The difficulty with CGNR is that an inaccurate solution can be returned
because A T DA can be ill-conditioned when D is ill-conditioned. To understand
the difficulty, consider the two-layered WLS problem, which is obtained
by subtituting (6) in the case
Observe that if sequence
A T Db;
constructed by CGNR is very close to
In other words, information about A 2 , D 2 and b 2 is lost when forming the
Krylov sequence. A different framework for interpreting this difficulty is
described in Section 5.
Another iterative method for least-squares problems is LSQR due to Paige
and Saunders [15]. This method shares the same difficulty with CGNR because
it works in the same Krylov space.
A standard technique for handling ill-conditioning in conjugate gradient
is reorthogonalization; see, for example, Paige [16] and Parlett and
Scott [17]. Reorthogonalization, however, cannot solve the difficulty with
ill-conditioning in (2) because even the act of forming the first Krylov vector
A T Db causes a loss of information.
Another technique for addressing ill-conditioned linear systems with iterative
methods is called "regularization"; a typical regularization technique
modifies the ill-conditioned system with additional terms. See Hanke [10].
Regularization does not appear to be a good approach for solving (1) because
(1) already has a well-defined solution (in particular, Theorem 1 implies that
solutions are not highly sensitive to perturbation of the data vector b). A
regularization technique would compute a completely different solution.
In our own previous work [3], we proposed an iterative method for (2)
based on "correcting" the standard CGNR search directions. We have since
dropped that approach because we found a case that seemingly could not be
handled or detected by that algorithm.
4 MINRES-L for Two Layers
In this section and the next we consider the two-layered case, that is,
in (6). We consider the two-layered case separately from the p-layered case
because the two-layered case contains all the main ideas of the general case
but is easier to write down and analyze. (In the our algorithm
reduces to MINRES applied to (2) and hence is not novel.) Furthermore, the
case is expected to occur commonly in practice. We mention also that
the two-layered WLS and LLS problems were considered in x22 of Lawson
and Hanson [13].
As noted in the preceding section, the two-layered WLS problem is written
in the form (7), in which the diagonal entries of D 1 ; D 2 on the order of 1 and
us introduce a new variable v such that
A T
Note that this equation always has a solution v because the right-hand side
is in the range of A T
1 . Multiplying (8) by ffi 2 and adding to (7) yields
A T
Putting (8) and (9) together, we get
A T
!/
x
A T
Our algorithm, which we call MINRES-L (for MINRES "layered"), is the
application of the MINRES iteration due to Paige and Saunders [14] to (10).
Note that (10) is a symmetric linear system.
In general, this linear system is rank deficient because if (x; v) is a solution
solution. Thus, (10) is
rank deficient whenever the rank of A 1 is less than n. This means we must
address existence and uniqueness of a solution. Existence follows because
the original WLS problem (7) is guaranteed to have a solution. Uniqueness
of x is established as follows: if we add times the first row of (10) to ffi 1
times the second row, we recover the original WLS problem (7). Since (7)
has a unique solution, (10) must uniquely determine x. Since x is uniquely
determined, so is A 1 v.
The question arises whether MINRES (in exact arithmetic) will find a solution
of (10). MINRES can find a solution only if it lies in the Krylov space,
which (because of rank deficiency) is not necessarily full dimensional. This
question was answered affirmatively by Theorem 2.4 of Brown and Walker [4].
(Their analysis concerns GMRES, but the same result applies to MINRES in
exact arithmetic.) Furthermore, their result states that, assuming the initial
guess is 0, the computed solution (x; v) will have minimum norm over all
possible solutions. Since x is uniquely determined, their result implies that
will have minimum norm.
Recall from Section 3 that the problem with applying conjugate gradient
directly to (7) is that the linear system may be ill-conditioned when
and hence conjugate gradient may return an inaccurate answer. Thus, it may
seem paradoxical that we remedy a problem caused by ill-conditioning with
an iterative method based on a truly rank-deficient system. One explanation
of this paradox concerns the limiting behavior as 1. In this case,
(7) tends to the linear system A T
This system will, in
general, not have a unique solution (because A 1 is not assumed to have rank
n), so CGNR will compute some solution that may have nothing to do with
. Thus, the CGNR solution is not expected to have the forward
accuracy that we demand.
On the other hand, as we see that (10) tends to
A T
!/
x
A T
This system is easily seen to be the Lagrange multiplier conditions for the
two-layered LLS problem: recall from Section 2 that the two-layered LLS
problem is
minimize kD 1=2
subject to A T
This is the correct limiting behavior: the WLS solution tends to the LLS
solution as explanation of MINRES-L's convergence
behavior follows.
Convergence Analysis for Two Layers
In this section we consider convergence of MINRES-L in the presence of
roundoff error for the case 2. As mentioned in the introduction, we make
a simplifying assumption concerning the effect of roundoff error in Krylov
space methods. The assumption concerns either CG or MINRES applied to
the symmetric linear system In our use of these algorithms, there
is no preconditioner, and the initial guess is x Further, in our use of
MINRES, c lies in the range-space of M (i.e., the system is consistent). In
our use of CG, M is positive definite. With these restrictions in mind, our
assumption about the effect of roundoff is that after a sufficient number of
iterations, either method will compute an iterate -
x satisfying
where C is a modest constant, ffl is machine epsilon, and x is the true solution.
(If multiple solutions exist, we take x to be the minimum-norm solution.)
As far as we know, this bound has not been rigorously proved, but it is
related to a bound proved by Greenbaum [9] in the case of conjugate gradient.
In particular, Greenbaum's result implies that (11) would hold for CG if we
were guaranteed that the recursively updated residual drops to well below
machine precision, which always happens in our test cases.
As for MINRES, less is known, but a bound like (11) is known to hold
for GMRES implemented with Householder transformations [5]. GMRES is
equivalent to MINRES augmented with a full reorthogonalization process.
We are content to assert (11) for MINRES, with evidence coming from our
computational experiments.
This bound sheds light on why MINRES-L can attain much better accuracy
than CGNR. For CGNR, the error bound (11) implies that kA T
A T DA- xk gets very small, where - x is the computed solution. This latter
quantity is the same as
x)k. But recall that we are seeking
a bound on the forward error, that is, on
xk. In this case, the factor
greatly skew the norm when is close to zero, so there is
no bound on kx \Gamma -
xk independent of ffi 1 =ffi 2 , that is, (5) is not expected to be
satisfied by CGNR. This is confirmed by our computational experiments.
In contrast, an analysis of MINRES-L starting from (11) does yield the
accuracy bound (5). We need the following preliminary lemma.
A be an m \Theta n matrix of rank n and -
A an r \Theta n submatrix
of A. Suppose the linear system -
consistent. Here, c is a
given vector, and -
D is a given diagonal positive definite matrix. Then for
any solution x,
-A \Delta kck (12)
and
Furthermore, there exists a solution x satisfying
-A
Proof. First, note the following preliminary result. Let H;K be two symmetric
n \Theta n matrices such that H is positive semidefinite and K is positive
definite. Let b be an n-vector in the range space of H. Then (H
converges to a solution of . This is proved by reducing to
the diagonal case using simultaneous diagonalization of H;K.
Let D be the extension of -
D to an m \Theta m diagonal matrix obtained by
filling in zeros, so that A T
A. Since A T
the limit of solution x of -
as noted in the preceding paragraph. Let M be an m \Theta m diagonal
matrix with 1's in diagonal positions corresponding to -
D and zeros elsewhere.
We have
ck
ck (15)
ffl?0
-A \Delta kck:
The last line was obtained by the transpose of (4). This proves (12). Note
that this holds for all x satisfying -
c, since this latter equation
uniquely determines -
Ax. Similarly, to demonstrate (13), we start from (15):
ck
ffl?0
ck
For the second part of the proof, observe by the first part that A T
and thus
ffl?0
Axk:
Combining this with (12) proves (14).
To resume the analysis of MINRES-L, we define
where (- x; -
v) is the solution computed by MINRES-L. Then (11) applied to
yields the bounds
In this formula, H 2 is shorthand for the coefficient matrix of (10).
We can extract another equation from (16) and (17); in particular, if we
multiply (16) by multiply (17) by ffi 1 and then add, we eliminate the terms
involving - v:
Let x be the exact solution to the WLS problem. The last two terms of this
equation can be replaced with terms involving x by using (7). Interchanging
the left- and right-hand sides yields
The goal is to derive an accuracy bound like (5) from (18) and (19). We
start by bounding the quantity on the right-hand side of (18). Note that
can be bounded by because the largest entries in D 1 ; D 2 are
bounded by -. We can bound kxk by -Akbk using Theorem 1. Next we turn
to bounding kvk in (18). Recall that, as mentioned in the preceding section,
v is not uniquely determined, but MINRES will find the minimum-norm v
satisfying (10). Recall that v is determined by the constraint
One way to pick such a v is to make it minimize kA 2 vk subject to the above
constraint. In this case, v is a layered least-squares solution with right-hand
side data (b yields the bound
for this choice of v. (The factor -
can be improved to -
-A by using
the analysis of Gonzaga and Lara [8].) Combining the x and v contributions
means that we have bounded the right-hand side of (18); let us rewrite (18)
with the new bound:
Next, we write new equations for r . Observe that r 1 lies in the range
of A Tand A T, so we can find h 1 satisfying
Similarly, by (17) there exists h 2 satisfying
By applying (13) to r 1 and r 2 separately, with "A T c" in the lemma taken to
be first r 1 and then r 2 , we conclude from (21) and (22) that
Substituting (21) and (22) into (19) yields
Notice (by analogy with (7)) that the preceding equation is exactly a weighted
least-squares computation where the "unknown" is -
and the right-hand
side data is Thus, by Theorem 1,
We now build a chain of inequalities: the right-hand side of the preceding
inequality is bounded by (23) and (24), and the right-hand side of (23) and
(24) is bounded by (20). Combining all of this yields
To obtain the preceding inequality, we used the facts that
assumption) and that kdiag(D \Gamma1
by assumption, since the
smallest entry in each D i is taken to be 1).
Thus, we have an error bound of the form (5) as desired; in particular,
there is no dependence of the error bound on ffi 2 =ffi 1 . Note that this bound depends
on -. Recall that - is defined to be the maximum entry in D
and is assumed to be small. Indeed, as noted in Section 2, we can always
assume that if we are willing to divide the problem into many layers.
6 MINRES-L for p Layers
In this section we present the MINRES-L algorithm for the p-layered WLS
problem. The algorithm is the application of MINRES to the symmetric
linear system H p is a square matrix of size (1
1)=2)n \Theta (1 +p(p \Gamma 1)=2)n, c p is a vector of that order, and w is the vector of
unknowns. Matrix H p is partitioned into (1
blocks each of size n \Theta n. Vectors c p and w are similarly partitioned. The
WLS solution vector is the first subvector of w.
In more detail, the vector w is composed of x concatenated with p(p\Gamma1)=2
n-vectors that we denote v i;j , where i lies in lies in
Recall that the p-layered WLS problem may be written
Let x be the solution to this equation. Then we see from this equation that
A T
lies in the span of [A T
Therefore, there exists
a solution [v to the equation
A T
This equation is the first block-row of H In other words, the first
block row of H p contains one copy of each of the matrices A T
and the
first block of c p is A T
In general, the (p 1)th block-row of H
the equation
A T
A T
A T
This completes the description of block-rows . We now
establish some properties of these block-rows, and we postpone the description
of block-rows
Lemma Suppose w is a solution to the linear equation (28) for each
denotes the concatenation of x and all of the v i;j 's. Then
x is the solution to the WLS problem (26).
Proof. For each i, multiply (28) by ffi i and then sum all p equations obtained
in this manner. Observe that all the v i;j terms cancel out and we end up
exactly with (26).
We also need the converse to be true.
Lemma 3 Suppose x is the solution to (26). Then there exist vectors v i;j
such that (28) is satisfied for each
Proof. The proof is by induction on (decreasing) We assume
that we have already determined v i;j for all
that (28) is satisfied for now we must
determine v k;j for for the particular value
k. The base case of the induction is that we can select v
to satisfy (28) in the case
lies in the range of [A T
because of (26).
Now for the induction case of k ! p. Rewrite (28) for the case
multiply through by
A T
Recall that our goal is to choose v k;j for to make this equation
valid.
Multiply (28) for each and add this to (29). After
rearranging and summations and cancelling common terms on the left-hand
side, we end up with
Dividing through by ffi k and separating out the v k;j terms from the second
summation yields:
A T
A T
A T
But from (26) we know that
lies in the range of
the rightmost summation of (31) also lies in the same
range. Therefore, there exist v k;j for
But then these same choices will make (29) valid because the algebraic steps
used to derive (31) from (29) can be reversed. This proves the lemma.
Note that the preceding proof actually demonstrates a strengthened version
of the lemma. The strengthened version states that if we are given x
satisfying (26) and, for some k, vectors v i;j for k - that satisfy
(28) for all then we can extend the given data to a solution of
(28) for all strengthened version is needed below.
We now explain the remaining p(p \Gamma 1)=2 block-rows of H p . These rows
exist solely for the purpose of making H p symmetric. First, we have to
order the variables and equations correctly. The variables will be listed in
the order (x; v The first
equations will be listed in the order (28) for 1. This means
that the first p rows of H p have the format [S
and T p is a p \Theta (p \Gamma 1)(p \Gamma 2)=2 matrix. Furthermore, it is easily checked
that S p is symmetric: its first block-row and first block-column both consist
of A T
listed in the order 1)st entry of its main
diagonal is \Gamma(ffi p =ffi i )A T
all its other blocks are
zeros. Then we define H p to be
We define c p as
A T
A T
where there are p(p \Gamma 1)=2 blocks of zeros. For example, the following linear
system is H 3
A T
A T
A T
A T
x
A T
A T
We now must consider whether has any solutions; in particular,
we must demonstrate that the new group of equations T T
with the first p rows. Here w 0 denotes the first p blocks of w, that is,
Studying the structure of T p , we see that there are
indexed by (i;
correspondence with the columns of T p , which correspond to variables v i;j for
in that range). The row indexed by (i; j) has exactly two nonzero block
entries that yield the equation
A T
A T
Our task is therefore to show that we can simultaneously satisfy (28) for
Our approach is to select the v p;j 's in the order v In
particular, assuming v are already selected, we define v p;j to
be any solution to
A T
The following lemma shows that this linear system is consistent.
Lemma 4 If the v p;j 's are chosen in reverse order to satisfy (33), then at
each step the linear system is consistent, and (32) is satisfied.
Proof. The proof is by reverse induction on j. The base case is
which case (33) has a solution because, as noted above, A T
lies in the span of [A T
In the case vacuously
true: there is no i in the specified range.
Now consider the case any i in the range
Start with the version of (33) satisfied by v p;i , which holds by the induction
hypothesis:
A T
Move the terms of the first summation to the right-hand
side:
A T
A T
A T
A T
A T
The second line was obtained from the first by applying (32) inductively
(with "j" in (32) taken to be k). The third line was obtained by merging the
two summations on the right.
But notice that the preceding equation means that v p;i satisfies the same
linear system as v p;j , that is (33), except with the right-hand side scaled by
This proves that (33) is consistent for the j case since we have constructed
a solution to it. Although this linear system does not necessarily
have a unique solution, a linear system of the form A T uniquely determines
Ax. Thus, we have also proved that
for all This result is actually a strengthening of (32) for j; for
that equation we need only the specific case of
The reader may have noticed that the preceding proof is apparently too
complicated and that we could establish the result more simply by solving
for v p;p\Gamma1 in (33) with setting v
2. This simpler approach does not yield the bounds on kv p;j k
needed in the next section.
This proof shows that the above method for selecting v
consistent and satisfies (32). We also see that (27) is satisfied; this follows
immediately from taking (33). To complete the proof that there is
a solution to H p need only verify (28) in the case
But recall from the proof of Lemma 3 that the remaining v i;j 's for
can be determined sequentially by using the construction in the
proof. Thus, the arguments of this section have established the following
theorem.
Theorem 2 There exists at least one solution w to H p
more, any such solution has as its first n entries the vector x that solves
(26).
7 Convergence Analysis for p Layers
The convergence analysis for p layers follows the same basic outline as the
convergence analysis for two layers. In particular, we use (11) as the starting
point for the error analysis. Observe that (11) has the norm of the true
solution on the right-hand side. Thus, to apply that bound, we must get a
norm bound on v i;j for all
We start with bounds on v p;j for Apply Lemma
1 to (33) in the case In the lemma, take -
As noted above, A T
lies in
the range of [A T
so (33) is consistent. The right-hand side of (33)
in the case has the form A T c with
Note that kD p (b 1)kbk. Thus, from (12),
-A
To derive the third line from the second, we used the facts that kD \Gamma1
for each i and
Now we use the same line of reasoning to get a bound on v p;p\Gamma2 based on
(33) for the case 2. In this case, the right-hand side of (33) has the
Thus, kck is bounded by ffi p\Gamma2 (- A
which is at
most
We continue this argument inductively. Each time the bound grows by a
factor 2- A to take into account the fact that v p;i appears on the right-hand
side for the equation determining v p;i\Gamma1 . In the end we conclude that
Next we must bound v i;j for 1 These vectors are
determined by (28). We can find a solution to (28) by first solving
for z i , where -
is already known to be consistent. Furthermore, in the preceding
equation. We set v Using (12), we conclude that
for each
We now claim that
p. This is proved by induction on decreasing i using
recurrence (36). The on the right-hand side of (36) is bounded
by (35), and the remaining terms are bounded by the induction hypothesis.
We omit the details.
For the right-hand side of (11) we need a bound on kv i;j k. Note that up
to now we have not uniquely determined v i;j itself. Recall that in each case
Lemma 1 was used to bound kA k v i;j k. We can force unique determination
by choosing the v i;j as in the proof of Lemma 1, yielding
by (14). Note that MINRES does not necessarily select this v i;j , but because
of its minimization property (that is, Theorem 2.4 of Brown and Walker [4]
described in Section 4), it will select v i;j whose norm is no larger than in the
preceding bound.
We now can apply (11). The other factor on the right-hand side, namely,
easily seen to be bounded by p
w be the solution
computed by MINRES-L, and let
substituting (37) on the right-hand side of (11) yields
Let r be the first p block-entries of r. Note that r j must lie in the
span of [A
in order for the equation H
to have a
solution, because it can be seen from (28) that the (p 1)st block-row of
us find h i that solves r
for each i. By (13) we know that k[A
Let - x be the first n entries of -
that is, the computed WLS solution. If
we multiply the (p \Gamma i+1)st block row of H
and add these p rows, we obtain
A T
The third line was obtained from the second by interchanging the order of
summation. Thus, we see from the third line above that -
solves a WLS
problem in which the ith entry of the data vector is A i
in this range, we conclude that the data vector is bounded
in norm by k. Then Theorem 1
implies that
A
Substituting (38) yields
A \Delta (4- A )
This is a bound of the form (5) as desired.
Computational Experiments
In this section we present computational experiments on MINRES-L and
CGNR to compare their accuracy and efficiency. The first few tests involve
a small node-arc adjacency matrix. The remaining tests are on matrices
arising in linear programming and boundary value problems. All tests were
conducted in Matlab 4.2 running on an Intel Pentium under Microsoft Windows
NT 4.0. Matlab is a software package and programming language for
numerical computation written by The Mathworks, Inc. All computations
are in IEEE double precision with machine epsilon approximately 2:2
Matlab sparse matrix operations were used in all tests.
Our implementation of CGNR is based on CGLS1 as in (3.2) of Bj-orck,
Elfving and Strako-s [2]. These authors conclude that CGLS1 is a good way
to organize CGNR. There are two matrix-vector products per CGLS1 it-
eration, one with matrix A T D 1=2 and one with D 1=2 A. In our implemen-
tation, the CGNR iteration terminates when the scaled computed residual
ks k k=kA T Dbk drops below 10 \Gamma13 . Our implementation of MINRES is
based on [14], except Givens rotations were used instead of 2 \Theta 2 Householder
matrices (so that there are some inconsequential sign differences).
The MINRES-L iteration terminates when the scaled computed residual
The first matrix A used in the following tests is the reduced node-arc
adjacency matrix of the graph depicted in Figure 1. A "node-arc adjacency"
matrix contains one column for each node of a graph and one row for each
edge. Each row contains exactly two nonzero entries, a +1 and a \Gamma1 in the
columns corresponding to the endpoints of the edge. (The choice of which
endpoint is assigned +1 and which is assigned \Gamma1 induces an orientation
on the edge, but often this orientation is irrelevant for the application.) A
reduced node-arc incidence (RNAI) matrix is obtained from a node-arc incidence
matrix by deleting one column. RNAI matrices arise in the analysis of
an electrical network with batteries and resistors; see [23]. They also arise in
network flow problems. In the case of Figure 1, the column corresponding to

Figure

1: An based on this graph was used for the first
group of tests. The column corresponding to the top node is deleted. Edges
marked with heavy lines are weighted 1, and edges marked with light lines
are weighted varies from test to test.
the top node was deleted. Thus, A is an 9 matrix. It is well known that
the RNAI matrix for a connected graph always has full rank. RNAI matrices
are known to have small values of -A and -
-A [23].
In all these tests, the weight matrix has two layers. We took
vary from experiment to experiment.
The rows of A in correspondence with D 2 are drawn as thinner lines in Figure
1. Finally, the right-hand side b was chosen to be the first prime numbers.
The results are displayed in Table 1, and the cases when
are plotted in Figure 2. The scaled error that is tabulated and
plotted in all cases is defined to be k- x \Gamma xk=kbk. We choose this particular
scaling for the error because our goal is to investigate stability bound
(5). The true solution x is computed using the COD method [12]. Note
that the accuracy of CGNR decays as ffi 2 gets smaller, whereas MINRES-L's
accuracy stays constant. MINRES-L requires many more flops than CGNR
because the system matrix is larger. The running time of CGNR is about
the same for the first four rows of the table as the ill-conditioning increases.
In the last two rows the running time of CGNR drops because the matrix
A T DA masquerades as a low-rank matrix for small values of ffi 2 , causing early
termination of the Lanczos process.
Besides returning an inaccurate solution, CGNR has the additional difficulty
that its residual (the quantity normally measured in practical use of
this algorithm) does not reflect the forward error, so there is no simple way

Table

1: Behavior of the two-layered MINRES-L algorithm compared to
CGNR for decreasing values of ffi 2 . The error reported is the scaled error
defined in the text. Note that the CG accuracy degrades while the MINRES-
accuracy stays about the same.
MINRES-L MINRES-L MINRES-L CGNR CGNR CGNR
Iterations Error Flops Iterations Error
to determine whether CGNR is computing good answers. In contrast, the
error and residual in MINRES-L are closely correlated. This correlation is
predicted by our theory.
The next computational test involved a larger matrix A taken from the
Netlib linear programming test set, namely, the matrix in problem AFIRO,
which is 51 \Theta 27. We used a matrix D with 1's in its first 27 diagonal positions
its remaining 24 positions (i.e., D
The right-hand side vector b was chosen to contain the first
primes. MINRES-L required 137 iterations and 250 kflops and yielded
a solution -
x with scaled error 3:0 with respect to the true solution
computed by the COD method. For this matrix, -A and -
-A are not known.
CGNR on this problem required 69 iterations and 61 kflops and returned an
answer with scaled error 2:2 . The convergence plots are depicted in

Figure

3.
The excessive number of iterations required by MINRES is apparently
caused by a loss of orthogonality in the Lanczos process. To verify this
hypothesis, we ran GMRES on the same layered matrix. GMRES [19] on
a symmetric matrix is equivalent to MINRES with full reorthogonalization.
(In exact arithmetic the two algorithms are identical.) We call this algorithm
GMRES-L. The same termination tests were used. The result is depicted in

Figure

4. In this case, GMRES-L ran for 50 iterations (fewer than (1
returned a more accurate answer, one with forward error
. However, the number of flops was higher, 350 k, because of the
scaled error
scaled residual
scaled error
scaled residual

Figure

2: Convergence behavior of CGNR and MINRES-L for the
RNAI test case. The plots are for In
these plots and all that follow, the x-axis is the iteration number. For both
algorithms the computed (i.e., recursively updated) residual is plotted rather
than the true residual. Other experiments (not reported here) indicate that
these are usually indistinguishable. The \Theta on the y-axis indicates the cutoff
below which the CGNR scaled residual must drop in order for (11) to be true
. The ffi on the y-axis is the analog for MINRES-L.

Figure

3: Convergence behavior of CGNR and MINRES-L for AFIRO. The
curves are labeled as in Figure 2.
Gram-Schmidt process in the GMRES main loop.
The next computational test involves a larger matrix A arising from finite-element
analysis. The application is the solution of the boundary value
problem r \Delta on the polygonal domain depicted in Figure 5 with
Dirichlet boundary conditions. The conductivity field c is 1 on the outer part
of the domain and is 10 12 on the darker triangles. As discussed in [24], this
type of problem gives rise to a weighted least-squares problem in which A
encodes information about the geometry and D encodes the ill-conditioned
conductivity field. The values of -A and -
-A for this matrix are not known,
although bounds are known for variants of these parameters. The particular
matrix A is 652 \Theta 136. The right-hand side vector b was chosen according to
the Dirichlet boundary conditions described in [24]. The MINRES-L method
for this problem gave scaled error of 1:3 iterations and 6.5
mflops. To compute the true solution, we used the NSHI method in [24].
In this case, surprisingly, CGNR gave almost as accurate an answer, but the
termination test was never activated. (We cut off CGNR after 10n iterations.)
The residual of CGNR is quite oscillatory as depicted in Figure 6. In the
finite-element literature, CGNR would be referred to as conjugate gradient
on the assembled stiffness matrix, which is A T DA.
A cause of this odd behavior of CGNR is as follows. Note that the region
of high conductivity is not incident on the boundary of the domain so

Figure

4: Convergence behavior of GMRES-L (-  and \Delta \Delta
Figure

5: Domain and finite element mesh used for the finite element exper-
iment. Conductivity in the dark triangles is 10 12 and in the light triangles is

Figure

Convergence of CGNR and MINRES-L for the finite element test
problem. The curves are labeled as in Figure 2.
Thus, A T
starts from a right-hand side that is already almost zero. Furthermore, this
right-hand side is nearly orthogonal to the span of A T
dominates
the stiffness matrix A T DA. Thus, CGNR has trouble making progress. The
surprisingly accurate answer from CGNR in this example is not so useful
in practice because there is no apparent way to detect that convergence is
underway.
The final test is a three-layered problem based on the matrix A from
ADLITTLE of the Netlib test set, a 138 \Theta 56 matrix. Matrix D has as its
first 28 diagonal entries 1, its next 28 diagonal entries 10 \Gamma8 and its last 82
entries . The right-hand side vector is the first 138 prime numbers.
The convergence is depicted in Figure 7. As expected, the scaled error of
MINRES-L decreased to while the scaled error of CGNR was 0:3.
Note the excessive number of iterations required by MINRES-L. Again, this
is apparently due to loss of orthogonality because the number of iterations
was only 118 for GMRES-L to achieve a scaled error of 9:4 In fact,
for this test GMRES-L was more efficient than MINRES-L in terms of flop
count.
In most cases we see that the MINRES-L algorithm performs essentially
as expected, except for the two cases in which a loss of orthogonality causes
many more iterations than expected. In every case, MINRES-L's running
time is higher than CGNR's, but CGNR can produce bad solutions as measured
by forward error.

Figure

7: Convergence of CGNR and MINRES-L for ADLITTLE. The curves
are labeled as in Figure 2. Note the excessive number of iterations for
MINRES-L caused by a loss of orthogonality.
9 An Issue for Interior-Point Methods
In this section we describe an issue that arises when using the MINRES-L
algorithm in an interior-point method for linear programming. Full consideration
of this matter is postponed to future work.
It is well known that the system of equations for the Newton step in an
interior-point method can be expressed as a weighted least-squares problem.
To be precise, consider the linear programming problem
subject to A T
whose dual is
subject to Ay
(which is standard form, except we have transposed A to be consistent with
least-squares notation). A primal-dual method starting at a feasible interior
point problem computes an update \Deltay to y satisfying
is an algorithm-dependent
parameter usually in [0; 1], - is the duality gap, and e is the vector of all 1's.
See Wright [26]. Since (40) has the form of a WLS problem, we can obtain
\Deltay using the MINRES-L algorithm.
One way to compute \Deltas is via \Deltas := \GammaA\Deltay. This method is not stable
because \Deltas has very small entries in positions where s has very small en-
these small entries must be computed accurately with respect to the
corresponding entry of s. In contrast, the error in all components of \Deltas
arising from the product A\Deltay is on the order of ffl \Delta ksk (where ffl is machine-
epsilon). A direct method for accurately computing all components of \Deltas
was proposed by Hough [11], who obtains a bound of the form
\Deltas
for each i. We will consider methods for extending MINRES-L to accurate
computation of \Deltas in future work. As noted by Hough, \Deltax is easily computed
from \Deltas with a similar accuracy bound assuming \Deltas satisfies (41).
Conclusions
We have presented an iterative algorithm MINRES-L for solving weighted
least squares. Theory and computational experiments indicate that the
method is more accurate than CGNR when the weight matrix is highly ill-
conditioned. This work raises a number of questions.
1. Is there an iterative method that does not require the layering assumption

2. If layering is indeed required, can we get a more parsimonious layered
linear system when p - 3? In particular, is there a 3n \Theta 3n system of
equations with all the desired properties for the 3-layered case (instead
of the 4n \Theta 4n system that we presented)?
3. What is the best way to handle loss of orthogonality in MINRES that
was observed in Section 8?
4. Can this work be extended to stable computation of \Deltax and \Deltas in an
interior-point method? (This question was raised in Section 9.)
5. What about preconditioning? In most of our computational tests, we
ran both MINRES and CG for more than n iterations because our aim
was to compute the solution vector as accurately as possible. In prac-
tice, one hopes for convergence in much fewer than n iterations. What
are techniques for preconditioning WLS problems? Note that the analysis
of MINRES-L's accuracy in Section 5 and Section 7 presupposes
that no preconditioner is used.

Acknowledgments

We had helpful discussions of this work with Anne Greenbaum and Mike
Overton of NYU; Roland Freund, David Gay, and Margaret Wright of Bell
Labs; Patty Hough of Sandia; Rich Lehoucq and Steve Wright of Argonne;
Homer Walker of Utah State; and Zden-ek Strako-s of the Czech Academy
of Sciences. We thank Patty Hough and Gail Pieper for carefully reading
an earlier draft of this paper. In addition, we received the Netlib linear
programming test cases in Matlab format from Patty Hough.



--R

Numerical methods for least squares problems.
Stability of conjugate gradient and Lanczos methods for linear least squares problems.
Iterative methods for weighted least squares.
GMRES on (nearly) singular systems.
Numerical stability of GMRES.
On linear least-squares problems with diagonally dominant weight matrices
Matrix Computations
A note on properties of condition numbers.
Estimating the attainable accuracy of recursively computed residual methods.
Conjugate gradient type methods for ill-posed problems
Stable computation of search directions for near-degenerate linear programming problems
Complete orthogonal decomposition for weighted least squares.
Solving Least Squares Problems.
Solution of sparse indefinite systems of linear equations.
LSQR: An algorithm for sparse linear equations and sparse least squares.
Practical use of the symmetric Lanczos process with re- orthogonalization
The Lanczos algorithm with selective reorthog- onalization
Iterative methods for sparse linear systems.
GMRES: A generalized minimum residual algorithm for solving nonsymmetric linear systems.
On scaled projections and pseudoinverses.
A framework for equilibrium equations.
A Dantzig-Wolfe-like variant of Karmarkar's interior-point linear programming algorithm
Stable numerical algorithms for equilibrium systems.
Stable finite elements for problems with wild coefficients.
A primal-dual interior point method whose running time depends only on the constraint matrix

--TR
krylov-space;iterative method;conjugate gradient;MINRES;weighted least squares;achievable accuracy