Example

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_row14, 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.