Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_lapack_dgelsy (f08ba)

## Purpose

nag_lapack_dgelsy (f08ba) computes the minimum norm solution to a real linear least squares problem
 min ‖b − Ax‖2 x
$minx ‖b-Ax‖2$
using a complete orthogonal factorization of A$A$. A$A$ is an m$m$ by n$n$ matrix which may be rank-deficient. Several right-hand side vectors b$b$ and solution vectors x$x$ can be handled in a single call.

## Syntax

[a, b, jpvt, rank, info] = f08ba(a, b, jpvt, rcond, 'm', m, 'n', n, 'nrhs_p', nrhs_p)
[a, b, jpvt, rank, info] = nag_lapack_dgelsy(a, b, jpvt, rcond, 'm', m, 'n', n, 'nrhs_p', nrhs_p)

## Description

The right-hand side vectors are stored as the columns of the m$m$ by r$r$ matrix B$B$ and the solution vectors in the n$n$ by r$r$ matrix X$X$.
nag_lapack_dgelsy (f08ba) first computes a QR$QR$ factorization with column pivoting
AP = Q
 ( R11 R12 ) 0 R22
,
$AP= Q R11 R12 0 R22 ,$
with R11${R}_{11}$ defined as the largest leading sub-matrix whose estimated condition number is less than 1 / rcond$1/{\mathbf{rcond}}$. The order of R11${R}_{11}$, rank, is the effective rank of A$A$.
Then, R22${R}_{22}$ is considered to be negligible, and R12${R}_{12}$ is annihilated by orthogonal transformations from the right, arriving at the complete orthogonal factorization
AP = Q
 ( T11 0 ) 0 0
Z .
$AP= Q T11 0 0 0 Z .$
The minimum norm solution is then
X = PZT
 ( T11 − 1 Q1T b ) 0
$X = PZT T11-1 Q1T b 0$
where Q1${Q}_{1}$ consists of the first rank columns of Q$Q$.

## References

Anderson E, Bai Z, Bischof C, Blackford S, Demmel J, Dongarra J J, Du Croz J J, Greenbaum A, Hammarling S, McKenney A and Sorensen D (1999) LAPACK Users' Guide (3rd Edition) SIAM, Philadelphia http://www.netlib.org/lapack/lug
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)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
The second dimension of the array must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
The m$m$ by n$n$ matrix A$A$.
2:     b(ldb, : $:$) – double array
The first dimension of the array b must be at least max (1,m,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}},{\mathbf{n}}\right)$
The second dimension of the array must be at least max (1,nrhs)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs}}\right)$
The m$m$ by r$r$ right-hand side matrix B$B$.
3:     jpvt( : $:$) – int64int32nag_int array
Note: the dimension of the array jpvt must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If jpvt(i)0${\mathbf{jpvt}}\left(i\right)\ne 0$, the i$i$th column of A$A$ is permuted to the front of AP$AP$, otherwise column i$i$ is a free column.
4:     rcond – double scalar
Used to determine the effective rank of A$A$, which is defined as the order of the largest leading triangular sub-matrix R11${R}_{11}$ in the QR$QR$ factorization of A$A$, whose estimated condition number is < 1 / rcond$\text{}<1/{\mathbf{rcond}}$.
Suggested value: if the condition number of a is not known then rcond = sqrt((ε) / 2)${\mathbf{rcond}}=\sqrt{\left(\epsilon \right)/2}$ (where ε$\epsilon$ is machine precision, see nag_machine_precision (x02aj)) is a good choice. Negative values or values less than machine precision should be avoided since this will cause a to have an effective rank = min (m,n)$\text{rank}=\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{m}},{\mathbf{n}}\right)$ that could be larger than its actual rank, leading to meaningless results.

### Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The first dimension of the array a.
m$m$, the number of rows of the matrix A$A$.
Constraint: m0${\mathbf{m}}\ge 0$.
2:     n – int64int32nag_int scalar
Default: The second dimension of the array a.
n$n$, the number of columns of the matrix A$A$.
Constraint: n0${\mathbf{n}}\ge 0$.
3:     nrhs_p – int64int32nag_int scalar
Default: The second dimension of the array b.
r$r$, the number of right-hand sides, i.e., the number of columns of the matrices B$B$ and X$X$.
Constraint: nrhs0${\mathbf{nrhs}}\ge 0$.

