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_opt_lsq_check_hessian (e04yb)

Purpose

nag_opt_lsq_check_hessian (e04yb) checks that a user-supplied function for evaluating the second derivative term of the Hessian matrix of a sum of squares is consistent with a user-supplied function for calculating the corresponding first derivatives.

Syntax

[fvec, fjac, b, user, ifail] = e04yb(m, lsqfun, lsqhes, x, lb, 'n', n, 'user', user)
[fvec, fjac, b, user, ifail] = nag_opt_lsq_check_hessian(m, lsqfun, lsqhes, x, lb, 'n', n, 'user', user)
Note: the interface to this routine has changed since earlier releases of the toolbox:
Mark 24: drop w and iw, introduce user in main routine and all call-backs
.

Description

Functions for minimizing a sum of squares of m$m$ nonlinear functions (or ‘residuals’), fi(x1,x2,,xn)${f}_{\mathit{i}}\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$ and mn$m\ge n$, may require you to supply a function to evaluate the quantities
 m bjk = ∑ fi(∂2fi)/( ∂ xj ∂ xk) i = 1
$bjk=∑i=1mfi ∂2fi ∂xj∂xk$
for j = 1,2,,n$j=1,2,\dots ,n$ and k = 1,2,,j$k=1,2,\dots ,j$. nag_opt_lsq_check_hessian (e04yb) is designed to check the bjk${b}_{jk}$ calculated by such user-supplied functions. As well as the function to be checked (lsqhes), you must supply a function (lsqfun) to evaluate the fi${f}_{i}$ and their first derivatives, and a point x = (x1,x2,,xn)T$x={\left({x}_{1},{x}_{2},\dots ,{x}_{n}\right)}^{\mathrm{T}}$ at which the checks will be made. Note that nag_opt_lsq_check_hessian (e04yb) checks functions of the form required by nag_opt_lsq_uncon_mod_deriv2_comp (e04he). nag_opt_lsq_check_hessian (e04yb) is essentially identical to CHKLSH in the NPL Algorithms Library.
nag_opt_lsq_check_hessian (e04yb) first calls user-supplied functions lsqfun and lsqhes to evaluate the first derivatives and the bjk${b}_{jk}$ at x$x$. Let J$J$ denote the m$m$ by n$n$ matrix of first derivatives of the residuals. The Hessian matrix of the sum of squares,
 G = JTJ + B, $G=JTJ+B,$
is calculated and projected onto two orthogonal vectors y$y$ and z$z$ to give the scalars yTGy${y}^{\mathrm{T}}Gy$ and zTGz${z}^{\mathrm{T}}Gz$ respectively. The same projections of the Hessian matrix are also estimated by finite differences, giving
 p = (yTg(x + hy) − yTg(x)) / h  and q = (zTg(x + hz) − zTg(x)) / h
$p=(yTg(x+hy)-yTg(x))/h and q=(zTg(x+hz)-zTg(x))/h$
respectively, where g()$g\left(\right)$ denotes the gradient vector of the sum of squares at the point in brackets and h$h$ is a small positive scalar. If the relative difference between p$p$ and yTGy${y}^{\mathrm{T}}Gy$ or between q$q$ and zTGz${z}^{\mathrm{T}}Gz$ is judged too large, an error indicator is set.

None.

Parameters

Compulsory Input Parameters

1:     m – int64int32nag_int scalar
The number m$m$ of residuals, fi(x)${f}_{i}\left(x\right)$, and the number n$n$ of variables, xj${x}_{j}$.
Constraint: 1nm$1\le {\mathbf{n}}\le {\mathbf{m}}$.
2:     lsqfun – function handle or string containing name of m-file
lsqfun must calculate the vector of values fi(x)${f}_{i}\left(x\right)$ and their first derivatives (fi)/(xj) $\frac{\partial {f}_{i}}{\partial {x}_{j}}$ at any point x$x$. (nag_opt_lsq_uncon_mod_deriv2_comp (e04he) gives you the option of resetting parameters of lsqfun to cause the minimization process to terminate immediately. nag_opt_lsq_check_hessian (e04yb) will also terminate immediately, without finishing the checking process, if the parameter in question is reset.)
[iflag, fvec, fjac, user] = lsqfun(iflag, m, n, xc, ldfjac, user)

