ORTVEC Call

Language Reference

ORTVEC Call

provides columnwise orthogonalization by the Gram-Schmidt process and stepwise QR decomposition by the Gram-Schmidt process

CALL ORTVEC( w, r, $\rho$ , lindep, v<, q>);

The ORTVEC subroutine returns the following values:

w

If the Gram-Schmidt process converges (lindep=0), w is the m ×1 vector w orthonormal to the columns of Q, which is assumed to have $n \leq m$ (nearly) orthonormal columns. If the Gram-Schmidt process does not converge (lindep=1), w is a vector of missing values. For stepwise QR decomposition, w is the (n+1)^th orthogonal column of the matrix Q. If there is no matrix Q, that is, if the q argument is not specified, w is the normalized value of the vector v,

$w= \frac{v}{\sqrt{v^' v}}$

r

If the Gram-Schmidt process converges (lindep=0), r specifies the n ×1 vector r of Fourier coefficients. If the Gram-Schmidt process does not converge (lindep=1), r is a vector of missing values. If the q argument is not specified, r is a vector with zero dimension. For stepwise QR decomposition, r contains the n upper triangular elements of the (n+1)^th column of R.

$\rho$

If the Gram-Schmidt process converges (lindep=0), $\rho$ specifies the distance from w to the range of Q. Even if the Gram-Schmidt process converges, if $\rho$ is sufficiently small, the vector v may be linearly dependent on the columns of Q. If the Gram-Schmidt process does not converge (lindep=1), $\rho$ is set to 0. For stepwise QR decomposition, $\rho$ contains the diagonal element of the (n+1)^th column of R.

lindep

returns a value of 1 if the Gram-Schmidt process does not converge in 10 iterations. In most cases, if lindep=1, the input vector v is linearly dependent on the n columns of the input matrix Q. In that case, $\rho$ is set to 0, and the results w and r contain missing values. If lindep=0, the Gram-Schmidt process did converge, and the results w, r, and $\rho$ are computed.

The inputs to the ORTVEC subroutine are as follows:

v: specifies an m ×1 vector v that is to be orthogonalized to the n columns of Q. For stepwise QR decomposition of a matrix, v is the (n+1)^th matrix column before its orthogonalization.
q: specifies an optional m ×n matrix Q that is assumed to have $n \leq m$ (nearly) orthonormal columns. Thus, the n ×n matrix Q' Q should approximate the identity matrix. The column orthonormality assumption is not tested in the ORTVEC call. If it is violated, the results are not predictable. The argument q can be omitted or can have zero rows and columns. For stepwise QR decomposition of a matrix, q contains the first n matrix columns that are already orthogonal.

The relevant formula for the ORTVEC subroutine is

$v= {Qr} + \rho w$

Assuming that the m ×n matrix Q has n (nearly) orthonormal columns, the ORTVEC subroutine orthogonalizes the vector v to the columns of Q. The vector r is the array of Fourier coefficients, and $\rho$ is the distance from w to the range of Q.

There are two special cases:

If m > n, ORTVEC normalizes the result w, so that w' w = 1.
If m = n, the output vector w is the null vector.

The case m < n is not possible since Q is assumed to have n (nearly) orthonormal columns.

To initialize a stepwise QR decomposition, ORTVEC can be called to normalize v only, that is, to compute $w= v/ \sqrt{v^' v}$ and $\rho = \sqrt{v^' v}$ only. There are two ways of using the ORTVEC call for this reason:

Omit the last argument q, as in call ortvec(w,r,rho,lindep,v);.
Provide a matrix Q with zero rows and columns, for example, by using the free q; command.

In both cases, r is a column vector with zero rows.

The ORTVEC subroutine is useful for the following applications:

performing stepwise QR decomposition. Compute Q and R, so that A = QR, where Q is column orthonormal, Q' Q = I, and R is upper triangular. The j^th step is applied to the j^th column, v, of A, and it computes the j^th column w of Q and the j^th column, $(r\; \rho \; 0 )^'$ , of R.
computing the m ×(m-n) null space matrix, Q₂, corresponding to an m ×n range space matrix, Q₁ (m>n), by the following stepwise process: set v = e_i (where e_i is the i^th unit vector) and try to make it orthogonal to all column vectors of Q₁ and the already generated Q₂, if the subroutine is successful, append w to Q₂; otherwise, try v = e_i+1.

The 4 ×3 matrix Q contains the unit vectors e₁, e₃, and e₄. The column vector v is pairwise linearly independent with the three columns of Q. As expected, the ORTVEC call computes the vector w as the unit vector e₂ with ${\mu}= (1,1,1)$ and $\rho = 1$ .

proc iml;
   q = { 1  0  0,
         0  0  0,
         0  1  0,
         0  0  1 };
   v = { 1, 1, 1, 1 };
   call ortvec(w,u,rho,lindep,v,q);
   print rho u w;

You can perform the QR decomposition of the linearly independent columns of an m ×n matrix A with the following statements:

proc iml;
   a = { . . . enter matrix A here . . . };
   nind = 0;  ndep = 0;  dmax = 0.;
   n = ncol(a);  m = nrow(a);
   free q;
   do j = 1 to n;
      v = a[ ,j];
      call ORTVEC(w,u,rho,lindep,v,q);
      aro = abs(rho);
      if aro > dmax then dmax = aro;
      if aro <= 1.e-10 * dmax then lindep = 1;
      if lindep = 0 then do;
         nind = nind + 1;
         q = q || w;
         if nind = n then r = r || (u // rho);
         else r = r || (u // rho // j(n-nind,1,0.));
      end;
      else do;
         print "Column " j " is linearly dependent.";
         ndep = ndep + 1; ind[ndep] = j;
      end;
   end;

Next, process the remaining columns of A:

   do j = 1 to ndep;
      k = ind[ndep-j+1];
      v = a[ ,k];
      call ORTVEC(w,u,rho,lindep,v,q);
      if lindep = 0 then do;
         nind = nind + 1;
         q = q || w;
         if nind = n then r = r || (u // rho);
         else r = r || (u // rho // j(n-nind,1,0.));
      end;
   end;

Now compute the null space in the last columns of Q:

   do i = 1 to m;
      if nind < m then do;
        v = j(m,1,0.); v[i] = 1.;
        call ORTVEC(w,u,rho,lindep,v,q);
        aro = abs(rho);
        if aro > dmax then dmax = aro;
        if aro <= 1.e-10 * dmax then lindep = 1;
        if lindep = 0 then do;
          nind = nind + 1;
          q = q || w;
        end;
        else print "Unit vector" i "linearly dependent.";
      end;
   end;
   if nind < m then do;
      print "This is theoretically not possible.";
   end;

Chapter Contents
Previous
Next
Top