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_linsys_real_gen_solve (f04jg)

 Contents

    1  Purpose
    2  Syntax
    7  Accuracy
    9  Example

Purpose

nag_linsys_real_gen_solve (f04jg) finds the solution of a linear least squares problem, Ax=b, where A is a real m by nmn matrix and b is an m element vector. If the matrix of observations is not of full rank, then the minimal least squares solution is returned.

Syntax

[a, b, svd, sigma, irank, work, ifail] = f04jg(a, b, tol, lwork, 'm', m, 'n', n)
[a, b, svd, sigma, irank, work, ifail] = nag_linsys_real_gen_solve(a, b, tol, lwork, 'm', m, 'n', n)

Description

The minimal least squares solution of the problem Ax=b is the vector x of minimum (Euclidean) length which minimizes the length of the residual vector r=b-Ax.
The real m by nmn matrix A is factorized as
A=Q R 0  
where Q is an m by m orthogonal matrix and R is an n by n upper triangular matrix. If R is of full rank, then the least squares solution is given by
x= R-1 0 QTb.  
If R is not of full rank, then the singular value decomposition of R is obtained so that R is factorized as
R=UDVT,  
where U and V are n by n orthogonal matrices and D is the n by n diagonal matrix
D=diagσ1,σ2,,σn,  
with σ1σ2σn0, these being the singular values of A. If the singular values σk+1,,σn are negligible, but σk is not negligible, relative to the data errors in A, then the rank of A is taken to be k and the minimal least squares solution is given by
x=V S-1 0 0 0 UT 0 0 I QTb,  
where S=diagσ1,σ2,,σk.
The function also returns the value of the standard error
σ = rTr m-k , if ​m>k, = 0, if ​m=k,  
rTr being the residual sum of squares and k the rank of A.

References

Lawson C L and Hanson R J (1974) Solving Least Squares Problems Prentice–Hall

Parameters

Compulsory Input Parameters

1:     aldan – double array
lda, the first dimension of the array, must satisfy the constraint ldam.
The m by n matrix A.
2:     bm – double array
The right-hand side vector b.
3:     tol – double scalar
A relative tolerance to be used to determine the rank of A. tol should be chosen as approximately the largest relative error in the elements of A. For example, if the elements of A are correct to about 4 significant figures then tol should be set to about 5×10-4. See Further Comments for a description of how tol is used to determine rank. If tol is outside the range ε,1.0, where ε is the machine precision, then the value ε is used in place of tol. For most problems this is unreasonably small.
4:     lwork int64int32nag_int scalar
The dimension of the array work.
Constraint: lwork4×n.

Optional Input Parameters

1:     m int64int32nag_int scalar
Default: the dimension of the array b and the first dimension of the array a. (An error is raised if these dimensions are not equal.)
m, the number of rows of a.
Constraint: mn.
2:     n int64int32nag_int scalar
Default: the second dimension of the array a.
n, the number of columns of a.
Constraint: 1nm.

Output Parameters

1:     aldan – double array
If svd is returned as false, a stores details of the QR factorization of A.
If svd is returned as true, the first n rows of a store the right-hand singular vectors, stored by rows; and the remaining rows of the array are used as workspace.
2:     bm – double array
The first n elements of b contain the minimal least squares solution vector x. The remaining m-n elements are used for workspace.
3:     svd – logical scalar
Is returned as false if the least squares solution has been obtained from the QR factorization of A. In this case A is of full rank. svd is returned as true if the least squares solution has been obtained from the singular value decomposition of A.
4:     sigma – double scalar
The standard error, i.e., the value rTr/m-k when m>k, and the value zero when m=k. Here r is the residual vector b-Ax and k is the rank of A.
5:     irank int64int32nag_int scalar
k, the rank of the matrix A. It should be noted that it is possible for irank to be returned as n and svd to be returned as true. This means that the matrix R only just failed the test for nonsingularity.
6:     worklwork – double array
If svd is returned as false, then the first n elements of work contain information on the QR factorization of A (see argument a above), and workn+1 contains the condition number RER-1E of the upper triangular matrix R.
If svd is returned as true, then the first n elements of work contain the singular values of A arranged in descending order and workn+1 contains the total number of iterations taken by the QR algorithm. The rest of work is used as workspace.
7:     ifail int64int32nag_int scalar
ifail=0 unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:
   ifail=1
Constraint: ldam.
Constraint: mn.
Constraint: n1.
On entry, lwork is too small.
   ifail=2
The QR algorithm has failed to converge to the singular values in 50×n iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.
   ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
   ifail=-399
Your licence key may have expired or may not have been installed correctly.
   ifail=-999
Dynamic memory allocation failed.

Accuracy

The computed factors Q, R, U, D and VT satisfy the relations
Q R 0 =A+E,  
Q U 0 0 I D 0 VT=A+F,  
where
E2c1ε A2,  
F2c2εA2,  
ε being the machine precision, and c1 and c2 being modest functions of m and n. Note that A2=σ1.
For a fuller discussion, covering the accuracy of the solution x see Lawson and Hanson (1974), especially pages 50 and 95.

Further Comments

If the least squares solution is obtained from the QR factorization then the time taken by the function is approximately proportional to n23m-n. If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to n23m+19n. The approximate proportionality factor is the same in each case.
This function is column biased and so is suitable for use in paged environments.
Following the QR factorization of A the condition number
cR=RE R-1E  
is determined and if cR is such that
cR×tol>1.0  
then R is regarded as singular and the singular values of A are computed. If this test is not satisfied, R is regarded as nonsingular and the rank of A is set to n. When the singular values are computed the rank of A, say k, is returned as the largest integer such that
σk>tol×σ1,  
unless σ1=0 in which case k is returned as zero. That is, singular values which satisfy σitol×σ1 are regarded as negligible because relative perturbations of order tol can make such singular values zero.

Example

This example obtains a least squares solution for r=b-Ax, where
A= 0.05 0.05 0.25 -0.25 0.25 0.25 0.05 -0.05 0.35 0.35 1.75 -1.75 1.75 1.75 0.35 -0.35 0.30 -0.30 0.30 0.30 0.40 -0.40 0.40 0.40 ,  b= 1 2 3 4 5 6  
and the value tol is to be taken as 5×10-4.
function f04jg_example


fprintf('f04jg example results\n\n');

% Solve linear least squares problem Ax = b for general A
a = [0.05, 0.05, 0.25, -0.25;
     0.25, 0.25, 0.05, -0.05;
     0.35, 0.35, 1.75, -1.75;
     1.75, 1.75, 0.35, -0.35;
     0.3, -0.3,  0.3,   0.3;
     0.4, -0.4,  0.4,   0.4];
b = [1;    2;     3;     4;     5;     6];

tol = 0.0005;
lwork = int64(32);
[a, x, svd, sigma, irank, work, ifail] = ...
  f04jg(a, b, tol, lwork);

disp('Solution');
disp(x(1:4)');
fprintf('Standard error = %6.3f, rank = %2d\n\n', sigma, irank);
if svd==1
  fprintf('solution obtained from SVD of A\n');
else
  fprintf('solution obtained from QU factorization of A\n');
end


f04jg example results

Solution
    4.9667   -2.8333    4.5667    3.2333

Standard error =  0.909, rank =  3

solution obtained from SVD of A

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–2015