Input Parameters

1:     iflag – int64int32nag_int scalar
To lsqfun, iflag will be set to 2$2$.
2:     m – int64int32nag_int scalar
The numbers m$m$ of residuals.
3:     n – int64int32nag_int scalar
The numbers n$n$ of variables.
4:     xc(n) – double array
The point x$x$ at which the values of the fi${f}_{i}$ and the (fi)/(xj) $\frac{\partial {f}_{i}}{\partial {x}_{j}}$ are required.
5:     ldfjac – int64int32nag_int scalar
The first dimension of the array fjac as declared in the (sub)program from which nag_opt_lsq_check_hessian (e04yb) is called.
6:     user – Any MATLAB object
lsqfun is called from nag_opt_lsq_check_hessian (e04yb) with the object supplied to nag_opt_lsq_check_hessian (e04yb).

Output Parameters

1:     iflag – int64int32nag_int scalar
If you reset iflag to some negative number in lsqfun and return control to nag_opt_lsq_check_hessian (e04yb), the function will terminate immediately with ifail set to your setting of iflag.
2:     fvec(m) – double array
Unless iflag is reset to a negative number, fvec(i)${\mathbf{fvec}}\left(\mathit{i}\right)$ must contain the value of fi${f}_{\mathit{i}}$ at the point x$x$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$.
3:     fjac(ldfjac,n) – double array
ldfjacm$\mathit{ldfjac}\ge {\mathbf{m}}$.
Unless iflag is reset to a negative number, fjac(i,j)${\mathbf{fjac}}\left(\mathit{i},\mathit{j}\right)$ must contain the value of (fi)/(xj) $\frac{\partial {f}_{\mathit{i}}}{\partial {x}_{\mathit{j}}}$ at the point x$x$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$ and j = 1,2,,n$\mathit{j}=1,2,\dots ,n$.
4:     user – Any MATLAB object
Note:  nag_opt_lsq_check_deriv (e04ya) should be used to check the first derivatives calculated by lsqfun before nag_opt_lsq_check_hessian (e04yb) is used to check the bjk${b}_{jk}$ since nag_opt_lsq_check_hessian (e04yb) assumes that the first derivatives are correct.
3:     lsqhes – function handle or string containing name of m-file
lsqhes must calculate the elements of the symmetric matrix
 m B(x) = ∑ fi(x)Gi(x), i = 1
$B(x)=∑i=1mfi(x)Gi(x),$
at any point x$x$, where Gi(x)${G}_{i}\left(x\right)$ is the Hessian matrix of fi(x)${f}_{i}\left(x\right)$. (As with lsqfun, a parameter can be set to cause immediate termination.)
[iflag, b, user] = lsqhes(iflag, m, n, fvec, xc, lb, user)

Input Parameters

1:     iflag – int64int32nag_int scalar
Is set to a non-negative number.
2:     m – int64int32nag_int scalar
The numbers m$m$ of residuals.
3:     n – int64int32nag_int scalar
The numbers n$n$ of variables.
4:     fvec(m) – double array
The value of the residual fi${f}_{\mathit{i}}$ at the point x$x$, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$, so that the values of the fi${f}_{\mathit{i}}$ can be used in the calculation of the elements of b.
5:     xc(n) – double array
The point x$x$ at which the elements of b are to be evaluated.
6:     lb – int64int32nag_int scalar
Gives the length of the array b.
7:     user – Any MATLAB object
lsqhes is called from nag_opt_lsq_check_hessian (e04yb) with the object supplied to nag_opt_lsq_check_hessian (e04yb).

Output Parameters

