Now it is time for a few illustrative examples ...
Suppose you are learning about normal-equations, orthogonal transformations, QR decompositions, etc.6 You have read the proper sections in your text(s), and you want to try your hand at it to see if you really understand it.
First you create an over-determined coefficient matrix, 3
parameters, and 5 equations (a
). Then you create an
observation matrix (b
):
> a = [3,4,1;0,2,2;0,0,7;zeros(2,3)]; > b = [14;10;21;6;2];
You've just read that the RLaB operator `\
' solves systems
of equations, so you try it out:
> x = a \ b x = 1 2 3
You check the answer (note that this is a contrived problem):
> b - a*x -7.11e-15 -1.78e-15 -1.42e-14 6 2
and it looks ``OK''. The residual in the first three rows is near the machine-epsilon7. Now you wish to follow the example in the text more closely, in an attempt to reinforce your reading. The text has stated that the ``normal equations'' are: (ATA)x = (ATb).
Not having read the chapter on Gaussian elimination, and matrix inverses yet you try:
> x = inv (a'*a) * (a'*b) x = 1 2 3
Fortunately, the problem we are working with is not ill-conditioned, otherwise we may have produced a terrible result with the above procedure. If you want to pursue the reasoning behind the previous statement I suggest you read the section ``Linear Systems of Equations''.
Well, this is all too easy, now you want to get dirty, so you move on to orthogonal transformations. You have read about the construction of Householder vectors and reflections; now you would like to try it first-hand. You know that:
Where v is the Householder vector used to form the reflection matrix. First you must construct the vector. Your text8 tells you that a good method for constructing the Householder vector is:
> a = rand(5,2); // Start out with a more difficult [A] > a a = 0.655 0.265 0.129 0.7 0.91 0.95 0.112 0.0918 0.299 0.902 > ac1 = a[;1]; // grab the 1st column of [a] to work with > u = norm (ac1, "2"); // compute the 2-norm of [ac1] > v[2:5] = ac1[2:5] / (ac1[1] + sign (ac1[1])*u) v = 0 0.0705 0.498 0.0611 0.164 > v[1] = 1; > v = v'; // make v a column vector
By using the matrix creation, and element referencing features you have generated the vector in 4 commands. We could have used a signal command
> v = [ 1 , a[;1][2:5]' / (a[1;1] + sign(a[1;1])*norm(a[;1],"2")) ]'
But, this is less than clear. Note that in this case, since we are working with vectors, we only use a single index when subscripting the variables.
Now that we have our Householder vector, we are ready to assemble the Householder reflection (matrix).
> P = eye (5,5) - 2*(v*v')/(v'*v) P = -0.558 -0.11 -0.776 -0.0952 -0.255 -0.11 0.992 -0.0547 -0.00671 -0.018 -0.776 -0.0547 0.614 -0.0474 -0.127 -0.0952 -0.00671 -0.0474 0.994 -0.0156 -0.255 -0.018 -0.127 -0.0156 0.958 > P*a -1.17 -1.2 -1.65e-17 0.596 -1.54e-16 0.22 -1.31e-17 0.00217 -5.39e-17 0.662
As you can see it worked out just like they said it would. All the elements of the first column of A, below the diagonal, have been zeroed out9. In this manner we can proceed to transform A into an upper triangular matrix.