RELATIVE ERRORS VERSUS RESIDUALS OF APPROXIMATE

SOLUTIONS TO LINEAR ALGEBRAIC EQUATIONS AND LEAST SQUARES PROBLEMS

Jiu Ding

Abstract. Lower and upper bounds are given to the ratio of the error of the perturbed solution and the distance of the unperturbed solution to the null space of the matrix for general consistent systems of linear equations and least squares problems.

1. Introduction. In this paper we consider the perturbation of the non– homogeneous system of linear algebraic equations

Ax = b,  (1)

where A ∈ ℝ^{m×n} is an m × n real matrix and b ∈ ℝ^{m} is an m dimensional real vector. Error estimates of numerically stable algorithms, such as Gaussian–typed direct methods for solving (1), are reduced to, via the backward error analysis of Wilkinson [4], the study of the difference of an exact solution of the perturbed system

(A + E)y = b + e  (2)

or its corresponding least squares problem to an exact solution of the original system (1) or its corresponding least squares problem. Hence, the perturbation analysis for solving (1) is important in numerical linear algebra [4].

When A in (1) is nonsingular, a well–known perturbation result says that if y approximates the exact solution x = A^{−1}b of (1), then

\frac{1}{κ} \frac{∥r_{y}∥}{∥b∥} ≤ \frac{∥y − x∥}{∥x∥} ≤ κ\frac{∥r_{y}∥}{∥b∥},  (3)

where r_{y} = Ay − b is the residual associated with y and κ = ∥A∥∥A^{−1}∥ is the condition number of A with ∥ · ∥ standing for the vector norm and the induced matrix norm.

The inequalities in (3) show that the relative error of the perturbed solution y for solving the nonsingular system (1) is equivalent to the norm of the residual r_{y} divided by ∥b∥, but the condition number κ determines how sensitive the perturbed solution is relative to the perturbation of the original system.

In this paper we extend (3) to general m × n matrices A. Specifically, we give lower and upper bounds to the ratio of the error of the perturbed solution and the distance of an unperturbed solution to the null space of A for general consistent systems of linear equations and least squares problems. Such bounds involve the residual of the approximate solution as in the classical case. Some ideas in this paper were originally due to [2] in which an optimal upper bound for the perturbation analysis of orthogonal projections has been obtained.

In the next section we consider the case when (1) is consistent, and the more general case of least squares problems will be dealt with in Section 3.

2. The Case When b ∈ R(A). Let N(A) and R(A) be the null space and the range of A, respectively. We need the concept of the generalized inverse A^{†} of A associated with two projections P: ℝ^{n} → N(A) and Q: ℝ^{m} → R(A), which is defined as the unique n × m matrix A^{†} satisfying

AA^{†}A = A, A^{†}AA^{†} = A^{†}, A^{†}A = I − P, AA^{†} = Q.

When the norms for ℝ^{n} and ℝ^{m} are induced from their respective inner products, e.g., the Euclidean norm ∥·∥_{2}, we further require that P and Q be self–adjoint (i.e., P^{∗} = P and Q^{∗} = Q), so that A^{†}A and AA^{†} are orthogonal projections with respect to their inner products. In this case A^{†} is called the Moore–Penrose generalized inverse of A. See [1] or [3] for more details on A^{†}.

Again let κ = ∥A∥∥A^{†}∥ be the condition number of A. Let d(x, N(A)) = inf{∥x − z∥ : z ∈ N(A)} be the distance of x to N(A) and let r_{y} = Ay − b be the residual associated with y.

\underline{Theorem 2.1}. Suppose b ∈ R(A). Then for any y ∈ ℝ^{n}, there is a solution x to (1) such that

\frac{1}{κ} \frac{∥r_{y}∥}{∥b∥} ≤ \frac{∥y − x∥}{∥A^{†}b∥} ≤ \frac{∥y − x∥}{d(x, N(A))} ≤ κ\frac{∥r_{y}∥}{∥b∥}.  (4)

\underline{Proof}. Let x = A^{†}b + (I − A^{†}A)y be the projection of y onto the solution set of (1) along N(A) [2]. Then