1:     iflag – int64int32nag_int scalar
If lsqhes resets iflag to some negative number, nag_opt_lsq_check_hessian (e04yb) will terminate immediately, with ifail set to your setting of iflag.
2:     b(lb) – double array
Unless iflag is reset to a negative number b must contain the lower triangle of the matrix B(x)$B\left(x\right)$, evaluated at the point in xc, stored by rows. (The upper triangle is not needed because the matrix is symmetric.) More precisely, b(j(j1) / 2 + k)${\mathbf{b}}\left(\mathit{j}\left(\mathit{j}-1\right)/2+\mathit{k}\right)$ must contain i = 1mfi(2fi)/(xjxk) $\sum _{\mathit{i}=1}^{m}{f}_{\mathit{i}}\frac{{\partial }^{2}{f}_{\mathit{i}}}{\partial {x}_{\mathit{j}}\partial {x}_{\mathit{k}}}$ evaluated at the point x$x$, for j = 1,2,,n$\mathit{j}=1,2,\dots ,n$ and k = 1,2,,j$\mathit{k}=1,2,\dots ,\mathit{j}$.
3:     user – Any MATLAB object
4:     x(n) – double array
n, the dimension of the array, must satisfy the constraint 1nm$1\le {\mathbf{n}}\le {\mathbf{m}}$.
x(j)${\mathbf{x}}\left(\mathit{j}\right)$, for j = 1,2,,n$\mathit{j}=1,2,\dots ,n$, must be set to the coordinates of a suitable point at which to check the bjk${b}_{jk}$ calculated by lsqhes. ‘Obvious’ settings, such as 0$0$ or 1$1$, should not be used since, at such particular points, incorrect terms may take correct values (particularly zero), so that errors could go undetected. For a similar reason, it is preferable that no two elements of x should have the same value.
5:     lb – int64int32nag_int scalar
The dimension of the array b as declared in the (sub)program from which nag_opt_lsq_check_hessian (e04yb) is called.
Constraint: lb(n + 1) × n / 2${\mathbf{lb}}\ge \left({\mathbf{n}}+1\right)×{\mathbf{n}}/2$.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: For n, the dimension of the array x.
The number m$m$ of residuals, fi(x)${f}_{i}\left(x\right)$, and the number n$n$ of variables, xj${x}_{j}$.
Constraint: 1nm$1\le {\mathbf{n}}\le {\mathbf{m}}$.
2:     user – Any MATLAB object
user is not used by nag_opt_lsq_check_hessian (e04yb), but is passed to lsqfun and lsqhes. Note that for large objects it may be more efficient to use a global variable which is accessible from the m-files than to use user.

Input Parameters Omitted from the MATLAB Interface

ldfjac iuser liw w lw

Output Parameters

1:     fvec(m) – double array
Unless you set iflag negative in the first call of lsqfun, fvec(i)${\mathbf{fvec}}\left(\mathit{i}\right)$ contains the value of fi${f}_{\mathit{i}}$ at the point supplied by you in x, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$.
2:     fjac(ldfjac,n) – double array
ldfjacm$\mathit{ldfjac}\ge {\mathbf{m}}$.
Unless you set iflag negative in the first call of lsqfun, fjac(i,j)${\mathbf{fjac}}\left(\mathit{i},\mathit{j}\right)$ contains the value of the first derivative (fi)/(xj) $\frac{\partial {f}_{\mathit{i}}}{\partial {x}_{\mathit{j}}}$ at the point given in x, as calculated by lsqfun, for i = 1,2,,m$\mathit{i}=1,2,\dots ,m$ and j = 1,2,,n$\mathit{j}=1,2,\dots ,n$.
3:     b(lb) – double array
Unless you set iflag negative in lsqhes, b(j × (j1) / 2 + k)${\mathbf{b}}\left(\mathit{j}×\left(\mathit{j}-1\right)/2+\mathit{k}\right)$ contains the value of bjk${b}_{\mathit{j}\mathit{k}}$ at the point given in x as calculated by lsqhes, for j = 1,2,,n$\mathit{j}=1,2,\dots ,n$ and k = 1,2,,j$\mathit{k}=1,2,\dots ,\mathit{j}$.
4:     user – Any MATLAB object
5:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Note: nag_opt_lsq_check_hessian (e04yb) may return useful information for one or more of the following detected errors or warnings.
Errors or warnings detected by the function:

Cases prefixed with W are classified as warnings and do not generate an error of type NAG:error_n. See nag_issue_warnings.

