Complete Orthogonal Decomposition for Weighted Least Squares.
This paper proposes a complete orthogonal decomposition (COD) algorithm for solving weighted least-squares problems. In applications, the weight matrix can be highly ill conditioned, and this can cause standard methods like QR factorization to return inaccurate answers in floating-point arithmetic. Stewart and Todd independently established a norm bound for the weighted least-squares problem that is independent of the weight matrix. Vavasis proposed a definition of a "stable" solution of weighted least squares based on this norm bound: The solution computed by a stable algorithm must satisfy an accuracy bound that is not affected by ill conditioning in the weight matrix. A forward error analysis shows that the COD algorithm is stable in this sense, but it is simpler and more efficient than the algorithm proposed by Vavasis. Our forward error bound is contrasted to the backward error analysis of other previous works on weighted least squares.
Introduction
We consider solving the problem
min
for y, where D is a symmetric positive definite m \Theta m matrix, A is an m \Theta n
matrix, y is an n-vector, and b is an m-vector. An equivalent way to write this
problem is "
x
y
b#
which is a special case of an equilibrium system [9]. Applications include optimization
involving a barrier function, finite elements, structural analysis, and
electrical networks. The following assumptions are made throughout the paper.
This work supported by an NSF Presidential Young Investigator grant, with matching funds received
from AT&T and Xerox Corp.
y Center for Applied Mathematics, Cornell University, Ithaca, New York 14853
z Department of Computer Science, Cornell University, Ithaca, New York 14853.
A1. A has rank n, i.e. full column rank.
A2. D is diagonal.
and A2 imply that (1) is a full-rank weighted least-squares problem with a
unique solution, and they allow us to use the norm bound obtained by Stewart.
That bound is given in the following theorem.
Theorem 1 [8] Let D denote the set of all positive definite m \Theta m real diagonal
matrices. Let A be an m\Thetan real matrix of rank n: Then there exist finite constants
-A and -
-A such that
a)
A
A
A
A
A similar result was obtained independently by Todd [10]. In this theorem, the
norm can be any matrix norm induced by a vector norm. In this paper,
Similarly, the condition number of a matrix M is the condition number of M in
the 2-norm, i.e. We make one more assumption.
A3. D is very ill-conditioned.
The ill-conditioning of D arises in certain classes of finite element problems [11],
electrical networks, and it always occurs in the barrier method for optimization
[14]. It also indicates that the coefficient matrix of the least-squares problem is
ill-conditioned. For this reason, the methods typically used to solve least-squares
problems can give highly inaccurate solutions y, as argued by Vavasis [12]. Since
D is the source of ill-conditioning, we use the definition of stability given by
Vavasis.
An algorithm for (1) is stable if, in the presence of finite-precision
arithmetic, an error bound of the form
is satisfied, where y is the true solution, -
y is the computed solution, f(A) is some
function of A not depending on D, and ffl ? 0 is machine roundoff.
Other standard terminology is modified analogously. For example, a well-conditioned
matrix is one for which there is an upper bound on the condition number that
does not depend on D: In order to show that our algorithm is stable, then, we
strive to obtain bounds on norms, condition numbers, and errors that do not
depend on D:
We now present the algorithm.
Algorithm: Orthogonal Decomposition
Step 1: QR factor, with column pivoting,
to get
where Q is an n \Theta n orthogonal matrix, R is an n \Theta m upper triangular matrix,
and P is an m \Theta m permutation matrix.
Step 2: Apply reduced QR factorization (without pivoting) to R T to get
where Z 1 is an m \Theta n matrix with orthonormal columns, and U 1 is an n \Theta n upper
triangular matrix.
Step 3: Solve the following system, via back substitution, for -
y:
Step 4: To get y, multiply the result of Step 3 by Q:
Solution of least-squares problems via QR factorization with column pivoting was
introduced by Golub [3]. Note that the QR factorization for the least-squares
problem occurs in Step 2. The QR factorization in Step 1 is to make the algorithm
stable.
First, we compare this algorithm to Vavasis's NSH method [12]. The NSH
method is the only algorithm in the literature known to stably (in the sense of
Definition 1) solve (1). The NSH algorithm employs nonstandard techniques,
particularly when choosing the nullspace basis for A T D In contrast, the complete
orthogonal decomposition algorithm uses standard techniques that are well
understood, namely QR decomposition and back substitution. Also, our algorithm
requires less work than the NSH algorithm. The NSH method requires
solving an m \Theta m system of equation and thus requires O(m 3 ) flops. The work
for the QR factorizations dominates the work required for the complete orthogonal
algorithm, so this algorithm requires O(mn 2 ) flops. Since
could be much smaller than m), the complete orthogonal decomposition algorithm
requires less work.
The rest of this paper is devoted to analysis of the stability of the orthogonal
decomposition algorithm. Before giving a rigorous stability analysis of the algo-
rithm, we offer an intuitive explanation of why this algorithm finds an accurate
solution. The first step is a QR factorization of a matrix that is well-conditioned
up to a scaling of the columns. So, the result is a computed upper triangular
matrix that is close to the exact upper triangular matrix. It would be useful
to know something about the condition number of this matrix as well. To
minimize confusion assume, without loss of generality, that A T D \Gamma1=2 has been
"pre-pivoted." This means that the columns of A T D \Gamma1=2 are ordered in such a
way that the norms of the first n columns are, loosely speaking, in decreasing
order. In addition, the norms of the first n columns are larger than those of the
last columns. One might suspect that this implies that the entries of D \Gamma1=2
are ordered in the same way. In other words, some inequality similar to
d \Gamma1=2
might hold. This ordering becomes significant in the second step of the algorithm.
Recall that
Notice the Q T A T is upper triangular. So let
Notice also that R is ill-conditioned and that the ill-conditioning arises from
We try to "offset" the effects of D \Gamma1=2 in the following naive way. Let
consider the following:
D
dm
r 1m
dm
It is clear that -
R is well-conditioned, so if the weights are indeed in the order
described above, then it is not difficult to show that there are upper bounds
on all entries. It can also be shown that there are upper bounds on the entries
of
. Using this information, it is not difficult to show that
D 1=2 R (and hence R T -
well-conditioned up
to a scaling of the columns. In the second step, then, we have a least-squares
problem with a coefficient matrix that is well-conditioned up to a scaling of the
columns, namely
min
In traditional analysis, R T being well-conditioned up to a scaling of the columns
indicates that the standard algorithms give an accurate solution. In this analysis,
therefore, one might expect a parallel result.
In the remainder of the paper, we present a rigorous stability analysis of
our proposed algorithm. The next section contains a discussion of numerical
tolerance. Section 3 examines the condition number and the accuracy of the
upper triangular matrix R computed in the first step of the algorithm. Section
4 does the same for U 1 , the upper triangular matrix computed in Step 2. An
analysis of Steps 3 and 4 of the algorithm is presented in Section 5.
A note on numerical tolerance
In the upcoming analysis we assume throughout that any occurrence of exact
linear dependence among the columns of A T is always determined correctly in
Step 1 of the algorithm (QR factorization with pivoting). This requires the use of
a numerical tolerance. To illustrate this point, consider applying the algorithm
when D
A T =B @
Observe that the third column of A T is dependent on the first two. If the QR
factorization of A T D \Gamma1=2 were done in exact arithmetic, this dependence would
be manifested as a '0' in the (3; 3) position of the factored A T D \Gamma1=2 after the first
two QR factorization steps, and column 4 would be chosen for the third pivot.
In finite-precision arithmetic, however, we would expect the (3; 3) entry to be
on the order of machine-epsilon rather than 0. Because column 4 is weighted by
\Gamma20 , the unwanted residual in the (3; position could cause column 3 to be
chosen for the next column pivot instead of column 4. Thus, without modifi-
cation, ordinary QR-factorization with column pivoting procedure has missed a
linear dependence.
We address this problem as follows: after the kth QR factorization step,
we check whether the residual portion (that is, positions n) of any
uneliminated column has become very small (according to some tolerance level)
with respect to the original norm of that column. If so, we change those entries
to zeros. Notice that this test requires very little additional work because the
usual QR factorization algorithm with column pivoting already monitors the
norms of the residual portions of the columns [3]. In this way, if there are exact
dependences among the rows of A, our algorithm does not miss them.
If this numerical test fails to detect exact dependence, then our stability analysis
no longer holds. It can be shown that the test we have proposed will fail only
in the case that there is near-dependence among the columns of A T . However,
in this case, the parameter -A given by Theorem 1 is large, and so our stability
bound (which depends on -A and -
-A ) is not practically applicable.
3 The First QR Factorization
The intuitive discussion in x1 asserted that the ordering of the weights produced
in Step 1 of the algorithm is important in stabilizing the algorithm. It is necessary,
then, to determine how the weights of the basis rows compare to those of other
basis rows and to those of the nonbasis rows. We begin by comparing the norms
of the rows A.
be an n \Theta n matrix whose columns are an arbitrary set of n
rows a i 1
of A that form a basis for the rowspace of A. Then
1-j-n
1-j-n
Proof: Let B be the basis defined above. Without loss of generality, suppose
that ka i Then, we can write
and partition B \Gamma1 as:
Then
This means v T a i n
1. By the Cauchy-Schwarz inequality, kvk - 1=ka i n k. Also,
The first inequality is because v is a row of B \Gamma1 , and the second is from [12].
Thus,
i.e.
min
1-j-n
Multiply both sides of this inequality by the inequality to
obtain
1-j-n
1-j-n
as required. 2
Notice that we have proved that there is a lower bound on any column of any
basis for the rowspace of A.
Suppose now that k ? 0 steps of the factorization have been completed.
Partition the resulting matrix by rows as follows:
a T
a T
Lemma 1 can be extended so that is applies to the residual portions of the rows
of A, i.e. - a T
be an set of n columns of A T such that the first k columns
of B are fa i 1
and B is a basis. Suppose that a i 1
are the first
columns chosen by the column pivoting in the QR factorization. As with -
A
above, write
. ff
Then
k- a
k- a i j k: (8)
Proof: Let defined as above. Then
where R is upper triangular. Comparing this matrix to the partitioned one above,
we see that the columns of -
are - a i k+1
To prove Lemma 2, then, we need
a lower bound on a typical column of -
B. Notice that
R
Therefore,
Now, as above, turn an upper bound on k -
into a lower bound on any column
of -
Using the previous two lemmas, we can determine the relationships between
the weights of the basis rows and those of the nonbasis rows. For the remainder
of the paper assume, as in the intuitive discussion in x1, that the columns of
A T D \Gamma1=2 have been reordered so that no pivoting is necessary. This implies that
not only do the first n rows of A form a basis for the rowspace of A; but there is
also an order that has been imposed on the rows of A: Let
denote the (reordered) entries of D.
Theorem 2 Suppose the first k - 0 steps of the QR factorization have been
completed. If d \Gamma1=2
k+1 is the weight assigned to a k+1 2 B; and d \Gamma1=2
j is the weight
assigned to a j 62 B; then
d k+1
provided a j is linearly independent of the first k basis vectors.
Proof: Let -
A be as defined before Lemma 2 for k - 0. (Notice that when
this is just the matrix A, partioned into rows.) Since B forms a basis for
the rowspace of A,
Assuming a j is linearly independent of the first k basis rows implies
c i a
which means c i 6= 0 for at least one i such that k us take
and
a
Notice that - a
for at least one i. Let l be such that k
- a l =c l@ -
a n g is a basis for f- a
and (8) hold for any basis for the rowspace of A,
k- a j k -
k- a
Recall that the columns of A T D \Gamma1=2 have been reordered so that no pivoting is
necessary. This means that there is an order imposed on the columns of A T D \Gamma1=2 :
More specifically, since we are at step k
k- a j k -
k- a k+1 k:
Thus,
d k+1
k- a k+1 k
k- a j k
k- a
-A
min k+1-i-n k- a i k
which is (9).We also need to know the relationships between the weights of the basis rows of
A. Suppose a i ; a It follows from (7), (8), and the implicit
order indicated by the absence of column pivoting that
Recall that the intuitive argument given in x1 relied on the weights being in the
following
d \Gamma1=2
What we have found, however, is that they are not ordered in exactly this way.
Instead, this ordering holds up to scaling by a constant, i.e.
d \Gamma1=2
d \Gamma1=2
This bound is sufficient for our arguments.
The second step of the algorithm performs a QR factorization on R T . In order
to analyze that step, then, it is necessary to know something about the condition
of R. The relationships between the weights of the rows of A are used in the
following theorem to show that R is well-conditioned up to a scaling of the rows
or the columns. Recall that for a rectangular matrix M , -(M) is the condition
number (in the 2-norm) of M , i.e.
Theorem 3 Let
D a RD 1=2\Gammaa ; where a - 0 and -
~
-( ~
for any 1 - k - n:
Proof: First, we must find an upper bound on k ~
Ck. Since ~
C is a submatrix
of C, then k ~
Therefore, it is sufficient show that there is an upper
bound on kCk. Write C as follows:
D a RD
D a Q T A T D \Gamma1=2 D
D a -
If the entries of C are written explicitly,
a
a
dm
a
r 1m
a
dm
a
Consider -
be the basis consisting of the first n rows
of A. Then -
r ii
i is the ith
row of -
R, then k-r These facts, (9), and (10) imply
the following:-A
fld a
Recall that Theorem 2, and thus (13), holds only when a j is linearly independent
of the first vectors. We must now consider the case not covered by
Theorem 2. Suppose that B and D are defined as before. For each nonbasis row
a j there is a 1 - k - n such that a j is linearly independent of the first
basis vectors, but is linearly dependent on the first k basis vectors. So,
where c k 6= 0: Now, suppose that k steps of the QR factorization have been
completed. Then
ff 1i
ff ki3
is as in Lemma 2). After this point, transformations act
only on - a j : This gives
ff 1i
ff ii3
telling us that - r these entries do not contribute to the norm
in (13). Thus, (13) holds even in the case where a j is not linearly independent
of the first vectors.
If c T
th row of C, then
1-i-n
1-i-n
The next step is to find an upper bound on
Recall that
If is easy to show that ~
1 is a submatrix of C \GammaT
To obtain
an upper bound on kC \GammaT
1 k, we use the following fact, which will be proved after
the current proof.
Fact: If C
kC \GammaT
and the bound on the condition number is
-( ~
Thus, we have proved our claim. 2
Recall that we must still prove the fact used in the above proof. We state it
in the form of the following lemma and give the proof below.
defined as in the previous theorem.
Then
kC \GammaT
Proof: Recall that
D a -
D a -
Notice that
where B is the rowspace basis consisting of the first n rows of A. If
then
D a
a
l 11
a
l 21
a
l 22
a
l n1
a
l
a
l
If c T
i is the i th row of C \GammaT, then
kC \GammaT
1-i-n
1-i-n
as claimed. 2
We now know that R is not only well-conditioned up to a scaling of the rows,
but it is well-conditioned up to a scaling of both the rows and the columns. This
result will be useful later in the analysis.
Standard results tell us that the differences between the columns of the computed
R and those of the exact matrix R are small. Since the second step
of the algorithm is a QR factorization of R T , we need to know that the same is
true of the rows. This is given by the following theorem.
Theorem 4 Let r T
j be the j th rows of R and -
R; respectively. Then
Proof: Recall that
According to Wilkinson [13],
where r i is the i th column of R; and - r i is the i th column of -
R. This gives us a
bound on the elementwise error, namely
We can now find a bound on the normwise error of the rows of R; namely
Consider kr
Suppose that k - 0 steps of the QR factorization have been
completed. Partition X in a manner similar to that of A for Lemma 2. Then
no pivoting is necessary, the columns of X are ordered such
that k- x
Thus,
Taking the square root of both sides gives the required result. 2
We have now established that the QR factorization in the first step of the
algorithm gives an upper triangular matrix R that is well-conditioned up to a
scaling of the rows or the columns. It also yields a computed matrix, -
R; whose
rows are close to those of R. With these results in hand, we move on to the
analysis of the second step of the algorithm.
4 The Second QR Factorization
Recall that in this step we use the "skinny" QR factorization,
where Z 1 is an m \Theta n matrix with orthonormal columns and U 1 is an n \Theta n upper
triangular matrix. The results of x3 tell us that R T is well-conditioned up to a
scaling of the columns. The QR factorization in this step, then, gives an upper
triangular matrix -
U 1 that is close to the exact upper triangular matrix U 1 . In
addition, the results of x3 can be used to prove similar results about U
let
is the coefficient
matrix of the system of equations in the back substitution step, the following
theorem concerning the condition of U 1 will be quite useful.
Theorem 5 Let U 1 and -
D be defined as above. Then
D a U 1
Proof: First, notice that
It follows from
the fact that U \GammaT
RR
~
R ~
where e k is the k th column of the k \Theta k identity matrix. So,
R ~
R T
z T
R ~
R T ~
z T
d 1=2
~
where z k is the kth column of Z is the kth column of R T and x is the last
column of
R ~
R T ~
. Multiplying both sides by d \Gammaa
~
d \Gammaa
~
z T
d 1=2\Gammaa
~
D a x:
We show that there is an upper bound on the right-hand side as follows:
z T
d 1=2\Gammaa
~
D a
~
D a xk
D a
R ~
R T ~
R ~
R T ~
R ~
R T ~
R ~
R T ~
RD \Gammaa D a ~
R T ~
RD \Gammaa
D a ~
R T ~
In the fifth line of the above inequality, - r 11 is the (1,1) entry of -
Notice that if \Gamma 1- a - 1; then Theorem 3 applies. So,
kd \Gammaa
~
RD \Gammaa
D a ~
R T ~
kd \Gammaa
1-i-n
kd \Gammaa
~
In order to find and upper bound on k -
D a U 1
D 1=2\Gammaa k; we do the
D a U 1
kd 1=2\Gammaa
D a u
1-i-n
kd 1=2\Gammaa
D a u
1-i-n
kd 1=2
where u i is the ith column of U 1 , r i is the ith column of R T , and -r i is the ith
column of -
R: The third line of the inequality is derived from the second using
(10), the fact that U 1 is upper triangular, and U last line is
obtained by applying Theorem 3. Thus for \Gamma 1- a - 1;
D a U 1
D a U 1
:Now that we know that -
D a U 1
D 1=2\Gammaa is well-conditioned for \Gamma 1- a - 1, we
move on to the analysis of the remainder of the algorithm.
5 Finding the Solution y
In analyzing the remainder of the algorithm, we first show that a small error
introduced in the back substitution step. In Step 3 of the algorithm, the upper
triangular system
is solved for -
y: (Note that this is slightly different from the system given in Step
3 of the algorithm as presented in x1 because we are assuming that the columns
of A T D \Gamma1=2 have been "pre-pivoted.") Instead of working with the system given
above, consider the following system:
n) as before. In working
through the steps of back substitution, we find that solving this system is equivalent
to solving the original one, even in floating point arithmetic. (In other words,
a rescaling of the rows does not change the numerical bounds.) Recall from the
last section that -
D a U 1
D 1=2\Gammaa is well-conditioned for \Gamma 1- a - 1. Therefore,
standard back substitution results apply when showing that the error at this
step is small. The following theorem states that error bound.
Theorem 6 Let -
y be the exact solution to -
y
be the computed solution. Then
Proof: Let -
y be the computed solution to the above system. Then -
y is the
exact solution to the nearby system of equations,
The matrix E accounts for errors during the back substitution and jEj - ffl \Delta
or
Substituting for - y on the right-hand side yields
E)
Thus,
E)
E)
as claimed. 2
In the theorem above, the errors in the computation of U 1 itself (which also
contribute to the error in -
y) are not included, but be could accounted for as
a somewhat larger perturbation matrix E. As we we have already argued at
the end of x3, the errors in computing U 1 are small. In fact, the errors made in
forming each row of U 1 are small with respect to the norm of that row. Therefore,
the perturbation matrix E is small with respect to -
Explicitly including
this analysis in the previous theorem would make the proof more complicated,
but the bound would be qualitatively the same.
The final step is to obtain y by multiplying - y by Q. Let - y be the computed
result. We assume that - y accounts for the errors during both this step and the
previous step. The error bound is obtained as follows:
Notice that the error bound is of the form
Thus, our algorithm satisfies the definition of stability.
6 Summary and Open Questions
The least-squares problem
min
where D is an m \Theta m diagonal, positive definite, ill-conditioned matrix; A is an
m \Theta n full-rank matrix; y and b are n-vectors, has a unique solution. Because of
the ill-conditioning of D, the standard methods for solving least-squares problems
do not find an accurate solution. In an attempt to find a standard algorithm
that will solve this problem accurately, we have employed complete orthogonal
decomposition. This involves four steps, given in x1. We then proceeded to show
that this algorithm is stable, as defined in x1.
The first step is a "stabilizing" QR factorization that gives an upper triangular
matrix that is well-conditioned up to a scaling of the rows or the columns. Using
this, we are able to show that the result of the second step is an upper triangular
matrix that is also well-conditioned up to a scaling of the rows or the columns. We
were then able to show that there is a bound on the error in the back substitution
step. It was then easy to show that there is a small error introduced in the
last step. Thus, this algorithm gives an accurate solution to this least-squares
problem.
Now that we know that the algorithm is stable, there are a couple of further
questions.
1. This algorithm has been implemented using dense methods. In many
applications, the matrix A is sparse. Can we implement this algorithm in such a
way that it takes advantage of that sparsity?
2. The results thus far are theoretical. This algorithm has not yet been
extensively tested in applications. The question, then, is whether or not this
algorithm is effective in applications. We are currently beginning tests of our
algorithm in interior point methods [6].
The problem of stably solving the ill-conditioned equilibrium system in barrier
methods for optimization has received a fair amount of attention [2]. In the
case of barrier methods for linear programming (that is, interior point methods),
the equibrium system reduces to weighted least squares, which is the problem
addressed by this paper. Other authors have recently looked at ill-conditioning
in barrier methods including Coleman and Liu [1], M. Wright [15], S. Wright [16],
Nash and Sofer [7], and Gould [5].
The differences between these other works and ours may be summarized as fol-
lows. These other works typically look at the more general problem minkH \Gamma1=2 (Ay\Gamma
b)k where H is symmetric, positive definite, but not necessarily diagonal, which
is a problem that we currently cannot address with our techniques. On the other
hand, these other authors all make an assumption that the large and small entries
on the diagonal of H have some correlation with the columns of A T . This
corresponds to a nondegeneracy assumption about the underlying optimization
problem. In contrast, our method does not involve any restrictions about where
"large" versus "small" entries of D can appear, and thus our method has no difficulty
when there is degeneracy or near-degeneracy in the underlying optimization
problem.



--R

An interior Newton method for quadratic pro- gramming
Practical Optimization.
Numerical methods for solving linear least squares problems.
Matrix Computations
On the accurate determination of search directions for simple differentiable penalty functions.
Complete orthogonal decomposition in interior point methods.
A barrier method for large-scale constrained opti- mization
On scaled projections and pseudoinverses.
A framework for equilibrium equations.
A Dantzig-Wolfe-like variant of Karmarkar's interior-point linear programming algorithm
Stable finite elements for problems with wild coefficients.
Stable numerical algorithms for equilibrium systems.
The Algebraic Eigenvalue Problem.
Interior methods for constrained optimization.
Determining subspace information from the Hessian of a barrier function.
Stability of linear algebra computations in interior-point methods for linear programming
--TR
numerical stability;interior-point methods;weighted least squares;equilibrium systems;QR factorization;forward error analysis