Language Reference: QR Call

QR Call

CALL QR( q, r, piv, lindep, a <, ord> <, b> ) ;

The QR subroutine produces the QR decomposition of a matrix by using Householder transformations.

The QR subroutine returns the following values:

q: specifies an orthogonal matrix $\text{[math]}$ that is the product of the Householder transformations applied to the $\text{[math]}$ matrix $\text{[math]}$ , if the b argument is not specified. In this case, the $\text{[math]}$ Householder transformations are applied, and q is an $\text{[math]}$ matrix. If the b argument is specified, q is the $\text{[math]}$ matrix $\text{[math]}$ that has the transposed Householder transformations $\text{[math]}$ applied on the $\text{[math]}$ columns of the argument matrix $\text{[math]}$ .
r: specifies a $\text{[math]}$ upper triangular matrix $\text{[math]}$ that is the upper part of the $\text{[math]}$ upper triangular matrix $\text{[math]}$ of the QR decomposition of the matrix $\text{[math]}$ . The matrix $\text{[math]}$ of the QR decomposition can be obtained by vertical concatenation (by using the operator //) of the $\text{[math]}$ zero matrix to the result matrix $\text{[math]}$ .
piv: specifies an $\text{[math]}$ vector of permutations of the columns of $\text{[math]}$ ; that is, on return, the QR decomposition is computed, not of $\text{[math]}$ , but of the permuted matrix whose columns are $\text{[math]}$ . The vector piv corresponds to an $\text{[math]}$ permutation matrix $\text{[math]}$ .
lindep: is the number of linearly dependent columns in matrix $\text{[math]}$ detected by applying the $\text{[math]}$ Householder transformations in the order specified by the argument vector piv.

The input arguments to the QR subroutine are as follows:

a

specifies an $\text{[math]}$ matrix $\text{[math]}$ that is to be decomposed into the product of the orthogonal matrix $\text{[math]}$ and the upper triangular matrix $\text{[math]}$ .

ord

specifies an optional $\text{[math]}$ vector that specifies the order of Householder transformations applied to matrix $\text{[math]}$ . When you specify the ord argument, the columns of $\text{[math]}$ can be divided into the following groups:

ord[j]>0: Column $\text{[math]}$ of $\text{[math]}$ is an initial column, meaning it has to be processed at the start in increasing order of ord[j]. This specification defines the first $\text{[math]}$ columns of $\text{[math]}$ that are to be processed.
ord[j]=0: Column $\text{[math]}$ of $\text{[math]}$ is a pivot column, meaning it is to be processed in order of decreasing residual Euclidean norms. The pivot columns of $\text{[math]}$ are processed after the $\text{[math]}$ initial columns and before the $\text{[math]}$ final columns.
ord[j]<0: Column $\text{[math]}$ of $\text{[math]}$ is a final column, meaning it has to be processed at the end in decreasing order of ord[j]. This specification defines the last $\text{[math]}$ columns of $\text{[math]}$ that are to be processed. If $\text{[math]}$ , some of these columns are not processed.

The default is ord[j]=j, in which case the Householder transformations are processed in the same order in which the columns are stored in matrix $\text{[math]}$ (without pivoting).

b

specifies an optional $\text{[math]}$ matrix $\text{[math]}$ that is to be multiplied by the transposed $\text{[math]}$ matrix $\text{[math]}$ . If b is specified, the result q contains the $\text{[math]}$ matrix $\text{[math]}$ . If b is not specified, the result q contains the $\text{[math]}$ matrix $\text{[math]}$ .

The QR subroutine decomposes an $\text{[math]}$ matrix $\text{[math]}$ into the product of an $\text{[math]}$ orthogonal matrix $\text{[math]}$ and an $\text{[math]}$ upper triangular matrix $\text{[math]}$ , so that

$\text{[math]}$

by means of $\text{[math]}$ Householder transformations.

The $\text{[math]}$ orthogonal matrix $\text{[math]}$ is computed only if the last argument b is not specified, as in the following example:

call qr(q, r, piv, lindep, a, ord);

In many applications, the number of rows, $\text{[math]}$ , is very large. In these cases, the explicit computation of the $\text{[math]}$ matrix $\text{[math]}$ might require too much memory or time.

In the usual case where $\text{[math]}$ ,

$\text{[math]}$	$\text{[math]}$	$\text{[math]}$
$\text{[math]}$	$\text{[math]}$	$\text{[math]}$
$\text{[math]}$	$\text{[math]}$	$\text{[math]}$

where $\text{[math]}$ is the result returned by the QR subroutine.

The $\text{[math]}$ columns of matrix $\text{[math]}$ provide an orthonormal basis for the $\text{[math]}$ columns of $\text{[math]}$ and are called the range space of $\text{[math]}$ . Since the $\text{[math]}$ columns of $\text{[math]}$ are orthogonal to the $\text{[math]}$ columns of $\text{[math]}$ , $\text{[math]}$ , they provide an orthonormal basis for the orthogonal complement of the columns of $\text{[math]}$ and are called the null space of $\text{[math]}$ .

In the case where $\text{[math]}$ ,

	$\text{[math]}$	$\text{[math]}$	$\text{[math]}$
	$\text{[math]}$

Specifying the argument ord as an $\text{[math]}$ vector lets you specify a special order of the columns in matrix $\text{[math]}$ on which the Householder transformations are applied. There are two special cases:

If you do not specify the ord argument, the default values ord $\text{[math]}$ are used. In this case, Householder transformations are done in the same order in which the columns are stored in $\text{[math]}$ (without pivoting).
If you set all components of ord to zero, the Householder transformations are done in order of decreasing Euclidean norms of the columns of $\text{[math]}$ .

To check the QR decomposition, use the following statements to compute the three residual sum of squares (represented by the variables SS0, SS1, and SS2), which should be close to zero:

a = shape(1:20, 5);
m = nrow(a); n = ncol(a);
ord = j(1, n, 0);
call qr(q, r, piv, lindep, a);
ss0 = ssq(a[ ,piv] - q[,1:n] * r);
ss1 = ssq(q * q` - i(m));
ss2 = ssq(q` * q - i(m));
print ss0 ss1 ss2;

Figure 23.238 Result of a QR Decomposition

ss0	ss1	ss2
6.231E-28	2.948E-31	2.862E-31

If the QR subroutine detects linearly dependent columns while processing matrix $\text{[math]}$ , the column order given in the result vector piv can differ from an explicitly specified order in the argument vector ord. If a column of $\text{[math]}$ is found to be linearly dependent on columns already processed, this column is swapped to the end of matrix $\text{[math]}$ . The order of columns in the result matrix $\text{[math]}$ corresponds to the order of columns processed in $\text{[math]}$ . The swapping of a linearly dependent column of $\text{[math]}$ to the end of the matrix corresponds to the swapping of the same column in $\text{[math]}$ and leads to a zero row at the end of the upper triangular matrix $\text{[math]}$ .

The scalar result lindep counts the number of linearly dependent columns that are detected in constructing the first $\text{[math]}$ Householder transformations in the order specified by the argument vector ord. The test of linear dependence depends on the singularity criterion, which is 1E $\text{[math]}$ 8 by default.

Solving the linear system $\text{[math]}$ with an upper triangular matrix $\text{[math]}$ whose columns are permuted corresponding to the result vector piv leads to a solution $\text{[math]}$ with permuted components. You can reorder the components of $\text{[math]}$ by using the index vector piv at the left-hand side of an expression, as follows:

a = {3  0  0 -1,
     0  1  2  0,
     4 -4 -1  1,
    -1  2  3  4};
b = {-1, 8, -3, 28};

n = ncol(a); p = ncol(b);
ord = j(1, n, 0);
call qr(qtb, r, piv, lindep, a, ord, b);
print piv;

x = j(n,1);
x[piv] = inv(r) * qtb[1:n, 1:p];
print x;

Figure 23.239 Solution to a Linear System

piv
1	4	2	3

x
1
2
3
4

The Full-Rank Linear Least Squares Problem

This example solves the full-rank linear least squares problem. Specify the argument b as an $\text{[math]}$ matrix $\text{[math]}$ , as follows:

call qr(q, r, piv, lindep, a, ord, b);

When you specify the b argument, the QR subroutine computes the matrix $\text{[math]}$ (instead of $\text{[math]}$ ) as the result q. Now you can compute the $\text{[math]}$ least squares solutions $\text{[math]}$ of an overdetermined linear system with an $\text{[math]}$ coefficient matrix $\text{[math]}$ , rank( $\text{[math]}$ ) = $\text{[math]}$ , and $\text{[math]}$ right-hand sides b $\text{[math]}$ stored as the columns of the $\text{[math]}$ matrix $\text{[math]}$ :

$\text{[math]}$

where $\text{[math]}$ is the Euclidean vector norm. This is accomplished by solving the $\text{[math]}$ upper triangular systems with back substitution:

$\text{[math]}$

For most applications, the number of rows of $\text{[math]}$ , $\text{[math]}$ , is much larger than $\text{[math]}$ , the number of columns of $\text{[math]}$ , or $\text{[math]}$ , the number of right-hand sides. In these cases, you are advised not to compute the large $\text{[math]}$ matrix $\text{[math]}$ (which can consume too much memory and time) if you can solve your problem by computing only the smaller $\text{[math]}$ matrix $\text{[math]}$ implicitly.

For example, use the first five columns of the $\text{[math]}$ Hilbert matrix $\text{[math]}$ , as follows:

a= {  36      -630      3360     -7560      7560     -2772,
    -630     14700    -88200    211680   -220500     83160,
    3360    -88200    564480  -1411200   1512000   -582120,
   -7560    211680  -1411200   3628800  -3969000   1552320,
    7560   -220500   1512000  -3969000   4410000  -1746360,
   -2772     83160   -582120   1552320  -1746360    698544 };
aa = a[, 1:5];
b= { 463, -13860, 97020, -258720, 291060, -116424};

m = nrow(aa); n = ncol(aa); p = ncol(b);
call qr(qtb, r, piv, lindep, aa, , b);

if lindep=0 then do;
   x=inv(r)*qtb[1:n];
   print x;  /* x solves aa*x=b */
end;
else /* handle linear dependence */;

Figure 23.240 Solution to Least Squares Problem

x
1
0.5
0.3333333
0.25
0.2

Note that you are using only the first $\text{[math]}$ rows, $\text{[math]}$ , of the qtb matrix. The IF-THEN statement of the preceding example can be replaced by the more efficient TRISOLV function:

if lindep=0 then
   x = trisolv(1, r, qtb[1:n], piv);

For information about solving rank-deficient linear least squares problems, see the RZLIND call.