hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_lapack_dgeqpf (f08be)

Purpose

nag_lapack_dgeqpf (f08be) computes the QRQR factorization, with column pivoting, of a real mm by nn matrix.

Syntax

[a, jpvt, tau, info] = f08be(a, jpvt, 'm', m, 'n', n)
[a, jpvt, tau, info] = nag_lapack_dgeqpf(a, jpvt, 'm', m, 'n', n)

Description

nag_lapack_dgeqpf (f08be) forms the QRQR factorization, with column pivoting, of an arbitrary rectangular real mm by nn matrix.
If mnmn, the factorization is given by:
AP = Q
(R)
0
,
AP= Q R 0 ,
where RR is an nn by nn upper triangular matrix, QQ is an mm by mm orthogonal matrix and PP is an nn by nn permutation matrix. It is sometimes more convenient to write the factorization as
AP =
(Q1Q2)
(R)
0
,
AP= Q1 Q2 R 0 ,
which reduces to
AP = Q1 R ,
AP= Q1 R ,
where Q1Q1 consists of the first nn columns of QQ, and Q2Q2 the remaining mnm-n columns.
If m < nm<n, RR is trapezoidal, and the factorization can be written
AP = Q
(R1R2)
,
AP= Q R1 R2 ,
where R1R1 is upper triangular and R2R2 is rectangular.
The matrix QQ is not formed explicitly but is represented as a product of min (m,n)min(m,n) elementary reflectors (see the F08 Chapter Introduction for details). Functions are provided to work with QQ in this representation (see Section [Further Comments]).
Note also that for any k < nk<n, the information returned in the first kk columns of the array a represents a QRQR factorization of the first kk columns of the permuted matrix APAP.
The function allows specified columns of AA to be moved to the leading columns of APAP at the start of the factorization and fixed there. The remaining columns are free to be interchanged so that at the iith stage the pivot column is chosen to be the column which maximizes the 22-norm of elements ii to mm over columns ii to nn.

References

Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore

Parameters

Compulsory Input Parameters

1:     a(lda, : :) – double array
The first dimension of the array a must be at least max (1,m)max(1,m)
The second dimension of the array must be at least max (1,n)max(1,n)
The mm by nn matrix AA.
2:     jpvt( : :) – int64int32nag_int array
Note: the dimension of the array jpvt must be at least max (1,n)max(1,n).
If jpvt(i)0jpvti0, then the ii th column of AA is moved to the beginning of APAP before the decomposition is computed and is fixed in place during the computation. Otherwise, the ii th column of AA is a free column (i.e., one which may be interchanged during the computation with any other free column).

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the array a.
mm, the number of rows of the matrix AA.
Constraint: m0m0.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
nn, the number of columns of the matrix AA.
Constraint: n0n0.

Input Parameters Omitted from the MATLAB Interface

lda work

Output Parameters

1:     a(lda, : :) – double array
The first dimension of the array a will be max (1,m)max(1,m)
The second dimension of the array will be max (1,n)max(1,n)
ldamax (1,m)ldamax(1,m).
If mnmn, the elements below the diagonal store details of the orthogonal matrix QQ and the upper triangle stores the corresponding elements of the nn by nn upper triangular matrix RR.
If m < nm<n, the strictly lower triangular part stores details of the orthogonal matrix QQ and the remaining elements store the corresponding elements of the mm by nn upper trapezoidal matrix RR.
2:     jpvt( : :) – int64int32nag_int array
Note: the dimension of the array jpvt must be at least max (1,n)max(1,n).
Details of the permutation matrix PP. More precisely, if jpvt(i) = kjpvti=k, then the kkth column of AA is moved to become the ii th column of APAP; in other words, the columns of APAP are the columns of AA in the order jpvt(1),jpvt(2),,jpvt(n)jpvt1,jpvt2,,jpvtn.
3:     tau(min (m,n)min(m,n)) – double array
Further details of the orthogonal matrix QQ.
4:     info – int64int32nag_int scalar
info = 0info=0 unless the function detects an error (see Section [Error Indicators and Warnings]).

Error Indicators and Warnings

  info = iinfo=-i
If info = iinfo=-i, parameter ii had an illegal value on entry. The parameters are numbered as follows:
1: m, 2: n, 3: a, 4: lda, 5: jpvt, 6: tau, 7: work, 8: info.
It is possible that info refers to a parameter that is omitted from the MATLAB interface. This usually indicates that an error in one of the other input parameters has caused an incorrect value to be inferred.

Accuracy

The computed factorization is the exact factorization of a nearby matrix (A + E)(A+E), where
E2 = O(ε) A2 ,
E2 = O(ε) A2 ,
and εε is the machine precision.

