Page 449 - Matrix Analysis & Applied Linear Algebra

5.13 Orthogonal Projection 445

5.13.22. Cimmino’s Reflection Method. In 1938 the Italian mathematician
Gianfranco Cimmino used the following elementary observation to con-
struct an iterative algorithm for solving linear systems. For a 2 × 2 sys-
tem Ax = b, let H 1 and H 2 be the two lines (hyperplanes) defined
by the two equations. For an arbitrary guess r 0 , let r 1 be the reflection
of r 0 about the line H 1 , and let r 2 be the reflection of r 0 about the
line H 2 . As illustrated in Figure 5.13.9, the three points r 0 , r 1 , and
r 2 lie on a circle whose center is H 1 ∩H 2 (the solution of the system).

Figure 5.13.9
The mean value m =(r 1 + r 2 )/2is strictly inside the circle, so m is a
better approximation to the solution than r 0 . It’s visually evident that
iteration produces a sequence that converges to the solution of Ax = b.
Prove this in general by using the following blueprint.
n
(a) For a scalar β and a vector u ∈ such that u =1,
2
T
consider the hyperplane H = {x | u x = β} (Exercise 5.13.17).
Use (5.6.8) to show that the reflection of a vector b about H
T
is r = b − 2(u b − β)u.
n×r
(b) For a system Ax = b in which the rows of A ∈ have been
scaled so that A i∗ =1 for each i, let H i = {x | A i∗ x = b i }
2
be the hyperplane defined by the i th equation. If r 0 ∈ r×1 is
an arbitrary vector, and if r i is the reflection of r 0 about H i ,
explain why the mean value of the reflections {r 1 , r 2 ,..., r n } is
T
m = r 0 − (2/n)A ε, where ε = Ar 0 − b.
T
(c) Iterating part (b) produces m k = m k−1 − (2/n)A ε k−1 , where
ε k−1 = Am k−1 − b. Show that if A is nonsingular, and if

T
x = A −1 b, then x−m k = I − (2/n)A A k (x−m 0 ). Note:
k
T
It can be proven that I − (2/n)A A → 0 as k →∞, so
m k → x for all m 0 . In fact, m k converges even if A is rank
deficient—if consistent, it converges to a solution, and, if incon-
sistent, the limit is a least squares solution. Cimmino’s method
also works with weighted means. If W = diag (w 1 ,w 2 ,...,w n ),
T
where w i > 0 and w i =1, then m k = m k−1 −ωA Wε k−1
is a convergent sequence in which 0 <ω < 2isa “relaxation
parameter” that can be adjusted to alter the rate of convergence.