| |
- If A is a n × m matrix of full column rank (hence m <= n)
then there exists a decomposition
A = U R (*)
with U n × m orthogonal and R m × m
upper triangular and invertible. If the length of the columns of U is normalized to one
and the diagonal elements of R are chosen real and positive, then this decomposition is
unique. U is an orthonormal basis of the range space of A then.
- Given such a decomposition a least squares problem
||Ax - b || = minx
with the euclidean norm ||.|| can be solved using the triangular solve
R x = UT b .
In order to see this choose W as the orthogonal complement of U,
multiply Ax - b by (U,W)T.
||WTb|| is the length of the optimal residual.
This technique avoids the explicit building of the normal equations and is
the preferred solution way.
- There exist several algorithms for computing the decomposition (*) above,
and we present here three of them: the classical Gram-Schmidt orthogonalization,
the so called modified Gram-Schmidt method, and Householders QR decomposition with
explicit construction of the range space.
In the following, A(.,i) denotes the i-th column of A and
similar for the other matrices.
- The classical Gram-Schmidt algorithm orthogonalizes one column after the
other w.r.t to all previously computed columns of U:
For i=1,2,...m:{
αj = U(.,j)T A(.,i) , j=1,..,i-1
Rj,i = αj , j=1,..,i-1
vi = A(.,i) - ∑j=1 to i-1 αj U(.,j)
Ri,i = || vi||
U(.,i) = vi / || vi||
}
- The modified Gram - Schmidt method orthogonalizes in step i all remaining columns of A
with respect to the just computed i-th column of U, and computes R rowwise here.
Here the original A is overwritten by these intermediate columns:
For i=1,2,...m:
{Ri,i = ||A(.,i)||
U(.,i) = A(.,i)/Ri,i
For j=i+1,...,m
{Ri,j = U(.,i)T A(.,j)
A(.,j) = A(.,j)-Ri,j U(.,i)
}}
With no roundoff both algorithms produce exactly the same result, but you
will see here that even in the high precision arithmetic IEEE double used here
this is far from being true. Thereas the classical method is practically useless
the modified algorithm can be shown (work done by Ake Bjoerck) to produce errors no larger
than machine precision times condition number of A (up to a small constant independent on A).
- Householders method works quite differently: Here m reflector matrices
Qi are applied from the left on A with the aim to transform it
to the upper triangular R with a zero matrix of dimension (n-m) × m
appended. For a least squares problem these are applied to b too, resulting in a vector
c say. Then there remains the linear system
R x = (c)1m
to be solved. In order to compute the orthogonal basis of the range space of A
one must apply these reflector matrices from the right to the unit matrix, i.e.
U = (Q1 × ... × Qm)(.,1:m)
This U will always be orthonormal to machine accuracy, as shown by J.H.Wilkinson.
For reasons of fairness we compute x here also using this U.
The reflector matrices are uniquely defined by the normal to the hyperplane
with respect to which the reflection is done. We use here the form
Qi = I - vi viT with ||vi || = √2
You have the choice to see these normals explicitly.
- This application uses these three solution processes in order to allow you
to assess their usefulness.
|