W ifail < 0${\mathbf{ifail}}<0$
A negative value of ifail indicates an exit from nag_opt_lsq_check_hessian (e04yb) because you have set iflag negative in user-supplied functions lsqfun or lsqhes. The setting of ifail will be the same as your setting of iflag. The check on lsqhes will not have been completed.
ifail = 1${\mathbf{ifail}}=1$
 On entry, m < n${\mathbf{m}}<{\mathbf{n}}$, or n < 1${\mathbf{n}}<1$, or ldfjac < m$\mathit{ldfjac}<{\mathbf{m}}$, or lb < (n + 1) × n / 2${\mathbf{lb}}<\left({\mathbf{n}}+1\right)×{\mathbf{n}}/2$, or liw < 1$\mathit{liw}<1$, or lw < 5 × n + m + m × n + n × (n − 1) / 2$\mathit{lw}<5×{\mathbf{n}}+{\mathbf{m}}+{\mathbf{m}}×{\mathbf{n}}+{\mathbf{n}}×\left({\mathbf{n}}-1\right)/2$, if n > 1${\mathbf{n}}>1$, or lw < 6 + 2 × m$\mathit{lw}<6+2×{\mathbf{m}}$, if n = 1${\mathbf{n}}=1$.
W ifail = 2${\mathbf{ifail}}=2$
You should check carefully the derivation and programming of expressions for the bjk${b}_{jk}$, because it is very unlikely that lsqhes is calculating them correctly.

Accuracy

ifail is set to 2$2$ if
 |yTGy − p| ≥ sqrt(h)(|yTGy| + 1.0) or |zTGz − q| ≥ sqrt(h)(|zTGz| + 1.0)
$|yTGy-p|≥h(|yTGy|+1.0) or |zTGz-q|≥h(|zTGz|+1.0)$
where h$h$ is set equal to sqrt(ε)$\sqrt{\epsilon }$ (ε$\epsilon$ being the machine precision as given by nag_machine_precision (x02aj)) and other quantities are defined as in Section [Description].

nag_opt_lsq_check_hessian (e04yb) calls lsqhes once and lsqfun three times.

Example

