Colin Champion
Summary: this page contains three C/C++ functions for multidimensional optimisation. Two of them are straight conversions of Michael Powell’s Fortran whereas quadprog is a simple implementation of an algorithm due to Philip Wolfe.
| function | description | |
doc : source : mdplib source : references : test : download
The prototypes are
double uoamin(double (*func)(double*),double *x,int n,
double xinc,double tol,int maxiter) ;
double uoamax(double (*func)(double*),double *x,int n,
double xinc,double tol,int maxiter) ;
int uoacond() ;
where
Choose between uoamin and uoamax depending on whether you want to minimise or maximise f.
On return x will hold the approximate location of the minimum/maximum and the return value will be the value of f as computed at that point. The algorithm terminates when one of the three following conditions is met:
A call to uoacond() returns the value of cond corresponding to the most recent call to newuoa.
If you’re looking for reliable software you’d be better off running Attractive Chaos’s C++ translation of Powell’s original Fortran code (or Steven Johnson’s), but if you want to look inside you may find my code more congenial.
My version used Attractive Chaos’s as a starting point. I tried to convert it into the sort of code I’d have written myself. At each step I checked correctness using the test program below. Unfortunately the rewriting was error-prone and I occasionally introduced bugs which weren’t shown up by the test, but which became evident as the code became clearer and faulty logic started to show through. I can’t be confident of having removed all the bugs I introduced, which is why I can’t recommend my code to anyone who wants something he or she can trust.
At first I kept my output bit-identical to Attractive Chaos’s on my test file: this required me to dispense with a couple of minor code reorganisations. Later, when I made a version of lincoa I found some duplicated code. I replaced this by shared code but this led to small changes in the output for both programs. At this point I allowed myself the minor code reorganisations.
Incidentally I changed every line of code while working through it.
• uoacond • biglag • bigden • sethd • trsapp • setrhodelta • setknewdistsq • uoamin • myfunc • uoamax
/* This is NEWUOA for unconstrained minimization. The codes were written by Powell in Fortran and then translated to C with f2c. I further modified the code to make it independent of libf2c and f2c.h. Please refer to "The NEWUOA software for unconstrained optimization without derivatives", which is available at www.damtp.cam.ac.uk, for more information. */ /* The original fortran codes are distributed without restrictions. The C++ codes are distributed under MIT license. */ /* The MIT License Copyright (c) 2004, by M.J.D. Powell <mjdp@cam.ac.uk> 2008, by Attractive Chaos <attractivechaos@aol.co.uk> 2015, by Colin Champion <colin.champion@masterlyinactivity.com> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include "memory.h" #include <math.h> #define pi 3.141592653589793 static int cond = -1 ; int uoacond() { return cond ; } int update(int,int,double **,double **,int,double *,double,int *) ; void shiftxbase(int n,int npt,double **xp,double *xopt,double *vlag,double **b, double **z,double *pq,double *hq,int idz) ; double calcvlagbeta(int n,int npt,double *vlag,double **z,double **b,double *w, double *xopt,double *d,double beta,double dsq,double xsq, int idz) ; int chooseknew(int n,int npt,double beta,double **z,double *vlag,double **xp, double *xopt,double detrat,double rhosq,int ktemp,int idz) ; double genpqw(int n,int npt,double *pqw,double **z,int knew,int idz) ; void updatehq(int n,int npt,double *hq,double pqknew,double **xpt,int knew) ; void updategopt(int n,double *gopt,double *hq,double *step) ; void updategopt2(int n,int npt,double *gopt,double *xopt,double *pq,double **) ; /* -------------------------------------------------------------------------- */ static double biglag(int n,int npt,double *xopt,double **xp,double **b, double **z,int idz,int knew,double delta,double *d) { /* N is the number of variables. NPT is the number of interpolation * equations. XOPT is the best interpolation point so far. XPT * contains the coordinates of the current interpolation * points. BMAT provides the last N columns of H. ZMAT and IDZ give * a factorization of the first NPT by NPT submatrix of H. NDIM is * the first dimension of BMAT and has the value NPT+N. KNEW is the * index of the interpolation point that is going to be moved. DELTA * is the current trust region bound. D will be set to the step from * XOPT to the new point. ALPHA will be set to the KNEW-th diagonal * element of the H matrix. HCOL, GC, GD, S and W will be used for * working space. */ /* The step D is calculated in a way that attempts to maximize the * modulus of LFUNC(XOPT+D), subject to the bound ||D|| <= DELTA, * where LFUNC is the KNEW-th Lagrange function. */ int i,j,k,iterc,isave ; double sp,ss,cf1,cf2,cf3,cf4,cf5,dhd,cth,tau,sth,sum,temp,step,angle,scale ; double denom,delsq,tempa,tempb,taubeg,tauold,taumax,dd,gg,alpha ; double *hcol=vector(npt),*w=vector(n),*s=vector(n),*gd=vector(n) ; double *gc=vector(n) ; delsq = delta * delta ; /* Set the first NPT components of HCOL to the leading elements of * the KNEW-th column of H. */ alpha = genpqw(n,npt,hcol,z,knew,idz) ; /* Set the unscaled initial direction D. Form the gradient of LFUNC * atXOPT, and multiply D by the second derivative matrix of LFUNC. */ for(dd=i=0;i<n;i++) { d[i] = xp[knew][i] - xopt[i] ; gc[i] = b[i][knew] ; dd += d[i]*d[i] ; } for(k=0;k<npt;k++) { for(temp=sum=j=0;j<n;j++) { temp += xp[k][j] * xopt[j] ; sum += xp[k][j] * d[j] ; } temp *= hcol[k] ; sum *= hcol[k] ; for(i=0;i<n;i++) { gc[i] += temp * xp[k][i] ; gd[i] += sum * xp[k][i] ; } } /* Scale D and GD, with a sign change if required. Set S to another * vector in the initial two dimensional subspace. */ for(gg=sp=dhd=i=0;i<n;i++) { gg += gc[i] * gc[i] ; sp += d[i] * gc[i] ; dhd += d[i] * gd[i] ; } scale = delta / sqrt(dd) ; if(sp*dhd<0) scale = -scale ; if(sp*sp>dd*.99*gg) temp = 1 ; else temp = 0 ; tau = scale * (fabs(sp)+0.5*scale*fabs(dhd)) ; if(gg*delsq<tau*.01*tau) temp = 1 ; for(i=0;i<n;i++) { d[i] *= scale ; gd[i] *= scale ; s[i] = gc[i] + temp*gd[i] ; } /* Begin the iteration by overwriting S with a vector that has the * required length and direction, except that termination occurs if * the given D and S are nearly parallel. */ for(iterc=0;iterc<n;iterc++) { for(dd=sp=ss=i=0;i<n;i++) { dd += d[i] * d[i] ; sp += d[i] * s[i] ; ss += s[i] * s[i] ; } temp = dd*ss - sp*sp ; if(temp<=dd*1e-8*ss) break ; denom = sqrt(temp); for(i=0;i<n;i++) { s[i] = (dd*s[i]-sp*d[i]) / denom ; w[i] = 0 ; } /* Calculate the coefficients of the objective function on the * circle, beginning with the multiplication of S by the second * derivative matrix. */ for(k=0;k<npt;k++) { for(sum=j=0;j<n;j++) sum += xp[k][j] * s[j] ; sum *= hcol[k] ; for(j=0;j<n;j++) w[j] += sum * xp[k][j] ; } for(cf1=cf2=cf3=cf4=cf5=i=0;i<n;i++) { cf1 += s[i] * w[i]; cf2 += d[i] * gc[i]; cf3 += s[i] * gc[i]; cf4 += d[i] * gd[i]; cf5 += s[i] * gd[i]; } cf1 /= 2 ; cf4 = cf4/2 - cf1 ; /* Seek the value of the angle that maximizes the modulus of TAU. */ taubeg = cf1 + cf2 + cf4 ; taumax = tauold = taubeg; temp = 2 * pi / 50.0 ; for(isave=0,i=1;i<50;tauold=tau,i++) { angle = (double) i*temp; cth = cos(angle) ; sth = sin(angle) ; tau = cf1 + (cf2 + cf4*cth) * cth + (cf3 + cf5*cth) * sth ; if(fabs(tau)>fabs(taumax)) { taumax = tau ; isave = i ; tempa = tauold ; } else if(i==isave+1) tempb = tau ; } if(isave==0) tempa = tau ; else if(isave==49) tempb = taubeg ; step = 0 ; if(tempa!=tempb) { tempa -= taumax ; tempb -= taumax ; step = 0.5 * (tempa-tempb) / (tempa+tempb) ; } angle = temp * (isave+step) ; /* Calculate the new D and GD. Then test for convergence. */ cth = cos(angle) ; sth = sin(angle) ; tau = cf1 + (cf2 + cf4*cth) * cth + (cf3 + cf5*cth) * sth ; for(i=0;i<n;i++) { d[i] = cth*d[i] + sth*s[i] ; gd[i] = cth*gd[i] + sth*w[i] ; s[i] = gc[i] + gd[i] ; } if(fabs(tau)<=fabs(taubeg)*1.1) break ; } free(hcol,w,s,gd,gc) ; return alpha ; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ static double bigden(int n,int npt,double *xopt,double **xp,double **b, double **z,int idz,int kopt,int knew, double *d,double *w, double *vlag) { /* N is the number of variables. * NPT is the number of interpolation equations. * XOPT is the best interpolation point so far. * XPT contains the coordinates of the current interpolation points. * BMAT provides the last N columns of H. * ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H. * NDIMis the first dimension of BMAT and has the value NPT+N. * KOPT is the index of the optimal interpolation point. * KNEW is the index of the interpolation point that is going to be moved. * D will be set to the step from XOPT to the new point, and on entry it * should be the D that was calculated by the last call of BIGLAG. The length of the initial D provides a trust region bound on the final D. * W will be set to Wcheck for the final choice of D. * VLAG will be set to Theta*Wcheck+e_b for the final choice of D. * D is calculated in a way that should provide a denominator with a large modulus in the updating formula when the KNEW-th interpolation point is shifted to the new position XOPT+D. * the RETURN value beta will be found in the updating formula when the KNEW-th interpolation point is moved to its new position. */ int i,j,k,isave,iterc,jc,nw,ksav,ndim=npt+n ; double dd,ds,ss,den[9],par[9],tau,sum,diff,temp,step,beta,alpha,angle ; double denex[9],tempa,tempb,tempc,ssden,dtest,xoptd,xopts,denold,denmax ; double densav,dstemp,sumold,sstemp,xoptsq ; double **prod=matrix(ndim,5),**wvec=matrix(ndim,5),*s=vector(n) ; /* Store the first NPT elements of the KNEW-th column of H in W(N+1) * to W(N+NPT). */ alpha = genpqw(n,npt,w+n,z,knew,idz) ; /* The initial search direction D is taken from the last call of * BIGBDLAG, and the initial S is set below, usually to the direction * from X_OPT to X_KNEW, but a different direction to an * interpolation point may be chosen, in order to prevent S from * being nearly parallel to D. */ for(dd=ds=ss=i=0;i<n;i++) { dd += d[i]*d[i] ; s[i] = xp[knew][i] - xopt[i] ; ds += d[i] * s[i] ; ss += s[i] * s[i] ; } if(ds*ds>dd*.99*ss) { ksav = knew ; dtest = ds * ds / ss ; for(k=0;k<npt;k++) if(k!=kopt) { for(dstemp=sstemp=i=0;i<n;i++) { diff = xp[k][i] - xopt[i] ; dstemp += d[i] * diff ; sstemp += diff*diff ; } if(dstemp*dstemp/sstemp<dtest) { ksav = k ; dtest = dstemp * dstemp / sstemp ; ds = dstemp ; ss = sstemp ; } } for(i=0;i<n;i++) s[i] = xp[ksav][i] - xopt[i] ; } ssden = dd*ss - ds*ds ; densav = 0 ; /* Begin the iteration by overwriting S with a vector that has the * required length and direction. */ /* ------------------------------------------------------------------------ */ for(iterc=1;;iterc++) { temp = 1.0 / sqrt(ssden) ; for(xoptsq=xoptd=xopts=i=0;i<n;i++) { s[i] = temp * (dd*s[i]-ds*d[i]) ; xoptd += xopt[i] * d[i] ; xopts += xopt[i] * s[i] ; xoptsq += xopt[i] * xopt[i] ; } /* Set the coefficients of the first 2 terms of BETA. */ tempa = 0.5 * xoptd * xoptd ; tempb = 0.5 * xopts * xopts ; den[0] = dd * (xoptsq + 0.5*dd) + tempa + tempb ; den[1] = 2 * xoptd * dd ; den[2] = 2 * xopts * dd ; den[3] = tempa - tempb ; den[4] = xoptd * xopts ; for(i=4;i<8;i++) den[i+1] = 0 ; /* Put the coefficients of Wcheck in WVEC. */ for(k=0;k<npt;k++) { for(tempa=tempb=tempc=i=0;i<n;i++) { tempa += xp[k][i] * d[i]; tempb += xp[k][i] * s[i]; tempc += xp[k][i] * xopt[i]; } wvec[k][0] = 0.25 * (tempa*tempa + tempb*tempb) ; wvec[k][1] = tempa * tempc ; wvec[k][2] = tempb * tempc ; wvec[k][3] = 0.25 * (tempa*tempa - tempb*tempb) ; wvec[k][4] = 0.5 * tempa * tempb ; } for(i=0;i<n;i++) { wvec[i+npt][0] = 0 ; wvec[i+npt][1] = d[i] ; wvec[i+npt][2] = s[i] ; wvec[i+npt][3] = 0 ; wvec[i+npt][4] = 0 ; } /* Put the coefficents of THETA*Wcheck in PROD. */ for(jc=0;jc<5;jc++) { if(jc==2||jc==3) nw = ndim ; else nw = npt ; for(k=0;k<npt;k++) prod[k][jc] = 0 ; for(j=0;j<npt-n-1;j++) { for(sum=k=0;k<npt;k++) sum += z[j][k] * wvec[k][jc] ; if(j<idz) sum = -sum ; for(k=0;k<npt;k++) prod[k][jc] += sum * z[j][k] ; } if(nw==ndim) for(k=0;k<npt;k++) { for(sum=j=0;j<n;j++) sum += b[j][k] * wvec[npt+j][jc] ; prod[k][jc] += sum ; } for(j=0;j<n;j++) { for(sum=i=0;i<nw;i++) sum += b[j][i] * wvec[i][jc] ; prod[npt+j][jc] = sum ; } } /* Include in DEN the part of BETA that depends on THETA. */ for(k=0;k<ndim;k++) { for(sum=i=0;i<5;i++) sum += ( par[i] = 0.5 * prod[k][i] * wvec[k][i] ) ; den[0] -= par[0] + sum ; tempa = prod[k][0] * wvec[k][1] + prod[k][1] * wvec[k][0] ; tempb = prod[k][1] * wvec[k][3] + prod[k][3] * wvec[k][1] ; tempc = prod[k][2] * wvec[k][4] + prod[k][4] * wvec[k][2] ; den[1] -= tempa + 0.5 * (tempb+tempc) ; den[5] -= 0.5 * (tempb-tempc) ; tempa = prod[k][0] * wvec[k][2] + prod[k][2] * wvec[k][0] ; tempb = prod[k][1] * wvec[k][4] + prod[k][4] * wvec[k][1] ; tempc = prod[k][2] * wvec[k][3] + prod[k][3] * wvec[k][2] ; den[2] -= tempa + 0.5 * (tempb-tempc) ; den[6] -= 0.5 * (tempb+tempc) ; tempa = prod[k][0] * wvec[k][3] + prod[k][3] * wvec[k][0] ; den[3] -= tempa + par[1] - par[2] ; tempa = prod[k][0] * wvec[k][4] + prod[k][4] * wvec[k][0] ; tempb = prod[k][1] * wvec[k][2] + prod[k][2] * wvec[k][1] ; den[4] -= tempa + 0.5 * tempb ; den[7] -= par[3] - par[4] ; tempa = prod[k][3] * wvec[k][4] + prod[k][4] * wvec[k][3] ; den[8] -= 0.5 * tempa ; } /* Extend DEN so that it holds all the coefficients of DENOM. */ for(sum=i=0;i<5;i++) { tempa = prod[knew][i] ; sum += ( par[i] = 0.5 * tempa * tempa ) ; } denex[0] = alpha*den[0] + par[0] + sum ; tempa = 2 * prod[knew][0] * prod[knew][1] ; tempb = prod[knew][1] * prod[knew][3] ; tempc = prod[knew][2] * prod[knew][4] ; denex[1] = alpha*den[1] + tempa + tempb + tempc ; denex[5] = alpha*den[5] + tempb - tempc ; tempa = 2 * prod[knew][0] * prod[knew][2] ; tempb = prod[knew][1] * prod[knew][4] ; tempc = prod[knew][2] * prod[knew][3] ; denex[2] = alpha*den[2] + tempa + tempb - tempc ; denex[6] = alpha*den[6] + tempb + tempc ; tempa = 2 * prod[knew][0] * prod[knew][3] ; denex[3] = alpha*den[3] + tempa + par[1] - par[2] ; tempa = 2 * prod[knew][0] * prod[knew][4] ; denex[4] = alpha*den[4] + tempa + prod[knew][1]*prod[knew][2] ; denex[7] = alpha*den[7] + par[3] - par[4] ; denex[8] = alpha*den[8] + prod[knew][3]*prod[knew][4] ; /* Seek the value of the angle that maximizes the modulus of DENOM. */ sum = denex[0] + denex[1] + denex[3] + denex[5] + denex[7]; denold = denmax = sum ; temp = 2 * pi / 50 ; par[0] = 1 ; for(isave=0,i=1;i<50;i++) { angle = i * temp ; par[1] = cos(angle) ; par[2] = sin(angle) ; for(j=4;j<9;j+=2) { par[j-1] = par[1]*par[j-3] - par[2]*par[j-2] ; par[j] = par[1]*par[j-2] + par[2]*par[j-3] ; } sumold = sum ; for(sum=j=0;j<9;j++) sum += denex[j]*par[j] ; if(fabs(sum)>fabs(denmax)) { denmax = sum ; isave = i ; tempa = sumold ; } else if(i==isave+1) tempb = sum ; } if(isave==0) tempa = sum ; else if(isave==49) tempb = denold ; step = 0 ; if(tempa!=tempb) { tempa -= denmax ; tempb -= denmax ; step = 0.5 * (tempa-tempb) / (tempa+tempb) ; } /* Calculate the new parameters of the denominator, the new VLAG * vector and the new D. Then test for convergence. */ angle = temp * (isave+step) ; par[1] = cos(angle) ; par[2] = sin(angle) ; for(j=4;j<9;j+=2) { par[j-1] = par[1]*par[j-3] - par[2]*par[j-2] ; par[j] = par[1]*par[j-2] + par[2]*par[j-3] ; } for(denmax=beta=j=0;j<9;j++) { beta += den[j]*par[j] ; denmax += denex[j]*par[j] ; } for(k=0;k<ndim;k++) for(vlag[k]=j=0;j<5;j++) vlag[k] += prod[k][j] * par[j] ; tau = vlag[knew] ; for(dd=tempa=tempb=i=0;i<n;i++) { d[i] = par[1]*d[i] + par[2]*s[i] ; w[i] = xopt[i] + d[i] ; dd += d[i] * d[i] ; tempa += d[i] * w[i] ; tempb += w[i] * w[i] ; } if(iterc>=n) break ; if(iterc>1) densav = fmax(densav, denold) ; if(fabs(denmax)<=fabs(densav)*1.1) break ; densav = denmax ; /* Set S to 0.5 the gradient of the denominator with respect to * D. Then branch for the next iteration. */ for(i=0;i<n;i++) { temp = tempa*xopt[i] + tempb*d[i] - vlag[npt+i] ; s[i] = tau*b[i][knew] + alpha*temp ; } for(k=0;k<npt;k++) { for(sum=j=0;j<n;j++) sum += xp[k][j] * w[j] ; temp = (tau*w[n+k]-alpha*vlag[k]) * sum ; for(i=0;i<n;i++) s[i] += temp * xp[k][i] ; } for(ss=ds=i=0;i<n;i++) { ss += s[i]*s[i] ; ds += d[i]*s[i] ; } ssden = dd*ss - ds*ds ; if(ssden<dd*1e-8*ss) break ; } /* Set the vector W before the RETURN from the subroutine. */ for(k=0;k<ndim;k++) for(w[k]=j=0;j<5;j++) w[k] += wvec[k][j] * par[j] ; vlag[kopt] += 1 ; freematrix(prod,wvec) ; free(s) ; return beta ; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* The following instructions set the vector HD to the vector D multiplied by the second derivative matrix of Q. */ static void sethd(double *hd,int n,double **xp,double *d,double *hq,double *pq,int npt) { int i,j,k ; double temp ; for(i=0;i<n;i++) hd[i] = 0 ; updategopt(n,hd,hq,d) ; updategopt2(n,npt,hd,d,pq,xp) ; } /* -------------------------------------------------------------------------- */ static double trsapp(int n, int npt, double *xopt, double **xp, double *gq, double *hq, double *pq,double delsq, double *step) { /* The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meanings, * in order to define the current quadratic model Q. * DELTA is the trust region radius, and has to be positive. STEP * will be set to the calculated trial step. CRVMIN will be returned as the * least curvature of H aint the conjugate directions that occur, * except that it is set to 0 ifSTEP goes all the way to the trust * region boundary. The calculation of STEP begins with the * truncated conjugate gradient method. If the boundary of the trust * region is reached, then further changes to STEP may be made, each * one being in the 2D space spanned by the current STEP and the * corresponding gradient of Q. Thus STEP should provide a * substantial reduction to Q within the trust region. */ int i,j,iterc,isave,itermax ; double dd,cf,dg,gg,ds,sg,ss,dhd,dhs,cth,sgk,shs,sth,qadd,qbeg,qred,qmin ; double temp,qsav,qnew,ggbeg,alpha,angle,reduc,ggsav,tempa,tempb,bstep ; double ratio,angtest,crvmin ; double *d=vector(n),*g=vector(n),*hd=vector(n),*hs=vector(n) ; itermax = n ; for(i=0;i<n;i++) d[i] = xopt[i] ; sethd(hd,n,xp,d,hq,pq,npt) ; /* Prepare for the first line search. */ for(qred=dd=i=0;i<n;i++) { step[i] = hs[i] = 0 ; g[i] = gq[i] + hd[i] ; d[i] = -g[i] ; dd += d[i]*d[i] ; } if(dd==0) { free(d,g,hd,hs) ; return 0 ; } ggbeg = gg = dd ; /* ------------------------------------------------------------------------ */ /* -- Calculate the step to the trust region boundary and the product HD -- */ for(ds=ss=0,iterc=1;ss<delsq;iterc++) { temp = delsq - ss; bstep = temp / (ds + sqrt(ds * ds + dd * temp)); sethd(hd,n,xp,d,hq,pq,npt) ; for(dhd=j=0;j<n;j++) dhd += d[j] * hd[j] ; /* Update CRVMIN and set the step-length ALPHA. */ alpha = bstep ; if(dhd>0) { if(iterc==1) crvmin = dhd/dd ; else crvmin = fmin(crvmin,dhd/dd) ; alpha = fmin(alpha,gg/dhd) ; } qred += ( qadd = alpha * (gg - 0.5*alpha*dhd) ) ; /* Update STEP and HS. */ for(ggsav=gg,gg=i=0;i<n;i++) { step[i] += alpha * d[i] ; hs[i] += alpha * hd[i] ; temp = g[i] + hs[i] ; gg += temp * temp ; } /* Begin another conjugate direction iteration if required. */ if(alpha>=bstep) break ; ds = 0 ; if(qadd>qred*.01&&gg>ggbeg*1e-4&&iterc<itermax) { temp = gg / ggsav ; for(dd=ss=i=0;i<n;i++) { d[i] = temp*d[i] - g[i] - hs[i] ; dd += d[i]*d[i] ; ds += d[i]*step[i] ; ss += step[i]*step[i] ; } } if(ds<=0) { free(d,g,hd,hs) ; return crvmin ; } } /* ------------------------------------------------------------------------ */ /* ------------ Test whether an alternative iteration is required --------- */ for(iterc++;gg>ggbeg*1e-4;iterc++) { for(sg=shs=i=0;i<n;i++) { sg += step[i]*g[i] ; shs += step[i]*hs[i] ; } sgk = sg + shs ; angtest = sgk / sqrt(gg*delsq) ; if(angtest<=-.99) break ; /* Begin the alternative iteration by calculating D and HD and some * scalar products. */ temp = sqrt(delsq*gg - sgk*sgk) ; tempa = delsq / temp ; tempb = sgk / temp ; for(i=0;i<n;i++) d[i] = tempa*(g[i]+hs[i]) - tempb*step[i] ; sethd(hd,n,xp,d,hq,pq,npt) ; for(dg=dhd=dhs=i=0;i<n;i++) { dg += d[i]*g[i] ; dhd += hd[i]*d[i] ; dhs += hd[i]*step[i] ; } /* Seek the value of the angle that minimizes Q. */ cf = 0.5 * (shs-dhd) ; qbeg = sg + cf ; qsav = qmin = qbeg ; temp = 2 * pi / 50 ; for(isave=0,i=1;i<50;i++) { angle = i * temp ; cth = cos(angle) ; sth = sin(angle) ; qnew = (sg+cf*cth)*cth + (dg+dhs*cth)*sth ; if(qnew<qmin) { qmin = qnew ; isave = i ; tempa = qsav ; } else if(i==isave+1) tempb = qnew ; qsav = qnew ; } if(isave==0) tempa = qnew ; else if(isave==49) tempb = qbeg ; angle = 0 ; if(tempa!=tempb) { tempa -= qmin ; tempb -= qmin ; angle = 0.5 * (tempa-tempb) / (tempa+tempb) ; } angle = temp * (isave+angle) ; /* Calculate the new STEP and HS. Then test for convergence. */ cth = cos(angle) ; sth = sin(angle) ; reduc = qbeg - (sg + cf * cth) * cth - (dg + dhs * cth) * sth ; for(gg=i=0;i<n;i++) { step[i] = cth*step[i] + sth*d[i] ; hs[i] = cth*hs[i] + sth*hd[i] ; temp = g[i] + hs[i] ; gg += temp*temp ; } qred += reduc ; ratio = reduc / qred ; if(iterc>=itermax||ratio<=0.01) break ; } free(d,g,hd,hs) ; return 0 ; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ static double setrhodelta(double *rho,double rhoend,double *delta) { double ratio=rho[0]/rhoend ; delta[0] = rho[0] / 2 ; if(ratio<=16) rho[0] = rhoend ; else if(ratio<=250) rho[0] = sqrt(ratio) * rhoend ; else rho[0] *= 0.1 ; delta[0] = fmax(delta[0],rho[0]) ; return ratio ; } static double setknewdistsq (int n,int npt,double delta,int *knew,double *xopt,double **xp) { int j,k ; double sum,temp,distsq=4*delta*delta ; for(k=0;k<npt;k++) { for(sum=j=0;j<n;j++) { temp = xp[k][j] - xopt[j] ; sum += temp*temp ; } if(sum>distsq) { knew[0] = k ; distsq = sum ; } } return distsq ; } /* -------------------------------------------------------------------------- */ double uoamin(double (*func)(double*),double *x,int n, double rhobeg,double rhoend,int maxfun) { /* This subroutine seeks the least value of a function of many * variables, by a trust region method that forms quadratic models * by interpolation. There can be some freedom in the interpolation * conditions, which is taken up by minimizing the Frobenius norm of * the change to the second derivative of the quadratic model, * beginning with a zero matrix. The arguments of the subroutine are * as follows. */ /* N must be set to the number of variables and must be at least * two. NPT is the number of interpolation conditions. Its value * must be in the interval [N+2,(N+1)(N+2)/2]. Initial values of the * variables must be set in X(1),X(2),...,X(N). They will be changed * to the values that give the least calculated F. xinc and tol * must be set to the initial and final values of a trust region * radius, so both must be positive with tol<=xinc. Typically * xinc should be about one tenth of the greatest expected change * to a variable, and tol should indicate the accuracy that is * required in the final values of the variables. MAXFUN must be set * to an upper bound on the number of calls of CALFUN. */ /* SUBROUTINE func(x) must be provided by the user. */ /* XBASE will hold a shift of origin that should reduce the contributions * from rounding errors to values of the model and Lagrange functions. * XOPT will be set to the displacement from XBASE of the vector of variables that provides the least calculated F so far. * XNEW will be set to the displacement from XBASE of the vector of variables for the current calculation of F. */ /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ /* XPT will contain the interpolation point coordinates relative to XBASE. * FVAL will hold the values of F at the interpolation points. * GQ will hold the gradient of the quadratic model at XBASE. * HQ will hold the explicit second derivatives of the quadratic model. * PQ will contain the parameters of the implicit second derivatives of the quadratic model. * BMAT will hold the last N columns of H. * ZMAT will hold the factorization of the leading NPT by NPT submatrix * of H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, * where the elements of DZ are plus or minus 1.0, asspecified by IDZ. * NDIM is the first dimension of BMAT and has the value NPT+N. * D is reserved for trial steps from XOPT. * VLAG will contain the values of the Lagrange functions at a new point X. * They are part of a product that requires VLAG to be of length NDIM. */ int i,j,k,ih,nf,nh,nfm,idz,ipt,jpt,nfmm,knew,kopt,nptm,ksave,nfsav,itemp ; int ktemp,itest,npt=2*n+1 ; double f,dsq,rho,sum,fbeg,diff,beta,gisq,temp,suma,sumb,fopt,gqsq ; double xipt,xjpt,diffa,diffb,diffc,alpha,delta,recip,reciq,fsave ; double dnorm,ratio,dstep,vquad,detrat,crvmin,distsq,xoptsq ; double *xbase=vector(n),*xopt=vector(n),*xnew=vector(n),*fval=vector(npt) ; double *gq=vector(n),*hq=vector((n*(n+1))/2),*pq=vector(npt),*d=vector(n) ; double *vlag=vector(npt+n),*w=vector(2*npt),*pqw=vector(npt) ; double **z=matrix(npt-(n+1),npt),**b=matrix(n,npt+n),**xp=matrix(npt,n) ; nh = (n*(n+1)) / 2 ; nptm = npt - (n+1) ; if(maxfun<1) maxfun = 1 ; for(j=0;j<n;j++) xbase[j] = x[j] ; /* Begin the initialization procedure. NF becomes 1 more than the * number of function values so far. The coordinates of the * displacement of the next initial interpolation point from XBASE * are set in XPT(NF,.). */ recip = 1 / (rhobeg*rhobeg) ; reciq = sqrt(0.5) * recip ; /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ for(nf=0;nf<npt;nf++) { nfm = nf-1 ; nfmm = nfm-n ; if(nfm<2*n) // always happens when npt = 2*n+1 { if(nfm>=0&&nfm<n) xp[nf][nfm] = rhobeg ; else if(nfm>=n) xp[nf][nfmm] = -rhobeg ; } else // the case which never happens { itemp = nfmm / n ; jpt = nfm - itemp * n - n ; ipt = jpt + itemp ; if(ipt>=n) { ipt -= n ; swap(ipt,jpt) ; } xipt = rhobeg ; if(fval[ipt+n+1] < fval[ipt+1]) xipt = -xipt ; xjpt = rhobeg ; if(fval[jpt+n+1] < fval[jpt+1]) xjpt = -xjpt ; xp[nf][ipt] = xipt ; xp[nf][jpt] = xjpt ; } /* Calculate the next value of F. The least function value so far and its * index are required. */ for(j=0;j<n;j++) x[j] = xp[nf][j] + xbase[j]; fval[nf] = f = func(x) ; if(nf==0) { fbeg = fopt = f ; kopt = 0 ; } else if(f<fopt) { fopt = f ; kopt = nf ; } /* Set the non0 initial elements of BMAT and the quadratic model * in the cases when NF is at most 2*N+1. */ if(nfm<2*n) { if(nfm>=0&&nfm<n) { gq[nfm] = (f-fbeg) / rhobeg ; if(npt<=nf+n) { b[nfm][0] = -1/rhobeg ; b[nfm][nf] = 1/rhobeg ; b[nfm][npt+nfm-1] = -(rhobeg*rhobeg)/2 ; } } else if(nfm>=n) { b[nfmm][nf-n] = 1/(2*rhobeg) ; b[nfmm][nf] = -1/(2*rhobeg) ; z[nfmm][0] = -(reciq+reciq) ; z[nfmm][nf] = z[nfmm][nf-n] = reciq ; ih = (nfmm * (nfmm+3)) / 2 ; temp = (fbeg-f) / rhobeg ; hq[ih] = (gq[nfmm]-temp) / rhobeg; gq[nfmm] = .5 * (gq[nfmm] + temp); } } /* Set the off-diagonal second derivatives of the Lagrange * functions and the initial quadratic model. */ else { ih = (ipt*(ipt+1))/2 + jpt ; if(xipt<0) ipt += n ; if(xjpt<0) jpt += n ; z[nfmm][0] = z[nfmm][nf] = recip ; z[nfmm][ipt] = z[nfmm][jpt] = -recip ; hq[ih] = (fbeg - fval[ipt+1] - fval[jpt+1] + f) / (xipt*xjpt) ; } } /* ------------------------------------------------------------------------ */ /* - Begin the iterative procedure, because the initial model is complete - */ diffa = diffb = itest = idz = 0 ; for(i=0;i<n;i++) xopt[i] = xp[kopt][i] ; /* Generate the next trust region step and test its length. Set KNEW * to -1 ifthe purpose of the next F will be to improve the * model. */ /* ------------------------------------------------------------------------ */ for(knew=-1,delta=rho=rhobeg,nfsav=nf;;) { if(knew<0) for(;;) { crvmin = trsapp(n,npt,xopt,xp,gq,hq,pq,delta*delta,d) ; for(dsq=i=0;i<n;i++) dsq += d[i]*d[i] ; dnorm = fmin(delta,sqrt(dsq)) ; if(dnorm>=rho/2) break ; // knew must be -1 delta /= 10 ; ratio = -1 ; if(delta<=rho*1.5) delta = rho ; temp = crvmin * .125 * rho * rho ; if(nf<=nfsav+2||temp<=diffa||temp<=diffb||temp<=diffc) { distsq = setknewdistsq(n,npt,delta,&knew,xopt,xp) ; /* If KNEW is positive, then set DSTEP, and branch back for the next * iteration, which will generate a "model step". */ if(knew>=0) { dstep = fmax(rho,fmin(sqrt(distsq)/10,delta/2)) ; dsq = dstep*dstep ; break ; } else if(delta>rho||dnorm>rho) continue ; } /* Return from the calculation, after another Newton-Raphson step, * if it is too short to have been tried before. */ if(rho<=rhoend) { for(i=0;i<n;i++) { xnew[i] = xopt[i] + d[i] ; x[i] = xbase[i] + xnew[i] ; } f = func(x) ; knew = -2 ; break ; } /* The calculations with the current value of RHO are complete. Pick * the next values of RHO and DELTA. */ ratio = setrhodelta(&rho,rhoend,&delta) ; nfsav = nf ; } if(knew==-2) { cond = 0 ; break ; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ /* Shift XBASE if XOPT may be too far from XBASE. First make the * changes to BMAT that do not depend on ZMAT. */ for(xoptsq=i=0;i<n;i++) xoptsq += xopt[i] * xopt[i] ; if(dsq<=xoptsq/1000) { updategopt(n,gq,hq,xopt) ; updategopt2(n,npt,gq,xopt,pq,xp) ; shiftxbase(n,npt,xp,xopt,vlag,b,z,pq,hq,idz) ; for(j=0;j<n;j++) { xbase[j] += xopt[j] ; xopt[j] = 0 ; } xoptsq = 0 ; } /* Pick the model step ifKNEW is positive. A different choice of D * may be made later, ifthe choice of D by BIGLAG causes * substantial cancellation in DENOM. */ if(knew>=0) alpha = biglag(n,npt,xopt,xp,b,z,idz,knew,dstep,d) ; /* Calculate VLAG and BETA for the current choice of D. The first * NPT components of W_check will be held in W. */ for(k=0;k<npt;k++) { for(vlag[k]=suma=sumb=j=0;j<n;j++) { suma += xp[k][j] * d[j] ; sumb += xp[k][j] * xopt[j] ; vlag[k] += b[j][k] * d[j] ; } w[k] = suma * (suma/2+sumb) ; } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ beta = calcvlagbeta(n,npt,vlag,z,b,w,xopt,d,beta,dsq,xoptsq,idz) ; vlag[kopt] += 1 ; /* If KNEW is positive and if the cancellation in DENOM is unacceptable, * then BIGDEN calculates an alternative model step. */ if(knew>=0&&fabs(1 + alpha*beta / (vlag[knew]*vlag[knew]))<=0.8) beta = bigden(n,npt,xopt,xp,b,z,idz,kopt,knew,d,w,vlag) ; /* Calculate the next value of the objective function. */ if(nf+1>=maxfun) { cond = 2 ; break ; } for(i=0;i<n;i++) { xnew[i] = xopt[i] + d[i] ; x[i] = xbase[i] + xnew[i] ; } nf += 1 ; f = func(x) ; /* Use the quadratic model to predict the change in F due to the * step D, and set DIFF to the error of this prediction. */ for(vquad=ih=j=0;j<n;j++) { vquad += d[j] * gq[j] ; for(i=0;i<=j;ih++,i++) { temp = d[i]*xnew[j] + d[j]*xopt[i] ; if(i==j) temp /= 2 ; vquad += temp * hq[ih] ; } } for(k=0;k<npt;k++) vquad += pq[k] * w[k] ; diff = f - fopt - vquad ; diffc = diffb ; diffb = diffa ; diffa = fabs(diff) ; if(dnorm>rho) nfsav = nf ; /* Update FOPT and XOPT ifthe new F is the least value of the * objective function so far. The branch when KNEW is positive * occurs ifD is not a trust region step. */ fsave = fopt ; if(f<fopt) { fopt = f ; for(i=0;i<n;i++) xopt[i] = xnew[i] ; } ksave = knew ; if(knew<0) { /* Pick the next value of DELTA after a trust region step. */ if(vquad>=0) { cond = 1 ; break ; } ratio = (f-fsave) / vquad ; if(ratio<=0.1) delta = dnorm/2 ; else if(ratio<=.7) delta = fmax(delta/2,dnorm) ; else delta = fmax(delta/2,2*dnorm) ; if(delta<=rho*1.5) delta = rho ; /* Set KNEW to the index of the next interpolation point to be * deleted. */ temp = fmax(delta/10,rho) ; if(f>=fsave) { ktemp = kopt ; detrat = 1 ; } else { ktemp = -1 ; detrat = 0 ; } knew = chooseknew(n,npt,beta,z,vlag,xp,xopt,detrat,temp*temp,ktemp,idz) ; } /* Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation * point can be moved. Begin the updating of the quadratic model, * starting with the explicit second derivative term. */ if(knew>=0) idz = update(n,npt,b,z,idz,vlag,beta,&knew) ; if(knew>=0) { /* Update the other second derivative parameters, and then the * gradient vector of the model. Also include the new interpolation * point. */ updatehq(n,npt,hq,pq[knew],xp,knew) ; pq[knew] = 0 ; genpqw(n,npt,pqw,z,knew,idz) ; for(k=0;k<npt;k++) pq[k] += diff * pqw[k] ; fval[knew] = f ; for(gqsq=i=0;i<n;i++) { gq[i] += diff * b[i][knew]; gqsq += gq[i] * gq[i]; xp[knew][i] = xnew[i]; } /* If a trust region step makes a small change to the objective * function, then calculate the gradient of the least Frobenius norm * interpolant at XBASE, and store it in W, using VLAG for a vector * of right hand sides. */ if(ksave==-1&&delta==rho&&fabs(ratio)>.01) itest = 0 ; else if(ksave==-1&&delta==rho) { for(k=0;k<npt;k++) vlag[k] = fval[k] - fval[kopt] ; for(gisq=i=0;i<n;i++) { for(w[i]=k=0;k<npt;k++) w[i] += b[i][k] * vlag[k] ; gisq += w[i]*w[i] ; } itest += 1 ; if(gqsq<gisq*100) itest = 0 ; if(itest>=3) { for(i=0;i<n;i++) gq[i] = w[i] ; for(ih=0;ih<nh;ih++) hq[ih] = 0 ; for(j=0;j<nptm;j++) { for(w[j]=k=0;k<npt;k++) w[j] += vlag[k] * z[j][k] ; if(j<idz) w[j] = -w[j] ; } for(k=0;k<npt;k++) for(pq[k]=j=0;j<nptm;j++) pq[k] += z[j][k] * w[j]; itest = 0 ; } } if(f<fsave) kopt = knew ; knew = -1 ; /* If a trust region step has provided a sufficient decrease in F, * then branch for another trust region calculation. The case * KSAVE>0 occurs when the new function value was calculated by a * model step. */ if(f<=fsave+vquad/10||ksave>=0) continue ; /* Alternatively, find out ifthe interpolation points are close * enough to the best point so far. */ } distsq = setknewdistsq(n,npt,delta,&knew,xopt,xp) ; /* If KNEW is positive, then set DSTEP, and branch back for the next * iteration, which will generate a "model step". */ if(knew>=0) { dstep = fmax(rho,fmin(sqrt(distsq)/10,delta/2)) ; dsq = dstep*dstep ; } else if(ratio<=0&&delta<=rho&&dnorm<=rho) /* The calculations with the current value of RHO are complete. Pick * the next values of RHO and DELTA. */ { if(rho<=rhoend) { cond = 0 ; break ; } ratio = setrhodelta(&rho,rhoend,&delta) ; nfsav = nf ; } } /* Return from the calculation, after another Newton-Raphson step, * if it is too short to have been tried before. */ if(fopt<=f) { for(i=0;i<n;i++) x[i] = xbase[i] + xopt[i] ; f = fopt ; } free(xbase,xopt,xnew,w,fval,gq) ; free(hq,pq,d,vlag,pqw) ; freematrix(z,b,xp) ; return f ; } /* -------------------------------------------------------------------------- */ static double(*f)(double*) ; static double myfunc(double *x) { return -f(x) ; } double uoamax(double (*func)(double*),double *x,int n, double xinc,double tol,int max_iter) { f = func ; return -uoamin(myfunc,x,n,xinc,tol,max_iter) ; }
• updategopt • updategopt2 • update • shiftxbase • calcvlagbeta • chooseknew • genpqw • updatehq
/* The original fortran codes are distributed without restrictions. The C++ codes are distributed under MIT license. */ /* The MIT License Copyright (c) 2004, by M.J.D. Powell <mjdp@cam.ac.uk> 2008, by Attractive Chaos <attractivechaos@aol.co.uk> 2015, by Colin Champion <colin.champion@masterlyinactivity.com> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include"memory.h" #include <math.h> #define max(a,b) ((a)>(b)?(a):(b)) #define min(a,b) ((a)<(b)?(a):(b)) void updategopt(int n,double *gopt,double *hq,double *step) { int i,j,ih ; for(ih=j=0;j<n;j++) for(i=0;i<=j;i++,ih++) { if(i<j) gopt[j] += hq[ih] * step[i] ; gopt[i] += hq[ih] * step[j] ; } } void updategopt2(int n,int npt,double *gopt,double *xopt,double *pq,double **xp) { int i,k ; double sum,temp ; for(k=0;k<npt;k++) { for(sum=i=0;i<n;i++) sum += xp[k][i] * xopt[i]; temp = pq[k] * sum ; for(i=0;i<n;i++) gopt[i] += temp * xp[k][i] ; } } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ int update(int n,int npt,double **b,double **z,int idz,double *vlag, double beta,int *kn) { /* The arrays BMAT and ZMAT with IDZ are updated, in order to shift * the interpolation point that has index KNEW. On entry, VLAG * contains the components of the vector Theta*Wcheck+e_b of the * updating formula (6.11), and BETA holds the value of the * parameter that has this name. */ int nptm=npt-n-1,i,j,ja,jb,jl,jp,iflag,knew=kn[0] ; double tau,temp,scala,scalb,alpha,denom,tempa,tempb,tausq,sqrtdn ; double *w=vector(npt+n) ; static int upcall=-1 ; upcall += 1 ; if(knew<0) throw up("knew=%d\n",knew) ; /* Apply the rotations that put zeros in the KNEW-th row of ZMAT. */ for(jl=0,j=1;j<npt-n-1;j++) { if(j==idz) jl = j ; else if(z[j][knew]!=0) { tempa = z[jl][knew]; tempb = z[j][knew]; temp = sqrt(tempa*tempa+tempb*tempb) ; tempa /= temp ; tempb /= temp ; for(i=0;i<npt;i++) { temp = tempa*z[jl][i] + tempb*z[j][i] ; z[j][i] = tempa*z[j][i] - tempb*z[jl][i] ; z[jl][i] = temp ; } z[j][knew] = 0; } } /* Put the first NPT components of the KNEW-th column of HLAG into * W, and calculate the parameters of the updating formula. */ if(idz<1) tempa = z[0][knew] ; else tempa = -z[0][knew] ; for(i=0;i<npt;i++) w[i] = tempa * z[0][i] ; if(jl>0) for(tempb=z[jl][knew],i=0;i<npt;i++) w[i] += tempb * z[jl][i] ; alpha = w[knew] ; tau = vlag[knew] ; tausq = tau * tau ; denom = alpha*beta + tausq ; vlag[knew] -= 1 ; if(denom==0) { kn[0] = -1 ; free(w) ; return idz ; } /* Complete the updating of ZMAT when there is only 1.0 nonzero * element in the KNEW-th row of the new matrix ZMAT, but, ifIFLAG * is set to 1.0, then the first column of ZMAT will be exchanged * with another 1.0 later. */ iflag = 0 ; sqrtdn = sqrt(fabs(denom)) ; /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ if(jl==0) { tempa = tau / sqrtdn ; tempb = z[0][knew] / sqrtdn ; for(i=0;i<npt;i++) z[0][i] = tempa*z[0][i] - tempb*vlag[i] ; if(denom<0) { if(idz==0) idz = 1 ; else iflag = 1 ; } } else { /* Complete the updating of ZMAT in the alternative case. */ if(beta>=0) ja = jl ; else ja = 0 ; jb = jl - ja ; temp = z[jb][knew] / denom ; tempa = temp * beta ; tempb = temp * tau ; temp = z[ja][knew] ; scala = 1 / sqrt(fabs(beta)*temp*temp+tausq) ; scalb = scala * sqrtdn ; for(i=0;i<npt;i++) { z[ja][i] = scala * (tau * z[ja][i] - temp * vlag[i]) ; z[jb][i] = scalb * (z[jb][i] - tempa * w[i] - tempb * vlag[i]) ; } if(denom<=0) { if(beta<0) idz += 1 ; else iflag = 1 ; } } if(denom==0) throw up("denom is 0") ; /* IDZ is reduced in the following case, and usually the first * column of ZMAT is exchanged with a later 1.0. */ if(iflag==1) for(idz-=1,i=0;i<npt;i++) swap(z[0][i],z[idz][i]) ; /* Finally, update the matrix BMAT. */ for(j=0;j<n;j++) { jp = npt + j ; w[jp] = b[j][knew] ; tempa = (alpha*vlag[jp] - tau*w[jp]) / denom ; tempb = -(beta*w[jp] + tau*vlag[jp]) / denom ; for(i=0;i<=jp;i++) { b[j][i] += tempa*vlag[i] + tempb*w[i] ; if(i>=npt) b[i-npt][jp] = b[j][i]; } } free(w) ; return idz ; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void shiftxbase(int n,int npt,double **xp,double *xopt,double *vlag,double **b, double **z,double *pq,double *hq,int idz) { int i,j,k,ih ; double tempq,temp,sum,xoptsq,sumz,*w=vector(2*npt) ; for(xoptsq=i=0;i<n;i++) xoptsq += xopt[i] * xopt[i] ; tempq = xoptsq/4 ; for(k=0;k<npt;k++) { for(sum=i=0;i<n;i++) sum += xp[k][i] * xopt[i]; sum -= xoptsq/2 ; w[npt+k] = sum ; for(i=0;i<n;i++) { xp[k][i] -= xopt[i]/2 ; vlag[i] = b[i][k]; w[i] = sum*xp[k][i] + tempq*xopt[i] ; for(j=0;j<=i;j++) b[j][npt+i] += vlag[i]*w[j] + w[i]*vlag[j] ; } } /* Then the revisions of BMAT that depend on ZMAT are calculated. */ for(k=0;k<npt-n-1;k++) { for(sumz=i=0;i<npt;i++) { sumz += z[k][i] ; w[i] = w[npt+i] * z[k][i] ; } for(j=0;j<n;j++) { sum = tempq * sumz * xopt[j]; for(i=0;i<npt;i++) sum += w[i] * xp[i][j]; vlag[j] = sum; if(k<idz) sum = -sum; for(i=0;i<npt;i++) b[j][i] += sum * z[k][i]; } for(i=0;i<n;i++) { if(k>=idz) temp = vlag[i] ; else temp = -vlag[i] ; for(j=0;j<=i;j++) b[j][npt+i] += temp * vlag[j] ; } } /* The following instructions complete the shift of XBASE, * including the changes to the parameters of the quadratic * model. */ for(ih=j=0;j<n;j++) { for(w[j]=k=0;k<npt;k++) { w[j] += pq[k] * xp[k][j] ; xp[k][j] -= xopt[j]/2 ; } for(i=0;i<=j;ih++,i++) { hq[ih] += w[i]*xopt[j] + xopt[i]*w[j] ; b[j][npt+i] = b[i][npt+j] ; } } free(w) ; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ //C Calculate VLAG and BETA for the current choice of STEP. The first NPT //C elements of VLAG are set to the values of the Lagrange functions at //C XPT(KOPT,.)+STEP(.). The first NPT components of W_check are held //C in W, where W_check is defined in a paper on the updating method. double calcvlagbeta(int n,int npt,double *vlag,double **z,double **b,double *w, double *xopt,double *d,double beta,double dsq,double xsq, int idz) { int i,j,k ; double sum,dx,bsum ; for(beta=k=0;k<npt-n-1;k++) { for(sum=i=0;i<npt;i++) sum += z[k][i] * w[i] ; if(k<idz) { beta += sum*sum ; sum = -sum ; } else beta -= sum*sum ; for(i=0;i<npt;i++) vlag[i] += sum * z[k][i]; } for(bsum=dx=j=0;j<n;j++) { for(sum=i=0;i<npt;i++) sum += w[i] * b[j][i]; bsum += sum * d[j] ; for(k=0;k<n;k++) sum += b[k][npt+j] * d[k]; vlag[npt+j] = sum ; bsum += sum * d[j] ; dx += d[j] * xopt[j] ; } // more numerical instability here? return dx*dx + dsq*((xsq+dx) + dx + dsq/2) + beta - bsum ; } /* -------------------------------------------------------------------------- */ int chooseknew(int n,int npt,double beta,double **z,double *vlag,double **xp, double *xopt,double detrat,double rhosq,int ktemp,int idz) { int i,j,k,knew ; double temp,qtemp,distsq,hdiag ; for(knew=-1,k=0;k<npt;k++) if(k!=ktemp) { for(hdiag=j=0;j<npt-n-1;j++) if(j<idz) hdiag -= z[j][k]*z[j][k] ; else hdiag += z[j][k]*z[j][k] ; temp = fabs(beta*hdiag+vlag[k]*vlag[k]) ; for(distsq=j=0;j<n;j++) { qtemp = xp[k][j] - xopt[j] ; distsq += qtemp*qtemp ; } if(distsq>=rhosq) { if(rhosq==0) temp *= distsq * distsq ; else { qtemp = distsq / rhosq ; temp *= qtemp*qtemp*qtemp ; } } if(temp>detrat) { detrat = temp ; knew = k ; } } return knew ; } /* -------------------------------------------------------------------------- */ double genpqw(int n,int npt,double *pqw,double **z,int knew,int idz) { int i,j,k ; double temp ; for(k=0;k<npt;k++) pqw[k] = 0 ; for(j=0;j<npt-n-1;j++) { if(j<idz) temp = -z[j][knew] ; else temp = z[j][knew] ; if(temp!=0) for(k=0;k<npt;k++) pqw[k] += temp * z[j][k] ; } return pqw[knew] ; } void updatehq(int n,int npt,double *hq,double pqknew,double **xpt,int knew) { int i,j,ih ; double temp ; for(ih=i=0;i<n;i++) { temp = pqknew * xpt[knew][i] ; for(j=0;j<=i;j++,ih++) hq[ih] += temp * xpt[knew][j] ; } }
The test program uoatest.c below calls uoamin to minimise the Chebyshev function used by Powell to illustrate his own Fortran code. It does this in 2, 4, 6, 8, 10 and 20 dimensions (Powell stopped at 8). The correct output is this:
2 31: 3.78846734e-19: 0.211325 0.788675 4 90: 3.52620610e-14: 0.102673 0.406204 0.593796 0.897327 6 193: 1.55535235e-12: 0.066876 0.288740 0.366683 0.633317 0.711260 ... 8 341: 3.51687373e-03: 0.043153 0.193091 0.266329 0.500000 0.500000 ... 10 563: 4.77271370e-03: 0.074748 0.171518 0.286434 0.359647 0.470750 ... 20 2418: 4.57295520e-03: 0.024600 0.070921 0.116574 0.176668 0.206808 ...
[Note that it is slightly changed since the bugfix of 22 Feb 2020.]
Attractive Chaos’s code gives the following:
2 31: 3.78847194e-19: 0.211325 0.788675 4 91: 2.19176132e-15: 0.102673 0.406204 0.593796 0.897327 6 190: 3.12702103e-14: 0.066877 0.288741 0.366682 0.633318 0.711259 ... 8 303: 3.51687373e-03: 0.043153 0.193091 0.266329 0.500000 0.500000 ... 10 668: 6.50395480e-03: 0.059620 0.166708 0.239170 0.398884 0.398885 ... 20 1934: 4.57295560e-03: 0.024600 0.070921 0.116573 0.176667 0.206807 ...
Notice that I perform noticeably more function evaluations for n=20 but get an appreciably smaller minimum for n=10.
Attractive Chaos’s results are not identical to Powell’s own, which are these:
At the return from NEWUOA Number of function values = 31
Least value of F = 3.788471046857957D-19 The corresponding X is:
2.113249D-01 7.886751D-01
At the return from NEWUOA Number of function values = 90
Least value of F = 3.526693206487107D-14 The corresponding X is:
1.026728D-01 4.062038D-01 5.937962D-01 8.973272D-01
At the return from NEWUOA Number of function values = 198
Least value of F = 4.343402133989936D-14 The corresponding X is:
6.687652D-02 2.887405D-01 3.666823D-01 6.333176D-01 7.112593D-01
9.331234D-01
At the return from NEWUOA Number of function values = 314
Least value of F = 3.516873725862449D-03 The corresponding X is:
4.315284D-02 1.930909D-01 2.663288D-01 5.000002D-01 4.999999D-01
7.336712D-01 8.069093D-01 9.568473D-01
Powell’s results are no better, just different, and the discrepancy isn’t surprising given the numerical delicacy of the algorithm. Powell remarked that “the results are still highly sensitive to computer rounding errors” and was in the habit of giving two results for each problem, differing through having the x-coordinates permuted (which makes no difference mathematically); his pairs of results typically differ more widely between themselves than ours do from his on the Chebyshev problem.
• initchebyfunc • closechebyfunc • chebyfunc • main
#include "memory.h" static int n=0,ncalls=0 ; static double **y=0 ; void initchebyfunc(int k) { n = k ; ncalls = 0 ; y = matrix(n+1,n) ; } int closechebyfunc() { int ret=ncalls ; freematrix(y) ; ncalls = n = 0 ; y = 0 ; return ret ; } double chebyfunc(double *x) { double q,qret ; int i,j ; for(j=0;j<n;j++) { y[0][j] = 1 ; y[1][j] = 2*x[j]-1 ; } for(i=2;i<=n;i++) for(j=0;j<n;j++) y[i][j] = 2*y[1][j]*y[i-1][j] - y[i-2][j] ; for(qret=i=0;i<=n;i++) { for(q=j=0;j<n;j++) q += y[i][j] ; q /= n ; if((i&1)==0) q += 1/((i-1.0)*(i+1.0)) ; qret += q*q ; } ncalls += 1 ; return qret ; } double uoamin(double (*)(double*),double*,int,double,double,int) ; double uoamax(double (*)(double*),double*,int,double,double,int) ; int main(int argc,char **argv) { int i,ind,n,ncalls,order[]={2,4,6,8,10,20,0,0} ; double *x,y ; for(ind=0;(n=order[ind]);ind++) { x = vector(n) ; for(i=0;i<n;i++) x[i] = (i+1.0) / (n+1.0) ; initchebyfunc(n) ; y = uoamin(chebyfunc,x,n,x[0]/5,1e-6,50000) ; ncalls = closechebyfunc() ; printf("%2d%6d: %.8e:",n,ncalls,y) ; for(i=0;i<n&&i<5;i++) printf(" %9.6f",x[i]) ; if(n>5) printf(" ...") ; printf("\n") ; free(x) ; } }
This link generates a perl script which you can download and execute to compile and run newuoa. It works in Firefox only (and needs popups enabled for this site).
newuoa.c mdplib.c uoatest.c utils: memory.h --- g++ -o uoatest (*.c) uoatest
doc : source : mdplib source : references : test : download
The prototypes are
double bfgsmin
(double *x,int n,double xinc,double (*f)(double *),double (*fdf)(double*,double*),
int (*decider)(int,double,double,double,double*,double*,double*,int)) ;
double bfgsmax
(double *x,int n,double step,double (*f)(double *),double (*fdf)(double*,double*),
int (*decider)(int,double,double,double,double*,double*,double*,int)) ;
where
Choose between bfgsmin and bfgsmax depending on whether you want to minimise or maximise f. The latter is just a wrapper for the former, and cannot be nested recursively.
On return x will hold the approximate location of the minimum/maximum and the return value will be the value of f as computed at that point. The decider function is invoked after each line search and returns a value to indicate whether further iterations should be performed (0 to continue, 1 to cease iterating).
The decider function should satisfy the prototype
int mydecide(int niter,double qorig,double qold,double q,
double *xold,double *x,double *grad,int n) ;
where
Note that bfgsmin does nothing to protect against an infinite sequence of line searches: you should do this yourself using the niter parameter. The following is a suitable simple decision function:
static int mydecide(int niter,double qorig,double qold,double q,
double *xold,double *x,double *grad,int n)
{ return niter>=500 || (qold-q)<q/200 ; }
The function bfgscond() returns a flag indicating the reason for termination on the most recent call. The meaning is as follows:
The case 3 may give cause for concern, since it is likely to arise if your differentiation is incorrect, or if the function is not differentiable at the current position.
• solvequad • wrap_f • wrap_df • interp_quad • check_extremum • interp_cubic • x3 • interpolate • minimise • bfgsmin • minusf • minusfdf • minuscriterion • bfgsmax
#include "memory.h" #include <math.h> # define EPS 2.2204460492503131e-16 int solvequad(double a,double b,double c,double *x0,double *x1) { double det,r,sgnb,temp,r1,r2 ; if(a==0) { if(b==0) return 0 ; else { *x0 = -c/b ; return 1 ; } } // linear if((det=b*b-4*a*c)>0) { if(b==0) { *x0 = -(*x1 = sqrt(-c/a)) ; return 2 ; } sgnb = (b > 0 ? 1 : -1); temp = -(b + sgnb * sqrt (det)) / 2 ; r1 = temp / a ; r2 = c / temp ; if(r1<r2) { *x0 = r1 ; *x1 = r2 ; } else { *x0 = r2 ; *x1 = r1 ; } return 2; } else if(det==0) { *x0 = *x1 = -b / (2*a) ; return 2 ; } else return 0; } static double wrap_f(double alpha,double *pt,double *x,double *p,int n,xy *key, double *xcache,double (*f)(double *)) { int i ; if(alpha==key->x) return key->y ; if(alpha!=xcache[0]) for(i=0;i<n;i++) { pt[i] = x[i] + alpha*p[i] ; } xcache[0] = key->x = alpha ; return key->y = (*f)(pt) ; } static xy wrap_df(double alpha,double *pt,double *x,double *p,double *g, int n,xy dkey,double *xcache,double *gkey, double (*df)(double *,double *)) { int i ; double q ; if(alpha==dkey.x) return dkey ; if(alpha!=xcache[0]) for(i=0;i<n;i++) { pt[i] = x[i] + alpha*p[i] ; } if(alpha!=gkey[0]) { (*df)(pt,g) ; gkey[0] = alpha ; } // calculate the slope for(q=i=0;i<n;i++) q += g[i] * p[i] ; return xy(xcache[0]=alpha,q) ; } /* Find a minimum in x=[0,1] of the interpolating quadratic through * (0,f0) (1,f1) with derivative fp0 at x=0. The interpolating * polynomial is q(x) = f0 + fp0 * z + (f1-f0-fp0) * z^2 */ static double interp_quad(double f0,double fp0,double f1,double zl,double zh) { double f , z , fh = f0 + zh*(fp0+zh*(f1-f0-fp0)) ; double c = 2*(f1-f0-fp0) ; /* curvature */ xy A = xy(zl,f0+zl*(fp0+zl*(f1-f0-fp0))) ; if(fh<A.y) A = xy(zh,f0+zh*(fp0+zh*(f1-f0-fp0))) ; if(c>0) /* positive curvature required for a minimum */ { z = -fp0 / c ; /* location of minimum */ if(z>zl&&z<zh) { f = f0 + z*(fp0+z*(f1-f0-fp0)) ; if(f<A.y) A = xy(z,f) ; } } return A.x ; } /* Find a minimum in x=[0,1] of the interpolating cubic through * (0,f0) (1,f1) with derivatives fp0 at x=0 and fp1 at x=1. * * The interpolating polynomial is: * * c(x) = f0 + fp0 * z + eta * z^2 + xi * z^3 * * where eta=3*(f1-f0)-2*fp0-fp1, xi=fp0+fp1-2*(f1-f0). */ static xy check_extremum(double c0,double c1,double c2,double c3,double z,xy A) { /* could make an early return by testing curvature >0 for minimum */ double y = c0 + z * (c1 + z * (c2 + z * c3)); if(y<A.y) return xy(z,y) ; else return A ; } static double interp_cubic(double c0,double c1,double f1,double fp1,double zl,double zh) { double c2=3*(f1-c0)-2*c1-fp1 , c3=c1+fp1-2*(f1-c0) , z0,z1 ; xy A = xy(zl,c0+zl*(c1+zl*(c2+zl*c3))) ; A = check_extremum(c0,c1,c2,c3,zh,A) ; if(solvequad(3*c3,2*c2,c1,&z0,&z1)==2&&z1>zl&&z1<zh) A = check_extremum(c0,c1,c2,c3,z1,A) ; else if(z0>zl&&z0<zh) A = check_extremum(c0,c1,c2,c3,z0,A) ; return A.x ; } struct x3 { double a,fa ; settable fpa ; x3() { a=fa=0 ; fpa = settable() ; } x3(double x,double y) { a = x ; fa = y ; fpa = settable() ; } x3(double x,double y,settable z) { a = x ; fa = y ; fpa = z ; } } ; static double interpolate(x3 a,x3 b,double xmin,double xmax) { /* Map [a,b] to [0,1] */ double z , q=b.a-a.a , zmin=(xmin-a.a)/q , zmax=(xmax-a.a)/q ; if(zmin>zmax) swap(zmin,zmax) ; if(!b.fpa.set) z = interp_quad(a.fa,a.fpa*q,b.fa,zmin,zmax) ; else z = interp_cubic(a.fa,a.fpa*q,b.fa,b.fpa*q,zmin,zmax) ; return a.a + z*q ; } /* -------------------------------------------------------------------------- */ static double minimise(double alpha1,double *x,double *p, double *xalpha,double *g,int n, double y,double gp,int *flag, double (*f)(double *),double (*df)(double *,double *)) { double rho=0.01,sigma=0.1,tau1=9,tau2=0.05,tau3=0.5 ;//recommended by Fletcher double falpha,fpalpha,delta,alpha_next,alpha=alpha1,alpha_new ; double xcache=0,gcache=0 ; xy cache_f(0,y),cache_df(0,gp) ; x3 a,b,alpha_prev ; int i ; a = x3(0,y,gp) ; b = x3(alpha,0,0) ; alpha_prev = x3(0,y,gp) ; /* Begin bracketing */ for(flag[0]=i=0;i<100;i++) { falpha = wrap_f(alpha,xalpha,x,p,n,&cache_f,&xcache,f) ; /* Fletcher's rho test */ if( falpha>y+alpha*rho*gp || falpha>=alpha_prev.fa ) { a = alpha_prev ; b = x3(alpha,falpha,settable()) ; break ; } cache_df = wrap_df(alpha,xalpha,x,p,g,n,cache_df,&xcache,&gcache,df) ; fpalpha = cache_df.y ; /* Fletcher's sigma test */ if(fabs(fpalpha)<=-sigma*gp) return cache_f.y ; if(fpalpha>=0) { a = x3(alpha,falpha,fpalpha) ; b = alpha_prev ; break ; } delta = alpha - alpha_prev.a ; alpha_next = interpolate(alpha_prev,x3(alpha,falpha,fpalpha), alpha+delta,alpha+tau1*delta) ; alpha_prev = x3(alpha,falpha,fpalpha) ; alpha = alpha_next ; } /* Sectioning of bracket [a,b] */ for(i=0;i<100;i++) { delta = b.a - a.a ; alpha = interpolate(a,b,a.a+tau2*delta,b.a-tau3*delta) ; falpha = wrap_f(alpha,xalpha,x,p,n,&cache_f,&xcache,f) ; if((a.a-alpha)*a.fpa<EPS) { flag[0] = 1 ; return cache_f.y ; }/* roundoff */ if( falpha>y+rho*alpha*gp || falpha>=a.fa ) b = x3(alpha,falpha) ; else { cache_df = wrap_df(alpha,xalpha,x,p,g,n,cache_df,&xcache,&gcache,df) ; fpalpha = cache_df.y ; if(fabs(fpalpha)<=-sigma*gp) return cache_f.y ; if( ((b.a-a.a)>=0&&fpalpha>=0) || ((b.a-a.a)<=0&&fpalpha<=0) ) b = a ; a = x3(alpha,falpha,fpalpha) ; } } return cache_f.y ; // can't happen } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ double bfgsmin (double *x,int n,double step,double (*f)(double *),double (*fdf)(double*,double*), int (*criterion)(int,double,double,double,double*,double*,double*,int)) { int i,iter,flag ; double q,u,alpha,gg,pp,dy,pg,oldy,dxg,dgg,dxdg,dgnorm,A,B,y,origy ; double *oldx=vector(n),*oldgrad=vector(n),*p=vector(n),*dx=vector(n) ; double *dg=vector(n),*grad=vector(n) ; origy = y = (*fdf)(x,grad) ; /* Use the gradient as the initial direction */ for(pp=i=0;i<n;i++) { p[i] = grad[i] ; pp += p[i]*p[i] ; } for(q=-1/sqrt(pp),i=0;i<n;i++) p[i] *= q ; for(iter=0;;iter++) // loop over iterations { for(pp=gg=i=0;i<n;i++) { gg += grad[i] * grad[i] ; pp += p[i] * p[i] ; oldgrad[i] = grad[i] ; oldx[i] = x[i] ; } if(iter==0) pg = -sqrt(gg) ; if(pp==0||pg==0) break ; if(iter>0&&dy<0) alpha = min(1,-2*max(-dy,10*EPS*fabs(y))/pg) ; else alpha = step?fabs(step):1 ; /* line minimisation with cubic interpolation; new point comes back as x */ oldy = y ; y = minimise(alpha,oldx,p,x,grad,n,y,pg,&flag,f,fdf) ; dy = y - oldy ; if(flag||criterion(iter+1,origy,oldy,y,oldx,x,grad,n)) break ; /* Choose a new direction for the next step */ /* This is the BFGS update: */ /* p' = g1 - A dx - B dg */ /* A = - (1+ dg.dg/dx.dg) B + dg.g/dx.dg */ /* B = dx.g/dx.dg */ for(i=0;i<n;i++) { dx[i] = x[i] - oldx[i] ; dg[i] = grad[i] - oldgrad[i] ; } for(dgnorm=dxdg=dgg=dxg=i=0;i<n;i++) { dxg += dx[i] * grad[i] ; dgg += dg[i] * grad[i] ; dxdg += dg[i] * dx[i] ; dgnorm += dg[i] * dg[i] ; } if(dxdg==0) A = B = 0 ; else { B = dxg/dxdg ; A = dgg/dxdg - (1+dgnorm/dxdg)*B ; } for(i=0;i<n;i++) p[i] = grad[i] - ( A*dx[i] + B*dg[i] ) ; for(pg=pp=i=0;i<n;i++) { pp += p[i] * p[i] ; pg += p[i] * grad[i] ; } pp = 1 / sqrt(pp) ; /* update direction and gp */ if(pg>=0) pp = -pp ; for(i=0;i<n;i++) p[i] *= pp ; pg *= pp ; } free(oldx,oldgrad,p,dx,dg,grad) ; return y ; } /* -------------------------------------------------------------------------- */ static int dim ; static int (*decision)(int,double,double,double,double*,double*,double*,int) ; static double (*g)(double *),(*gdg)(double*,double*) ; static double minusf(double *x) { return -g(x) ; } static double minusfdf(double *x,double *dx) { double q = gdg(x,dx) ; for(int i=0;i<dim;i++) dx[i] = -dx[i] ; return -q ; } static int minuscriterion(int niter,double origy,double oldy,double y, double *oldx,double *x,double *grad,int n) { int i,choice ; for(i=0;i<dim;i++) grad[i] = -grad[i] ; choice = decision(niter,-origy,-oldy,-y,oldx,x,grad,n) ; for(i=0;i<dim;i++) grad[i] = -grad[i] ; return choice ; } double bfgsmax (double *x,int n,double step,double (*f)(double *),double (*fdf)(double*,double*), int (*criterion)(int,double,double,double,double*,double*,double*,int)) { g = f ; gdg = fdf ; decision = criterion ; return bfgsmin(x,n,step,minusf,minusfdf,minuscriterion) ; }
bfgsmin.c utils: memory.h
doc : algorithm : references : source : quadprogtest : quadprogtest2 : download
quadprog maximises a quadratic function of
int quadprog(double **Q,double *A,double *x,int n0,int n1,
double **le,double *lerhs,int nle,
double **ge,double *gerhs,int nge,
double **eq,double *eqrhs,int neq) ;
or
int quadprog(double **Q,double *A,double *x,int n0) ;
where x0... xn0-1 are required to be non-negative and xn0... xn-1 may be positive or negative, and in which the function to maximise is
xTQx + 2A.x
and where
The return value is
The short calling sequence maximises the quadratic function subject only to x being non-negative. In the long calling sequence the equation and inequality arrays won’t be looked at if their counts are 0, so you can supply null pointers.
To be precise, the constraints are
∑le[i][j]×x[j] ≤ lerhs[i] for i from 0 to nle–1,
∑ge[i][j]×x[j] ≥ gerhs[i] for i from 0 to nge–1, and
∑eq[i][j]×x[j] = eqrhs[i] for i from 0 to neq–1.
It is not required that the RHSs of the constraints should be positive. If they are negative, the inequalities will be converted between ‘≤’ and ‘≥’ and the equations will have their signs flipped.
Negative return codes seem to arise from cycling or ill conditioning. Often the right solution has been found but hasn’t been recognised as correct. So you should probably report the error but proceed on the assumption that the value returned in x is correct.
void quadprogprefs(int debug,int rescale,int perturb) ; void quadprogeps(double eps) ; double quadprogresid() ;
quadprogprefs may be used to modify the behaviour of quadprog.
If debug is 0 no debug prints will be made; if it is 1 then tableaux will be printed at each stage without perturbation terms; if it is 2 then tableaux will be printed with perturbation terms. The default is 0.
If rescale is non-zero then a rough and ready rescaling of the tableau will be performed; if it is 0 then your input will be taken as given. The default is non-zero.
If perturb is non-zero then the perturbation method will be used to protect against cycling on degenerate problems; if it is 0 then perturbations are not used. The default is non-zero. [I disabled this option. I may reinstate it on maturer consideration or else provide some alternative remedy for cycling. I wonder whether brute force would be more effective; i.e. whenever the best pivot leaves the objective function unchanged, choose the first pivot with this property in a cyclic lexicographic order, and when I get back to the starting point admit failure. This is costly whenever the algorithm stalls but adds no work at all the rest of the time.]
quadprogeps may be used to adjust the floating-point tolerance; e.g. quadprogeps(1e-9) tells quadprog to treat two numbers as equal to within floating-point error if their difference is no greater than 10-9. A negative argument restores the default (currently 10-10).
quadprogresid returns the tableau residue from the most recent call to quadprog. This is the value of the objective function on termination. The objective function is initially negative and increases on each iteration (or perhaps stays the same) until no further increase is possible. If the equations are consistent the residue will be zero; otherwise it will be negative.
quadprog iterates for as long as pivoting operations are available, and then considers itself to have found a consistent answer (return code 0) if it is left with a residue that is 0 to within floating-point error; otherwise it considers itself to have been given an insoluble problem (return code -2). However you are given the final value of x in any case. If quadprog reports failure you may want to check the residue: if it is close to 0, the recovery may in fact be valid but the accumulated error larger than quadprog allowed for (though I don’t think this is likely).
The first thing to say is that quadratic programming is a very large subject – a web search will bring up lots of resources for devotees. But it is also frustrating that there is little aimed at readers who have a quadratic programming problem to solve on the way to some other objective and who want to dispatch it with the minimum of fuss. It seems to be a characteristic of operations research that it focuses on the infrequent solution of large resource allocation problems. In the case of linear programming this may not do any harm: I’ve never been in the position of having a linear programming problem to solve. But quadratic programming is different: as Wolfe himself remarks, it can be applied “to find the best least-squares fit to given data, where certain parameters are known a priori to satisfy inequality constraints”. This is a fundamental problem in pattern-matching, yet I’ve never seen anyone discuss it there. Moreover, as Roger Fletcher mentions, quadratic programming has a natural place in more general constrained optimisation problems.
So I wrote quadprog because I wanted a subroutine that did its job and that I could call from deep inside a C program. It will not be time-critical and I cannot claim that my code is efficient (and as for its being correct I can only cross my fingers).
I used Wolfe’s algorithm because it is simple and popular and because Wolfe’s description is clear and easy to follow. My implementation of the simplex algorithm started from Numerical Recipes.
People who are serious about operations research consider the tableau algorithm to be inherently inefficient (see Robert Fourer’s blog). I think that this is true in the sense in which Gaussian elimination is inefficient: for certain classes of problem (which for some people are the important classes) big speedups are available.
Wolfe says that his quadratic program required changes to only 11 instructions in his linear programming code. I feel that his linear program must have been ideally suited to the purpose: my own quadprog differs quite a lot from what I would write for linear programming.
Anyone working from Wolfe’s paper should be warned of a small error: a minus sign is needed in front of A'31 and A'32 in the bottom 2 rows of his (34).
I look at all columns and for each legal column find the best row to pivot on. For linear programming people usually do something simpler (and cheaper): they find the column with the largest coefficient in the objective function and then choose a row. I don’t know if this is guaranteed to work in quadratic programming: the chosen column might have only one positive cell, and a complementarity condition may prohibit its selection. But a less exhaustive search would be possible (e.g. sort on column coefficients and then search in decreasing order of their magnitude, terminating as soon as a valid pivot is available).
I use the perturbation method to protect against degeneracy. This is clumsy and inefficient: as I implemented it, it appreciably adds to the cost even when no degeneracy is present. The received wisdom (except for Prof. Fletcher’s) is that degeneracy is more a theoretical than a practical risk, but really the simplex algorithm isn’t sound unless it’s dealt with and my conscience made itself heard. If you resent the cpu cycles it burns you can turn the method off with quadprogprefs.
The perturbation method is essentially Wolfe’s own recommendation. In fact he refers to the Dantzig, Orden and Wolfe paper which seems to me a masterpiece of obfuscation – a simple computational expedient dressed up as heavyweight algebra – but Robert Fourer says that the algebra is actually the point. For my part I haven’t understood it well enough to know why the constant vector is m+2 long rather than m+1. A clearer account of the perturbation method is available from Vincent Conitzer’s web tutorial (so long as you don’t interpret the description of the εi as ‘positive abstract symbols’ to imply that their values must be non-negative).
Bland’s method is a great deal slicker but doesn’t carry over from linear to quadratic programming. In fact my web searches for ways of handling degeneracy in quadratic programming didn’t bring much up except for Fletcher’s papers, which sound from their abstracts to be the sort of thing I was looking for but whose contents are hidden from me. That he pays close attention to rounding error is another point in their favour.
Debug prints do not show perturbations unless you specifically ask for them, but it’s worth seeing the results of a pivoting operation at a point at which the Simplex algorithm stalls on Wolfe’s example:
====|==============================|========|========================
| x0 x1 v1 v2 U0 | | perturbations
----|------------------------------|--------|------------------------
x2 |-1.00 1.00 0.00 0.00 0.00 | 1.00 | 1.00 0.00 0.00 0.00
v0 | 1.00 -0.00 -0.00 -0.00 1.00 | 1.00 | 0.00 1.00 0.00 0.00
z1 | 0.00 -1.00 1.00 0.00 1.00 | 0.00 | 0.00 0.00 1.00 0.00
z2 | 1.00 -1.00 0.00 1.00 -1.00 | 1.00 |-1.00 0.00 0.00 1.00
----|------------------------------|--------|------------------------
obj |-1.00 2.00 -1.00 -1.00 0.00 | -1.00 | 3.00 2.00 1.00 1.00
====|========================|========|========================
| x0 v1 v2 U0 | | perturbations
----|------------------------|--------|------------------------
x2 |-1.00 1.00 0.00 1.00 | 1.00 | 1.00 0.00 1.00 0.00
v0 | 1.00 -0.00 -0.00 1.00 | 1.00 | 0.00 1.00 0.00 0.00
x1 | 0.00 1.00 0.00 1.00 | 0.00 | 0.00 0.00 1.00 0.00
z2 | 1.00 -1.00 1.00 -2.00 | 1.00 |-1.00 0.00 -1.00 1.00
----|------------------------|--------|------------------------
obj |-1.00 1.00 -1.00 2.00 | -1.00 | 3.00 2.00 3.00 1.00
The pink cell is the chosen pivot. No improvement is made to the objective function nor to the first two perturbation terms, and it is the third perturbation which shows the pivot to be advantageous. The artificial variable z1 is moved to the right and eliminated, and in this case the pivot cell would have been accepted without the perturbation rule (or even if considered by it to be disadvantageous) because elimination of cells is irreversible and therefore brings no risk of cycling.
I consider two values to be equal up to floating-point error if they differ by no more than 10-10. For this criterion to be reasonable I rescale the problem to have terms of order around 1, but the rescaling is very rough and ready. I test for equality ‘up to rounding error’ much more often than Numerical Recipes, which often tests for equality tout court – I hope what I’ve done is right.
• quadprogprefs • quadprogeps • var • varvector • swap • release • delrow • delcol • quadprint • tableau • pivot • findpivot • simplexalg • quadprog • quadprogresid
/* Copyright (c) 2016, by Colin Champion <colin.champion@masterlyinactivity.com> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include "memory.h" #include <math.h> #define EPS 1e-10 static int debug=0,rescale=1,perturb=1,debugper=0 ; static double qresid=0,eps=EPS ; void quadprogprefs(int d,int r,int p) { rescale = r ; perturb = p ; debug = (d>0) ; debugper = p&&(d>1) ; } void quadprogeps(double e) { if(e>=0) eps = e ; else eps = EPS ; } struct var { char type ; int ind ; var() { type = 0 ; ind = -1 ; } var(char c,int i) { type = c ; ind = i ; } } ; inline bool operator==(var v1,var v2) { return v1.type==v2.type && v1.ind==v2.ind ; } static var *varvector(int n) { return (var *) cjcalloc(n,sizeof(var)) ; } static void swap(var &a,var &b) { var c=b ; b = a ; a = c ; } struct tableau { double **a,*scale,rheps ; int nrow,ncol,nle,nge,neq,n,n0,n1 ; var *lh,*rh ; void release() { free(lh,rh,scale) ; freematrix(a) ; } void delrow(int row) { int j ; if(row<nrow-1) { for(j=0;j<=ncol;j++) a[row][j] = a[nrow-1][j] ; lh[row] = lh[nrow-1] ; } for(j=0;j<=ncol;j++) a[nrow-1][j] = a[nrow][j] ; nrow -= 1 ; } void delcol(int col) { int i ; if(col<ncol-1) { for(i=0;i<=nrow;i++) a[i][col] = a[i][ncol-1] ; rh[col] = rh[ncol-1] ; } for(i=0;i<=nrow;i++) a[i][ncol-1] = a[i][ncol] ; ncol -= 1 ; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // tableau format: // the total number of columns in the tableau is ncol+1 (the extra '1' // is for the constant terms). The columns // include artificial variables and KKT multipliers. // Variables are coded as follows: // x[j] <-> true variables // y[i] <-> slack variables for inequalities, J&B's v // w[j] <-> artificial variables for constraints, NR's z // v[j] <-> KKT multipliers, complementary with x, J&B's y // u[i] <-> " " , " " y, J&B's mu // z[j] <-> artificial variables for stationarity condition // upper case 'X' are unrestricted; upper case 'U' come from equations rather // than inequalities and are therefore also unrestricted; upper case 'W' come // from equations rather than inequalities and can be eliminated as soon as // they get to the RHS. X and U never move from the LHS. void quadprint(char *s,int phase) { int i,j ; char c,ltype ; printf("\n%s",s) ; for(i=0;i<=nrow;i++) { if(i==0) { printf("====|") ; for(j=0;j<ncol;j++) printf("======") ; printf("|=======") ; printf("\n |") ; for(j=0;j<ncol;j++) printf(" %c%-3d",rh[j].type,rh[j].ind) ; printf("|") ; if(debugper) printf(" | perturbations") ; printf("\n----|") ; for(j=0;j<ncol;j++) printf("------") ; printf("|-------") ; printf("\n") ; } ltype = lh[i].type ; c = ' ' ; if(i<nrow) { if((ltype=='z'||ltype=='w')&&a[i][ncol]>rheps) c = '>' ; else if(ltype=='W'&&fabs(a[i][ncol])>rheps) c = '>' ; if(phase==0&<ype=='z') c = ' ' ; printf("%c%c%-2d|",c,ltype,lh[i].ind) ; } else printf("obj |") ; for(j=0;j<ncol;j++) { if(fabs(a[i][j])<0.01) printf(" - ") ; else printf("%5.2f ",a[i][j]) ; } if(fabs(a[i][ncol])<0.01) printf("| -") ; else printf("|%7.2f",a[i][ncol]) ; if(debugper) { printf(" |") ; for(j=ncol+1;j<=ncol;j++) printf("%5.2f ",a[i][j]) ; } else if(i==nrow&&fabs(a[i][ncol])<0.02) printf(" (%.1e)",a[i][ncol]) ; else if(i<nrow&&(ltype!='U'&<ype!='X')&&a[i][ncol]<-rheps) printf(" -> 0 ***") ; else if(c=='>'&&fabs(a[i][ncol])<0.02) printf(" (%.1e)",a[i][ncol]) ; printf("\n") ; if(i==nrow-1) { printf("----|") ; for(j=0;j<ncol;j++) printf("------") ; printf("|-------") ; printf("\n") ; } } } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ tableau(double **Q,double *A,int na,int nb,double **le,double *lerhs,int m1, double **ge,double *gerhs,int m2,double **eq,double *eqrhs,int m3) { int i,j,ind,m=m1+m2+m3 ; double q,sqrtn ; n0 = na ; n1 = nb ; n = n0 + n1 ; nrow = m + n ; ncol = m + 2*n + n0 ; nle = nge = neq = 0 ; a = matrix(nrow+1,ncol+1) ; scale = vector(n) ; // variable ids lh = varvector(nrow) ; rh = varvector(ncol) ; // first m rows of the tableau state primal feasibility // first true >= for(ind=i=0;i<m1;i++) if(lerhs[i]>=0) { for(j=0;j<n;j++) a[ind][j] = -le[i][j] ; a[ind++][ncol] = lerhs[i] ; } for(i=0;i<m2;i++) if(gerhs[i]<0) { for(j=0;j<n;j++) a[ind][j] = ge[i][j] ; a[ind++][ncol] = -gerhs[i] ; } nle = ind ; nge = m1 + m2 - nle ; neq = m3 ; // now true >= for(i=0;i<m1;i++) if(lerhs[i]<0) { for(j=0;j<n;j++) a[ind][j] = le[i][j] ; a[ind++][ncol] = -lerhs[i] ; } for(i=0;i<m2;i++) if(gerhs[i]>=0) { for(j=0;j<n;j++) a[ind][j] = -ge[i][j] ; a[ind++][ncol] = gerhs[i] ; } // and finally equalities for(i=0;i<m3;ind++,i++) if(eqrhs[i]>=0) { for(j=0;j<n;j++) a[ind][j] = -eq[i][j] ; a[ind][ncol] = eqrhs[i] ; } else { for(j=0;j<n;j++) a[ind][j] = eq[i][j] ; a[ind][ncol] = -eqrhs[i] ; } for(i=0;i<n;i++) scale[i] = 1 ; if(rescale) { sqrtn = sqrt(n) ; for(j=0;j<n;j++) if(Q[j][j]<0) scale[j] = 1 / sqrt(-Q[j][j]) ; for(i=0;i<m;i++) for(j=0;j<n;j++) a[i][j] *= scale[j] ; for(i=0;i<m;i++) { for(q=j=0;j<n;j++) q += a[i][j] * a[i][j] ; a[i][ncol] *= ( q = 1 / sqrt(q/sqrtn) ) ; for(j=0;j<n;j++) a[i][j] *= q ; } } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ // remaining n rows state the KKT stationarity condition // note that we don't really need the artificial z[i] when A[i] is // negative - we could put the v[i] into the basis directly for(i=0;i<n;i++) { for(j=0;j<n;j++) a[m+i][j] = Q[i][j] * scale[i] * scale[j] ; if(i<n0) a[m+i][n+i] = 1 ; for(j=0;j<nle;j++) a[m+i][n+n0+j] = a[j][i] ; for(;j<nle+nge;j++) a[m+i][n+n0+j] = -a[j][i] ; // correcting Wolfe for(;j<m;j++) a[m+i][n+n0+j] = a[j][i] ; a[m+i][ncol] = A[i] * scale[i] ; lh[m+i] = var('z',i) ; rh[m+n+n0+i] = var('Z',i) ; if(A[i]<0) { for(j=0;j<=ncol;j++) a[m+i][j] = -a[m+i][j] ; swap(lh[m+i],rh[m+n+n0+i]) ; } a[m+i][m+n+n0+i] = 1 ; } for(i=0;i<n0;i++) rh[i] = var('x',i) ; for(i=n0;i<n;i++) rh[i] = var('X',i) ; for(i=0;i<n0;i++) rh[n+i] = var('v',i) ; for(i=0;i<nle+nge;i++) rh[n+n0+i] = var('u',i) ; for(i=nle+nge;i<m;i++) rh[n+n0+i] = var('U',i) ; for(i=0;i<nle;i++) lh[i] = var('y',i) ; // y for <= for(i=nle;i<nle+nge;i++) lh[i] = var('w',i) ; // w for >= for(i=nle+nge;i<m;i++) lh[i] = var('W',i) ; // w for == qresid = 0 ; for(q=eps,i=0;i<nrow;i++) q += a[i][ncol] * a[i][ncol] ; rheps = eps * sqrt(q/nrow) ; } /* ------------------------------------------------------------------------ */ void pivot(int pivrow,int pivcol) { int i,j ; double piv=1/a[pivrow][pivcol],q ; for(i=0;i<=nrow;i++) if(i!=pivrow) { for(q=a[i][pivcol]*piv,j=0;j<=ncol;j++) a[i][j] -= a[pivrow][j] * q ; a[i][pivcol] = q ; } for(j=0;j<=ncol;j++) a[pivrow][j] *= -piv ; a[pivrow][pivcol] = piv ; swap(rh[pivcol],lh[pivrow]) ; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ // return -2 if no pivot found, -1 if obj fn can be increased without limit. // Note that we mustn't multiply qobj into r until the end of the loop because // if qobj was 0 we wouldn't find the row variable which first became -ve. int findpivot(int phase) { int i,pivrow,pivcol,ret,rtype,ltype,lforce ; double rmax,r,rbest,qobj,qcell,rforce ; for(ret=-2,pivcol=0;pivcol<ncol;pivcol++) { rtype = rh[pivcol].type ; if(phase==0&&(rtype=='u'||rtype=='U'||rtype=='v')) continue ; qobj = a[nrow][pivcol] ; pivrow = -1 ; if(phase) // see if complementarity forces a pivot for this column { if(rtype=='u') ltype = 'y' ; else if(rtype=='y') ltype = 'u' ; else if(rtype=='v') ltype = 'x' ; else if(rtype=='x') ltype = 'v' ; else ltype = '$' ; // dummy if(ltype!='$') for(i=0;pivrow<0&&i<nrow;i++) if(lh[i].type==ltype&&lh[i].ind==rh[pivcol].ind) pivrow = i ; } if(pivrow>=0) // a variable whose complement is in the basis { if(qobj>eps||(a[nrow][ncol]<-rheps&&qobj>=0)) { if((qcell=a[pivrow][pivcol])>=0) continue ; } else continue ; rbest = -a[pivrow][ncol] / qcell ; for(i=0;pivrow>=0&&i<nrow;i++) if(lh[i].type!='U'&&lh[i].type!='X') if((qcell=a[i][pivcol])<0&&-a[i][ncol]/qcell<rbest-eps) pivrow = -1 ; if(pivrow<0) continue ; } else if(qobj>eps||(a[nrow][ncol]<-rheps&&qobj>=0)) // restricted variable { for(i=0;i<nrow;i++) if(lh[i].type!='U'&&lh[i].type!='X') { if((qcell=a[i][pivcol])<0) { r = -a[i][ncol] / qcell ; // r is positive (up to fp error) if(pivrow<0||r<rbest) { pivrow = i ; rbest = r ; } } else if(qcell<eps&&pivrow<0) pivrow = -2 ; // cell 0 up to fp error } if(qobj>eps&&pivrow==-1) return -1 ; } else if(qobj<-eps&&(rtype=='U'||rtype=='X')) // unrestricted variable { for(i=0;i<nrow;i++) if((qcell=a[i][pivcol])>0&&lh[i].type!='U'&&lh[i].type!='X') { r = - a[i][ncol] / qcell ; if(pivrow<0||r>rbest) { pivrow = i ; rbest = r ; } } if(qobj<-eps&&pivrow==-1) return -1 ; } else continue ; if(qobj<=eps) rbest = -1 ; // only take a column with obj 0 as last resort if(pivrow>=0&&(ret<0||rbest*qobj>rmax)) { rmax = rbest*qobj ; ret = pivrow + nrow*pivcol ; } } // end loop over pivcol return ret ; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ int simplexalg(double *x) { int i,j,row,col,ret,iter,maxiter,phase ; double qmax,q ; char empty[100] ; for(i=0;i<n;i++) x[i] = 0 ; // objective function for(ret=phase=0;ret>=0&&phase<2;phase++) { if(phase) for(i=0;i<ncol;i++) // discard ... the unused cpts of z1, z2 if(rh[i].type=='z'||rh[i].type=='Z') { delcol(i) ; i -= 1 ; } for(j=0;j<=ncol;j++) a[nrow][j] = 0 ; if(phase==0) { for(i=0;i<nrow;i++) if(lh[i].type=='w'||lh[i].type=='W') for(j=0;j<=ncol;j++) a[nrow][j] -= a[i][j] ; } else for(i=0;i<nrow;i++) { // some w rows may still be present if their constant terms are 0 if(lh[i].type=='w'||lh[i].type=='W') { delrow(i) ; i -= 1 ; } else if(lh[i].type=='z'||lh[i].type=='Z') for(j=0;j<=ncol;j++) a[nrow][j] -= a[i][j] ; } for(q=1,i=0;i<nrow;i++) q *= (nrow+ncol-i) / (1.0+i) ; maxiter = min(q,20000) ; if(debug) { printf("maxiter= %d; rheps=%.1e\n",maxiter,rheps) ; quadprint("",phase) ; } for(iter=0;;iter++) { if(phase==0) { for(ret=i=0;ret==0&&i<nrow;i++) if( (lh[i].type=='w'&&a[i][ncol]>rheps) || (lh[i].type=='W'&&fabs(a[i][ncol])>rheps) ) ret = 1 ; } else for(ret=i=0;ret==0&&i<nrow;i++) if((lh[i].type=='z'||lh[i].type=='Z')&&a[i][ncol]>rheps) ret = 1 ; if(ret!=0&&a[nrow][ncol]>=0) { ret = -4 ; break ; } // don't go haywire if(ret!=0&&iter==maxiter) { ret = -3 ; break ; } // too many iters if(ret==0||(ret=findpivot(phase))<0) break ; row = ret % nrow ; col = ret / nrow ; if(debug==0) empty[0] = 0 ; else sprintf(empty,"iter %d: %c%d <-> %c%d (%.1e <-> %.1e)\n", iter,lh[row].type,lh[row].ind,rh[col].type,rh[col].ind, a[row][ncol],a[nrow][col]) ; pivot(row,col) ; // unrestricted LH Lagrangians U and artificial RH vbles can be zapped if(debug==0&&lh[row].type=='U') delrow(row) ; if( rh[col].type=='W' || (phase==1&&(rh[col].type=='z'||rh[col].type=='Z')) ) delcol(col) ; // the artificial variable for a '>=' needs to become a slack variable else if(rh[col].type=='w') { rh[col].type = 'y' ; a[nrow][col] = -(1+a[nrow][col]) ; for(i=0;i<nrow;i++) a[i][col] = -a[i][col] ; } if(debug) quadprint(empty,phase) ; for(row=0;row<nrow;row++) if(lh[row].type!= 'X'&&lh[row].type!= 'U') { if(a[row][ncol]<-10*rheps) return -4 ; else if(a[row][ncol]<0) a[row][ncol] = 0 ; } } // end loop over iter qresid = a[nrow][ncol] ; // there's an auxiliary variable stuck on the LHS and not equal to 0 // but the residue is small anyway, perhaps we're alright if(ret==-2&&qresid>=-rheps) ret = 0 ; } // end loop over phase for(i=0;i<nrow;i++) if(lh[i].type=='x'||lh[i].type=='X') x[lh[i].ind] = a[i][ncol] * scale[lh[i].ind] ; return ret ; } } ; /* -------------------------------------------------------------------------- */ int quadprog(double **Q,double *A,double *x,int n0,int n1, double **le,double *lerhs,int nle, double **ge,double *gerhs,int nge, double **eq,double *eqrhs,int neq) { tableau t(Q,A,n0,n1,le,lerhs,nle,ge,gerhs,nge,eq,eqrhs,neq) ; int ret = t.simplexalg(x) ; t.release() ; return ret ; } int quadprog(double **Q,double *A,double *x,int n) { return quadprog(Q,A,x,n,0, 0,0,0, 0,0,0, 0,0,0) ; } double quadprogresid() { return qresid ; }
The test program quadprogtest solves 6 problems. The first four find the closest point to (-1,-1) subject to various constraints (partly chosen to verify the handling of unconstrained variables); the last 2 are the examples given by Wolfe and by J&B.
The correct output comes after the source code.
#include "memory.h"
int quadprog(double **qmat,double *qadd,double *x,int n0,int n1,
double **le,double *lerhs,int nle,
double **ge,double *gerhs,int nge,
double **eq,double *eqrhs,int neq) ;
void quadprogprefs(int d,int r,int p) ;
static int ngood=0,ndone=0 ;
void printsoln(int ret,double *x,int n,double *ans)
{ int t ;
double q ;
ndone += 1 ;
if(ret==0)
{ printf("---- solution: ") ;
for(t=0;t<n;t++) printf("%.2f ",x[t]) ; printf("\n") ;
for(q=t=0;t<n;t++) q += (x[t]-ans[t]) * (x[t]-ans[t]) ;
if(q<1e-12) ngood += 1 ;
}
else printf("---- no solution: ret=%d\n",ret) ;
}
int main(int argc,char **argv)
{ int i,j,ret,probno=0 ;
double x[3],ans[3],**Q=matrix(3,3),A[3],**eq=matrix(3,3),eqrhs[3] ;
if(argc>=2) probno = atoi(argv[1]) ;
quadprogprefs(1,0,0) ;
// the first four problems find the closest point to (-1,-1)
// subject to various linear constraints
for(i=0;i<3;i++) { Q[i][i] = -1 ; A[i] = -1 ; }
eq[0][0] = eq[0][1] = 1 ;
eqrhs[0] = 2 ;
// x + y = 2
if(probno==1||probno<=0)
{ printf("---- problem 1 [right ans = (1.00,1.00)]\n") ;
ans[1] = ans[0] = 1 ;
ret = quadprog(Q,A,x,2,0, 0,0,0, 0,0,0, eq,eqrhs,1 ) ;
printsoln(ret,x,2,ans) ;
}
// x + y >= 2
if(probno==2||probno<=0)
{ printf("\n---- problem 2 [right ans = (1.00,1.00)]\n") ;
ans[1] = ans[0] = 1 ;
ret = quadprog(Q,A,x,2,0, 0,0,0, eq,eqrhs,1, 0,0,0 ) ;
printsoln(ret,x,2,ans) ;
}
// x + 2y >= -1
if(probno==3||probno<=0)
{ printf("\n---- problem 3 [right ans = (0.00,-0.50)]\n") ;
ans[0] = 0 ;
ans[1] = -0.5 ;
eq[0][0] = 1 ;
eq[0][1] = 2 ;
eqrhs[0] = -1 ;
ret = quadprog(Q,A,x,1,1, 0,0,0, eq,eqrhs,1, 0,0,0 ) ;
printsoln(ret,x,2,ans) ;
}
// 3x + 4y <= -32
if(probno==4||probno<=0)
{ printf("\n---- problem 4 [right ans = (-4.00,-5.00)]\n") ;
ans[0] = -4 ;
ans[1] = -5 ;
eq[0][0] = 3 ;
eq[0][1] = 4 ;
eqrhs[0] = -32 ;
ret = quadprog(Q,A,x,0,2, eq,eqrhs,1, 0,0,0, 0,0,0 ) ;
printsoln(ret,x,2,ans) ;
}
if(probno==5||probno<=0)
{ printf("\n---- problem 5 (J&B) [right ans = (3.00,2.00)]\n") ;
ans[0] = 3 ;
ans[1] = 2 ;
Q[0][0] = -1 ;
Q[1][1] = -4 ;
A[0] = 4 ;
A[1] = 8 ;
eq[0][0] = eq[0][1] = 1 ;
eqrhs[0] = 5 ;
eq[1][0] = 1 ;
eq[1][1] = 0 ;
eqrhs[1] = 3 ;
ret = quadprog(Q,A,x,2,0, eq,eqrhs,2, 0,0,0, 0,0,0 ) ;
printsoln(ret,x,2,ans) ;
}
if(probno==6||probno<=0)
{ printf("\n---- problem 6 (Wolfe) [right ans = (0.00,0.50,1.50)]\n") ;
ans[0] = 0 ;
ans[1] = 0.5 ;
ans[2] = 1.5 ;
Q[0][0] = Q[1][1] = Q[2][2] = -1 ;
A[0] = -1 ;
A[1] = 0 ;
A[2] = 2 ;
eq[0][0] = eq[0][2] = 1 ;
eq[0][1] = -1 ;
eqrhs[0] = 1 ;
ret = quadprog(Q,A,x,3,0, 0,0,0, 0,0,0, eq,eqrhs,1 ) ;
printsoln(ret,x,3,ans) ;
}
// closest point to (1,2) s.t. x >= 1, y >= 2, y-x >= 2
if(probno==7||probno<=0)
{ printf("\n---- problem 7 [right ans = (1.00,3.00)]\n") ;
ans[0] = 1 ;
ans[1] = 3 ;
Q[0][0] = Q[1][1] = -1 ;
A[0] = 1 ;
A[1] = 2 ;
eq[0][0] = 1 ;
eq[0][1] = eq[1][0] = 0 ;
eq[2][1] = eq[1][1] = 1 ;
eq[2][0] = -1 ;
eqrhs[0] = 1 ;
eqrhs[2] = eqrhs[1] = 2 ;
ret = quadprog(Q,A,x,2,0, 0,0,0, eq,eqrhs,3, 0,0,0 ) ;
printsoln(ret,x,2,ans) ;
}
if(ngood<ndone) printf("*** %d failures\n",ndone-ngood) ;
else printf("congratulations - all problems solved!\n") ;
}
---- problem 1 [right ans = (1.00,1.00)]
maxiter= 120; rheps=1.4e-10
====|==========================================|=======
| x0 x1 v0 v1 U0 z0 z1 |
----|------------------------------------------|-------
>W0 |-1.00 -1.00 - - - - - | 2.00
Z0 | 1.00 - -1.00 - 1.00 1.00 - | 1.00
Z1 | - 1.00 - -1.00 1.00 - 1.00 | 1.00
----|------------------------------------------|-------
obj | 1.00 1.00 - - - - - | -2.00
iter 0: W0 <-> x0 (2.0e+00 <-> 1.0e+00)
====|====================================|=======
| x1 v0 v1 U0 z0 z1 |
----|------------------------------------|-------
x0 |-1.00 - - - - - | 2.00
Z0 |-1.00 -1.00 - 1.00 1.00 - | 3.00
Z1 | 1.00 - -1.00 1.00 - 1.00 | 1.00
----|------------------------------------|-------
obj | - - - - - - | - (0.0e+00)
maxiter= 35; rheps=1.4e-10
====|========================|=======
| x1 v0 v1 U0 |
----|------------------------|-------
x0 |-1.00 - - - | 2.00
>Z0 |-1.00 -1.00 - 1.00 | 3.00
>Z1 | 1.00 - -1.00 1.00 | 1.00
----|------------------------|-------
obj | - 1.00 1.00 -2.00 | -4.00
iter 0: Z1 <-> U0 (1.0e+00 <-> -2.0e+00)
====|==================|=======
| x1 v0 v1 |
----|------------------|-------
x0 |-1.00 - - | 2.00
>Z0 |-2.00 -1.00 1.00 | 2.00
U0 |-1.00 - 1.00 | -1.00
----|------------------|-------
obj | 2.00 1.00 -1.00 | -2.00
iter 1: Z0 <-> x1 (2.0e+00 <-> 2.0e+00)
====|============|=======
| v0 v1 |
----|------------|-------
x0 | 0.50 -0.50 | 1.00
x1 |-0.50 0.50 | 1.00
U0 | 0.50 0.50 | -2.00
----|------------|-------
obj | - - | - (0.0e+00)
---- solution: 1.00 1.00
---- problem 2 [right ans = (1.00,1.00)]
maxiter= 120; rheps=1.4e-10
====|==========================================|=======
| x0 x1 v0 v1 u0 z0 z1 |
----|------------------------------------------|-------
>w0 |-1.00 -1.00 - - - - - | 2.00
Z0 | 1.00 - -1.00 - -1.00 1.00 - | 1.00
Z1 | - 1.00 - -1.00 -1.00 - 1.00 | 1.00
----|------------------------------------------|-------
obj | 1.00 1.00 - - - - - | -2.00
iter 0: w0 <-> x0 (2.0e+00 <-> 1.0e+00)
====|==========================================|=======
| y0 x1 v0 v1 u0 z0 z1 |
----|------------------------------------------|-------
x0 | 1.00 -1.00 - - - - - | 2.00
Z0 | 1.00 -1.00 -1.00 - -1.00 1.00 - | 3.00
Z1 | - 1.00 - -1.00 -1.00 - 1.00 | 1.00
----|------------------------------------------|-------
obj | - - - - - - - | - (0.0e+00)
maxiter= 56; rheps=1.4e-10
====|==============================|=======
| y0 x1 v0 v1 u0 |
----|------------------------------|-------
x0 | 1.00 -1.00 - - - | 2.00
>Z0 | 1.00 -1.00 -1.00 - -1.00 | 3.00
>Z1 | - 1.00 - -1.00 -1.00 | 1.00
----|------------------------------|-------
obj |-1.00 - 1.00 1.00 2.00 | -4.00
iter 0: Z1 <-> u0 (1.0e+00 <-> 2.0e+00)
====|========================|=======
| y0 x1 v0 v1 |
----|------------------------|-------
x0 | 1.00 -1.00 - - | 2.00
>Z0 | 1.00 -2.00 -1.00 1.00 | 2.00
u0 | - 1.00 - -1.00 | 1.00
----|------------------------|-------
obj |-1.00 2.00 1.00 -1.00 | -2.00
iter 1: Z0 <-> x1 (2.0e+00 <-> 2.0e+00)
====|==================|=======
| y0 v0 v1 |
----|------------------|-------
x0 | 0.50 0.50 -0.50 | 1.00
x1 | 0.50 -0.50 0.50 | 1.00
u0 | 0.50 -0.50 -0.50 | 2.00
----|------------------|-------
obj | - - - | - (0.0e+00)
---- solution: 1.00 1.00
---- problem 3 [right ans = (0.00,-0.50)]
maxiter= 84; rheps=1.0e-10
====|====================================|=======
| x0 X1 v0 u0 z0 z1 |
----|------------------------------------|-------
y0 | 1.00 2.00 - - - - | 1.00
Z0 | 1.00 - -1.00 -1.00 1.00 - | 1.00
Z1 | - 1.00 - -2.00 - 1.00 | 1.00
----|------------------------------------|-------
obj | - - - - - - | - (0.0e+00)
maxiter= 35; rheps=1.0e-10
====|========================|=======
| x0 X1 v0 u0 |
----|------------------------|-------
y0 | 1.00 2.00 - - | 1.00
>Z0 | 1.00 - -1.00 -1.00 | 1.00
>Z1 | - 1.00 - -2.00 | 1.00
----|------------------------|-------
obj |-1.00 -1.00 1.00 3.00 | -2.00
iter 0: y0 <-> X1 (1.0e+00 <-> -1.0e+00)
====|========================|=======
| x0 y0 v0 u0 |
----|------------------------|-------
X1 |-0.50 0.50 - - | -0.50
>Z0 | 1.00 - -1.00 -1.00 | 1.00
>Z1 |-0.50 0.50 - -2.00 | 0.50
----|------------------------|-------
obj |-0.50 -0.50 1.00 3.00 | -1.50
iter 1: Z0 <-> v0 (1.0e+00 <-> 1.0e+00)
====|==================|=======
| x0 y0 u0 |
----|------------------|-------
X1 |-0.50 0.50 - | -0.50
v0 | 1.00 - -1.00 | 1.00
>Z1 |-0.50 0.50 -2.00 | 0.50
----|------------------|-------
obj | 0.50 -0.50 2.00 | -0.50
iter 2: Z1 <-> u0 (5.0e-01 <-> 2.0e+00)
====|============|=======
| x0 y0 |
----|------------|-------
X1 |-0.50 0.50 | -0.50
v0 | 1.25 -0.25 | 0.75
u0 |-0.25 0.25 | 0.25
----|------------|-------
obj | - - | - (0.0e+00)
---- solution: 0.00 -0.50
---- problem 4 [right ans = (-4.00,-5.00)]
maxiter= 56; rheps=1.8e-09
====|==============================|=======
| X0 X1 u0 z0 z1 |
----|------------------------------|-------
>w0 | 3.00 4.00 - - - | 32.00
Z0 | 1.00 - 3.00 1.00 - | 1.00
Z1 | - 1.00 4.00 - 1.00 | 1.00
----|------------------------------|-------
obj |-3.00 -4.00 - - - | -32.00
iter 0: Z1 <-> X1 (1.0e+00 <-> -4.0e+00)
====|==============================|=======
| X0 Z1 u0 z0 z1 |
----|------------------------------|-------
>w0 | 3.00 4.00 -16.00 - -4.00 | 28.00
Z0 | 1.00 - 3.00 1.00 - | 1.00
X1 | - 1.00 -4.00 - -1.00 | -1.00
----|------------------------------|-------
obj |-3.00 -4.00 16.00 - 4.00 | -28.00
iter 1: w0 <-> z1 (2.8e+01 <-> 4.0e+00)
====|==============================|=======
| X0 Z1 u0 z0 y0 |
----|------------------------------|-------
z1 | 0.75 1.00 -4.00 - 0.25 | 7.00
Z0 | 1.00 - 3.00 1.00 - | 1.00
X1 |-0.75 - - - -0.25 | -8.00
----|------------------------------|-------
obj | - - - - - | - (0.0e+00)
maxiter= 20; rheps=1.8e-09
====|==================|=======
| X0 u0 y0 |
----|------------------|-------
>z1 | 0.75 -4.00 0.25 | 7.00
>Z0 | 1.00 3.00 - | 1.00
X1 |-0.75 - -0.25 | -8.00
----|------------------|-------
obj |-1.75 1.00 -0.25 | -8.00
iter 0: Z0 <-> X0 (1.0e+00 <-> -1.8e+00)
====|============|=======
| u0 y0 |
----|------------|-------
>z1 |-6.25 0.25 | 6.25
X0 |-3.00 - | -1.00
X1 | 2.25 -0.25 | -7.25
----|------------|-------
obj | 6.25 -0.25 | -6.25
iter 1: z1 <-> u0 (6.2e+00 <-> 6.2e+00)
====|======|=======
| y0 |
----|------|-------
u0 | 0.04 | 1.00
X0 |-0.12 | -4.00
X1 |-0.16 | -5.00
----|------|-------
obj | - | - (0.0e+00)
---- solution: -4.00 -5.00
---- problem 5 (J&B) [right ans = (3.00,2.00)]
maxiter= 495; rheps=5.3e-10
====|================================================|=======
| x0 x1 v0 v1 u0 u1 Z0 Z1 |
----|------------------------------------------------|-------
y0 |-1.00 -1.00 - - - - - - | 5.00
y1 |-1.00 - - - - - - - | 3.00
z0 |-1.00 - 1.00 - -1.00 -1.00 1.00 - | 4.00
z1 | - -4.00 - 1.00 -1.00 - - 1.00 | 8.00
----|------------------------------------------------|-------
obj | - - - - - - - - | - (0.0e+00)
maxiter= 210; rheps=5.3e-10
====|====================================|=======
| x0 x1 v0 v1 u0 u1 |
----|------------------------------------|-------
y0 |-1.00 -1.00 - - - - | 5.00
y1 |-1.00 - - - - - | 3.00
>z0 |-1.00 - 1.00 - -1.00 -1.00 | 4.00
>z1 | - -4.00 - 1.00 -1.00 - | 8.00
----|------------------------------------|-------
obj | 1.00 4.00 -1.00 -1.00 2.00 1.00 | -12.00
iter 0: z1 <-> x1 (8.0e+00 <-> 4.0e+00)
====|==============================|=======
| x0 v0 v1 u0 u1 |
----|------------------------------|-------
y0 |-1.00 - -0.25 0.25 - | 3.00
y1 |-1.00 - - - - | 3.00
>z0 |-1.00 1.00 - -1.00 -1.00 | 4.00
x1 | - - 0.25 -0.25 - | 2.00
----|------------------------------|-------
obj | 1.00 -1.00 - 1.00 1.00 | -4.00
iter 1: y0 <-> x0 (3.0e+00 <-> 1.0e+00)
====|==============================|=======
| y0 v0 v1 u0 u1 |
----|------------------------------|-------
x0 |-1.00 - -0.25 0.25 - | 3.00
y1 | 1.00 - 0.25 -0.25 - | -
>z0 | 1.00 1.00 0.25 -1.25 -1.00 | 1.00
x1 | - - 0.25 -0.25 - | 2.00
----|------------------------------|-------
obj |-1.00 -1.00 -0.25 1.25 1.00 | -1.00
iter 2: y1 <-> u0 (0.0e+00 <-> 1.2e+00)
====|==============================|=======
| y0 v0 v1 y1 u1 |
----|------------------------------|-------
x0 | - - - -1.00 - | 3.00
u0 | 4.00 - 1.00 -4.00 - | -
>z0 |-4.00 1.00 -1.00 5.00 -1.00 | 1.00
x1 |-1.00 - - 1.00 - | 2.00
----|------------------------------|-------
obj | 4.00 -1.00 1.00 -5.00 1.00 | -1.00
iter 3: z0 <-> u1 (1.0e+00 <-> 1.0e+00)
====|========================|=======
| y0 v0 v1 y1 |
----|------------------------|-------
x0 | - - - -1.00 | 3.00
u0 | 4.00 - 1.00 -4.00 | -
u1 |-4.00 1.00 -1.00 5.00 | 1.00
x1 |-1.00 - - 1.00 | 2.00
----|------------------------|-------
obj | - - - - | - (0.0e+00)
---- solution: 3.00 2.00
---- problem 6 (Wolfe) [right ans = (0.00,0.50,1.50)]
maxiter= 1001; rheps=1.2e-10
====|============================================================|=======
| x0 x1 x2 v0 v1 v2 U0 z0 Z1 Z2 |
----|------------------------------------------------------------|-------
>W0 |-1.00 1.00 -1.00 - - - - - - - | 1.00
Z0 | 1.00 - - -1.00 - - 1.00 1.00 - - | 1.00
z1 | - -1.00 - - 1.00 - 1.00 - 1.00 - | -
z2 | - - -1.00 - - 1.00 -1.00 - - 1.00 | 2.00
----|------------------------------------------------------------|-------
obj | 1.00 -1.00 1.00 - - - - - - - | -1.00
iter 0: W0 <-> x0 (1.0e+00 <-> 1.0e+00)
====|======================================================|=======
| x1 x2 v0 v1 v2 U0 z0 Z1 Z2 |
----|------------------------------------------------------|-------
x0 | 1.00 -1.00 - - - - - - - | 1.00
Z0 | 1.00 -1.00 -1.00 - - 1.00 1.00 - - | 2.00
z1 |-1.00 - - 1.00 - 1.00 - 1.00 - | -
z2 | - -1.00 - - 1.00 -1.00 - - 1.00 | 2.00
----|------------------------------------------------------|-------
obj | - - - - - - - - - | - (0.0e+00)
maxiter= 210; rheps=1.2e-10
====|====================================|=======
| x1 x2 v0 v1 v2 U0 |
----|------------------------------------|-------
x0 | 1.00 -1.00 - - - - | 1.00
>Z0 | 1.00 -1.00 -1.00 - - 1.00 | 2.00
z1 |-1.00 - - 1.00 - 1.00 | -
>z2 | - -1.00 - - 1.00 -1.00 | 2.00
----|------------------------------------|-------
obj | - 2.00 1.00 -1.00 -1.00 -1.00 | -4.00
iter 0: x0 <-> x2 (1.0e+00 <-> 2.0e+00)
====|====================================|=======
| x1 x0 v0 v1 v2 U0 |
----|------------------------------------|-------
x2 | 1.00 -1.00 - - - - | 1.00
>Z0 | - 1.00 -1.00 - - 1.00 | 1.00
z1 |-1.00 - - 1.00 - 1.00 | -
>z2 |-1.00 1.00 - - 1.00 -1.00 | 1.00
----|------------------------------------|-------
obj | 2.00 -2.00 1.00 -1.00 -1.00 -1.00 | -2.00
iter 1: Z0 <-> v0 (1.0e+00 <-> 1.0e+00)
====|==============================|=======
| x1 x0 v1 v2 U0 |
----|------------------------------|-------
x2 | 1.00 -1.00 - - - | 1.00
v0 | - 1.00 - - 1.00 | 1.00
z1 |-1.00 - 1.00 - 1.00 | -
>z2 |-1.00 1.00 - 1.00 -1.00 | 1.00
----|------------------------------|-------
obj | 2.00 -1.00 -1.00 -1.00 - | -1.00
iter 2: z1 <-> x1 (0.0e+00 <-> 2.0e+00)
====|========================|=======
| x0 v1 v2 U0 |
----|------------------------|-------
x2 |-1.00 1.00 - 1.00 | 1.00
v0 | 1.00 - - 1.00 | 1.00
x1 | - 1.00 - 1.00 | -
>z2 | 1.00 -1.00 1.00 -2.00 | 1.00
----|------------------------|-------
obj |-1.00 1.00 -1.00 2.00 | -1.00
iter 3: z2 <-> U0 (1.0e+00 <-> 2.0e+00)
====|==================|=======
| x0 v1 v2 |
----|------------------|-------
x2 |-0.50 0.50 0.50 | 1.50
v0 | 1.50 -0.50 0.50 | 1.50
x1 | 0.50 0.50 0.50 | 0.50
U0 | 0.50 -0.50 0.50 | 0.50
----|------------------|-------
obj | - - - | - (0.0e+00)
---- solution: 0.00 0.50 1.50
---- problem 7 [right ans = (1.00,3.00)]
maxiter= 2002; rheps=1.7e-10
====|======================================================|=======
| x0 x1 v0 v1 u0 u1 u2 Z0 Z1 |
----|------------------------------------------------------|-------
>w0 |-1.00 - - - - - - - - | 1.00
>w1 | - -1.00 - - - - - - - | 2.00
>w2 | 1.00 -1.00 - - - - - - - | 2.00
z0 |-1.00 - 1.00 - 1.00 - -1.00 1.00 - | 1.00
z1 | - -1.00 - 1.00 - 1.00 1.00 - 1.00 | 2.00
----|------------------------------------------------------|-------
obj | - 2.00 - - - - - - - | -5.00
iter 0: w1 <-> x1 (2.0e+00 <-> 2.0e+00)
====|======================================================|=======
| x0 y1 v0 v1 u0 u1 u2 Z0 Z1 |
----|------------------------------------------------------|-------
>w0 |-1.00 - - - - - - - - | 1.00
x1 | - 1.00 - - - - - - - | 2.00
w2 | 1.00 -1.00 - - - - - - - | -
z0 |-1.00 - 1.00 - 1.00 - -1.00 1.00 - | 1.00
z1 | - -1.00 - 1.00 - 1.00 1.00 - 1.00 | -
----|------------------------------------------------------|-------
obj | - 1.00 - - - - - - - | -1.00
iter 1: w0 <-> x0 (1.0e+00 <-> 0.0e+00)
====|======================================================|=======
| y0 y1 v0 v1 u0 u1 u2 Z0 Z1 |
----|------------------------------------------------------|-------
x0 | 1.00 - - - - - - - - | 1.00
x1 | - 1.00 - - - - - - - | 2.00
>w2 | 1.00 -1.00 - - - - - - - | 1.00
z0 |-1.00 - 1.00 - 1.00 - -1.00 1.00 - | -
z1 | - -1.00 - 1.00 - 1.00 1.00 - 1.00 | -
----|------------------------------------------------------|-------
obj |-1.00 1.00 - - - - - - - | -1.00
iter 2: z1 <-> y1 (0.0e+00 <-> 1.0e+00)
====|======================================================|=======
| y0 z1 v0 v1 u0 u1 u2 Z0 Z1 |
----|------------------------------------------------------|-------
x0 | 1.00 - - - - - - - - | 1.00
x1 | - -1.00 - 1.00 - 1.00 1.00 - 1.00 | 2.00
>w2 | 1.00 1.00 - -1.00 - -1.00 -1.00 - -1.00 | 1.00
z0 |-1.00 - 1.00 - 1.00 - -1.00 1.00 - | -
y1 | - -1.00 - 1.00 - 1.00 1.00 - 1.00 | -
----|------------------------------------------------------|-------
obj |-1.00 -1.00 - 1.00 - 1.00 1.00 - 1.00 | -1.00
iter 3: w2 <-> Z1 (1.0e+00 <-> 1.0e+00)
====|======================================================|=======
| y0 z1 v0 v1 u0 u1 u2 Z0 y2 |
----|------------------------------------------------------|-------
x0 | 1.00 - - - - - - - - | 1.00
x1 | 1.00 - - - - - - - 1.00 | 3.00
Z1 | 1.00 1.00 - -1.00 - -1.00 -1.00 - 1.00 | 1.00
z0 |-1.00 - 1.00 - 1.00 - -1.00 1.00 - | -
y1 | 1.00 - - - - - - - 1.00 | 1.00
----|------------------------------------------------------|-------
obj | - - - - - - - - - | - (0.0e+00)
maxiter= 792; rheps=1.7e-10
====|==========================================|=======
| y0 v0 v1 u0 u1 u2 y2 |
----|------------------------------------------|-------
x0 | 1.00 - - - - - - | 1.00
x1 | 1.00 - - - - - 1.00 | 3.00
>Z1 | 1.00 - -1.00 - -1.00 -1.00 1.00 | 1.00
z0 |-1.00 1.00 - 1.00 - -1.00 - | -
y1 | 1.00 - - - - - 1.00 | 1.00
----|------------------------------------------|-------
obj | - -1.00 1.00 -1.00 1.00 2.00 -1.00 | -1.00
iter 0: z0 <-> y0 (0.0e+00 <-> 0.0e+00)
====|====================================|=======
| v0 v1 u0 u1 u2 y2 |
----|------------------------------------|-------
x0 | 1.00 - 1.00 - -1.00 - | 1.00
x1 | 1.00 - 1.00 - -1.00 1.00 | 3.00
>Z1 | 1.00 -1.00 1.00 -1.00 -2.00 1.00 | 1.00
y0 | 1.00 - 1.00 - -1.00 - | -
y1 | 1.00 - 1.00 - -1.00 1.00 | 1.00
----|------------------------------------|-------
obj |-1.00 1.00 -1.00 1.00 2.00 -1.00 | -1.00
iter 1: y0 <-> u2 (0.0e+00 <-> 2.0e+00)
====|====================================|=======
| v0 v1 u0 u1 y0 y2 |
----|------------------------------------|-------
x0 | - - - - 1.00 - | 1.00
x1 | - - - - 1.00 1.00 | 3.00
>Z1 |-1.00 -1.00 -1.00 -1.00 2.00 1.00 | 1.00
u2 | 1.00 - 1.00 - -1.00 - | -
y1 | - - - - 1.00 1.00 | 1.00
----|------------------------------------|-------
obj | 1.00 1.00 1.00 1.00 -2.00 -1.00 | -1.00
iter 2: Z1 <-> u0 (1.0e+00 <-> 1.0e+00)
====|==============================|=======
| v0 v1 u1 y0 y2 |
----|------------------------------|-------
x0 | - - - 1.00 - | 1.00
x1 | - - - 1.00 1.00 | 3.00
u0 |-1.00 -1.00 -1.00 2.00 1.00 | 1.00
u2 | - -1.00 -1.00 1.00 1.00 | 1.00
y1 | - - - 1.00 1.00 | 1.00
----|------------------------------|-------
obj | - - - - - | - (0.0e+00)
---- solution: 1.00 3.00
congratulations - all problems solved!
quadprogtest2 gives quadprog a more strenuous workout. It is based on the test program for Michael Powell’s lincoa. It uses the same set of 200 inequalities in 12 variables. It generates 1000 random points xt and for each point uses quadratic programming to find the point x't closest to xt while satisfying the constraints. It checks in each case that it has received a return code of 0 and that the solution does indeed satisfy the constraints. It prints out the largest residue from any solution, which I used as a guide for setting ε in quadprog (but which I found led to unduly tight tolerances).
I get the following output:
largest residue in absolute value is -2.5e-15 (iter 987)
#include "memory.h" #include "random.h" #include <math.h> int quadprog(double **Q,double *A,double *x,int n0,int n1, double **le,double *lerhs,int nle, double **ge,double *gerhs,int nge, double **eq,double *eqrhs,int neq) ; double quadprogresid() ; #define pi 3.141592653589793 int main() { int np=50,ia=12,n=12,m=200,i,j,k,iter,ret,epsiter ; double xp[50],yp[50],zp[50],theta,sumx,sumy,sumz,sigma=0.5,maxerr,maxeps,q ; double b[200],x[12],**a=matrix(m,n),**Q=matrix(12,12),A[12] ; for(sumx=sumy=sumz=j=0;j<np;j++) { theta = j * pi / (np-1) ; sumx += xp[j] = cos(theta) * cos(2*theta) ; sumy += yp[j] = sin(theta) * cos(2*theta) ; sumz += zp[j] = sin(2*theta) ; } sumx /= np ; sumy /= np ; sumz /= np ; for(j=0;j<np;j++) { xp[j] -= sumx ; yp[j] -= sumy ; zp[j] -= sumz ; } for(k=0;k<m;k++) b[k] = 1 ; for(i=0;i<n;i++) for(k=0;k<m;k++) a[k][i] = 0 ; for(j=0;j<np;j++) for(i=0;i<4;i++) { a[4*j+i][3*i ] = xp[j] ; a[4*j+i][3*i+1] = yp[j] ; a[4*j+i][3*i+2] = zp[j] ; } for(i=0;i<n;i++) Q[i][i] = -1 ; for(maxeps=iter=0;iter<1000;iter++) { ranset(iter) ; for(i=0;i<n;i++) A[i] = sigma * gaussv() ; ret = quadprog(Q,A,x,0,n, a,b,m, 0,0,0, 0,0,0) ; for(maxerr=i=0;i<m;i++) { for(sumx=-b[i],j=0;j<n;j++) sumx += a[i][j] * x[j] ; if(i==0||sumx>maxerr) maxerr = sumx ; } q = quadprogresid() ; if(maxerr>1e-10||ret!=0) printf("ret=%d, err=%.1e\n",ret,maxerr) ; else if(fabs(q)>fabs(maxeps)) { maxeps = q ; epsiter = iter ; } } printf("largest residue in absolute value is %.1e (iter %d)\n", maxeps,epsiter) ; freematrix(a,Q) ; }
This link generates a perl script which you can download and execute to compile and run quadprog. It works in Firefox only (and needs popups enabled for this site).
quadprog.c quadprogtest.c quadprogtest2.c utils: memory.h random.h ranf.c --- g++ -o quadprogtest -O2 -g quadprog.c quadprogtest.c g++ -o quadprogtest2 -O2 -g quadprog.c quadprogtest2.c ranf.c ./quadprogtest ./quadprogtest2
doc : algorithm : references : source : quadprogtest : quadprogtest2 : download
doc :
lincoa minimises or maximises a given function subject to the m constraints
ai · x ≤ bi
The prototype is
/* minimisers */ double lincoa(double (*f)(double *),double *x,int n,double **a,double *b, int m,double rhobeg,double rhoend,int maxfunc) ; double lincoa(double (*g)(double *,double *),double *x,int n, double *aux,int naux,double **a,double *b,int m, double rhobeg,double rhoend,int maxfunc) ; /* maximisers */ double lincoamax(double (*f)(double *),double *x,int n,double **a,double *b, int m,double rhobeg,double rhoend,int maxfunc) ; double lincoamax(double (*g)(double *,double *),double *x,int n, double *aux,int naux,double **a,double *b,int m, double rhobeg,double rhoend,int maxfunc) ; /* return condition */ int lincoacond() ;
where
The constraints in C notation are
∑j a[i][j]∗x[j] ≤ b[i]
The start point is not required to satisfy the constraints: if it doesn’t, it will be adjusted to do so before the hillclimb begins.
If the function has side values (i.e. is g rather than f) then they work as follows.
The values in the auxiliary array on input to lincoa will be the initial contents of the second argument to g on each call.
The values in the auxiliary array on return from lincoa will be the values of the second argument to g on return from the call which optimises g.
Side values may be used to supply additional implicit input parameters or to collect values computed at the optimum besides the main input and output. The former task can otherwise be performed by using variables whose scope is common to g and your calling program; the second can only be achieved otherwise by an additional call to g using variables of shared scope. This is inelegant as well as slightly inefficient.
On return x will hold the approximate location of the optimum and the return value will be the value of f/g as computed at that point. The algorithm terminates when one of the five following conditions is met:
A call to lincoacond() returns the value of cond corresponding to the most recent call to lincoa. Negative values of cond are set if you call lincoa with illegal parameters as follows:
References will be found in some of the same places as those for newuoa.
I made the conversion by first making a literal C++ translation using fable, and then by rewriting the functions one by one to look like manually produced code. In particular I allocate work arrays where they are used, I transposed Powell’s matrices, and I use C vector-of-vector notation for matrices. Most subroutines are now functions. goto’s have been abolished. (The fable translation has been preserved for posterity in the owl archive.)
Initially I was able to keep results which were bit-identical to those from the fable translation (except for the last place of decimals, where I suspect the difference lies in the behaviour of the I/O libraries invoked). I temporarily lost this equivalence and it later reappeared; I don’t regard it as essential.
Powell’s implementation stretches the constraints to include the initial point if the point supplied doesn’t satisfy them. I don’t imagine he thought that this was the right design for a general purpose subroutine, but he wasn’t interested in providing friendly software. I used my own quadratic program for this purpose although Powell had his own implementation of the the Goldfarb/Idnani method and says that this is embedded in lincoa (without it being visible there to the naked eye).
In the course of making this conversion I noticed some code shared with newuoa. I ripped this out as a separate library (mdplib). Powell’s duplicates weren’t identical; I’ve tried to take the better versions. In particular his update subroutine checks against β being zero in lincoa but not in newuoa. I retained the check and slightly modified the newuoa logic to take account of the possibility (with no guarantees, obvs.).
The return condition cond=4 was added by CJC under the following reasoning (based on an incomplete understanding of lincoa).
In an application I found lincoa getting into an infinite loop. I traced this down to the decision explained in the comment
Set VQUAD to the change to the quadratic model when the move STEP is made from XOPT. If STEP is a trust region step, then VQUAD should be negative. If it is nonnegative due to rounding errors in this case, there is a branch to label 530 to try to improve the model.
The trouble is that the attempt to improve the model may lead to VQUAD again being non-negative, and thus ad infinitum. But if VQUAD is evaluated as non-negative the quadratic model must be flat in the search direction. The jump to 530 generates a new search direction, but in the worst case the quadratic model will be flat in all directions and an infinite loop will ensue.
Now the behaviour of a quadratic function is determined by its behaviour in just n independent directions. Any reasonable algorithm for finding search directions (such as conjugate gradients which underlies lincoa) will find a direction which isn’t flat in at most n attempts if any such direction exists. Therefore limiting the number of directions tried to n will obviate the risk of infinite loops without sacrificing results for any soluble problem; so that is what I have done.
Powell had obviously observed that his algorithm didn’t work if he limited himself to a single search direction. I looked to see what happened if I limited the search to 2. In my application the limit was often reached but overall behaviour was unchanged (which probably doesn’t mean very much). My test program using Powell’s example function usually failed to find the true minimum.
When I limited the number of search directions to 3 I found that all results were identical to those obtained with no limit except that I broke out of the the one problem which otherwise led to an infinite loop. So I hope that a limit of n will indeed be safe.
At the same time as making this change I put in an explicit limit on the number of iterations (which is not the same as the number of function evaluations). This is because I think it better to have such a check than to rely on deep logic to guarantee termination (which even Powell hadn’t fully worked through). lincoa throws up if the number of iterations is more than n×maxfunc.
• rearrange • getact • trstep • qmstep • lincupdate • prelim • lincoacond • lincob • lincoastub • lincoa • lincoamax
/* The original fortran codes are distributed without restrictions. The C++ codes are distributed under MIT license. */ /* The MIT License Copyright (c) nd., by M.J.D. Powell <mjdp@cam.ac.uk> 2016, by Colin Champion <colin.champion@masterlyinactivity.com> Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. */ #include"memory.h" #include <math.h> #define max(a,b) ((a)>(b)?(a):(b)) #define min(a,b) ((a)<(b)?(a):(b)) static int visit=-1 ; int update(int,int,double **,double **,int,double *,double,int *) ; void shiftxbase(int n,int npt,double **xp,double *xopt,double *vlag,double **b, double **z,double *pq,double *hq,int idz) ; double calcvlagbeta(int n,int npt,double *vlag,double **z,double **b,double *w, double *xopt,double *d,double beta,double dsq,double xsq, int idz) ; int chooseknew(int n,int npt,double beta,double **z,double *vlag,double **xp, double *xopt,double detrat,double rhosq,int ktemp,int idz) ; void genpqw(int n,int npt,double *pqw,double **z,int knew,int idz) ; void updatehq(int n,int npt,double *hq,double pqknew,double **xpt,int knew) ; void updategopt(int n,double *gopt,double *hq,double *step) ; void updategopt2(int n,int npt,double *gopt,double *xopt,double *pq,double **) ; //C These instructions rearrange the active constraints so that the new //C value of IACT(NACT) is the old value of IACT(IC). A sequence of //C Givens rotations is applied to the current QFAC and RFAC. Then NACT //C is reduced by one. static int rearrange(int ic,int nact,int n,int *iact,double *resact, double *rfac,double **qfac,double *vlam,double *resnew) { int jc,i,j,idiag,jw,jdiag ; double temp,cval,sval ; if((resnew[iact[ic]]=resact[ic])<1e-60) resnew[iact[ic]] = 1e-60 ; for(jc=ic;jc<nact-1;jc++) { jdiag = (jc*(jc+1)) / 2 ; idiag = ((jc+1)*(jc+2)) / 2 ; jw = idiag + jc ; temp = sqrt(rfac[jw]*rfac[jw]+rfac[jw+1]*rfac[jw+1]) ; cval = rfac[jw+1] / temp ; sval = rfac[jw] / temp ; rfac[jw] = sval * rfac[idiag-1] ; rfac[jw+1] = cval * rfac[idiag-1] ; rfac[idiag-1] = temp ; for(jw+=2,j=jc+2;j<nact;j++,jw+=j) { temp = sval*rfac[jw+jc] + cval*rfac[jw+jc+1] ; rfac[jw+jc+1] = cval*rfac[jw+jc] - sval*rfac[jw+jc+1] ; rfac[jw+jc] = temp ; } for(i=0;i<jc;i++) swap(rfac[idiag+i],rfac[jdiag+i]) ; for(i=0;i<n;i++) { temp = sval*qfac[jc][i] + cval*qfac[jc+1][i] ; qfac[jc+1][i] = cval*qfac[jc][i] - sval*qfac[jc+1][i] ; qfac[jc][i] = temp ; } iact[jc] = iact[jc+1] ; resact[jc] = resact[jc+1] ; vlam[jc] = vlam[jc+1] ; } return nact-1 ; } /* -------------------------------------------------------------------------- */ double getact(int n,int m,double **amat, double *b,int& nact, int *iact,double **qfac,double *rfac,double & snorm, double *resnew,double *resact,double *g,double *dw) { //C N, M, AMAT, B, NACT, IACT, QFAC and RFAC are the same as the terms //C with these names in SUBROUTINE LINCOB. The current values must be //C set on entry. NACT, IACT, QFAC and RFAC are kept up to date when //C GETACT changes the current active set. //C SNORM, RESNEW, RESACT, G and DW are the same as the terms with these //C names in SUBROUTINE TRSTEP. The elements of RESNEW and RESACT are //C also kept up to date. //C VLAM and W are used for working space, the vector VLAM being reserved //C for the Lagrange multipliers of the calculation. Their lengths must //C be at least N. //C The main purpose of GETACT is to pick the current active set. It is //C defined by the property that the projection of -G into the space //C orthogonal to the active constraint normals is as large as possible, //C subject to this projected steepest descent direction moving no closer //C to the boundary of every constraint whose current residual is at most //C 0.2*SNORM. On return, the settings in NACT, IACT, QFAC and RFAC are //C all appropriate to this choice of active set. //C Occasionally this projected direction is zero, and then the final value //C of W(1) is set to zero. Otherwise, the direction itself is returned //C in DW, and W(1) is set to the square of the length of the direction. //C //C Set some constants and a temporary VLAM. int i,j,k,ic,idiag,jw,l ; double tdel=0.2*snorm,ddsav,violmx,rdiag,sprod,temp,ctol ; double cosv,sinv,sum,dd,test,dnorm,vmult ; double *vlam=vector(n),*w=vector(n) ; for(ddsav=i=0;i<n;i++) { ddsav += g[i]*g[i] ; vlam[i] = 0 ; } ddsav *= 2 ; //C Set the initial QFAC to the identity matrix in the case NACT=0. if(nact==0) for(i=0;i<n;i++) { for(j=0;j<i;j++) qfac[j][i] = qfac[i][j] = 0 ; qfac[i][i] = 1 ; } else //C Remove any constraints from the initial active set whose residuals //C exceed TDEL. { for(ic=nact-1;ic>=0;ic--) if(resact[ic]>tdel) nact = rearrange(ic,nact,n,iact,resact,rfac,qfac,vlam,resnew) ; //C Remove any constraints from the initial active set whose Lagrange //C multipliers are nonnegative, and set the surviving multipliers. for(ic=nact-1;ic>=0;ic--) { for(temp=i=0;i<n;i++) temp += qfac[ic][i] * g[i] ; idiag = ((ic+1)*(ic+2))/2 ; for(jw=idiag+ic,j=ic+1;j<nact;j++,jw+=j) temp -= rfac[jw] * vlam[j] ; if(temp>=0) ic = nact = rearrange(ic,nact,n,iact,resact,rfac,qfac,vlam,resnew) ; else vlam[ic] = temp / rfac[idiag-1] ; } } //C Set the new search direction D. Terminate if the 2-norm of D is zero //C or does not decrease, or if NACT=N holds. The situation NACT=N //C occurs for sufficiently large SNORM if the origin is in the convex //C hull of the constraint gradients. while(nact<n) { for(j=nact;j<n;j++) for(w[j]=i=0;i<n;i++) w[j] += qfac[j][i] * g[i] ; for(dd=i=0;i<n;i++) { for(dw[i]=0,j=nact;j<n;j++) dw[i] -= w[j] * qfac[j][i] ; dd += dw[i] * dw[i] ; } if(dd>=ddsav||dd==0) { free(vlam,w) ; return 0 ; } ddsav = dd ; dnorm = sqrt(dd) ; //C Pick the next integer L or terminate, a positive value of L being //C the index of the most violated constraint. The purpose of CTOL //C below is to estimate whether a positive value of VIOLMX may be //C due to computer rounding errors. l = -1 ; if(m>0) { test = dnorm / snorm ; for(violmx=j=0;j<m;j++) if(resnew[j]>0&&resnew[j]<=tdel) { for(sum=i=0;i<n;i++) sum += amat[j][i] * dw[i] ; if(sum>test*resnew[j]&&sum>violmx) { l = j ; violmx = sum ; } } ctol= 0 ; temp = 0.01 * dnorm ; if(violmx>0&&violmx<temp) for(k=0;k<nact;k++) { for(j=iact[k],sum=i=0;i<n;i++) sum += dw[i] * amat[j][i] ; if(fabs(sum)>ctol) ctol = fabs(sum) ; } } if(l<0||violmx<=10*ctol) { free(vlam,w) ; return dd ; } //C Apply Givens rotations to the last (N-NACT) columns of QFAC so that //C the first (NACT+1) columns of QFAC are the ones required for the //C addition of the L-th constraint, and add the appropriate column //C to RFAC. idiag = (nact*(nact+1))/2 ; for(rdiag=0,j=n-1;j>=0;j--) { for(sprod=i=0;i<n;i++) sprod += qfac[j][i] * amat[l][i] ; if(j<nact) rfac[idiag+j] = sprod ; else if(fabs(rdiag)<=1e-20*fabs(sprod)) rdiag = sprod ; else { temp = sqrt(sprod*sprod+rdiag*rdiag) ; cosv = sprod / temp ; sinv = rdiag / temp ; rdiag = temp ; for(i=0;i<n;i++) { temp = cosv*qfac[j][i] + sinv*qfac[j+1][i] ; qfac[j+1][i] = -sinv*qfac[j][i] + cosv*qfac[j+1][i] ; qfac[j][i] = temp ; } } } if(rdiag<0) for(i=0;i<n;i++) qfac[nact][i] = -qfac[nact][i] ; rfac[idiag+nact] = fabs(rdiag) ; iact[nact] = l ; resact[nact] = resnew[l] ; vlam[nact] = resnew[l] = 0 ; nact += 1 ; //C Set the components of the vector VMU in W. for(ic=0;violmx>0||ic>=0;) // while(violmx>0) but with at least 1 iteration { temp = rfac[(nact*(nact+1))/2-1] ; w[nact-1] = 1/(temp*temp) ; for(i=nact-2;i>=0;i--) { idiag = ((i+1)*(i+2))/2 ; for(sum=0,jw=idiag+i,j=i+1;j<nact;j++,jw+=j) sum -= rfac[jw] * w[j] ; w[i] = sum / rfac[idiag-1] ; } //C Calculate the multiple of VMU to subtract from VLAM, and update VLAM. vmult = violmx ; for(ic=-1,j=0;j<nact-1;j++) if(vlam[j]>=vmult*w[j]) { ic = j ; vmult = vlam[j] / w[j] ; } for(j=0;j<nact;j++) vlam[j] -= vmult * w[j] ; if(ic>=0) vlam[ic] = 0 ; violmx -= vmult ; if(ic<0||violmx<0) violmx = 0 ; //C Reduce the active set if necessary, so that all components of the //C new VLAM are negative, with resetting of the residuals of the //C constraints that become inactive. for(ic=nact-1;ic>=0;ic--) if(vlam[ic]>=0) nact = rearrange(ic,nact,n,iact,resact,rfac,qfac,vlam,resnew) ; //C Calculate the next VMU if VIOLMX is positive. Return if NACT=N holds, //C as then the active constraints imply D=0. Otherwise, go to label //C 100, to calculate the new D and to test for termination. } // end while(violmx>0) } // end while(nact<n) free(vlam,w) ; return 0 ; } /* -------------------------------------------------------------------------- */ int trstep(int n,int npt,int m,double **amat,double *b,double **xpt,double *hq, double *pq,int& nact,int *iact,double *rescon,double **qfac, double *rfac,double& snorm,double *step,double *g) { //C N, NPT, M, AMAT, B, XPT, HQ, PQ, NACT, IACT, RESCON, QFAC and RFAC //C are the same as the terms with these names in LINCOB. If RESCON(J) //C is negative, then |RESCON(J)| must be no less than the trust region //C radius, so that the J-th constraint can be ignored. //C SNORM is set to the trust region radius DELTA initially. On the //C return, however, it is the length of the calculated STEP, which is //C set to zero if the constraints do not allow a long enough step. //C STEP is the total calculated step so far from the trust region centre, //C its final value being given by the sequence of CG iterations, which //C terminate if the trust region boundary is reached. //C G must be set on entry to the gradient of the quadratic model at the //C trust region centre. It is used as working space, however, and is //C always the gradient of the model at the current STEP, except that //C on return the value of G(1) is set to ONE instead of to ZERO if //C and only if GETACT is called more than once. //C RESNEW, RESACT, D, DW and W are used for working space. A negative //C value of RESNEW(J) indicates that the J-th constraint does not //C restrict the CG steps of the current trust region calculation, a //C zero value of RESNEW(J) indicates that the J-th constraint is active, //C and otherwise RESNEW(J) is set to the greater of TINY and the actual //C residual of the J-th constraint for the current STEP. RESACT holds //C the residuals of the active constraints, which may be positive. //C D is the search direction of each line search. DW is either another //C search direction or the change in gradient along D. The length of W //C must be at least MAX[M,2*N]. int i,j,k,continuing,ncall,ir,jsav,icount,ih ; double ctest=0.01,snsq=snorm*snorm,tiny=1e-60,dnorm,scale,resmax,ss,sum,rhs ; double temp,gamma,sumrhs,dd,ds,ad,adw,dg,dgd,alpha,alphm,alpht,alpbd,reduct ; double wgd,beta ; double *resnew=vector(m),*resact=vector(n),*d=vector(n),*dw=vector(n) ; double *w=vector(max(m,2*n)) ; static int thingy=0 ; //C Set the initial elements of RESNEW, RESACT and STEP. for(j=0;j<m;j++) { resnew[j] = rescon[j] ; if(resnew[j]>=snorm) resnew[j] = -1 ; else if(resnew[j]>=0&&resnew[j]<tiny) resnew[j] = tiny ; } if(m>0) for(k=0;k<nact;k++) { resact[k] = rescon[iact[k]] ; resnew[iact[k]] = 0 ; } for(i=0;i<n;i++) step[i] = 0 ; //C GETACT picks the active set for the current STEP. It also sets DW to //C the vector closest to -G that is orthogonal to the normals of the //C active constraints. DW is scaled to have length 0.2*SNORM, as then //C a move of DW from STEP is allowed by the linear constraints. for(ss=reduct=ncall=0,continuing=1;continuing;) { dnorm = getact(n,m,amat,b,nact,iact,qfac,rfac,snorm,resnew,resact,g,dw) ; ncall += 1 ; if(dnorm==0) break ; scale = 0.2*snorm/sqrt(dnorm) ; for(i=0;i<n;i++) dw[i] *= scale ; //C If the modulus of the residual of an active constraint is substantial, //C then set D to the shortest move from STEP to the boundaries of the //C active constraints. for(resmax=k=0;k<nact;k++) if(resact[k]>resmax) resmax = resact[k] ; gamma = 0 ; if(resmax>snorm*1e-4) { for(ir=k=0;k<nact;k++,ir++) { for(temp=resact[k],i=0;i<k;ir++,i++) temp -= rfac[ir] * w[i] ; w[k] = temp / rfac[ir] ; } for(i=0;i<n;i++) for(d[i]=k=0;k<nact;k++) d[i] += w[k] * qfac[k][i] ; //C The vector D that has just been calculated is also the shortest move //C from STEP+DW to the boundaries of the active constraints. Set GAMMA //C to the greatest steplength of this move that satisfies the trust //C region bound. for(rhs=snsq,ds=dd=i=0;i<n;i++) { sum = step[i] + dw[i] ; rhs -= sum * sum ; ds += d[i] * sum ; dd += d[i] * d[i] ; } if(rhs>0) { temp = sqrt(ds*ds+dd*rhs) ; if(ds<=0) gamma = (temp-ds) / dd ; else gamma = rhs / (temp+ds) ; } //C Reduce the steplength GAMMA if necessary so that the move along D //C also satisfies the linear constraints. for(j=0;j<m&&gamma>0;j++) if(resnew[j]>0) { for(ad=adw=i=0;i<n;i++) { ad += amat[j][i] * d[i] ; adw += amat[j][i] * dw[i] ; } if(ad>0) { if(0>(temp=(resnew[j]-adw)/ad)) temp = 0 ; if(temp<gamma) gamma = temp ; } } if(gamma>1) gamma = 1 ; } // end if(resmax>snorm*1e-4) //C Set the next direction for seeking a reduction in the model function //C subject to the trust region bound and the linear constraints. if(gamma<=0) { for(i=0;i<n;i++) d[i] = dw[i] ; icount = nact ; } else { for(i=0;i<n;i++) d[i] = dw[i] + gamma*d[i] ; icount = nact-1 ; } //C Set ALPHA to the steplength from STEP along D to the trust region //C boundary. Return if the first derivative term of this step is //C sufficiently small or if no further progress is possible. for(alpbd=1;;alpbd=0) { icount += 1 ; rhs = snsq - ss ; if(rhs<=0) { continuing = 0 ; break ; } // break both loops for(dg=ds=dd=i=0;i<n;i++) { dg += d[i] * g[i] ; ds += d[i] * step[i] ; dd += d[i] * d[i] ; } if(dg>=0) { continuing = 0 ; break ; } temp = sqrt(rhs*dd+ds*ds) ; if(ds<=0) alpha = (temp-ds)/dd ; else alpha = rhs/(temp+ds) ; if(-alpha*dg<=ctest*reduct) { continuing = 0 ; break ; } //C Set DW to the change in gradient along D. for(j=0;j<n;j++) dw[j] = 0 ; updategopt(n,dw,hq,d) ; updategopt2(n,npt,dw,d,pq,xpt) ; //C Set DGD to the curvature of the model along D. Then reduce ALPHA if //C necessary to the value that minimizes the model. for(dgd=i=0;i<n;i++) dgd += d[i] * dw[i] ; alpht = alpha ; if(dg+alpha*dgd>0) alpha = -dg/dgd ; //C Make a further reduction in ALPHA if necessary to preserve feasibility, //C and put some scalar products of D with constraint gradients in W. alphm = alpha ; for(jsav=-1,j=0;j<m;j++) { w[j] = 0 ; if(resnew[j]>0) { for(i=0;i<n;i++) w[j] += amat[j][i] * d[i] ; if(alpha*w[j]>resnew[j]) { alpha = resnew[j] / w[j] ; jsav = j ; } } } if(alpha<alpbd) alpha = alpbd ; if(alpha>alphm) alpha = alphm ; if(icount==nact&&alpha>1) alpha = 1 ; //C Update STEP, G, RESNEW, RESACT and REDUCT. for(ss=i=0;i<n;i++) { step[i] += alpha*d[i] ; ss += step[i]*step[i] ; g[i] += alpha*dw[i] ; } for(j=0;j<m;j++) if(resnew[j]>0) if((resnew[j]-=alpha*w[j])<tiny) resnew[j] = tiny ; if(icount==nact) for(k=0;k<nact;k++) resact[k] *= (1-gamma) ; reduct -= alpha * (dg+0.5*alpha*dgd) ; //C Test for termination. Branch to label 40 if there is a new active //C constraint and if the distance from STEP to the trust region //C boundary is at least 0.2*SNORM. if( alpha==alpht || -alphm*(dg+0.5*alphm*dgd)<=ctest*reduct ) { continuing = 0 ; break ; } if(jsav>=0) { continuing = (ss<=0.64*snsq) ; break ; } if(icount==n) { continuing = 0 ; break ; } //C Calculate the next search direction, which is conjugate to the //C previous one except in the case ICOUNT=NACT. if(nact>0) { for(j=nact;j<n;j++) for(w[j]=i=0;i<n;i++) w[j] += g[i] * qfac[j][i] ; for(i=0;i<n;i++) { for(temp=0,j=nact;j<n;j++) temp += qfac[j][i]*w[j] ; w[n+i] = temp ; } } else for(i=0;i<n;i++) w[n+i] = g[i] ; if(icount==nact) beta = 0 ; else { for(wgd=i=0;i<n;i++) wgd += w[n+i]*dw[i] ; beta = wgd / dgd ; } for(i=0;i<n;i++) d[i] = -w[n+i] + beta*d[i] ; } } //C Return from the subroutine. if(reduct>0) snorm = sqrt(ss) ; else snorm = 0 ; free(resnew,resact,dw,d,w) ; return ncall ; } /* -------------------------------------------------------------------------- */ int qmstep(int n,int npt,int m,double **amat,double **xpt,double *xopt,int nact, int *iact,double *rescon,double **qfac,int kopt,int knew,double del, double *step,double *gl,double *pqw) { //C N, NPT, M, AMAT, B, XPT, XOPT, NACT, IACT, RESCON, QFAC, KOPT are the //C same as the terms with these names in SUBROUTINE LINCOB. //C KNEW is the index of the interpolation point that is going to be moved. //C DEL is the current restriction on the length of STEP, which is never //C greater than the current trust region radius DELTA. //C STEP will be set to the required step from XOPT to the new point. //C GL must be set on entry to the gradient of LFUNC at XBASE, where LFUNC //C is the KNEW-th Lagrange function. It is used also for some other //C gradients of LFUNC. //C PQW provides the second derivative parameters of LFUNC. //C RSTAT and W are used for working space. Their lengths must be at least //C M and N, respectively. RSTAT(J) is set to -1.0, 0.0, or 1.0 if the //C J-th constraint is irrelevant, active, or both inactive and relevant, //C respectively. //C IFEAS will be set to 0 or 1 if XOPT+STEP is infeasible or feasible. //C //C STEP is chosen to provide a relatively large value of the modulus of //C LFUNC(XOPT+STEP), subject to ||STEP|| .LE. DEL. A projected STEP is //C calculated too, within the trust region, that does not alter the //C residuals of the active constraints. The projected step is preferred //C if its value of | LFUNC(XOPT+STEP) | is at least one fifth of the //C original one, but the greatest violation of a linear constraint must //C be at least 0.2*DEL, in order to keep the interpolation points apart. //C The remedy when the maximum constraint violation is too small is to //C restore the original step, which is perturbed if necessary so that //C its maximum constraint violation becomes 0.2*DEL. int i,j,k,ksav,jsav,ifeas,*istat=ivector(m) ; double test=0.2*del,temp,vbig,ss,sp,stp,vlag,stpsav,gg,vgrad,ghg,ww,ctol ; double sum,vnew,resmax,bigv,*w=vector(n) ; //C Replace GL by the gradient of LFUNC at the trust region centre, and //C set the elements of RSTAT. for(k=0;k<npt;k++) { for(temp=j=0;j<n;j++) temp += xpt[k][j] * xopt[j] ; temp *= pqw[k] ; for(i=0;i<n;i++) gl[i] += temp*xpt[k][i] ; } for(j=0;j<m;j++) { if(fabs(rescon[j])>=del) istat[j] = -1 ; else istat[j] = 1 ; } if(m>0) for(k=0;k<nact;k++) istat[iact[k]] = 0 ; //C Find the greatest modulus of LFUNC on a line through XOPT and //C another interpolation point within the trust region. for(ksav=-1,vbig=k=0;k<npt;k++) if(k!=kopt) { for(sp=ss=i=0;i<n;i++) { temp = xpt[k][i] - xopt[i] ; ss += temp * temp ; sp += gl[i] * temp ; } stp = -del / sqrt(ss) ; if(k==knew) { if(sp*(sp-1)<0) stp = -stp ; vlag = fabs(stp*sp) + stp*stp*fabs(sp-1) ; } else vlag = fabs(stp*(1-stp)*sp) ; if(vlag>vbig) { ksav = k ; stpsav = stp ; vbig = vlag ; } } if(ksav<0) // this was an error condition in Powell's code { printf("*** lincoa %d: qmstep didn\'t find a nonzero modulus\n",visit) ; free(istat,w) ; return -1 ; } //C Set STEP to the move that gives the greatest modulus calculated above. //C This move may be replaced by a steepest ascent step from XOPT. for(gg=i=0;i<n;i++) { gg += gl[i]*gl[i] ; step[i] = stpsav*(xpt[ksav][i]-xopt[i]) ; } vgrad = del * sqrt(gg) ; if(vgrad>vbig/10) { //C Make the replacement if it provides a larger value of VBIG. for(ghg=k=0;k<npt;k++) { for(temp=j=0;j<n;j++) temp += xpt[k][j] * gl[j] ; ghg += pqw[k] * temp * temp ; } vnew = vgrad + fabs(0.5*del*del*ghg/gg) ; if(vnew>vbig) { vbig = vnew ; stp = del / sqrt(gg) ; if(ghg<0) stp = -stp ; for(i=0;i<n;i++) step[i] = stp * gl[i] ; } if(nact!=0&&nact!=n) { //C Overwrite GL by its projection. Then set VNEW to the greatest //C value of |LFUNC| on the projected gradient from XOPT subject to //C the trust region bound. If VNEW is sufficiently large, then STEP //C may be changed to a move along the projected gradient. for(k=nact;k<n;k++) for(w[k]=i=0;i<n;i++) w[k] += gl[i] * qfac[k][i] ; for(gg=i=0;i<n;i++) { for(gl[i]=0,k=nact;k<n;k++) gl[i] += qfac[k][i] * w[k] ; gg += gl[i] * gl[i] ; } vgrad = del * sqrt(gg) ; if(vgrad>vbig/10) { for(ghg=k=0;k<npt;k++) { for(temp=j=0;j<n;j++) temp += xpt[k][j] * gl[j] ; ghg += pqw[k] * temp * temp ; } vnew = vgrad + fabs(0.5*del*del*ghg/gg) ; //C Set W to the possible move along the projected gradient. stp = del / sqrt(gg) ; if(ghg<0) stp = -stp ; for(ww=i=0;i<n;i++) { w[i] = stp * gl[i] ; ww += w[i]*w[i] ; } //C Set STEP to W if W gives a sufficiently large value of the modulus //C of the Lagrange function, and if W either preserves feasibility //C or gives a constraint violation of at least 0.2*DEL. The purpose //C of CTOL below is to provide a check on feasibility that includes //C a tolerance for contributions from computer rounding errors. if(vnew/vbig>=0.2) { for(ifeas=1,bigv=j=0;j<m;j++) { if(istat[j]==1) { for(temp=-rescon[j],i=0;i<n;i++) temp += w[i] * amat[j][i] ; if(temp>bigv) bigv = temp ; } if(bigv>=test) { ifeas = 0 ; break ; } } ctol = 0 ; temp = sqrt(ww) / 100 ; if(bigv>0&&bigv<temp) for(k=0;k<nact;k++) { for(j=iact[k],sum=i=0;i<n;i++) sum += w[i] * amat[j][i] ; if(fabs(sum)>ctol) ctol = fabs(sum) ; } if(bigv<=10*ctol||bigv>=test) { for(i=0;i<n;i++) step[i] = w[i] ; free(istat,w) ; return ifeas ; } } } } } //C Calculate the greatest constraint violation at XOPT+STEP with STEP at //C its original value. Modify STEP if this violation is unacceptable. for(ifeas=1,bigv=resmax=jsav=j=0;j<m;j++) if(istat[j]>=0) { for(temp=-rescon[j],i=0;i<n;i++) temp += step[i] * amat[j][i] ; if(resmax<temp) resmax = temp ; if(temp>=test) { ifeas = 0 ; break ; } if(temp>bigv) { bigv = temp ; jsav = j ; ifeas = -1 ; } } if(ifeas<0) { for(i=0;i<n;i++) step[i] += (test-bigv) * amat[jsav][i] ; ifeas = 0 ; } free(istat,w) ; return ifeas ; } /* -------------------------------------------------------------------------- */ void lincupdate(int n,int npt,double **xpt,double **bmat,double **zmat,int &idz, int ndim,double *sp,double *step,int kopt,int &knew) { int i,j,k,nptm=npt-n-1,jp ; double sum,ssq,temp,tempa,tempb,distsq,hdiag,denom,denabs,denmax ; double beta,*vlag=vector(npt+n),*w=vector(npt) ; //C The arguments N, NPT, XPT, BMAT, ZMAT, IDZ, NDIM ,SP and STEP are //C identical to the corresponding arguments in SUBROUTINE LINCOB. //C KOPT is such that XPT(KOPT,.) is the current trust region centre. //C KNEW on exit is usually positive, and then it is the index of an //C interpolation point to be moved to the position XPT(KOPT,.)+STEP(.). //C It is set on entry either to its final value or to 0. In the latter //C case, the final value of KNEW is chosen to maximize the denominator //C of the matrix updating formula times a weighting factor. //C VLAG and W are used for working space, the first NPT+N elements of //C both of these vectors being required. //C //C The arrays BMAT and ZMAT with IDZ are updated, the new matrices being //C the ones that are suitable after the shift of the KNEW-th point to //C the new position XPT(KOPT,.)+STEP(.). A return with KNEW set to zero //C occurs if the calculation fails due to a zero denominator in the //C updating formula, which should never happen. for(k=0;k<npt;k++) { w[k] = sp[npt+k] * (sp[npt+k]/2+sp[k]) ; for(vlag[k]=j=0;j<n;j++) vlag[k] += bmat[j][k] * step[j] ; } for(ssq=i=0;i<n;i++) ssq += step[i] * step[i] ; beta = calcvlagbeta(n,npt,vlag,zmat,bmat,w,xpt[kopt],step,beta,ssq,sp[kopt],idz) ; vlag[kopt] += 1 ; //C If KNEW is zero initially, then pick the index of the interpolation //C point to be deleted, by maximizing the absolute value of the //C denominator of the updating formula times a weighting factor. if(knew<0) knew = chooseknew(n,npt,beta,zmat,vlag,xpt,xpt[kopt],-1,0,-1,idz) ; if(knew<0) throw up("knew=%d in lincupdate; visit=%d\n",knew,visit) ; //C and call the NEWUOA update to do the rest idz = update(n,npt,bmat,zmat,idz,vlag,beta,&knew) ; free(vlag,w) ; } /* -------------------------------------------------------------------------- */ int prelim(int n,int npt,int m,double **amat,double *b,double *x,double rhobeg, double *xbase,double **xpt,double *fval,double *xsav,double *gopt, int &kopt,double *pq,double **bmat,double **zmat,int ndim, double *sp,double *rescon,int maxflag,double (*func)(double *), double (*gunc)(double *,double *),double *aux,int naux) { int nptm=npt-n-1,kbase=0,i,j,k,itemp,ipt,jpt,jsav,jp,nf,idz ; double rhosq=rhobeg*rhobeg,recip=1/rhosq,test=0.2*rhobeg,feas,bigv,resid,f ; double temp,reciq=sqrt(0.5)/rhosq ; double *step=vector(n),*w=vector(n+npt),*auxtemp ; if(func) auxtemp = 0 ; else auxtemp = vector(naux) ; //C The arguments N, NPT, M, AMAT, B, X, RHOBEG, IPRINT, XBASE, XPT, FVAL, //C XSAV, XOPT, GOPT, HQ, PQ, BMAT, ZMAT, NDIM, SP and RESCON are the //C same as the corresponding arguments in SUBROUTINE LINCOB. //C KOPT is set to the integer such that XPT(KOPT,.) is the initial trust //C region centre. //C IDZ is going to be set to one, so that every element of Diag(DZ) is //C one in the product ZMAT times Diag(DZ) times ZMAT^T, which is the //C factorization of the leading NPT by NPT submatrix of H. //C STEP, PQW and W are used for working space, the arrays STEP and PQW //C being taken from LINCOB. The length of W must be at least N+NPT. //C //C SUBROUTINE PRELIM provides the elements of XBASE, XPT, BMAT and ZMAT //C for the first iteration, an important feature being that, if any of //C of the columns of XPT is an infeasible point, then the largest of //C the constraint violations there is at least 0.2*RHOBEG. It also sets //C the initial elements of FVAL, XOPT, GOPT, HQ, PQ, SP and RESCON. //C Set the initial elements of XPT, BMAT, SP and ZMAT to zero. idz = 0 ; for(j=0;j<n;j++) { xbase[j] = x[j] ; for(k=0;k<npt;k++) xpt[k][j] = 0 ; for(i=0;i<ndim;i++) bmat[j][i] = 0 ; } for(k=0;k<npt;k++) { sp[k] = 0 ; for(j=0;j<nptm;j++) zmat[j][k] = 0 ; } //C Set the nonzero coordinates of XPT(K,.), K=1,2,...,min[2*N+1,NPT], //C but they may be altered later to make a constraint violation //C sufficiently large. The initial nonzero elements of BMAT and of //C the first min[N,NPT-N-1] columns of ZMAT are set also. //C for(j=0;j<n;j++) { xpt[j+1][j] = rhobeg ; if(j<nptm) { jp = n + j + 1 ; xpt[jp][j] = -rhobeg ; bmat[j][j+1] = rhobeg/2 ; bmat[j][jp] = -rhobeg/2 ; zmat[j][0] = -2*reciq ; zmat[j][j+1] = zmat[j][jp] = reciq ; } else { bmat[j][0] = -(bmat[j][j+1] = 1/rhobeg) ; bmat[j][npt+j] = -rhosq/2 ; } } //C Set the remaining initial nonzero elements of XPT and ZMAT when the //C number of interpolation points exceeds 2*N+1. for(k=n;k<nptm;k++) { itemp = k / n ; ipt = k - itemp*n ; jpt = ipt + itemp ; if(jpt>=n) jpt -= n ; xpt[n+k+1][ipt] = xpt[n+k+1][jpt] = rhobeg ; zmat[k][0] = zmat[k][n+k+1] = recip ; zmat[k][ipt+1] = zmat[k][jpt+1] = -recip ; } //C Update the constraint right hand sides to allow for the shift XBASE. for(j=0;j<m;j++) { for(temp=i=0;i<n;i++) temp += amat[j][i] * xbase[i] ; b[j] -= temp ; } //C Go through the initial points, shifting every infeasible point if //C necessary so that its constraint violation is at least test (=0.2*RHOBEG) for(nf=0;nf<npt;nf++) // n = no of variables, m = no of equations { feas = 1 ; bigv = 0 ; if(nf>0) for(j=0;j<m;j++) { for(resid=-b[j],i=0;i<n;i++) resid += xpt[nf][i] * amat[j][i] ; if(resid<=bigv) continue ; bigv = resid ; jsav = j ; if(resid>test) { feas = 0 ; break ; } else feas = -1 ; } if(feas<0) { for(i=0;i<n;i++) step[i] = xpt[nf][i] + (test-bigv)*amat[jsav][i] ; for(k=0;k<npt;k++) for(sp[npt+k]=j=0;j<n;j++) sp[npt+k] += xpt[k][j]*step[j] ; lincupdate(n,npt,xpt,bmat,zmat,idz,ndim,sp,step,kbase,nf) ; for(i=0;i<n;i++) xpt[nf][i] = step[i] ; } //C Calculate the objective function at the current interpolation point, //C and set KOPT to the index of the first trust region centre. for(j=0;j<n;j++) x[j] = xbase[j] + xpt[nf][j] ; for(j=0;j<n;j++) if(abnormal(x[j])) throw up("1: %.1e\n",x[j]) ; if(func) f = func(x) ; else { for(j=0;j<naux;j++) auxtemp[j] = aux[j] ; f = gunc(x,naux?auxtemp:0) ; } if(isnan(f)||isinf(f)) throw up("function returns %.1e",f) ; if(maxflag) f = -f ; if(nf==0||(f<fval[kopt]&&feas>0)) { kopt = nf ; for(j=0;j<naux;j++) xsav[n+j] = auxtemp[j] ; } fval[nf] = f ; } //C Set PQ for the first quadratic model. for(j=0;j<nptm;j++) for(w[j]=k=0;k<npt;k++) w[j] += zmat[j][k]*fval[k] ; for(k=0;k<npt;k++) for(pq[k]=j=0;j<nptm;j++) pq[k] += zmat[j][k] * w[j] ; //C Set XOPT, SP, GOPT and HQ for the first quadratic model. for(j=0;j<n;j++) { xsav[j] = xbase[j] + xpt[kopt][j] ; gopt[j] = 0 ; } for(k=0;k<npt;k++) { for(sp[k]=j=0;j<n;j++) sp[k] += xpt[k][j] * xpt[kopt][j] ; temp = pq[k] * sp[k] ; for(j=0;j<n;j++) gopt[j] += fval[k]*bmat[j][k] + temp*xpt[k][j] ; } //C Set the initial elements of RESCON. for(j=0;j<m;j++) { for(rescon[j]=b[j],i=0;i<n;i++) rescon[j] -= xpt[kopt][i] * amat[j][i] ; if(rescon[j]<0) rescon[j] = 0 ; if(rescon[j]>=rhobeg) rescon[j] = -rescon[j] ; } free(step,w,auxtemp) ; return idz ; } /* -------------------------------------------------------------------------- */ static int cond=-7 ; int lincoacond() { return cond ; } double lincob(int n,int npt,int m,double **amat,double *b,double *x, double rhobeg,double rhoend,int maxfun,int ndim,int maxflag, double (*func)(double *), double (*gunc)(double *,double *),double *aux,int naux) { //C The arguments N, NPT, M, X, RHOBEG, RHOEND, IPRINT and MAXFUN are //C identical to the corresponding arguments in SUBROUTINE LINCOA. //C AMAT is a matrix whose columns are the constraint gradients, scaled //C so that they have unit length. //C B contains on entry the right hand sides of the constraints, scaled //C as above, but later B is modified for variables relative to XBASE. //C XBASE holds a shift of origin that should reduce the contributions //C from rounding errors to values of the model and Lagrange functions. //C XPT contains the interpolation point coordinates relative to XBASE. //C FVAL holds the values of F at the interpolation points. //C XSAV holds the best feasible vector of variables so far, without any //C shift of origin. //C XOPT is set to XSAV-XBASE, which is the displacement from XBASE of //C the feasible vector of variables that provides the least calculated //C F so far, this vector being the current trust region centre. //C GOPT holds the gradient of the quadratic model at XSAV = XBASE+XOPT. //C HQ holds the explicit second derivatives of the quadratic model. //C PQ contains the parameters of the implicit second derivatives of the //C quadratic model. //C BMAT holds the last N columns of the big inverse matrix H. //C ZMAT holds the factorization of the leading NPT by NPT submatrix //C of H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, //C where the elements of DZ are plus or minus one, as specified by IDZ. //C NDIM is the first dimension of BMAT and has the value NPT+N. //C STEP is employed for trial steps from XOPT. It is also used for working //C space when XBASE is shifted and in PRELIM. //C SP is reserved for the scalar products XOPT^T XPT(K,.), K=1,2,...,NPT, //C followed by STEP^T XPT(K,.), K=1,2,...,NPT. //C XNEW is the displacement from XBASE of the vector of variables for //C the current calculation of F, except that SUBROUTINE TRSTEP uses it //C for working space. //C IACT is an integer array for the indices of the active constraints. //C RESCON holds useful information about the constraint residuals. Every //C nonnegative RESCON(J) is the residual of the J-th constraint at the //C current trust region centre. Otherwise, if RESCON(J) is negative, the //C J-th constraint holds as a strict inequality at the trust region //C centre, its residual being at least |RESCON(J)|; further, the value //C of |RESCON(J)| is at least the current trust region radius DELTA. //C QFAC is the orthogonal part of the QR factorization of the matrix of //C active constraint gradients, these gradients being ordered in //C accordance with IACT. When NACT is less than N, columns are added //C to QFAC to complete an N by N orthogonal matrix, which is important //C for keeping calculated steps sufficiently close to the boundaries //C of the active constraints. //C RFAC is the upper triangular part of this QR factorization, beginning //C with the first diagonal element, followed by the two elements in the //C upper triangular part of the second column and so on. //C PQW is used for working space, mainly for storing second derivative //C coefficients of quadratic functions. Its length is NPT+N. //C The array W is also used for working space. The required number of //C elements, namely MAX[M+3*N,2*M+N,2*NPT], is set in LINCOA. int i,j,k,np=n+1,nptm=npt-n-1,ncall,nf,nvala,nvalb,knew,iflag,iter,loopflag ; int itest,kopt,ksave,ifeas,nact,idz=-1,ih,wlen=max(m+3*n,max(2*m+n,2*npt)) ; double delta,fopt,rho,fsave,ratio,xoptsq,sum,sumz,qoptsq,temp,delsav,snorm ; double vquad,xdiff,del,vqalt,diff,f,ssq,dffalt,distsq ; double *xbase=vector(n),**xpt=matrix(npt,n),*fval=vector(npt) ; double *xsav=vector(n+naux),*xopt=vector(n),*gopt=vector(n),*auxtemp ; double *hq=vector((n*(n+1))/2),*pq=vector(npt),**bmat=matrix(n,ndim) ; double **zmat=matrix(nptm,npt),*step=vector(n),*sp=vector(2*npt) ; double *xnew=vector(n),*rescon=vector(m),**qfac=matrix(n,n) ; double *rfac=vector((n*(n+1))/2),*pqw=vector(npt),*w=vector(wlen) ; int *iact=ivector(n) ; if(func) auxtemp = 0 ; else auxtemp = vector(naux) ; //C Set the elements of XBASE, XPT, FVAL, XSAV, XOPT, GOPT, HQ, PQ, BMAT, //C ZMAT and SP for the first iteration. An important feature is that, //C if the interpolation point XPT(K,.) is not feasible, where K is any //C integer from [1,NPT], then a change is made to XPT(K,.) if necessary //C so that the constraint violation is at least 0.2*RHOBEG. Also KOPT //C is set so that XPT(KOPT,.) is the initial trust region centre. idz = prelim(n,npt,m,amat,b,x,rhobeg,xbase,xpt,fval,xsav,gopt,kopt,pq,bmat, zmat,ndim,sp,rescon,maxflag,func,gunc,aux,naux) ; for(j=0;j<n;j++) xopt[j] = xpt[kopt][j] ; //C Begin the iterative procedure. fopt = fval[kopt] ; delta = rho = rhobeg ; itest = 3 ; ifeas = nact = nvala = nvalb = loopflag = 0 ; cond = knew = -1 ; nf = npt ; for(iter=0;;iter++) // begin main loop { if(iter>=n*maxfun) throw up("max iters %d*%d exceeded by lincoa",n,maxfun) ; fsave = fopt ; //C Shift XBASE if XOPT may be too far from XBASE. First make the changes //C to BMAT that do not depend on ZMAT. for(xoptsq=i=0;i<n;i++) xoptsq += xopt[i] * xopt[i] ; if(xoptsq>=1e4*delta*delta) { shiftxbase(n,npt,xpt,xopt,step,bmat,zmat,pq,hq,idz) ; for(k=0;k<npt;k++) sp[k] = 0 ; //C Update the right hand sides of the constraints. for(j=0;j<m;j++) { for(temp=i=0;i<n;i++) temp += amat[j][i] * xopt[i] ; b[j] -= temp ; } for(j=0;j<n;j++) { xbase[j] += xopt[j] ; xopt[j] = xpt[kopt][j] = 0 ; } } //C In the case KNEW=0, generate the next trust region step by calling //C TRSTEP, where SNORM is the current trust region radius initially. //C The final value of SNORM is the length of the calculated step, //C except that SNORM is zero on return if the projected gradient is //C unsuitable for starting the conjugate gradient iterations. delsav = delta ; ksave = knew ; iflag = 0 ; if(knew==-1) { snorm = delta ; for(i=0;i<n;i++) xnew[i] = gopt[i] ; ncall = trstep(n,npt,m,amat,b,xpt,hq,pq,nact,iact,rescon,qfac,rfac, snorm,step,xnew) ; //C A trust region step is applied whenever its length, namely SNORM, is at //C least HALF*DELTA. It is also applied if its length is at least 0.1999 //C times DELTA and if a line search of TRSTEP has caused a change to the //C active set. Otherwise there is a branch below to label 530 or 560. // [coded here by setting iflag according to the branch chosen - cjc] if(ncall>1) temp = 0.1999 * delta ; else temp = delta / 2 ; if(snorm<=temp) { if(delta<=2.8*rho) delta = rho ; else delta /= 2 ; nvala += 1 ; nvalb += 1 ; temp = snorm / rho ; if(delsav>rho) temp = 1 ; if(temp>=0.5) nvala = 0 ; if(temp>=0.1) nvalb = 0 ; if(delsav>rho||(nvala<5&&nvalb<3)) iflag = 1 ; else if(snorm>0) iflag = 2 ; else iflag = 3 ; } else nvala = nvalb = 0 ; } else { //C Alternatively, KNEW is positive. Then the model step is calculated //C within a trust region of radius DEL, after setting the gradient at //C XBASE and the second derivative parameters of the KNEW-th Lagrange //C function in W(1) to W(N) and in PQW(1) to PQW(NPT), respectively. if(0.1*delta>rho) del = 0.1*delta ; else del = rho ; for(i=0;i<n;i++) w[i] = bmat[i][knew] ; genpqw(n,npt,pqw,zmat,knew,idz) ; ifeas = qmstep(n,npt,m,amat,xpt,xopt,nact,iact,rescon,qfac,kopt,knew,del, step,w,pqw) ; if(ifeas<0) { ifeas = 0 ; cond = 4 ; break ; } } //C Set VQUAD to the change to the quadratic model when the move STEP is //C made from XOPT. If STEP is a trust region step, then VQUAD should be //C negative. If it is nonnegative due to rounding errors in this case, //C there is a branch to label 530 to try to improve the model. if(iflag==0) { for(vquad=ih=j=0;j<n;j++) for(vquad+=step[j]*gopt[j],i=0;i<=j;i++,ih++) { temp = step[i] * step[j] ; if(i==j) temp /= 2 ; vquad += temp * hq[ih] ; } for(k=0;k<npt;k++) { for(temp=j=0;j<n;j++) temp += xpt[k][j] * step[j] ; sp[npt+k] = temp ; vquad += 0.5 * pq[k] * temp * temp ; } if(ksave==-1&&vquad>=0) iflag = 1 ; } //C Calculate the next value of the objective function. The difference //C between the actual new value of F and the value predicted by the //C model is recorded in DIFF. if(iflag==0) { nf += 1 ; if(nf>maxfun) { cond = 2 ; break ; } for(xdiff=i=0;i<n;i++) { xnew[i] = xopt[i] + step[i] ; x[i] = xbase[i] + xnew[i] ; temp = x[i] - xsav[i] ; xdiff += temp * temp ; } if(ksave==-2) xdiff = rho ; else xdiff = sqrt(xdiff) ; if(xdiff<=0.1*rho||xdiff>=2*delta) { ifeas = 0 ; cond = 1 ; break ; } if(ksave<0) ifeas = 1 ; for(i=0;i<n;i++) if(abnormal(x[i])) throw up("2: %.1e\n",x[i]) ; if(func) f = func(x) ; else { for(i=0;i<naux;i++) auxtemp[i] = aux[i] ; f = gunc(x,naux?auxtemp:0) ; } if(isnan(f)||isinf(f)) throw up("function returns %.1e",f) ; if(maxflag) f = -f ; if(ksave==-2) { cond = 0 ; break ; } diff = f - fopt - vquad ; //C If X is feasible, then set DFFALT to the difference between the new //C value of F and the value predicted by the alternative model. if(ifeas==1&&itest<3) { for(k=0;k<npt;k++) { pqw[k] = 0 ; w[k] = fval[k] - fval[kopt] ; } for(j=0;j<nptm;j++) { for(sum=k=0;k<npt;k++) sum += w[k] * zmat[j][k] ; if(j<idz) sum = -sum ; for(k=0;k<npt;k++) pqw[k] += sum * zmat[j][k] ; } for(vqalt=k=0;k<npt;k++) { for(sum=j=0;j<n;j++) sum += bmat[j][k] * step[j] ; vqalt += sum * w[k] ; vqalt += pqw[k] * sp[npt+k] * (sp[k]+sp[npt+k]/2) ; } dffalt = f - fopt - vqalt ; } else if(itest==3) { dffalt = diff ; itest = 0 ; } //C Pick the next value of DELTA after a trust region step. if(ksave==-1) { ratio = (f-fopt) / vquad ; if(ratio<0.1) delta /= 2 ; else if(ratio<0.7) { if(delta/2>snorm) delta /= 2 ; else delta = snorm ; } else { temp = sqrt(2)*delta ; if(delta/2>2*snorm) delta /= 2 ; else delta = 2*snorm ; if(delta>temp) delta = temp ; } if(delta<=1.4*rho) delta = rho ; } //C Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point //C can be moved. If STEP is a trust region step, then KNEW is zero at //C present, but a positive value is picked by subroutine UPDATE. lincupdate(n,npt,xpt,bmat,zmat,idz,ndim,sp,step,kopt,knew) ; if(knew==-1) { cond = 3 ; break ; } //C If ITEST is increased to 3, then the next quadratic model is the //C one whose second derivative matrix is least subject to the new //C interpolation conditions. Otherwise the new model is constructed //C by the symmetric Broyden method in the usual way. if(ifeas==1) { if(fabs(dffalt)>=0.1*fabs(diff)) itest = 0 ; else itest += 1 ; } //C Update the second derivatives of the model by the symmetric Broyden //C method, using PQW for the second derivative parameters of the new //C KNEW-th Lagrange function. The contribution from the old parameter //C PQ(KNEW) is included in the second derivative matrix HQ. W is used //C later for the gradient of the new KNEW-th Lagrange function. if(itest<3) { updatehq(n,npt,hq,pq[knew],xpt,knew) ; pq[knew] = 0 ; genpqw(n,npt,pqw,zmat,knew,idz) ; for(k=0;k<npt;k++) pq[k] += diff * pqw[k] ; } //C Include the new interpolation point with the corresponding updates of //C SP. Also make the changes of the symmetric Broyden method to GOPT at //C the old XOPT if ITEST is less than 3. fval[knew] = f ; sp[knew] = sp[kopt] + sp[npt+kopt] ; for(ssq=i=0;i<n;i++) { xpt[knew][i] = xnew[i] ; ssq += step[i]*step[i] ; } sp[npt+knew] = sp[npt+kopt] + ssq ; if(itest<3) { for(i=0;i<n;i++) w[i] = bmat[i][knew] ; for(k=0;k<npt;k++) for(temp=pqw[k]*sp[k],i=0;i<n;i++) w[i] += temp * xpt[k][i] ; for(i=0;i<n;i++) gopt[i] += diff * w[i] ; } //C Update FOPT, XSAV, XOPT, KOPT, RESCON and SP if the new F is the //C least calculated value so far with a feasible vector of variables. if(f<fopt&&ifeas==1) { fopt = f ; for(j=0;j<n;j++) { xsav[j] = x[j] ; xopt[j] = xnew[j] ; } for(j=0;j<naux;j++) xsav[n+j] = auxtemp[j] ; kopt = knew ; snorm = sqrt(ssq) ; for(j=0;j<m;j++) if(rescon[j]>=delta+snorm) rescon[j] = snorm - rescon[j] ; else { rescon[j] += snorm ; if(rescon[j]+delta>0) { for(temp=b[j],i=0;i<n;i++) temp -= xopt[i] * amat[j][i] ; if(temp<0) temp = 0 ; if(temp>=delta) rescon[j] = -temp ; else rescon[j] = temp ; } } for(k=0;k<npt;k++) sp[k] += sp[npt+k] ; //C Also revise GOPT when symmetric Broyden updating is applied. if(itest<3) { updategopt(n,gopt,hq,step) ; for(k=0;k<npt;k++) for(temp=pq[k]*sp[npt+k],i=0;i<n;i++) gopt[i] += temp * xpt[k][i] ; } } // end if(f<fopt&&ifeas==1) if(itest==3) { //C Replace the current model by the least Frobenius norm interpolant if //C this interpolant gives substantial reductions in the predictions //C of values of F at feasible points. for(k=0;k<npt;k++) { pq[k] = 0 ; w[k] = fval[k] - fval[kopt] ; } for(j=0;j<nptm;j++) { for(sum=k=0;k<npt;k++) sum += w[k] * zmat[j][k] ; if(j<idz) sum = -sum ; for(k=0;k<npt;k++) pq[k] += sum * zmat[j][k] ; } for(j=0;j<n;j++) for(gopt[j]=i=0;i<npt;i++) gopt[j] += w[i] * bmat[j][i] ; for(k=0;k<npt;k++) for(temp=pq[k]*sp[k],i=0;i<n;i++) gopt[i] += temp * xpt[k][i] ; for(i=0;i<(n*(n+1))/2;i++) hq[i] = 0 ; } //C If a trust region step has provided a sufficient decrease in F, then //C branch for another trust region calculation. Every iteration that //C takes a model step is followed by an attempt to take a trust region //C step. knew = -1 ; if(ksave>=0||ratio>=0.1) { loopflag = 0 ; continue ; } } // end if(iflag==0) //C Alternatively, find out if the interpolation points are close enough //C to the best point so far (old label 530) if(iflag<2) { distsq = max(delta*delta,4*rho*rho) ; for(k=0;k<npt;k++) { for(sum=j=0;j<n;j++) { temp = xpt[k][j] - xopt[j] ; sum += temp * temp ; } if(sum>distsq) { knew = k ; distsq = sum ; } } //C If KNEW is positive, then branch back for the next iteration, which //C will generate a "model step". Otherwise, if the current iteration //C has reduced F, or if DELTA was above its lower bound when the last //C trust region step was calculated, then try a "trust region" step //C instead. if(knew>=0||fopt<fsave) { loopflag = 0 ; continue ; } else if(loopflag==n) { cond = 4 ; break ; } else if(delsav>rho) { loopflag += 1 ; continue ; } } // end if(iflag<2) loopflag = 0 ; //C The calculations with the current value of RHO are complete. //C Pick the next value of RHO. if(rho<=rhoend) { cond = 0 ; break ; } delta = rho /2 ; if(rho>250*rhoend) rho *= 0.1 ; else if(rho<16*rhoend) rho = rhoend ; else rho = sqrt(rho*rhoend) ; delta = max(delta,rho) ; knew = -1 ; nvala = nvalb = 0 ; } // end for(iter=0;;iter++) //C Return from the calculation, after branching to label 220 for another //C Newton-Raphson step if it has not been tried before. if(iflag==2) { for(nf++,i=0;i<n;i++) x[i] = xbase[i] + xopt[i] + step[i] ; for(j=0;j<n;j++) if(abnormal(x[j])) throw up("3: %.1e\n",x[j]) ; if(func) f = func(x) ; else { for(i=0;i<naux;i++) auxtemp[i] = aux[i] ; f = gunc(x,naux?auxtemp:0) ; } if(isnan(f)||isinf(f)) throw up("function returns %.1e",f) ; if(maxflag) f = -f ; } if(fopt<=f||(ifeas==0&&iflag!=2)) { for(f=fopt,i=0;i<n;i++) x[i] = xsav[i] ; for(i=0;i<naux;i++) aux[i] = xsav[n+i] ; } else for(i=0;i<naux;i++) aux[i] = auxtemp[i] ; free(xbase,fval,xsav,xopt,gopt,hq) ; free(pq,step,sp,xnew,rescon,rfac) ; free(pqw,w,iact,auxtemp) ; freematrix(bmat,zmat,qfac,xpt) ; if(maxflag) return -f ; else return f ; } /* -------------------------------------------------------------------------- */ int quadprog(double **Q,double *A,double *x,int n0,int n1, double **le,double *lerhs,int nle, double **ge,double *gerhs,int nge, double **eq,double *eqrhs,int neq) ; double quadprogresid() ; void quadprogprefs(int debug,int rescale,int perturb) ; double lincoastub(double (*func)(double *),double *x,int n, double (*gunc)(double *,double *),double *aux,int naux, double **a,double *b,int m, double rhobeg,double rhoend,int maxfun,int npt,int maxflag) { //C This subroutine seeks the least value of a function of many variables, //C subject to general linear inequality constraints, by a trust region //C method that forms quadratic models by interpolation. Usually there //C is much freedom in each new model after satisfying the interpolation //C conditions, which is taken up by minimizing the Frobenius norm of //C the change to the second derivative matrix of the model. One new //C function value is calculated on each iteration, usually at a point //C where the current model predicts a reduction in the least value so //C far of the objective function subject to the linear constraints. //C Alternatively, a new vector of variables may be chosen to replace //C an interpolation point that may be too far away for reliability, and //C then the new point does not have to satisfy the linear constraints. //C The arguments of the subroutine are as follows. //C //C N must be set to the number of variables and must be at least two. //C NPT must be set to the number of interpolation conditions, which is //C required to be in the interval [N+2,(N+1)(N+2)/2]. Typical choices //C of the author are NPT=N+6 and NPT=2*N+1. Larger values tend to be //C highly inefficent when the number of variables is substantial, due //C to the amount of work and extra difficulty of adjusting more points. //C M must be set to the number of linear inequality constraints. //C A is a matrix whose columns are the constraint gradients, which are //C required to be nonzero. //C IA is the first dimension of the array A, which must be at least N. //C B is the vector of right hand sides of the constraints, the J-th //C constraint being that the scalar product of A(.,J) with X(.) is at //C most B(J). The initial vector X(.) is made feasible by increasing //C the value of B(J) if necessary. //C X is the vector of variables. Initial values of X(1),X(2),...,X(N) //C must be supplied. If they do not satisfy the constraints, then B //C is increased as mentioned above. X contains on return the variables //C that have given the least calculated F subject to the constraints. //C RHOBEG and RHOEND must be set to the initial and final values of a //C trust region radius, so both must be positive with RHOEND<=RHOBEG. //C Typically, RHOBEG should be about one tenth of the greatest expected //C change to a variable, and RHOEND should indicate the accuracy that //C is required in the final values of the variables. //C The value of IPRINT should be set to 0, 1, 2 or 3, which controls the //C amount of printing. Specifically, there is no output if IPRINT=0 and //C there is output only at the return if IPRINT=1. Otherwise, the best //C feasible vector of variables so far and the corresponding value of //C the objective function are printed whenever RHO is reduced, where //C RHO is the current lower bound on the trust region radius. Further, //C each new value of F with its variables are output if IPRINT=3. //C MAXFUN must be set to an upper bound on the number of calls of CALFUN, //C its value being at least NPT+1. //C W is an array used for working space. Its length must be at least //C M*(2+N) + NPT*(4+N+NPT) + N*(9+3*N) + MAX [ M+3*N, 2*M+N, 2*NPT ]. //C On return, W(1) is set to the final value of F, and W(2) is set to //C the total number of function evaluations plus 0.5. //C //C SUBROUTINE CALFUN (N,X,F) has to be provided by the user. It must set //C F to the value of the objective function for the variables X(1), //C X(2),...,X(N). The value of the argument F is positive when CALFUN //C is called if and only if the current X satisfies the constraints //C to working accuracy. //C //C Check that N, NPT and MAXFUN are acceptable. double **amat,*bnorm,temp,sum,res,**Q,*A ; int i,j,iflag,u ; if(n<2) { cond = -1 ; return 0 ; } if(npt<n+2||npt>((n+2)*(n+1))/2) { cond = -2 ; return 0 ; } if(maxfun<npt) { cond = -3 ; return 0 ; } if(rhobeg<=0) throw up("lincoa call with rhobeg=%.1e",rhobeg) ; visit += 1 ; amat = matrix(m,n) ; bnorm = vector(m) ; //C Normalize the constraints, and copy the resultant constraint matrix //C and right hand sides into working space, after increasing the right //C hand sides if necessary so that the starting point is feasible. for(iflag=i=0;i<m;i++) { for(sum=temp=j=0;j<n;j++) { sum += a[i][j]*x[j] ; temp += a[i][j]*a[i][j] ; } if(temp==0) { cond = -4 ; freematrix(amat) ; free(bnorm) ; return 0 ; } if(sum>b[i]) iflag = 1 ; temp = sqrt(temp) ; bnorm[i] = b[i] / temp ; for(j=0;j<n;j++) amat[i][j] = a[i][j] / temp ; } if(iflag) { Q = matrix(n,n) ; A = vector(n) ; for(i=0;i<n;i++) { Q[i][i] = -1 ; A[i] = x[i] ; } u = quadprog(Q,A,x,0,n, a,b,m, 0,0,0, 0,0,0) ; freematrix(Q) ; free(A) ; if(u<0) printf("*** lincoa %d: return %d from quadprog (resid=%.1e)\n", visit,u,quadprogresid()) ; } res = lincob(n,npt,m,amat,bnorm,x,rhobeg,rhoend,maxfun,npt+n,maxflag, func,gunc,aux,naux) ; freematrix(amat) ; free(bnorm) ; return res ; } /* -------------------------------------------------------------------------- */ double lincoa(double (*func)(double *),double *x,int n,double **a,double *b, int m,double rhobeg,double rhoend,int maxfun) { return lincoastub(func,x,n,0,0,0,a,b,m,rhobeg,rhoend,maxfun,2*n+1,0) ; } double lincoa(double (*gunc)(double *,double *),double *x,int n, double *aux,int naux,double **a,double *b,int m, double rhobeg,double rhoend,int maxfun) { if(aux==0) naux = 0 ; return lincoastub(0,x,n,gunc,aux,naux,a,b,m,rhobeg,rhoend,maxfun,2*n+1,0) ; } double lincoamax(double (*func)(double *),double *x,int n,double **a,double *b, int m,double rhobeg,double rhoend,int maxfun) { return lincoastub(func,x,n,0,0,0,a,b,m,rhobeg,rhoend,maxfun,2*n+1,1) ; } double lincoamax(double (*gunc)(double *,double *),double *x,int n, double *aux,int naux,double **a,double *b,int m, double rhobeg,double rhoend,int maxfun) { if(aux==0) naux = 0 ; return lincoastub(0,x,n,gunc,aux,naux,a,b,m,rhobeg,rhoend,maxfun,2*n+1,1) ; }
The test program lincoatest.c below calls lincoa to minimise Powell's test function in 12 dimensions with interpolations differing in the number of points they use (2n+1 are recommended). The current output is this:
optimum = 2.761312711773343e+00 from 120 calls; cond=0 optimum = 2.761312522306842e+00 from 156 calls; cond=0 optimum = 2.751575558247701e+00 from 177 calls; cond=0 optimum = 2.761312521923594e+00 from 218 calls; cond=0 optimum = 2.761312522794534e+00 from 172 calls; cond=0 optimum = 2.761312539377540e+00 from 188 calls; cond=0
For comparison the fable translation of Powell’s original Fortran gives the following output:
Output from LINCOA with NPT = 15 and RHOEND = 1.0000D-06
At the return from LINCOA Number of function values = 115
Least value of F = 2.761312522760941D+00 The corresponding X is:
-1.233791D+00 -1.231348D+00 -8.870616D-01 1.129888D+00 -1.128285D+00
1.084034D+00 8.431307D-01 7.535403D-01 -8.461757D-01 -6.585848D-01
1.556805D+00 5.853412D-01
Output from LINCOA with NPT = 20 and RHOEND = 1.0000D-06
At the return from LINCOA Number of function values = 156
Least value of F = 2.761312522306842D+00 The corresponding X is:
-1.233791D+00 -1.231348D+00 -8.870616D-01 1.129888D+00 -1.128285D+00
1.084034D+00 8.431308D-01 7.535398D-01 -8.461758D-01 -6.585848D-01
1.556805D+00 5.853412D-01
Output from LINCOA with NPT = 25 and RHOEND = 1.0000D-06
At the return from LINCOA Number of function values = 177
Least value of F = 2.751575558247701D+00 The corresponding X is:
-1.210293D+00 -1.212501D+00 -9.512774D-01 1.157576D+00 -1.159364D+00
1.048821D+00 7.788128D-01 1.062499D+00 -7.734505D-01 -7.068258D-01
1.408297D+00 6.554270D-01
Output from LINCOA with NPT = 30 and RHOEND = 1.0000D-06
At the return from LINCOA Number of function values = 218
Least value of F = 2.761312521923595D+00 The corresponding X is:
-1.233791D+00 -1.231348D+00 -8.870616D-01 1.129888D+00 -1.128285D+00
1.084034D+00 8.431309D-01 7.535394D-01 -8.461759D-01 -6.585848D-01
1.556805D+00 5.853412D-01
Output from LINCOA with NPT = 35 and RHOEND = 1.0000D-06
At the return from LINCOA Number of function values = 172
Least value of F = 2.761312522794534D+00 The corresponding X is:
-8.431307D-01 7.535403D-01 8.461757D-01 1.233791D+00 -1.231348D+00
8.870616D-01 -1.129888D+00 -1.128285D+00 -1.084034D+00 6.585848D-01
1.556805D+00 -5.853412D-01
Output from LINCOA with NPT = 40 and RHOEND = 1.0000D-06
At the return from LINCOA Number of function values = 203
Least value of F = 2.761312522731237D+00 The corresponding X is:
-8.431307D-01 7.535403D-01 8.461757D-01 1.233791D+00 -1.231348D+00
8.870616D-01 -1.129888D+00 -1.128285D+00 -1.084034D+00 6.585848D-01
1.556805D+00 -5.853412D-01
#include "memory.h"
#include <math.h>
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
double lincoastub(double (*func)(double *),double *x,int n,
double (*gunc)(double *,double *),double *aux,int naux,
double **a,double *b,int m,
double rhobeg,double rhoend,int maxfun,int npt,int maxflag) ;
int lincoacond() ;
#define pi 3.141592653589793
static double flim ;
static int ncall ;
double func(double *x)
{ double v12,v13,v14,v23,v24,v34,del1,del2,del3,del4,f ;
ncall += 1 ;
v12 = x[0] * x[4] - x[3] * x[1] ;
v13 = x[0] * x[7] - x[6] * x[1] ;
v14 = x[0] * x[10] - x[9] * x[1] ;
v23 = x[3] * x[7] - x[6] * x[4] ;
v24 = x[3] * x[10] - x[9] * x[4] ;
v34 = x[6] * x[10] - x[9] * x[7] ;
del1 = v23 * x[11] - v24 * x[8] + v34 * x[5] ;
del2 = -v34 * x[2] - v13 * x[11] + v14 * x[8] ;
del3 = -v14 * x[5] + v24 * x[2] + v12 * x[11] ;
del4 = -v12 * x[8] + v13 * x[5] - v23 * x[2] ;
if(del1>0&&del2>0&&del3>0&&del4>0)
{ f = del1 + del2 + del3 + del4 ;
f = (f*f*f) / (del1*del2*del3*del4) ;
return min(f/6,flim) ;
}
else return flim ;
}
double gunc(double *x,double *aux)
{ int t ;
for(t=0;t<12;t++) if(aux[t]!=t) throw up("aux[%d]=%.1f",t,aux[t]) ;
double q = func(x) ;
for(t=0;t<12;t++) aux[t] = x[11-t] ;
return q ;
}
//C Calculate the tetrahedron of least volume that encloses the points
//C (XP(J),YP(J),ZP(J)), J=1,2,...,NP. Our method requires the origin
//C to be strictly inside the convex hull of these points. There are
//C twelve variables that define the four faces of each tetrahedron
//C that is considered. Each face has the form ALPHA*X+BETA*Y+GAMMA*Z=1,
//C the variables X(3K-2), X(3K-1) and X(3K) being the values of ALPHA,
//C BETA and GAMMA for the K-th face, K=1,2,3,4. Let the set T contain
//C all points in three dimensions that can be reached from the origin
//C without crossing a face. Because the volume of T may be infinite,
//C the objective function is the smaller of FMAX and the volume of T,
//C where FMAX is set to an upper bound on the final volume initially.
//C There are 4*NP linear constraints on the variables, namely that each
//C of the given points (XP(J),YP(J),ZP(J)) shall be in T. Let XS = min
//C XP(J), YS = min YP(J), ZS = min ZP(J) and SS = max XP(J)+YP(J)+ZP(J),
//C where J runs from 1 to NP. The initial values of the variables are
//C X(1)=1/XS, X(5)=1/YS, X(9)=1/ZS, X(2)=X(3)=X(4)=X(6)=X(7) =X(8)=0
//C and X(10)=X(11)=X(12)=1/SS, which satisfy the linear constraints,
//C and which provide the bound FMAX=(SS-XS-YS-ZS)**3/6. Other details
//C of the test calculation are given below, including the choice of
//C the data points (XP(J),YP(J),ZP(J)), J=1,2,...,NP. The smaller final
//C value of the objective function in the case NPT=35 shows that the
//C problem has local minima.
//C
int main(int argc,char const* argv[])
{ int np=50,ia=12,n=12,m=200,i,j,k ;
double xp[50],yp[50],zp[50],theta,sumx,sumy,sumz,xs,ys,zs,ss,res,res2 ;
double b[200],xval[12],x[12],aux[12],**a=matrix(m,n) ;
for(sumx=sumy=sumz=j=0;j<np;j++)
{ theta = j * pi / (np-1) ;
sumx += xp[j] = cos(theta) * cos(2*theta) ;
sumy += yp[j] = sin(theta) * cos(2*theta) ;
sumz += zp[j] = sin(2*theta) ;
}
sumx /= np ;
sumy /= np ;
sumz /= np ;
for(j=0;j<np;j++) { xp[j] -= sumx ; yp[j] -= sumy ; zp[j] -= sumz ; }
for(k=0;k<m;k++) b[k] = 1 ;
for(i=0;i<n;i++) for(k=0;k<m;k++) a[k][i] = 0 ;
for(j=0;j<np;j++) for(i=0;i<4;i++)
{ a[4*j+i][3*i ] = xp[j] ;
a[4*j+i][3*i+1] = yp[j] ;
a[4*j+i][3*i+2] = zp[j] ;
}
//C Set the initial vector of variables. The JCASE=1,6 loop gives six
//C different choices of NPT when LINCOA is called.
for(xs=ys=zs=ss=j=0;j<np;j++)
{ xs = min(xs,xp[j]) ;
ys = min(ys,yp[j]) ;
zs = min(zs,zp[j]) ;
ss = max(ss,xp[j]+yp[j]+zp[j]) ;
}
flim = ss - xs - ys - zs ;
flim *= flim * flim / 6 ;
for(i=0;i<12;i++) xval[i] = 0 ;
xval[0] = 1 / xs ;
xval[4] = 1 / ys ;
xval[8] = 1 / zs ;
xval[9] = xval[10] = xval[11] = 1 / ss ;
for(j=0;j<6;j++)
{ for(i=0;i<12;i++) x[i] = xval[i] ;
ncall = 0 ;
res = lincoastub(&func,x,n,0,0,0,a,b,m,1,1e-6,10000,5*j+15,0) ;
printf("optimum = %.15e from %d calls; cond=%d\n",res,ncall,lincoacond()) ;
for(i=0;i<12;i++) { x[i] = xval[i] ; aux[i] = i ; }
res2 = lincoastub(0,x,n,&gunc,aux,12,a,b,m,1,1e-6,10000,5*j+15,0) ;
if(res2!=res) throw up("discrepant returns %.15f/%.15f lincoa",res,res2) ;
for(i=0;i<12;i++) if(x[i]!=aux[11-i])
throw up("x[%d]=%.15f; aux[%d]=%.15f",i,x[i],11-i,aux[11-i]) ;
}
freematrix(a) ;
}
In the course of converting lincoa I had to transpose some of the C matrices. I wrote a little program trans for this purpose. I include the source code purely for my own benefit.
#include "memory.h" #include <string.h> #include <ctype.h> int main(int argc,char **argv) { int fileno,n,i,j,len,dep,nchange ; FILE *ifl,*ofl ; char mat[100],buf[50000],arg1[100] ; strcpy(mat,argv[1]) ; len = strlen(mat) ; mat[len] = '[' ; mat[len+1] = 0 ; for(nchange=0,fileno=2;fileno<argc;fileno++) { ifl = fopenread(argv[fileno]) ; n = fread(buf,1,50000,ifl) ; fclose(ifl) ; ofl = fopenwrite(argv[fileno]) ; for(i=0;i<n;i=j+1) { for(j=i;j<n&&((j>0&&isalnum(buf[j-1]))||strncmp(buf+j,mat,len+1));j++) ; fwrite(buf+i,1,j-i,ofl) ; if(j==n) break ; fwrite(mat,1,len+1,ofl) ; for(i=j+len+1,dep=0,j=i;j<n&&(buf[j]!=']'||buf[j+1]!='['||dep>0);j++) if(buf[j]=='[') dep += 1 ; else if(buf[j]==']') dep -= 1 ; if(j==n) break ; strncpy(arg1,buf+i,j-i) ; arg1[j-i] = 0 ; for(i=j+2,dep=0,j=i;j<n&&(buf[j]!=']'||dep>0);j++) if(buf[j]=='[') dep += 1 ; else if(buf[j]==']') dep -= 1 ; fwrite(buf+i,1,j-i,ofl) ; fwrite("][",1,2,ofl) ; fwrite(arg1,1,strlen(arg1),ofl) ; fwrite("]",1,1,ofl) ; nchange += 1 ; } fclose(ofl) ; } printf("%d changes\n",nchange) ; }
This link generates a perl script which you can download and execute to compile and run lincoa. It works in Firefox only (and needs popups enabled for this site).
lincoa.c mdplib.c quadprog.c lincoatest.c utils: memory.h --- g++ -o lincoatest -g -O2 (*.c) lincoatest
doc :