### Input Parameters Omitted from the MATLAB Interface

lda ldb work lwork

### Output Parameters

1:     a(lda, : $:$) – double array
The first dimension of the array a will be max (1,m)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$
The second dimension of the array will be max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$
ldamax (1,m)$\mathit{lda}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}}\right)$.
a stores details of its complete orthogonal factorization.
2:     b(ldb, : $:$) – double array
The first dimension of the array b will be max (1,m,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}},{\mathbf{n}}\right)$
The second dimension of the array will be max (1,nrhs)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{nrhs}}\right)$
ldbmax (1,m,n)$\mathit{ldb}\ge \mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{m}},{\mathbf{n}}\right)$.
The n$n$ by r$r$ solution matrix X$X$.
3:     jpvt( : $:$) – int64int32nag_int array
Note: the dimension of the array jpvt must be at least max (1,n)$\mathrm{max}\phantom{\rule{0.125em}{0ex}}\left(1,{\mathbf{n}}\right)$.
If jpvt(i) = k${\mathbf{jpvt}}\left(i\right)=k$, then the i$i$th column of AP$AP$ was the k$k$th column of A$A$.
4:     rank – int64int32nag_int scalar
The effective rank of A$A$, i.e., the order of the sub-matrix R11${R}_{11}$. This is the same as the order of the sub-matrix T11${T}_{11}$ in the complete orthogonal factorization of A$A$.
5:     info – int64int32nag_int scalar
info = 0${\mathbf{info}}=0$ unless the function detects an error (see Section [Error Indicators and Warnings]).

## Error Indicators and Warnings

info = i${\mathbf{info}}=-i$
If info = i${\mathbf{info}}=-i$, parameter i$i$ had an illegal value on entry. The parameters are numbered as follows:
1: m, 2: n, 3: nrhs_p, 4: a, 5: lda, 6: b, 7: ldb, 8: jpvt, 9: rcond, 10: rank, 11: work, 12: lwork, 13: 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

See Section 4.5 of Anderson et al. (1999) for details of error bounds.

The complex analogue of this function is nag_lapack_zgelsy (f08bn).

## Example

```function nag_lapack_dgelsy_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 = [7.4;
4.2;
-8.3;
1.8;
8.6;
2.1];
jpvt = [int64(0);0;0;0;0];
rcond = 0.01;
[aOut, bOut, jpvtOut, rank, info] = nag_lapack_dgelsy(a, b, jpvt, rcond)
```
```

aOut =

-2.9366   -0.7273    1.9134   -1.0317   -0.0890
0.5234    2.8681    0.1721    0.7329   -0.1692
0.4966   -0.4778   -2.3028   -0.6896    0.1052
0.3657   -0.5520   -0.3306    1.2344   -0.3173
-0.0268    0.6259   -0.0197    0.1934   -0.0034
0.5335   -0.0259    0.0087   -0.3961   -0.6352

bOut =

0.6344
0.9699
-1.4402
3.3678
3.3992
-0.0035

jpvtOut =

1
5
4
2
3

rank =

4

info =

0

```
```function f08ba_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 = [7.4;
4.2;
-8.3;
1.8;
8.6;
2.1];
jpvt = [int64(0);0;0;0;0];
rcond = 0.01;
[aOut, bOut, jpvtOut, rank, info] = f08ba(a, b, jpvt, rcond)
```
```

aOut =

-2.9366   -0.7273    1.9134   -1.0317   -0.0890
0.5234    2.8681    0.1721    0.7329   -0.1692
0.4966   -0.4778   -2.3028   -0.6896    0.1052
0.3657   -0.5520   -0.3306    1.2344   -0.3173
-0.0268    0.6259   -0.0197    0.1934   -0.0034
0.5335   -0.0259    0.0087   -0.3961   -0.6352

bOut =

0.6344
0.9699
-1.4402
3.3678
3.3992
-0.0035

jpvtOut =

1
5
4
2
3

rank =

4

info =

0

```