```function nag_opt_lsq_check_hessian_example
m = int64(15);
x = [0.19; -1.34; 0.88];
lb = int64(6);

y=[0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,0.37,0.58,0.73,0.96,1.34,2.10,4.39];

t = [1.0, 15.0, 1.0;
2.0, 14.0, 2.0;
3.0, 13.0, 3.0;
4.0, 12.0, 4.0;
5.0, 11.0, 5.0;
6.0, 10.0, 6.0;
7.0, 9.0, 7.0;
8.0, 8.0, 8.0;
9.0, 7.0, 7.0;
10.0, 6.0, 6.0;
11.0, 5.0, 5.0;
12.0, 4.0, 4.0;
13.0, 3.0, 3.0;
14.0, 2.0, 2.0;
15.0, 1.0, 1.0];

user = {y; t};

[fvec, fjac, b, user, ifail] = ...
nag_opt_lsq_check_hessian(m, @lsqfun, @lsqhes, x, lb, 'user', user)

function [iflag, fvecc, fjacc, user] = lsqfun(iflag, m, n, xc, ljc, user)
y = user{1};
t = user{2};

fvecc = zeros(m, 1);
fjacc = zeros(ljc, n);

for i = 1:double(m)
denom = xc(2)*t(i,2) + xc(3)*t(i,3);
fvecc(i) = xc(1) + t(i,1)/denom - y(i);
if (iflag ~= 0)
fjacc(i,1) = 1;
dummy = -1/(denom*denom);
fjacc(i,2) = t(i,1)*t(i,2)*dummy;
fjacc(i,3) = t(i,1)*t(i,3)*dummy;
end
end

function [iflag, b, user] = lsqhes(iflag, m, n, fvecc, xc, lb, user)
t = user{2};
b = zeros(lb, 1);

sum22 = 0;
sum32 = 0;
sum33 = 0;
for i = 1:double(m)
dummy = 2*t(i,1)/(xc(2)*t(i,2)+xc(3)*t(i,3))^3;
sum22 = sum22 + fvecc(i)*dummy*t(i,2)^2;
sum32 = sum32 + fvecc(i)*dummy*t(i,2)*t(i,3);
sum33 = sum33 + fvecc(i)*dummy*t(i,3)^2;
end
b(3) = sum22;
b(5) = sum32;
b(6) = sum33;
```
```

fvec =

-0.0020
-0.1076
-0.2330
-0.3785
-0.5836
-0.8689
-1.3464
-2.3739
-2.9750
-4.0132
-5.3226
-7.2917
-10.5703
-17.1274
-36.8087

fjac =

1.0000   -0.0406   -0.0027
1.0000   -0.0969   -0.0138
1.0000   -0.1785   -0.0412
1.0000   -0.3043   -0.1014
1.0000   -0.5144   -0.2338
1.0000   -0.9100   -0.5460
1.0000   -1.8098   -1.4076
1.0000   -4.7259   -4.7259
1.0000   -6.0762   -6.0762
1.0000   -7.8765   -7.8765
1.0000  -10.3970  -10.3970
1.0000  -14.1777  -14.1777
1.0000  -20.4789  -20.4789
1.0000  -33.0813  -33.0813
1.0000  -70.8885  -70.8885

b =

1.0e+04 *

0
0
1.5715
0
1.5712
1.5710

user =

[ 1x15 double]
[15x3  double]

ifail =

0

```
```function e04yb_example
m = int64(15);
x = [0.19; -1.34; 0.88];
lb = int64(6);

y=[0.14,0.18,0.22,0.25,0.29,0.32,0.35,0.39,0.37,0.58,0.73,0.96,1.34,2.10,4.39];

t = [1.0, 15.0, 1.0;
2.0, 14.0, 2.0;
3.0, 13.0, 3.0;
4.0, 12.0, 4.0;
5.0, 11.0, 5.0;
6.0, 10.0, 6.0;
7.0, 9.0, 7.0;
8.0, 8.0, 8.0;
9.0, 7.0, 7.0;
10.0, 6.0, 6.0;
11.0, 5.0, 5.0;
12.0, 4.0, 4.0;
13.0, 3.0, 3.0;
14.0, 2.0, 2.0;
15.0, 1.0, 1.0];

user = {y; t};

[fvec, fjac, b, user, ifail] = ...
e04yb(m, @lsqfun, @lsqhes, x, lb, 'user', user)

function [iflag, fvecc, fjacc, user] = lsqfun(iflag, m, n, xc, ljc, user)
y = user{1};
t = user{2};

fvecc = zeros(m, 1);
fjacc = zeros(ljc, n);

for i = 1:double(m)
denom = xc(2)*t(i,2) + xc(3)*t(i,3);
fvecc(i) = xc(1) + t(i,1)/denom - y(i);
if (iflag ~= 0)
fjacc(i,1) = 1;
dummy = -1/(denom*denom);
fjacc(i,2) = t(i,1)*t(i,2)*dummy;
fjacc(i,3) = t(i,1)*t(i,3)*dummy;
end
end

function [iflag, b, user] = lsqhes(iflag, m, n, fvecc, xc, lb, user)
t = user{2};
b = zeros(lb, 1);

sum22 = 0;
sum32 = 0;
sum33 = 0;
for i = 1:double(m)
dummy = 2*t(i,1)/(xc(2)*t(i,2)+xc(3)*t(i,3))^3;
sum22 = sum22 + fvecc(i)*dummy*t(i,2)^2;
sum32 = sum32 + fvecc(i)*dummy*t(i,2)*t(i,3);
sum33 = sum33 + fvecc(i)*dummy*t(i,3)^2;
end
b(3) = sum22;
b(5) = sum32;
b(6) = sum33;
```
```

fvec =

-0.0020
-0.1076
-0.2330
-0.3785
-0.5836
-0.8689
-1.3464
-2.3739
-2.9750
-4.0132
-5.3226
-7.2917
-10.5703
-17.1274
-36.8087

fjac =

1.0000   -0.0406   -0.0027
1.0000   -0.0969   -0.0138
1.0000   -0.1785   -0.0412
1.0000   -0.3043   -0.1014
1.0000   -0.5144   -0.2338
1.0000   -0.9100   -0.5460
1.0000   -1.8098   -1.4076
1.0000   -4.7259   -4.7259
1.0000   -6.0762   -6.0762
1.0000   -7.8765   -7.8765
1.0000  -10.3970  -10.3970
1.0000  -14.1777  -14.1777
1.0000  -20.4789  -20.4789
1.0000  -33.0813  -33.0813
1.0000  -70.8885  -70.8885

b =

1.0e+04 *

0
0
1.5715
0
1.5712
1.5710

user =

[ 1x15 double]
[15x3  double]

ifail =

0

```