hide long namesshow long 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

3 Description

4 References

▸▿ 5 Parameters

5.1 Compulsory Input Parameters

5.2 Optional Input Parameters

5.3 Output Parameters

6 Error Indicators and Warnings

7 Accuracy

8 Further Comments

9 Example

Purpose

nag_linsys_real_gen_solve (f04jg) finds the solution of a linear least squares problem,

A x = b

, where

A

is a real

m

n (m \geq n)

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

A x = b

is the vector

x

of minimum (Euclidean) length which minimizes the length of the residual vector

r = b - A x

The real

m

n (m \geq n)

matrix

A

is factorized as

A = Q (\begin{array}{l} R \\ 0 \end{array})

where

Q

is an

m

m

orthogonal matrix and

R

is an

n

n

upper triangular matrix. If

R

is of full rank, then the least squares solution is given by

x = (\begin{array}{l} R^{- 1} & 0 \end{array}) Q^{T} b .

R

is not of full rank, then the singular value decomposition of

R

is obtained so that

R

is factorized as

R = U D V^{T},

where

U

and

V

are

n

n

orthogonal matrices and

D

is the

n

n

diagonal matrix

D = diag (σ_{1}, σ_{2}, \dots, σ_{n}),

with

σ_{1} \geq σ_{2} \geq \dots \geq σ_{n} \geq 0

, these being the singular values of

A

. If the singular values

σ_{k + 1}, \dots, σ_{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 (\begin{array}{l} S^{- 1} & 0 \\ 0 & 0 \end{array}) (\begin{array}{l} U^{T} & 0 \\ 0 & I \end{array}) Q^{T} b,

where

S = diag (σ_{1}, σ_{2}, \dots, σ_{k})

The function also returns the value of the standard error

\begin{array}{ccl} σ & = & \sqrt{\frac{r^{T} r}{m - k}}, & if ​ m > k, \\ = & 0, & if ​ m = k, \end{array}

r^{T} r

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: $a (lda, n)$ – double array: lda, the first dimension of the array, must satisfy the constraint $lda \geq m$ .
The $m$ by $n$ matrix $A$ .
2: $b (m)$ – 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 \times 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: $lwork \geq 4 \times 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: $m \geq n$ .
2: $n$ – int64int32nag_int scalar: Default: the second dimension of the array a.
$n$ , the number of columns of $a$ .

Constraint: $1 \leq n \leq m$ .

Output Parameters

1: $a (lda, n)$ – double array: If svd is returned as false, a stores details of the $Q R$ 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: $b (m)$ – 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 $Q R$ 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 $\sqrt{r^{T} r / (m - k)}$ when $m > k$ , and the value zero when $m = k$ . Here $r$ is the residual vector $b - A x$ 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: $work (lwork)$ – double array: If svd is returned as false, then the first $n$ elements of work contain information on the $Q R$ factorization of $A$ (see argument a above), and $work (n + 1)$ contains the condition number ${‖R‖}_{E} {‖R^{- 1}‖}_{E}$ 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 $work (n + 1)$ contains the total number of iterations taken by the $Q R$ 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: $lda \geq m$ .

Constraint: $m \geq n$ .

Constraint: $n \geq 1$ .

On entry, lwork is too small.

$ifail = 2$: The $Q R$ algorithm has failed to converge to the singular values in $50 \times 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

V^{T}

satisfy the relations

Q (\begin{array}{l} R \\ 0 \end{array}) = A + E,

Q (\begin{array}{l} U & 0 \\ 0 & I \end{array}) (\begin{array}{l} D \\ 0 \end{array}) V^{T} = A + F,

where

{‖E‖}_{2} \leq c_{1} ε {‖A‖}_{2},

{‖F‖}_{2} \leq c_{2} ε {‖A‖}_{2},

ε

being the machine precision, and

c_{1}

and

c_{2}

being modest functions of

m

and

n

. Note that

{‖A‖}_{2} = σ_{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

Q R

factorization then the time taken by the function is approximately proportional to

n^{2} (3 m - n)

. If the least squares solution is obtained from the singular value decomposition then the time taken is approximately proportional to

n^{2} (3 m + 19 n)

. 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

Q R

factorization of

A

the condition number

c (R) = {‖R‖}_{E} {‖R^{- 1}‖}_{E}

is determined and if

c (R)

is such that

c (R) \times 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 \times σ_{1},

unless

σ_{1} = 0

in which case

k

is returned as zero. That is, singular values which satisfy

σ_{i} \leq tol \times σ_{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 - A x

, where

A = (\begin{array}{r} 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 \end{array}), b = (\begin{array}{l} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \end{array})

and the value tol is to be taken as

5 \times 10^{- 4}

Open in the MATLAB editor: f04jg_example

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