∥A^{†}b∥ = ∥x − (I − A^{†}A)y∥ ≥ d(x, N(A)),

so the middle inequality of (4) is valid. Since Ax = b, we have

A(y − x) = Ay − b.  (5)

On the other hand, by the construction of x,

y − x = A^{†}(Ay − b).  (6)

Let z ∈ N(A) be arbitrary. Then (6) gives that

\begin{aligned}\frac{∥y − x∥}{∥x − z∥} & ≤ \frac{∥A^{†}∥∥Ay − b∥}{∥x − z∥} = \frac{∥A∥∥A^{†}∥∥r_{y}∥}{∥A∥∥x − z∥} \\ & ≤ \frac{κ∥r_{y}∥}{∥Ax∥} = κ\frac{∥r_{y}∥}{∥b∥}, \end{aligned}

which gives the right inequality of (4). Now (5) implies that

\begin{aligned}\frac{∥r_{y}∥}{∥b∥} & = \frac{∥Ay − b∥}{∥b∥} ≤ \frac{∥A∥∥y − x∥}{∥b∥} \\ & ≤ \frac{∥A^{†}∥∥A∥∥y − x∥}{∥A^{†}∥∥b∥} ≤ κ\frac{∥y − x∥}{∥A^{†}b∥}. \end{aligned}

This gives the left inequality of (4).

\underline{Remark 2.1}. When the norms ∥ · ∥ are given via the inner products, ∥y − x∥ is exactly the minimal distance of y to the affine set of all solutions to (1), and in this case (4) gives lower and upper bounds of this distance with respect to the distance of the solution x to the null space of A.

Since RankA = n implies that N(A) = {0} and x = A^{†}b, we immediately have the following.

\underline{Corollary 2.1}. If in addition RankA = n, then

\frac{1}{κ} \frac{∥r_{y}∥}{∥b∥} ≤ \frac{∥y − A^{†}b∥}{∥A^{†}b∥} ≤ κ\frac{∥r_{y}∥}{∥b∥}.  (7)

In particular, if A is nonsingular, then (4) is reduced to (3).

The residual r_{y} depends on E and e if y turns out to be a solution of (2). In this case, subtracting (1) from (2) gives that

r_{y} = A(y − x) = e − Ey,

from which it follows that

∥e − Ey∥ = ∥r_{y}∥ ≤ ∥e∥ + ∥E∥∥y∥.  (8)

Therefore we have the following corollary.

\underline{Corollary 2.2}. If in addition b + e ∈ R(A + E), then for any solution y to (2), there is a solution x to (1) such that

\frac{1}{κ} \frac{∥e − Ey∥}{∥b∥} ≤ \frac{∥y − x∥}{∥A^{†}b∥} ≤ \frac{∥y − x∥}{d(x, N(A))} ≤ κ\frac{∥e∥ + ∥E∥∥y∥}{∥b∥}.  (9)

3. Least Squares Problems. Now we consider the least squares problem

∥Ax − b∥ = min_{z∈ℝ^{n}} ∥Az − b∥  (10)

and its perturbed form

∥(A + E)y − (b + e)∥ = min_{z∈ℝ^{n}} ∥(A + E)z − (b + e)∥.  (11)

Here we assume that the norms on ℝ^{n} and ℝ^{m} are determined by inner products, or more generally the norms are such that A^{†}b is a solution of (10).

\underline{Theorem 3.1}. For any y ∈ ℝ^{n} there is a solution x to (10) such that

\frac{1}{κ} \frac{∥r_{y}∥}{∥AA^{†}b∥} ≤ \frac{∥y − x∥}{∥A^{†}b∥} ≤ \frac{∥y − x∥}{d(x, N(A))} ≤ κ\frac{∥r_{y}∥}{∥AA^{†}b∥},  (12)

where r_{y} = Ay − AA^{†}b.

\underline{Proof}. Let x = A^{†}b + (I − A^{†}A)y be the (orthogonal) projection of y onto the solution set of (10) along N(A) (see, e.g., [2]). Then the middle inequality of (12) is obvious. Now, as in the previous section,

