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)

Purpose

nag_linsys_real_gen_solve (f04jg) finds the solution of a linear least squares problem, Ax = bAx=b, where AA is a real mm by n(mn)n(mn) matrix and bb is an mm 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 = bAx=b is the vector xx of minimum (Euclidean) length which minimizes the length of the residual vector r = bAxr=b-Ax.
The real mm by n(mn)n(mn) matrix AA is factorized as
A = Q
(U)
0
A=Q U 0
where QQ is an mm by mm orthogonal matrix and UU is an nn by nn upper triangular matrix. If UU is of full rank, then the least squares solution is given by
x = (U10)QTb.
x=(U-10)QTb.
If UU is not of full rank, then the singular value decomposition of UU is obtained so that UU is factorized as
U = RDPT,
U=RDPT,
where RR and PP are nn by nn orthogonal matrices and DD is the nn by nn diagonal matrix
D = diag(σ1,σ2,,σn),
D=diag(σ1,σ2,,σn),
with σ1σ2σn0σ1σ2σn0, these being the singular values of AA. If the singular values σk + 1,,σnσk+1,,σn are negligible, but σkσk is not negligible, relative to the data errors in AA, then the rank of AA is taken to be kk and the minimal least squares solution is given by
x = P
(S10)
0 0
(RT0)
0 I
QTb,
x=P S-1 0 0 0 RT 0 0 I QTb,
where S = diag(σ1,σ2,,σk)S=diag(σ1,σ2,,σk).
This function obtains the factorizations by a call to nag_eigen_real_gen_qu_svd (f02wd).
The function also returns the value of the standard error
σ = sqrt((rTr)/(mk)), if ​m > k,
= 0, if ​m = k,
σ = rTr m-k , if ​m>k, = 0, if ​m=k,
rTrrTr being the residual sum of squares and kk the rank of AA.

References

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

Parameters

Compulsory Input Parameters

1:     a(lda,n) – double array
lda, the first dimension of the array, must satisfy the constraint ldamldam.
The mm by nn matrix AA.
2:     b(m) – double array
m, the dimension of the array, must satisfy the constraint mnmn.
The right-hand side vector bb.
3:     tol – double scalar
A relative tolerance to be used to determine the rank of AA. tol should be chosen as approximately the largest relative error in the elements of AA. For example, if the elements of AA are correct to about 44 significant figures then tol should be set to about 5 × 1045×10-4. See Section [Further Comments] for a description of how tol is used to determine rank. If tol is outside the range (ε,1.0)(ε,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 as declared in the (sub)program from which nag_linsys_real_gen_solve (f04jg) is called.
Constraint: lwork4 × nlwork4×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.)
mm, the number of rows of a.
Constraint: mnmn.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
nn, the number of columns of aa.
Constraint: 1nm1nm.

Input Parameters Omitted from the MATLAB Interface

lda

Output Parameters

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

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,n < 1n<1,
orm < nm<n,
orlda < mlda<m,
orlwork < 4 × nlwork<4×n.
  ifail = 2ifail=2
The QRQR algorithm has failed to converge to the singular values in 50 × n50×n iterations. This failure can only happen when the singular value decomposition is employed, but even then it is not likely to occur.

Accuracy

The computed factors QQ, UU, RR, DD and PTPT satisfy the relations
Q
(U)
0
= A + E,Q
(R0)
0 I
(D)
0
PT = A + F,
Q U 0 =A+E,Q R 0 0 I D 0 PT=A+F,
where
E2c1ε A2,
E2c1ε A2,
F2c2εA2,
F2c2εA2,
εε being the machine precision, and c1c1 and c2c2 being modest functions of mm and nn. Note that A2 = σ1A2=σ1.
For a fuller discussion, covering the accuracy of the solution xx see Lawson and Hanson (1974), especially pages 50 and 95.

Further Comments

If the least squares solution is obtained from the QUQU factorization then the time taken by the function is approximately proportional to n2(3mn)n2(3m-n). If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to n2(3m + 19n)n2(3m+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 QUQU factorization of AA the condition number
c(U) = UE U1E
c(U)=UE U-1E
is determined and if c(U)c(U) is such that
c(U) × tol > 1.0
c(U)×tol>1.0
then UU is regarded as singular and the singular values of AA are computed. If this test is not satisfied, UU is regarded as nonsingular and the rank of AA is set to nn. When the singular values are computed the rank of AA, say kk, is returned as the largest integer such that
σk > tol × σ1,
σk>tol×σ1,
unless σ1 = 0σ1=0 in which case kk is returned as zero. That is, singular values which satisfy σitol × σ1σitol×σ1 are regarded as negligible because relative perturbations of order tol can make such singular values zero.

Example

function nag_linsys_real_gen_solve_example
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);
[aOut, bOut, svd, sigma, irank, work, ifail] = nag_linsys_real_gen_solve(a, b, tol, lwork)
 

aOut =

   -0.5000   -0.5000   -0.5000    0.5000
    0.5000    0.5000   -0.5000    0.5000
   -0.5000    0.5000   -0.5000   -0.5000
    0.5000   -0.5000   -0.5000   -0.5000
    0.1604   -0.5793   -0.0473    0.1714
    0.2138   -0.7724   -0.0630    0.2228


bOut =

    4.9667
   -2.8333
    4.5667
    3.2333
    1.1258
    0.8292


svd =

     1


sigma =

    0.9092


irank =

                    3


work =

    3.0000
    2.0000
    1.0000
    0.0000
    1.0000
   -1.7333
    0.4000
   -7.8000
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0


function f04jg_example
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);
[aOut, bOut, svd, sigma, irank, work, ifail] = f04jg(a, b, tol, lwork)
 

aOut =

   -0.5000   -0.5000   -0.5000    0.5000
    0.5000    0.5000   -0.5000    0.5000
   -0.5000    0.5000   -0.5000   -0.5000
    0.5000   -0.5000   -0.5000   -0.5000
    0.1604   -0.5793   -0.0473    0.1714
    0.2138   -0.7724   -0.0630    0.2228


bOut =

    4.9667
   -2.8333
    4.5667
    3.2333
    1.1258
    0.8292


svd =

     1


sigma =

    0.9092


irank =

                    3


work =

    3.0000
    2.0000
    1.0000
    0.0000
    1.0000
   -1.7333
    0.4000
   -7.8000
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0



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