The ORVEC subroutine provides columnwise orthogonalization and stepwise QR decomposition by using the Gram-Schmidt process.

The ORTVEC subroutine returns the following values:

w

is an vector. If the Gram-Schmidt process converges (lindep=0), w is orthonormal to the columns of , which is assumed to have (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 th orthogonal column of the matrix . If the q argument is not specified, w is the normalized value of the vector v,

r

is a vector. If the Gram-Schmidt process converges (lindep=0), r contains 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 upper triangular elements of the th column of .

rho

is a scalar value. If the Gram-Schmidt process converges (lindep=0), rho specifies the distance from w to the range of . Even if the Gram-Schmidt process converges, if rho is sufficiently small, the vector v can be linearly dependent on the columns of . 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 th column of . In formulas, the value rho is denoted by .

lindep

returns a value of 1 if the Gram-Schmidt process does not converge in 10 iterations. A value of 1 often indicates that the input vector v is linearly dependent on the columns of the input matrix . 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 input arguments to the ORTVEC subroutine are as follows:

v

specifies an vector v that is to be orthogonalized to the columns of . For stepwise QR decomposition of a matrix, v is the th matrix column before its orthogonalization.

q

specifies an optional matrix that is assumed to have (nearly) orthonormal columns. Thus, the matrix 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 matrix columns that are already orthogonal.

The relevant formula for the ORTVEC subroutine is

In the formula, if the matrix has (nearly) orthonormal columns, the vector v is orthogonal to the columns of and is the distance from w to the range of .

There are two special cases:

• If , ORTVEC normalizes the result w, so that ww.

• If , the output vector w is the null vector.

The case is not possible since is assumed to have (nearly) orthonormal columns.

To initialize a stepwise QR decomposition, the ORTVEC subroutine can be called to normalize v only (that is, to compute and ). There are two ways of using the ORTVEC subroutine 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 and , so that , where is column orthonormal, , and is upper triangular. The th step is applied to the th column, v, of , and it computes the th column w of and the th column, , of .

• computing the null space matrix, , that corresponds to an range space matrix, , by the following stepwise process:

  1. Set (where is the th unit vector) and try to make it orthogonal to all column vectors of and the already generated .

  2. If the subroutine is successful, append w to ; otherwise, try .

The matrix contains the unit vectors , and . The column vector v is pairwise linearly independent with the three columns of . As expected, the ORTVEC subroutine computes the vector w as the unit vector with and .

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;

a = {1  2  1,
     2  4  2,
     1  4 -1,
     1  0  3}; /* use any matrix A */
nind = 0;  ndep = 0;  dmax = 0.;
n = ncol(a);  m = nrow(a); ind = j(1,n,0);
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;
print q r;

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;
print q r;

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;
print q;