y − x = A^{†}(Ay − b) = A^{†}(Ay − AA^{†}b),  (13)

from which it follows that

A(y − x) = AA^{†}(Ay − AA^{†}b) = Ay − AA^{†}b.  (14)

The proof from here on is exactly the same as that for Theorem 2.1 except that b is replaced with AA^{†}b, so it will be omitted.

\underline{Remark 3.1}. Since the least squares problem (10) is equivalent to the consistent system of linear equations Ax = AA^{†}b, (12) is also a direct consequence of (4) by replacing b with AA^{†}b.

Now we estimate ∥r_{y}∥ if y solves (11). Let

Ax − b = r, (A + E)y − (b + e) = r′,  (15)

and subtract the first equality from the second one, we have

r_{y} = A(y − x) = e − Ey + r′ − r.  (16)

Since A^{†}r = 0 and (A + E)^{†}r′ = 0, (13) and (16) imply that

y − x = A^{†}r_{y} = A^{†}(e − Ey + r′ − r)
= A^{†}(e − Ey + r′) = A^{†}(e − Ey) + [A^{†} − (A + E)^{†}]r′.  (17)

By the decomposition formula [3] that

\begin{aligned}A^{†} − (A + E)^{†} & = A^{†}E(A + E)^{†} − A^{†}(A^{†})^{∗}E^{∗}[I − (A + E)(A + E)^{†}] \\ & − (I − A^{†}A)E^{∗}[(A + E)^{†}]^{∗}(A + E)^{†}, \end{aligned}

we obtain

[A^{†} − (A + E)^{†}]r′ = −A^{†}(A^{†})^{∗}E^{∗}r′,  (18)

from which (17) can be written as

y − x = A^{†}(e − Ey) − A^{†}(A^{†})^{∗}E^{∗}r′.  (19)

It follows from (14) and (19) that

r_{y} = A(y − x) = AA^{†}(e − Ey) − AA^{†}(A^{†})^{∗}E^{∗}r′
= AA^{†}(e − Ey) − (EA^{†})^{∗}r′.  (20)

Therefore,

∥AA^{†}(e − Ey) − (EA^{†})^{∗}r′∥ = ∥r_{y}∥ ≤ ∥e∥ + ∥E∥∥y∥ + ∥EA^{†}∥∥r′∥. (21)

Thus, we have the following corollary.

\underline{Corollary 3.1}. For any solution y to (11), there is a solution x to (10) such that

\frac{1}{κ} \frac{∥AA^{†}(e − Ey) − (EA^{†})^{∗}r′∥}{∥AA^{†}b∥} ≤ \frac{∥y − x∥}{∥A^{†}b∥} ≤ \frac{∥y − x∥}{d(x, N(A))}

≤ κ\frac{∥e∥ + ∥E∥∥y∥ + ∥EA^{†}∥∥r′∥}{∥AA^{†}b∥}.  (22)

Finally it is noted that Theorem 2.1 is still true for bounded linear operators with closed range between general Banach spaces such that the generalized inverse is well–defined and Theorem 3.1 can be easily extended to the case of general Hilbert spaces.

\underline{Acknowledgement}. Research was supported in part by an Aubrey K. Lucas Endowment award from the University of Southern Mississippi.

\underline{References}

1. S. Campbell and C. Meyer, Generalized Inverses of Linear Transformations, Pitman, London, England, 1979.

2. J. Ding, “Perturbation Results for Projecting a Point Onto a Linear Manifold,” SIAM J. Matrix Anal. Appl., 19 (1998), 696–700.

3. C. L. Lawson and R. J. Hanson, Solving Least Squares Problems, Prentice–Hall, Englewood Cliffs, New Jersey, 1974.

4. D. S. Watkins, Fundamentals of Matrix Computations, John Wiley & Sons, New York, 1991.

Jiu Ding Department of Mathematics The University of Southern Mississippi Hattiesburg, MS 39406–5045 email: jding@yizhi.st.usm.edu