} /* Next i */ return mG; /* Return matrix of individual gradients */ }

MaxBHHH(const fnF, const avX, const adFunc, const amInvHess, const fNumDer) { decl vX0,vX1, /* Parameter vectors */ vF0,vF1,vF2, /* Function values (n*1)-vectors */ mG0,vG0,mG1,mG2, /* Gradients in matrix and vector form, Hessian, direction */ mH,mHi,vS, dEps1,dEps2, /* Tolerance for convergence tests */ fConv=MAX_NOCONV,fStepOK, /* Flags for indicating convergence etc*/ dLambda; /* Line search parameter, step length */ adFunc[0] = M_NAN; amInvHess[0] = M_NAN*unit(sizerc(avX[0]));

decl iMaxIter,iPrint,bCompact; [iMaxIter,iPrint,bCompact]=GetMaxControl(); // println("/* INITIALIZATION */"); vX0 = avX[0]; [dEps1,dEps2] = GetMaxControlEps();

// println("/* Check initial function values */"); decl dStart = timer(); fnF(vX0,&vF0,0,0); decl sTimeSpan = sprint((timer() - dStart) / 100); // println("\n\t- One function evaluation ",fNumDer ? "":"(excluding analytical gradients) ", // "took ", sTimeSpan," sec.");

dStart = timer(); decl fInitFunc = fnF(vX0,&vF0,fNumDer ? 0 : &mG0,0); // println("\n\t- One function evaluation ",fNumDer ? "":"(including analytical gradients) ", // "took ", sprint((timer() - dStart) / 100)," sec.");

if (!fInitFunc){ format(90); println("\nInitial function evaluation failed!", "\nInitial position:",vX0', "Initial function values:",vF0' ); return MAX_FUNC_FAIL; }

for (decl i=0;i<iMaxIter;++i){ if (isnan(vF0)) print("Function evaluation failed! vX0'=",vX0'); decl dIterStartTime = timer();

// println("\n/* STEP 1: CALCULATE FUNCTION VALUE, GRADIENT AND HESSIAN */"); // println("vF0 = ",double(sumc(vF0))); // mG = fnGradient(fnF,vX0,vF0); if (fNumDer) mG0 = fnGradient(fnF,vX0,vF0); if (isnan(mG0)){ println(vX0',mG0); return MAX_FUNC_FAIL; } vG0 = sumc(mG0); mH = mG0'mG0;

// println("/* STEP 2: CALCULATE DIRECTION */"); // mHi = invertsym(mH); /* This should/could be updated to LU-decomp solution */ // vS = mHi*vG0';

decl vscale = diagonal(mH)', /* determine new scale factors */ ml,vdiag,fsingular; vscale = vscale .> 0 .? 1 ./ sqrt(vscale) .: 1; // Override scale // println(vscale'); /* solve the system q=Qb as Sq=SQSa, Sa=b */ if (decldl(vscale' .* mH .* vscale, &ml, &vdiag) == 1) { /* solve system for scaled score; solution stored in vdelta */ vS = solveldl(ml, vdiag, vG0' .* vscale) .* vscale; fsingular = FALSE; } else { ml = mH; /* remember for return value */ println(mH); vS = -vG0'; /* steepest decent */ // println("vX0' = ",vX0',"vF0' = ",vF0',"mH = ",mH,"vS' = ",vS'); fsingular = TRUE; } if (vG0*vS<0){ println("Changing direction");vS=-vS; }

// println("/* STEP 2B: LINE SEARCH, Contraction or extension of step length */"); // println("vF0 = ","%12.8f",double(sumc(vF0))); dLambda = 1; /* Was 1 initialize step length dLambda */ decl fStepOK=FALSE, bContraction, bF1 = fnF(vX0+dLambda*vS,&vF1,0,0), // We don't need the derivatives in the line search bF2 = fnF(vX0+2*dLambda*vS,&vF2,0,0);

if (!bF1 || !bF2 || sumc(vF2)<sumc(vF1) ) // Full step lengt implies function fail or lower likelihood bContraction = TRUE; // Reduce step lengt else bContraction = FALSE; // otherwise try to increase step lengt

while (dLambda >1e-15 && dLambda < 1e5){ // Min/max step lengths dLambda = bContraction ? 0.5*dLambda : 2*dLambda; // New potential dLambda (step length)

decl bF2 = fnF(vX0+dLambda*vS,&vF2,0,0);

// println( "bF1 = ",bF1," sumc(vF1)=","%12.8f",double(sumc(vF1)), // "\nbF2 = ",bF2," sumc(vF2)=","%12.8f",double(sumc(vF2))," dLambda=",double(dLambda)); // Case 1. dLambda imples missing value => Continue if we are in contraction phase, // other wise stop. if (bF2==FALSE || bF1==FALSE) if (bContraction){ // println("Case 1a. (bF2==FALSE or bF1==FALSE and contraction => Continue to contract) ",dLambda); bF1=bF2; vF1=vF2; } else{ // println("Case 1b. (bF2==FALSE and expansion => Stop expanding) ",dLambda); dLambda = bContraction ? 2*dLambda : 0.5*dLambda; // Go one step back. fStepOK = TRUE; fnF(vX0+dLambda*vS,&vF2,fNumDer ? 0 : &mG1,0); // We need gradient at chosen step length break; } // Case 2. dLambda implies lower likelihood => Stop else if (sumc(vF2)<sumc(vF1)){ // println("Case 2. (sumc(vF2)<sumc(vF1) ",double(sumc(vF2))," ",double(sumc(vF1)), " dLambda=",dLambda); dLambda = bContraction ? 2*dLambda : 0.5*dLambda; // Go one step back. fStepOK = TRUE; fnF(vX0+dLambda*vS,&vF2,fNumDer ? 0 : &mG1,0); // We need gradient at chosen step length break; }

// Case 3. dLambda implies higher likelihood => Continue the search // in the same direction. else if (sumc(vF2)>sumc(vF1)){ // println("Case 3. (sumc(vF2)>sumc(vF1) ",double(sumc(vF2))," ",double(sumc(vF1))," dLambda=",dLambda); bF1=bF2; vF1=vF2; } } // vF1 and mG1 now contains the function values evaluated at vX1=vX0+dLambda*vS

// println("/* STEP 3: UPDATE PARAMETER VECTOR */"); vX1 = vX0 + dLambda*vS;

// println("/* STEP 4: TEST STRONG CONVERGENCE AT vX0, rel grad.<eps1 and rel. diff. in vX<eps2 */"); decl vpp = vX0 .? 1.0 / sqr(sizerc(vX0)) + fabs(vX0) .: 1; fConv = (fabs(vG0' .* vpp) < dEps1 && fabs((vX1-vX0) ./ vpp) < 10 * dEps1) ? MAX_CONV : MAX_NOCONV;

if (iPrint > 0 && imod(i, iPrint) == 0){ MaxMonitor("BHHH", i, vX0, vG0', vS, double(sumc(vF0)), double(sumc(vF0)), double(dLambda), 0, fConv, bCompact); // println("Iteration timer = ",timespan(dIterStartTime)); }

if (fConv==MAX_CONV || !fStepOK){ break; } else if (i==iMaxIter-1){ return MAX_MAXIT; } else{ // Update old parameter and function values vX0 = vX1; vF0 = vF1; if (fNumDer==FALSE) mG0=mG1; } } // println("// Update input arguments"); avX[0] = vX0; adFunc[0] = double(sumc(vF0)); mHi = invertsym(mH); amInvHess[0] = mHi ? mHi : M_NAN*mH;

if (fConv==MAX_CONV) return MAX_CONV; else{ decl vpp = vX0 .? 1.0 / sqr(sizerc(vX0)) + fabs(vX0) .: 1; fConv = fabs(vG0' .* vpp) < dEps2 ? MAX_WEAK_CONV : MAX_NOCONV; return fConv; } }

[此贴子已经被作者于2005-8-14 18:03:54编辑过]