Further Comments

The total number of floating point operations is approximately (2/3) n2 (3mn) 23 n2 (3m-n)  if mnmn or (2/3) m2 (3nm) 23 m2 (3n-m)  if m < nm<n.
To form the orthogonal matrix QQ nag_lapack_dgeqpf (f08be) may be followed by a call to nag_lapack_dorgqr (f08af):
[a, info] = f08af(a(:,1:m), tau);
but note that the second dimension of the array a must be at least m, which may be larger than was required by nag_lapack_dgeqpf (f08be).
When mnmn, it is often only the first nn columns of QQ that are required, and they may be formed by the call:
[a, info] = f08af(a, tau);
To apply QQ to an arbitrary real rectangular matrix CC, nag_lapack_dgeqpf (f08be) may be followed by a call to nag_lapack_dormqr (f08ag). For example,
[c, info] = f08ag('Left','Transpose', a(:,1:min(m,n)), tau, c);
forms C = QTCC=QTC, where CC is mm by pp.
To compute a QRQR factorization without column pivoting, use nag_lapack_dgeqrf (f08ae).
The complex analogue of this function is nag_lapack_zgeqpf (f08bs).

Example

function nag_lapack_dgeqpf_example
a = [-0.09,  0.14, -0.46,  0.68,  1.29;
     -1.56,  0.2,   0.29,  1.09,  0.51;
     -1.48, -0.43,  0.89, -0.71, -0.96;
     -1.09,  0.84,  0.77,  2.11, -1.27;
      0.08,  0.55, -1.13,  0.14,  1.74;
     -1.59, -0.72,  1.06,  1.24,  0.34];
b = [-0.01, -0.04;
      0.04, -0.03;
      0.05,  0.01;
     -0.03, -0.02;
      0.02,  0.05;
     -0.06,  0.07];
jpvt = [int64(0);0;0;0;0];
% Compute the QR factorization of a
[a, jpvt, tau, info] = nag_lapack_dgeqpf(a, jpvt);

% Choose tol to reflect the relative accuracy of the input data
tol = 0.01;

% Determine which columns of R to use
k = find(abs(diag(a)) <= tol*abs(a(1,1)));
if numel(k) == 0
  k = numel(diag(a));
else
  k = k(1)-1;
end

% Compute c = (q^t)*b,
[c, info] = nag_lapack_dormqr('Left', 'Transpose', a, tau, b);

% Compute least-squares solution by backsubstitution in r*b = c
c(1:k, :) = inv(triu(a(1:k,1:k)))*c(1:k,:);
c(k+1:5, :) = 0;

% Unscramble the least-squares solution stored in c
x = zeros(5, 2);
for i=1:5
  x(jpvt(i), :) = c(i, :);
end

fprintf('\nLeast-squares solution\n');
disp(x);
 

Least-squares solution
   -0.0370   -0.0044
    0.0647   -0.0335
         0         0
   -0.0515    0.0018
    0.0066    0.0102


function f08be_example
a = [-0.09,  0.14, -0.46,  0.68,  1.29;
     -1.56,  0.2,   0.29,  1.09,  0.51;
     -1.48, -0.43,  0.89, -0.71, -0.96;
     -1.09,  0.84,  0.77,  2.11, -1.27;
      0.08,  0.55, -1.13,  0.14,  1.74;
     -1.59, -0.72,  1.06,  1.24,  0.34];
b = [-0.01, -0.04;
      0.04, -0.03;
      0.05,  0.01;
     -0.03, -0.02;
      0.02,  0.05;
     -0.06,  0.07];
jpvt = [int64(0);0;0;0;0];
% Compute the QR factorization of a
[a, jpvt, tau, info] = f08be(a, jpvt);

% Choose tol to reflect the relative accuracy of the input data
tol = 0.01;

% Determine which columns of R to use
k = find(abs(diag(a)) <= tol*abs(a(1,1)));
if numel(k) == 0
  k = numel(diag(a));
else
  k = k(1)-1;
end

% Compute c = (q^t)*b,
[c, info] = f08ag('Left', 'Transpose', a, tau, b);

% Compute least-squares solution by backsubstitution in r*b = c
c(1:k, :) = inv(triu(a(1:k,1:k)))*c(1:k,:);
c(k+1:5, :) = 0;

% Unscramble the least-squares solution stored in c
x = zeros(5, 2);
for i=1:5
  x(jpvt(i), :) = c(i, :);
end

fprintf('\nLeast-squares solution\n');
disp(x);
 

Least-squares solution
   -0.0370   -0.0044
    0.0647   -0.0335
         0         0
   -0.0515    0.0018
    0.0066    0.0102



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013