Hello all,
A while ago I asked for a means of trapping .NaNs during the course of
maximization iterations. Since there were no suggestions, I hacked
together a crude solution of my own. The attached code is the MaxBFGS
source distributed with Ox, but with a few lines added by me. This
checks the function value, the parameter vector, and the score vector
for .NaNs, and exits the iterations if any are found.
By way of explanation, and in my experience, when MaxBFGS encounters
.NaNs, Ox 2.10 under Linux starts eating up memory, then swap space,
until the X server crashes. Your mileage may vary.
To use the modified version, compile the source, maximize.ox using
oxl -c maximize.ox
Then move the maximize.oxo file to the /../ox/include directory
Or, just copy the attached compiled file maximize.oxo into your
/../ox/include directory.
I suggest backing up the original files first. Cheers.
--
Michael Creel [log in to unmask]
Tel. 93-581-1696 Fax. 93-581-2012
Dept. of Economics and Economic History
Universitat Autonoma de Barcelona
"Even paranoids have real enemies."
/*--------------------------------------------------------------------------
* maximize.ox - code for numerical differentiation and function maximization
*
* (C) Jurgen Doornik 1994-98
*
*--------------------------------------------------------------------------*/
#include <oxstd.h>
#include <oxfloat.h>
#include <maximize.h>
/*========================= numerical derivatives ==========================*/
const decl SQRT_EPS =1E-8; /* appr. square root of machine precision */
const decl DIFF_EPS =1E-8; /* Rice's formula: log(DIFF_EPS)=log(MACH_EPS)/2 */
const decl DIFF_EPS1=5E-6; /* Rice's formula: log(DIFF_EPS)=log(MACH_EPS)/3 */
const decl DIFF_EPS2=1E-4; /* Rice's formula: log(DIFF_EPS)=log(MACH_EPS)/4 */
static dFiniteDiff0(const x)
{
return max( (fabs(x) + SQRT_EPS) * SQRT_EPS, DIFF_EPS);
}
static dFiniteDiff1(const x)
{
return max( (fabs(x) + SQRT_EPS) * SQRT_EPS, DIFF_EPS1);
}
static dFiniteDiff2(const x)
{
return max( (fabs(x) + SQRT_EPS) * SQRT_EPS, DIFF_EPS2);
}
/*----------------------------- Num1Derivative -----------------------------*/
Num1Derivative(const func, vP, const avScore)
{
decl i, cp = rows(vP), left, right, fknowf = FALSE, p, h, f, fm, fp, v;
v = new matrix[cp][1];
for (i = 0; i < cp; i++) /* get 1st derivative by central difference */
{
p = double(vP[i][0]);
h = dFiniteDiff1(p);
vP[i][0] = p + h;
right = func(vP, &fp, 0, 0);
vP[i][0] = p - h;
left = func(vP, &fm, 0, 0);
vP[i][0] = p; /* restore original parameter */
if (left && right)
v[i][0] = (fp - fm) / (2 * h); /* take central difference */
else if (left)
{
if (!fknowf) /* see if we already know f */
{ fknowf = func(vP, &f, 0, 0); /* not: try to get it */
if (!fknowf)
return FALSE;
}
v[i][0] = (f - fm) / h; /* take left difference */
}
else if (right)
{
if (!fknowf)
{ fknowf = func(vP, &f, 0, 0);
if (!fknowf)
return FALSE;
}
v[i][0] = (fp - f) / h; /* take right difference */
}
else
return FALSE;
}
avScore[0] = v;
return TRUE;
}
/*--------------------------- END Num1Derivative ---------------------------*/
/*----------------------------- Num2Derivative -----------------------------*/
Num2Derivative(const func, vP, const amHessian)
{
decl i, j, ret_val, pi, pj, hi, hj, f, fpp, fmm, fpm, fmp, cp = rows(vP), m;
m = new matrix[cp][cp];
ret_val = FALSE;
func(vP, &f, 0, 0);/* get function at vP, around which to differentiate */
for (i = 0; i < cp; i++)/* approximate 2nd deriv. by central difference */
{
pi = double(vP[i][0]);
hi = dFiniteDiff2(pi);
ret_val = TRUE;
for (j = 0; j < i; j++) /* off-diagonal elements */
{
pj = double(vP[j][0]);
hj = dFiniteDiff2(pj);
vP[i][0] = pi + hi; vP[j][0] = pj + hj; /* +1 +1 */
ret_val = ret_val && func(vP, &fpp, 0, 0);
vP[i][0] = pi - hi; vP[j][0] = pj - hj; /* -1 -1 */
ret_val = ret_val && func(vP, &fmm, 0, 0);
vP[i][0] = pi + hi; vP[j][0] = pj - hj; /* +1 -1 */
ret_val = ret_val && func(vP, &fpm, 0, 0);
vP[i][0] = pi - hi; vP[j][0] = pj + hj; /* -1 +1 */
ret_val = ret_val && func(vP, &fmp, 0, 0);
if (!ret_val)
return FALSE;
m[i][j] = m[j][i] = ((fpp + fmm) - (fpm + fmp)) / (4 * hi * hj);
vP[j][0] = pj;
}
/* diagonal elements */
vP[i][0] = pi + 2 * hi; /* +1 +1 */
ret_val = ret_val && func(vP, &fpp, 0, 0);
vP[i][0] = pi - 2 * hi; /* -1 -1 */
ret_val = ret_val && func(vP, &fmm, 0, 0);
if (!ret_val)
return FALSE;
m[i][i] = ((fpp + fmm) - (2 * f)) / (4 * hi * hi);
vP[i][0] = pi;
}
amHessian[0] = m;
return ret_val;
}
/*---------------------------- END Num2Derivative --------------------------*/
/*------------------------------- NumJacobian ------------------------------*/
NumJacobian(const func, vU, const amJacobian)
{
decl i, ret_val = TRUE, p, h, fp, fm, cu = sizerc(vU), mj = 0;
/* approximate first derivative by taking central difference */
for (i = 0; i < cu && ret_val; i++) /* compute the gradient delta */
{
p = double(vU[i]);
h = dFiniteDiff1(p);
vU[i] = p + h;
ret_val = ret_val && func(&fp, vU);
vU[i] = p - h;
ret_val = ret_val && func(&fm, vU);
vU[i] = p; /* restore original parameter */
if (i == 0)
mj = zeros(rows(vec(fp)), cu);
mj[][i] = vec(fp - fm) / (2 * h); /* central difference */
}
amJacobian[0] = mj;
return ret_val;
}
/*----------------------------- END NumJacobian ----------------------------*/
/*========================= function maximization ==========================*/
static decl s_mxIter = 1000;
static decl s_dEps1 = 1e-4, s_dEps2 = 5e-3; /* convergence criteria */
static decl s_iPrint = 0; /* print results every s_iPrint'th iteration */
static decl s_bCompact = 0; /* TRUE: print iterations in compact format */
static decl s_dSteptolLinear = 1e-14;
/*------------------- fIsConvergence/iConvergenceCode ----------------------*/
static fIsStrongConvergence(vP, const vScore, const vDelta)
{
vP = vP .? vP .: 1; /* replace 0 by 1 */
if (fabs(vScore .* vP) < s_dEps1 && fabs(vDelta) < 10 * s_dEps1 * fabs(vP))
return TRUE;
return FALSE;
}
static fIsWeakConvergence(vP, const vScore)
{
vP = vP .? vP .: 1; /* replace 0 by 1 */
if (fabs(vScore .* vP) < s_dEps2)
return TRUE;
return FALSE;
}
static iConvergenceCode(const fConv, const fStepOK, const fItMax)
{
if (!fStepOK)
return (fConv) ? MAX_WEAK_CONV : MAX_LINE_FAIL;
else if (fConv)
return MAX_CONV;
else
return MAX_MAXIT;
}
MaxConvergenceMsg(const iCode)
{
if (iCode == MAX_CONV)
return "Strong convergence";
else if (iCode == MAX_WEAK_CONV)
return "Weak convergence (no improvement in line search)";
else if (iCode == MAX_MAXIT)
return "No convergence (maximum no of iterations reached)";
else if (iCode == MAX_LINE_FAIL)
return "No convergence (no improvement in line search)";
else if (iCode == MAX_FUNC_FAIL)
return "No convergence (function evaluation failed)";
else
return "No convergence";
}
MaxControlEps(const dEps1, const dEps2)
{
if (dEps1 > 0) ::s_dEps1 = dEps1;
if (dEps2 > 0) ::s_dEps2 = dEps2;
}
MaxControl(const mxIter, const iPrint, ...)
{
if (mxIter >= 0) ::s_mxIter = mxIter;
if (iPrint >= 0) ::s_iPrint = iPrint;
decl va = va_arglist();
if (sizeof(va) > 0 && va[0] >= 0)
s_bCompact = va[0];
}
GetMaxControl()
{
return { s_mxIter, s_iPrint, s_bCompact };
}
GetMaxControlEps()
{
return { s_dEps1, s_dEps2 };
}
/*----------------- END fIsConvergence/iConvergenceCode --------------------*/
/*----------------------------- monitor ------------------------------------*/
static monitor(const sMethod, const cIter, const vP, const vScore, const vDelta,
const dFuncInit, const dFunc, const steplen, const condhes, const ret_val)
{
if (s_bCompact)
{
if (ismatrix(vScore))
{ decl dgmax = max(fabs(vScore));
decl vpp = vP .? vP .: 1; /* replace 0 by 1 */
decl tol1 = max(fabs(vScore .* vpp));
decl tol2 = max(fabs(vDelta) ./ (10 * fabs(vpp)));
println("it", "%-4d", cIter, " f=", "%#14.7g", dFunc,
" df=", "%#10.4g", dgmax, " e1=", "%#10.4g", tol1,
" e2=", "%#10.4g",tol2, " step=", "%.4g", steplen);
}
else
println("it", "%-4d", cIter, " f=", "%#14.7g", dFunc,
" step=", "%.4g", steplen);
if (ret_val != MAX_NOCONV)
println(MaxConvergenceMsg(ret_val));
return;
}
if (cIter <= 0)
println("\nStarting values");
else
{ println("\nPosition after ", cIter, " ", sMethod, " iterations");
if (ret_val != MAX_NOCONV)
println("Status: ", MaxConvergenceMsg(ret_val));
}
print("parameters", vP');
if (columns(vScore))
print("gradients", vScore');
if (cIter <= 0)
print("Initial function = ", "%19.12g", dFuncInit);
else
{ print("function value = ", "%19.12g", dFunc);
if (condhes != 0)
print(" condhes = ", condhes);
if (steplen != 1)
print(" steplen = ", steplen);
}
print("\n");
}
/*----------------------------- END monitor --------------------------------*/
/*----------------------------- fLineSearch --------------------------------*/
static fLineSearch(const func, const avP, const adFunc, const vDelta,
const adStep, const dFuncBase)
{
decl dstep = adStep[0], f;
do /* move towards zero with linear line search */
{ dstep /= 2;
avP[0] -= dstep * vDelta;
if (func(avP[0], adFunc, 0, 0) && adFunc[0] > dFuncBase)
{
/* function & better function value; continue while improvement */
while (func(avP[0] - dstep * vDelta / 2, &f, 0, 0) && f > adFunc[0])
{ dstep /= 2;
avP[0] -= dstep * vDelta;
adFunc[0] = f;
}
adStep[0] = dstep;
return 1;
}
} while (dstep >= s_dSteptolLinear);
/* if we're here, the line search failed */
/* reset function value and parameters to that at steplength 0 */
avP[0] -= dstep * vDelta; /* this removes the remaining step length */
adFunc[0] = dFuncBase; /* dFuncBase was function value of previous iter */
adStep[0] = dstep;
return 0; /* failed to improve in line search */
}
/*------------------------- END fLineSearch --------------------------------*/
/*------------------------------ MaxBFGS -----------------------------------*/
MaxBFGS(const func, const avP, const adFunc, const amHessian, const fNumDer)
{
decl itno = -1, cp = sizer(avP[0]);
decl fstepok, fconv, fok, ret_val = MAX_NOCONV, fhesreset = FALSE;
decl dsteplen, dmxstep, dfprev, dfinit, doldstep;
decl d_2, d_3;
decl vscore, vdelta, vgamma, vhgamm, minvhess;
vscore = vdelta = vgamma = vhgamm = zeros(cp, 1);
dsteplen = dmxstep = 1.0;
if (isarray(amHessian)) /* process amHessian argument */
minvhess = amHessian[0];
else if (ismatrix(amHessian))
minvhess = amHessian;
else
minvhess = unit(cp); /* default initial Hessian */
if (!isarray(avP) || !isarray(adFunc))
{ eprint("MaxBFGS(): arguments 2 and 3 must be an address\n");
return MAX_FUNC_FAIL;
}
if (columns(avP[0]) > 1)
{ eprint("MaxBFGS(): expecting parameters in column vector\n");
return MAX_FUNC_FAIL;
}
if (rows(minvhess) != cp || columns(minvhess) != cp)
{ eprint("MaxBFGS(): initial Hessian dimension error, using I\n");
minvhess = unit(cp); /* default initial Hessian */
}
if (!func(avP[0], adFunc, 0, 0))
adFunc[0] = -DBL_MAX;
dfinit = adFunc[0];
if (cp == 0)
return MAX_CONV;
do /* start iterating, first loop is initialization */
{ itno++; fstepok = TRUE; fconv = FALSE;
doldstep = dsteplen;
dsteplen = dmxstep;
dfprev = adFunc[0]; /* keep old dfunc; compute new parameters */
avP[0] += dsteplen * vdelta; /* first time: add 0 */
vgamma = vscore; /* save old score */
/* evaluate function: 1 is success */
fok = func(avP[0], adFunc, fNumDer ? 0 : &vscore, 0);
/* now have analytical score if used, else get numerical later */
if (isnan(adFunc[0]))
{
print("\nERROR during iterations!\nFunction value is a .NAN! Better luck next time\n");
return;
}
if (isnan(avP[0]))
{
print("\nERROR during iterations!\nParameter vector includes .NANs! Better luck next time\n");
return;
}
if (isnan(vscore[0]))
{
print("\nERROR during iterations!\nScore vector includes .NANs! Better luck next time\n");
return;
}
if (!fok) /* yes: function evaluation failed */
{ if (itno == 0)
{ eprint("MaxBFGS(): function evaluation failed at starting values\n");
return MAX_FUNC_FAIL;
}
adFunc[0] = -DBL_MAX; /* very large negative number */
}
if (itno > 0 && (adFunc[0] < dfprev || !fok) )
{ /* failure to evaluate or improve: determine steplength in (0,1] */
fstepok = fLineSearch(func, avP, adFunc, vdelta, &dsteplen, dfprev);
if (!fstepok) /* no improvement in line search */
{ /* loglik/param have been reset to that at steplength 0 */
vscore = vgamma; /* also reset old score */
if (!fhesreset) /* don't do twice in a row */
{
fhesreset = TRUE; /* try once more, with a unit Hessian */
vdelta = vscore; /* vdelta = q */
continue; /* try again */
}
else
{ fconv = fIsWeakConvergence(avP[0], vscore);
break; /* stop iterating */
}
}
if (!fok) /* yes: failed to evaluate at steplength 1 */
dmxstep = max(dsteplen, 0.1);/*switch to restr. step method */
/* success in line search: now need analytical score */
if (!fNumDer && !func(avP[0], adFunc, &vscore, 0))
return MAX_FUNC_FAIL;
vdelta *= dsteplen; /* actual delta after line search */
}
fhesreset = FALSE; /* was a succesful iteration */
if (itno > 0) /* adjust restricted step method */
{ if (dsteplen >= dmxstep) /* grow maximum steplength */
{ dmxstep *= 2;
if (dmxstep > 1) dmxstep = 1.0;
}
else if (dmxstep < 1 && dsteplen < dmxstep / 2) /* shrink */
{ dmxstep /= 2;
if (dmxstep < 1e-6) dmxstep = 1e-6;
}
}
/* yes: get numerical score; else already have analytical score */
if (fNumDer && !Num1Derivative(func, avP[0], &vscore))
{
eprint("MaxBFGS(): numerical score failed\n");
return MAX_FUNC_FAIL;
}
vgamma -= vscore; /* ç = change in score */
vhgamm = minvhess * vgamma;
d_2 = double(vdelta'vgamma);
d_3 = double(vgamma'vhgamm);
if (d_2 < 0) /* do not change in == 0, else H always reset at start */
{ d_2 = 0;
minvhess = unit(cp); /* reset Hessian */
}
else if (d_2 != 0)
d_2 = 1 / d_2; /* d_2 = 1/d'g */
d_3 = (1 + d_3 * d_2) * d_2; /* (1 + g'Hg/d'g)/d'g */
/* Hnew = Hold + [(1 + g'Hg/d'g)/d'g]dd' - (dg'H + Hgd')/d'g */
/* first time: d_2 == d_3 == 0, so Hessian doesn't change */
minvhess += d_3 * vdelta * vdelta' -
(vdelta * vhgamm' + vhgamm * vdelta') * d_2;
vdelta = minvhess * vscore; /* vdelta = Hq */
fconv = fIsStrongConvergence(avP[0], vscore, vdelta);
if (!fconv && s_iPrint > 0 && imod(itno, s_iPrint) == 0)
monitor("BFGS", itno, avP[0], vscore, vdelta, dfinit, adFunc[0],
dsteplen, 0, ret_val);
} while (!(fconv || itno >= s_mxIter));
ret_val = iConvergenceCode(fconv, fstepok, itno >= s_mxIter);
if (s_iPrint > 0)
monitor("BFGS", itno, avP[0], vscore, vdelta, dfinit, adFunc[0],
dsteplen, 0, ret_val);
if (isarray(amHessian))
amHessian[0] = minvhess;
return ret_val;
}
/*----------------------------- END MaxBFGS --------------------------------*/
/*------------------------------ MaxSimplex --------------------------------*/
static mInitSimplex(const func, const vP, const dSteplen, const vDelta)
{
decl i, cp = rows(vP), msimplex, df;
msimplex = new matrix[cp+1][cp+1]; /* construct initial simplex */
/* each column has parameter point; function value in bottom row */
msimplex[][] = unit(cp) * vP + diag(dSteplen * vDelta);
msimplex[][cp] = vP;
for (i = 0; i <= cp; i++) /* evaluate function at simplex */
{
if (!func(msimplex[0:cp-1][i], &df, 0, 0))
df = -DBL_MAX;
msimplex[cp][i] = df;
}
return msimplex;
}
static fCheckSimplex(const func, vP, const dFunc, const avDelta, const avScore,
const dEps)
{
decl i, ret_val = TRUE, p, h, fm, fp, cp = rows(vP), v;
avDelta[0] = new matrix[cp][1];
v = ones(cp,1);
for (i = 0; i < cp; i++)
{
p = vP[i][0];
avDelta[0][i][0] = h = dFiniteDiff1(p);
vP[i][0] = p + h;
ret_val = ret_val && func(vP, &fp, 0, 0);
vP[i][0] = p - h;
ret_val = ret_val && func(vP, &fm, 0, 0);
ret_val = ret_val && (fp <= dFunc || fp - dFunc < dEps)
&& (fm <= dFunc || fm - dFunc < dEps);
v[i][0] = (fp - fm) / (2 * h); /* take central difference */
vP[i][0] = p; /* restore original parameter */
}
avScore[0] = v;
return ret_val;
}
MaxSimplex(const func, const avP, const adFunc, vDelta)
{
decl i, itno = -1, cp = sizer(avP[0]), cp1, iconv;
decl fcontract, fconv, fok, ret_val = MAX_NOCONV, fweak;
decl dsteplen, dfinit, dfmin, dfmax, dfnew, dfext;
decl vpmin, vpmax, vpnew, vpext, msimplex, vpcentroid, vscore;
decl alpha = 1.0, beta = 0.5, gamma = 2.0, algam = alpha * gamma;
dsteplen = 1.0;
cp1 = cp - 1;
iconv = 6;
if (!isarray(avP) || !isarray(adFunc))
{ eprint("MaxSimplex(): avP and adFunc must be an address\n");
return MAX_FUNC_FAIL;
}
if (rows(vDelta) == 0 && !Num1Derivative(func, avP[0], &vDelta))
{ eprint("MaxSimplex(): numerical derivatives failed\n");
return MAX_FUNC_FAIL;
}
else if (rows(vDelta) != cp)
{ eprint("MaxSimplex(): input dimension error\n");
return MAX_FUNC_FAIL;
}
msimplex = mInitSimplex(func, avP[0], dsteplen, vDelta);
dfinit = dfmax = msimplex[cp][cp];
vpmax = msimplex[:cp1][cp];
do /* start iterating */
{ itno++; iconv--; fcontract = fconv = fweak = FALSE;
msimplex = sortbyr(msimplex, cp);
if (s_iPrint > 0 && imod(itno, s_iPrint) == 0)
monitor("Simplex", itno, vpmax, 0, 0, dfinit, dfmax, dsteplen, 0, ret_val);
dfmin = msimplex[cp][0]; vpmin = msimplex[:cp1][0];
dfmax = msimplex[cp][cp]; vpmax = msimplex[:cp1][cp];
/* check for convergence */
if (iconv < 0 && (dfmax - dfmin) < s_dEps1 * 1e-8)
{ /* yes, now check for minimum */
if (!fCheckSimplex(func, vpmax, dfmax, &vDelta, &vscore, s_dEps1 * 1e-8) )
{
msimplex = mInitSimplex(func, vpmax, 1.0, vDelta);
if (itno < s_mxIter)
{ iconv = 6;
continue;
}
else
{ break;
}
}
fconv = fIsStrongConvergence(vpmax, vscore, 0);
if (fconv)
break;
}
vpcentroid = sumr(msimplex[:cp1][1:]) / cp; /* omit lowest point */
/* reflection */
vpnew = (1 + alpha) * vpcentroid - alpha * vpmin;
if (!func(vpnew, &dfnew, 0, 0))
dfnew = -DBL_MAX;
if (dfnew > dfmax) /* is an improvement on previous best */
{ /* try an extension step */
vpext = (1 + algam) * vpcentroid - algam * vpmin;
if (!func(vpext, &dfext, 0, 0))
dfext = -DBL_MAX;
if (dfext > dfnew) /* is an improvement on the reflection */
vpnew = vpext, dfnew = dfext; /* so keep the extension */
}
else if (dfnew < msimplex[cp][1]) /* not even better than 2nd worst */
{ /* contract the simplex */
if (dfnew < dfmin) /* yes: worse than previous worst */
vpnew = vpmin, dfnew = dfmin;
vpext = (1 - beta) * vpcentroid + beta * vpnew; /* contract */
if (!func(vpext, &dfext, 0, 0))
dfext = -DBL_MAX;
if (dfext > dfnew) /* contraction is an improvement */
vpnew = vpext, dfnew = dfext; /* so keep the contraction */
else
{ /* do a general contraction around the current best */
for (i = 0; i < cp; i++)
{
vpnew = msimplex[:cp1][i] = (msimplex[:cp1][i] + vpmax) * beta;
if (!func(vpnew, &dfnew, 0, 0))
dfnew = -DBL_MAX;
msimplex[cp][i] = dfnew;
}
fcontract = TRUE;
}
}
if (!fcontract)
msimplex[][0] = vpnew, msimplex[cp][0] = dfnew;
} while (!(fconv || itno >= s_mxIter));
ret_val = iConvergenceCode(fconv, TRUE, itno >= s_mxIter);
avP[0] = vpmax; adFunc[0] = dfmax;
if (s_iPrint > 0)
monitor("Simplex", itno, vpmax, vscore, 0, dfinit, dfmax, 1.0, 0, ret_val);
return ret_val;
}
/*------------------------------- END MaxSimplex ---------------------------*/
/*-------------------------------- MaxNewton -------------------------------*/
MaxNewton(const func, const avP, const adFunc, const amHessian, const fNumDer)
{
decl i, itno = -1, cp = sizer(avP[0]);
decl fstepok, fconv, fsingular = FALSE;
decl p, h, dfprev, dfinit, dsteplen, ret_val = MAX_NOCONV;
decl vdelta, vscale, vdiag, vscore, vscoreprev, ml, mhess;
vscore = vdelta = zeros(cp, 1);
mhess = unit(cp);
if (!isarray(avP) || !isarray(adFunc))
{ eprint("MaxNewton(): arguments 2 and 3 must be an address\n");
return MAX_FUNC_FAIL;
}
if (columns(avP[0]) > 1)
{ eprint("MaxNewton(): expecting parameters in column vector\n");
return MAX_FUNC_FAIL;
}
adFunc[0] = -DBL_MAX; /* to avoid (null) value */
/* first time (itno == 0): evaluate starting values */
do /* start iterating */
{ itno++; dsteplen = 1; fstepok = TRUE; fconv = FALSE;
dfprev = adFunc[0]; /* keep old dfunc; compute new parameters */
vscoreprev = vscore;
avP[0] += vdelta; /* first time: add 0 */
/* evaluate function: 1 is success */
if (!func(avP[0], adFunc, &vscore, (fNumDer) ? 0 : &mhess))
{
if (itno == 0)
{ eprint("MaxNewton(): function evaluation failed at starting values\n");
return MAX_FUNC_FAIL;
}
adFunc[0] = -DBL_MAX; /* very large negative number */
}
/* now have analytical Hessian if used, else get numerical later */
if (itno == 0)
dfinit = adFunc[0];
else if (adFunc[0] < dfprev || fsingular)
{ /* failure to evaluate or improve: determine steplength in (0,1] */
fstepok = fLineSearch(func, avP, adFunc, vdelta, &dsteplen, dfprev);
if (!fstepok && !fsingular) /* not better, try steepest decent */
{ /* loglik/param have been reset to that at steplength 0 */
avP[0] += vscoreprev;
vdelta = vscoreprev;
if (!func(avP[0], adFunc, 0, 0))
adFunc[0] = -DBL_MAX;
dsteplen = 1;
fstepok = fLineSearch(func, avP, adFunc, vdelta, &dsteplen, dfprev);
}
if (!fstepok) /* no improvement in line search */
{ /* loglik/param have been reset to that at steplength 0 */
vscore = vscoreprev;/* ml, vscale are still of steplength 0 */
fconv = fIsWeakConvergence(avP[0], vscore);
break; /* stop iterating */
}
/* get derivatives at new point */
if (!func(avP[0], adFunc, &vscore, (fNumDer) ? 0 : &mhess))
return MAX_FUNC_FAIL;
}
if (fNumDer)
{
vdiag = zeros(cp, 1);
for (i = 0; i < cp; i++) /* numerical approx.to minus hessian */
{
p = avP[0][i];
h = dFiniteDiff1(p);
avP[0][i] = p + h;
if (!func(avP[0], adFunc, &vdiag, 0))
return MAX_FUNC_FAIL;
mhess[][i] = (vscore - vdiag) / h; /* forward difference */
avP[0][i] = p; /* restore original parameter */
}
mhess = (mhess + mhess') / 2; /* make symmetric */
}
else
mhess = -mhess; /* use minus 2nd der., Q = -H */
vscale = diagonal(mhess)'; /* determine new scale factors */
vscale = vscale .> 0 .? 1 ./ sqrt(vscale) .: 1;
/* solve the system q=Qb as Sq=SQSa, Sa=b */
if (decldl(vscale' .* mhess .* vscale, &ml, &vdiag) == 1)
{ /* solve system for scaled score; solution stored in vdelta */
vdelta = solveldl(ml, vdiag, vscore .* vscale) .* vscale;
fsingular = FALSE;
}
else
{ ml = mhess; /* remember for return value */
vdelta = vscore; /* steepest descent */
fsingular = TRUE;
}
fconv = fIsStrongConvergence(avP[0], vscore, vdelta);
if (!fconv && s_iPrint > 0 && imod(itno, s_iPrint) == 0)
monitor("Newton", itno, avP[0], vscore, vdelta, dfinit, adFunc[0],
dsteplen, 0, ret_val);
} while (!(fconv || itno >= s_mxIter));
ret_val = iConvergenceCode(fconv, fstepok, itno >= s_mxIter);
if (s_iPrint > 0)
monitor("Newton", itno, avP[0], vscore, vdelta, dfinit, adFunc[0],
dsteplen, 0, ret_val);
if (isarray(amHessian))
{ if (fsingular)
amHessian[0] = -invertgen(ml, 1); /* use generalized inverse */
else
amHessian[0] = -vscale' .* solveldl(ml, vdiag, unit(cp)) .* vscale;
}
return ret_val;
}
/*----------------------------- END MaxNewton ------------------------------*/
|
