Many functions are included with the RLaB source distribution.
Functions can be found in the distribution subdirectories
./rlib
, ./toolbox
, and ./examples
. These
directories are normally installed under /usr/local/lib/rlab
.
We will continue with some simple examples demonstrating function creation and usage. We will carry on with the exercise of learning least-squares techniques.
In some earlier examples we played with solving a set of normal equations, and tried a simple experiment with Householder reflections. Now we want to try out this technique, and decompose a matrix into two matrices: Q and R; such that A = QR.13
We are going to decompose an entire matrix, so we will want to automate the procedures we used in previous examples. The first was creating a Householder vector. Instead of typing in our function at the command-line, we will use a text-editor to create the function in a file so that we can correct our mistakes without retyping the entire function.
// // house_v(): Given an N-vector V, generate an n-vector V // with V[1] = 1, such that (eye(n,n) - 2*(V*V')/(V'*V))*X // is zero in all but the 1st component. // static (sign) house_v = function(v) { local(v) n = max( size(v) ); u = norm(v, "2"); if(u != 0) { b = v[1] + sign(v[1])*u; if(n > 1) { v[2:n] = v[2:n]/b; } } v[1] = 1; return v; }; sign = function ( X ) { if (X >= 0) { return 1; else return -1; } }
Note that our new function is more complicated than our earlier
``one-liner''. This is due to the fact that the function is more
efficient, and does some input-checking. Notice that the variables
b
, n
, u
, and v
are local; these
local variables will never be seen by the user, and will not
interfere with any pre-existing variables by the same name in the
global-workspace.
Also note that the function argument, v
is copied to the
function local variable v
. This prevents the function from
changing the values in the input argument, and thus destroying the
contents of the caller's variable.
One other important feature of the new function is the usage of
the sign
function. house_v
requires that the sign
function return either 1 or -1. The RLaB built-in sign
function will return zero when its argument is zero; this
behavior is unacceptable, so we have written our own sign
function. Declaring our new sign function to be static means that
it will only affect statements within the file house_v.r
.
To use our new function type:
> a = rand (4,2); > x = house_v (a[;1]) x = 1 0.434 0.042 0.412
Now that we can generate a Householder vector, we need to automatically form the Householder reflection, and use it to reduce A to upper triangular form.
// P.r // P: Generate P matrix P = function ( V ) { m = max( size(V) ); return [ eye(m,m) - 2*(V*V')./(V'*V) ]; };
This function is a small one, and simply implements the formula we demonstrated earlier. We can test out our two new functions like so:
> ( p1 = P (house_v (a[;1])) ) * a -1.49 -1.11 -4.09e-17 0.0174 -2.87e-18 0.354 -8.89e-17 -0.779 > p1' * p1 1 4.79e-17 2.75e-18 1.8e-17 4.79e-17 1 -4.78e-18 -2.51e-18 2.75e-18 -4.78e-18 1 9e-18 1.8e-17 -2.51e-18 9e-18 1
Our new function seems to be working as expected. The computed
Householder reflection, p1
, zeros out all of the elements
below the first in column one of a
. Additionally p1
is an orthogonal transformation, as demonstrated by computing
P1TP1. It is usually more efficient to build up
programs as a collection of simple functions, like we are doing
here, testing each as it is written, and making the appropriate
fixes.
Function debugging can be easily accomplished by simply removing key
`;
' to turn expression printing on. Additionally, one can
comment out local
statements so that a function's variables
can be examined after function execution.
Now there are two more pieces, a better implementation of P
called house_row
14, and the final function (house_qr
) that
will apply the transformations in sequence to A to produce the
upper-triangular R, and the orthogonal Q.
// // house_row(): Given an MxN matrix A and a non-zero M-vector V // with V[1] = 1, the following algorithm overwrites A with // P*A, where P = eye(m,m) - 2*(V*V')/(V'*V) // house_row = function(A, v) { local(A) b = -2/(v'*v); w = b*A'*v; A = A + v*w'; return A; };
// house_qr.r // Given A, with M >= N, the following function finds Householder // matrices H1,...Hn, such that if Q = H1*...Hn, then Q'*A = R is // upper triangular. // House.qr returns a list containing `q', and `r' rfile house_row rfile house_v rfile P house_qr = function ( A ) { local (A) m = A.nr; n = A.nc; v = zeros(m,1); q = eye (m, m); for(j in 1:n) { // Generate the Householder vector v[j:m] = house_v( A[j:m;j] ); // Apply the Householder reflector to A A[j:m;j:n] = house_row( A[j:m;j:n], v[j:m] ); // Create Q if(j < m) { q = P ([ zeros (j-1,1); 1; v[j+1:m] ]) * q; } } return << q = q'; r = A >>; };
Notice the three rfile
statements near the top of the file.
These statements ensure that the user-function dependencies are
resolved before we try and execute the function. Also note
how house_qr
returns two matrices in a list, with element
names q
, and r
. We can use, and test, this new
function like:
> x = house_qr ( a ) q r > x.q * x.r 0.7 0.96 0.95 0.915 0.0918 0.441 0.902 0.0735 > a a = 0.7 0.96 0.95 0.915 0.0918 0.441 0.902 0.0735
A visual comparison shows that our function does indeed work.
Now we wish to use this factorization to compute a solution to our
original least-squares problem. Since we have decomposed A into
two matrices, one of which is upper triangular, we can reformulate
the problem as a simple back-substitution. The RLaB built-in
function solve
will do this for us, since sr
is
already upper-triangular, all solve
will do is the
back-substitution.
> a = [3,4,1;0,2,2;0,0,7;zeros(2,3)]; > b = [14;10;21;6;2]; > x = house_qr (a); > sq = x.q[;1:3]; > sr = x.r[1:3;]; > z = b'*sq; > sol = solve (sr, z') sol = 1 2 3 > b - a*sol 0 0 0 6 2
Now that we have built our own qr()
I should tell you that
RLaB has a built-in qr()
that is much more robust, and
significantly faster.