#include <minimizer.h>
Inheritance diagram for matvec::Minimizer:
Public Methods | |
| Minimizer (void) | |
| virtual | ~Minimizer (void) |
| virtual double | minfun (const Vector< double > &x, const int n)=0 |
| void | min_prtlevel (const int pl) |
| double | praxis (Vector< double > &x, const int n, int &maxfun, const double tol=1.0e-16, const double epsilon=1.0e-8, const int pl=2) |
Protected Attributes | |
| int | minfun_indx |
Private Methods | |
| void | min_first_dir (const int j, const int nits, double *d2, double *x1, double f1, const int fk, Vector< double > &x, const int n, const Matrix< double > &v, Vector< double > &w, Matrix< double > &q01) |
| void | quadratic (Vector< double > &x, const int n, const Matrix< double > &v, Vector< double > &w, Matrix< double > &q01) |
Private Attributes | |
| int | prtlevel |
|
|
Definition at line 46 of file minimizer.h. References minfun_indx, and prtlevel.
00046 {prtlevel = 2; minfun_indx = 0;}
|
|
|
Definition at line 47 of file minimizer.h.
00047 {;}
|
|
||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 406 of file min_praxis.cpp. References matvec::dmin, matvec::flin(), matvec::fx, matvec::h, matvec::ldt, matvec::m2, matvec::m4, matvec::macheps, minfun(), matvec::nf, matvec::nl, matvec::SmaLL, and matvec::t. Referenced by praxis(), and quadratic().
00410 {
00411 int k, i, dz;
00412 double x2, xm, f0, f2, fm, d1, t2, s, sf1, sx1;
00413
00414 sf1 = f1; sx1 = *x1;
00415 k = 0; xm = 0.0; fm = f0 = fx; dz = *d2 < macheps;
00416 ///////////////////////////////////
00417 // find step size //
00418 //////////////////////////////////
00419 s = 0;
00420 for (i=0; i<n; i++) s += x[i]*x[i];
00421 s = std::sqrt(s);
00422 if (dz) { t2 = m4*std::sqrt(std::fabs(fx)/dmin + s*ldt) + m2*ldt; }
00423 else { t2 = m4*std::sqrt(std::fabs(fx)/(*d2) + s*ldt) + m2*ldt; }
00424 s = s*m4 + t;
00425 if (dz && t2 > s) t2 = s;
00426 if (t2 < SmaLL) t2 = SmaLL;
00427 if (t2 > 0.01*h) t2 = 0.01 * h;
00428 if (fk && f1 <= fm) {
00429 xm = *x1;
00430 fm = f1;
00431 }
00432 if (!fk || std::fabs(*x1) < t2) {
00433 *x1 = (*x1 > 0 ? t2 : -t2);
00434 flin(*x1, j,x,n,v,w,q01);
00435 nf++;
00436 f1 = minfun(w,n);
00437 }
00438 if (f1 <= fm) {
00439 xm = *x1;
00440 fm = f1;
00441 }
00442 next:
00443 if (dz) {
00444 x2 = (f0 < f1 ? -(*x1) : 2*(*x1));
00445 flin(x2, j,x,n,v,w,q01);
00446 nf++;
00447 f2 = minfun(w,n);
00448 if (f2 <= fm) {
00449 xm = x2;
00450 fm = f2;
00451 }
00452 *d2 = (x2*(f1-f0) - (*x1)*(f2-f0))/((*x1)*x2*((*x1)-x2));
00453 }
00454 d1 = (f1-f0)/(*x1) - *x1**d2; dz = 1;
00455 if (*d2 <= SmaLL) { x2 = (d1 < 0 ? h : -h); }
00456 else { x2 = - 0.5*d1/(*d2); }
00457 if (std::fabs(x2) > h) x2 = (x2 > 0 ? h : -h);
00458 tryit:
00459 flin(x2, j,x,n,v,w,q01);
00460 nf++;
00461 f2 = minfun(w,n);
00462 if ((k < nits) && (f2 > f0)) {
00463 k++;
00464 if ((f0 < f1) && (*x1*x2 > 0.0)) goto next;
00465 x2 *= 0.5;
00466 goto tryit;
00467 }
00468 nl++;
00469 if (f2 > fm) x2 = xm; else fm = f2;
00470 if (std::fabs(x2*(x2-*x1)) > SmaLL) {
00471 *d2 = (x2*(f1-f0) - *x1*(fm-f0))/(*x1*x2*(*x1-x2));
00472 }
00473 else {
00474 if (k > 0) *d2 = 0;
00475 }
00476 if (*d2 <= SmaLL) *d2 = SmaLL;
00477 *x1 = x2; fx = fm;
00478 if (sf1 < fx) {
00479 fx = sf1;
00480 *x1 = sx1;
00481 }
00482 if (j != -1) for (i=0; i<n; i++) x[i] += (*x1)*v[i][j];
00483 }
|
|
|
Definition at line 50 of file minimizer.h. References prtlevel.
00050 {prtlevel = pl;}
|
|
||||||||||||
|
Implemented in matvec::Model. Referenced by min_first_dir(), and praxis(). |
|
||||||||||||||||||||||||||||
|
Definition at line 88 of file min_praxis.cpp. References matvec::Vector< T >::begin(), matvec::Matrix< T >::begin(), matvec::dmin, matvec::fx, matvec::h, matvec::large, matvec::ldt, matvec::m2, matvec::m4, matvec::macheps, min_first_dir(), minfun(), matvec::nf, matvec::nl, matvec::Matrix< T >::print(), matvec::Vector< T >::print(), matvec::print(), prtlevel, matvec::qd, matvec::qf1, matvec::ranf(), matvec::SmaLL, matvec::sort(), matvec::svd(), matvec::t, matvec::vlarge, and matvec::vsmall. Referenced by matvec::Model::multipoint_ml_estimate(), and matvec::Model::vce_dfreml().
00090 {
00091 ////////////////////////////////////////////////////////////////////
00092 // references:
00093 // - powell, m.j.d., 1964. an efficient method for finding
00094 // the minimum of a function in several variables without
00095 // calculating derivatives, computer journal, 7, 155-162
00096 // - brent, r.p., 1973. algorithms for minimization without
00097 // derivatives, prentice hall, englewood cliffs.
00098 //
00099 // problems, suggestions or improvements are always wellcome
00100 // karl gegenfurtner 07/08/87
00101 // c - version
00102 // On Thu Oct 21 19:56:35 CDT 1993, Tianlin made it into C++ version
00103 // note that lots of globel variables have been removed.
00104 //
00105 //
00106 // usage: min = praxis(fun, x, n);
00107 //
00108 // fun the function to be minimized. fun is called from
00109 // praxis with x and n as arguments
00110 // x a double array containing the initial guesses for
00111 // the minimum, which will contain the solution on
00112 // return
00113 // n an integer specifying the number of unknown
00114 // parameters
00115 // min praxis returns the least calculated value of fun
00116 //
00117 // some additional global variables control some more aspects of
00118 // the inner workings of praxis. setting them is optional, they
00119 // are all set to some reasonable default values given below.
00120 //
00121 // prin controls the printed output from the routine.
00122 // 0 -> no output
00123 // 1 -> print only starting and final values
00124 // 2 -> detailed map of the minimization process
00125 // 3 -> print also eigenvalues and vectors of the
00126 // search directions
00127 // the default value is 1
00128 // tol is the tolerance allowed for the precision of the
00129 // solution. praxis returns if the criterion
00130 // 2 * ||x[k]-x[k-1]|| <= sqrt(macheps) * ||x[k]|| + tol
00131 // is fulfilled more than ktm times.
00132 // the default value 1.0e-16
00133 // tol is usually equal to EPISILON*EPISILON
00134 // where
00135 // ktm see just above. default is 1, and a value of 4 leads
00136 // to a very(!) cautious stopping criterion.
00137 // step is a steplength parameter and should be set equal
00138 // to the expected distance from the solution.
00139 // exceptionally small or large values of step lead to
00140 // slower convergence on the first few iterations
00141 // the default value for step is 1.0
00142 // scbd is a scaling parameter. 1.0 is the default and
00143 // indicates no scaling. if the scales for the different
00144 // parameters are very different, scbd should be set to
00145 // a value of about 10.0.
00146 // illc should be set to true (1) if the problem is known to
00147 // be ill-conditioned. the default is false (0). this
00148 // variable is automatically set, when praxis finds
00149 // the problem to be ill-conditioned during iterations.
00150 // maxfun is the maximum number of calls to fun allowed. praxis
00151 // will return after maxfun calls to fun even when the
00152 // minimum is not yet found. the default value of 0
00153 // indicates no limit on the number of calls.
00154 // this return condition is only checked every n
00155 // iterations.
00156 //
00157 //////////////////////////////////////////////////////////////////////
00158 int i,j,k,k2, W = 6 + 6;
00159
00160 int ktm = 1;
00161 int illc = 0;
00162 double scbd = 1.0;
00163 double step = 10.0; // after the first n*10 loop, step reset to 1
00164 Vector<double> y(n);
00165 Vector<double> z(n);
00166 Vector<double> d(n);
00167 // d[0] = 0.0;
00168
00169 Vector<double> w(n); // working space for flin
00170 Matrix<double> v(n,n);
00171 Matrix<double> q01(2,n);
00172
00173 qd[0] = 0.0;
00174 prtlevel = pl;
00175 macheps = epsilon;
00176 SmaLL = macheps*macheps;
00177 vsmall = SmaLL*SmaLL;
00178 large = 1.0/SmaLL;
00179 vlarge = 1.0/vsmall;
00180 m2 = std::sqrt(macheps); m4 = std::sqrt(m2);
00181
00182 nl = 0;
00183 nf = 1;
00184 double t2,sl,dn,df,sf,lds,f1,s;
00185 h = step;
00186
00187 dmin = SmaLL;
00188 t2 = SmaLL + std::fabs(tol);
00189 t = t2;
00190 double ldfac = (illc ? 0.1 : 0.01);
00191 int kt = 0;
00192 int kl;
00193
00194 fx = minfun(x, n);
00195 qf1 = fx;
00196
00197 if (h < 100.0*t) h = 100.0*t;
00198 ldt = h;
00199 for (i=0; i<n; i++) for (j=0; j<n; j++) v[i][j] = (i == j ? 1.0 : 0.0);
00200 for (i=0; i<n; i++) q01[1][i] = x[i];
00201 if (prtlevel > 1) {
00202 std::cout << "\n------------- enter function praxis -----------\n";
00203 std::cout << "-- current parameter settings\n";
00204 std::cout << "-- scaling: " << std::setw(W) << scbd << "\n";
00205 std::cout << "-- tol : " << std::setw(W) << t << "\n";
00206 std::cout << "-- maxstep: " << std::setw(W) << h << "\n";
00207 std::cout << "-- illc : " << std::setw(W) << illc << "\n";
00208 std::cout << "-- ktm : " << std::setw(W) << ktm << "\n";
00209 std::cout << "-- maxite : " << std::setw(W) << maxfun << "\n";
00210 }
00211 if (prtlevel) print(n,x);
00212
00213 mloop:
00214 sf = d[0];
00215 s = d[0] = 0.0;
00216
00217 /* minimize along first direction */
00218 min_first_dir(0,2,&d[0],&s,fx,0,x,n,v,w,q01);
00219 if (s <= 0.0) for (i=0; i < n; i++) v[i][0] = -v[i][0];
00220 if ((sf <= (0.9 * d[0])) || ((0.9 * sf) >= d[0])) {
00221 for (i=1; i<n; i++) d[i] = 0.0;
00222 }
00223 for (k=1; k<n; k++) {
00224 for (i=0; i<n; i++) y[i] = x[i];
00225 sf = fx;
00226 illc = illc || (kt > 0);
00227 next:
00228 kl = k;
00229 df = 0.0;
00230 if (illc) { /* random step to get off resolution valley */
00231 for (i=0; i<n; i++) {
00232 z[i] = (0.1*ldt + t2*pow(10.0,(double)kt))*(ranf()-0.5);
00233 s = z[i];
00234 for (j=0; j < n; j++) x[j] += s * v[j][i];
00235 }
00236 fx = minfun(x, n);
00237 nf++;
00238 }
00239 /* minimize along non-conjugate directions */
00240 for (k2=k; k2<n; k2++) {
00241 sl = fx;
00242 s = 0.0;
00243 min_first_dir(k2,2,&d[k2],&s,fx,0,x,n,v,w,q01);
00244 if (illc) {
00245 s = d[k2] * (s + z[k2]) * (s + z[k2]);
00246 }
00247 else
00248 s = sl - fx;
00249 if (df < s) {
00250 df = s;
00251 kl = k2;
00252 }
00253 }
00254
00255 if (!illc && (df < std::fabs(100.0 * macheps * fx))) {
00256 illc = 1;
00257 goto next;
00258 }
00259 if (k == 1 && prtlevel > 1) {
00260 std::cout << "\n-- New Direction:\n";
00261 d.print();
00262 }
00263 // minimize along conjugate directions
00264 for (k2=0; k2<=k-1; k2++) {
00265 s = 0.0;
00266 min_first_dir(k2,2,&d[k2],&s,fx,0,x,n,v,w,q01);
00267 }
00268 f1 = fx;
00269 fx = sf;
00270 lds = 0.0;
00271 for (i=0; i<n; i++) {
00272 sl = x[i];
00273 x[i] = y[i];
00274 y[i] = sl - y[i];
00275 sl = y[i];
00276 lds = lds + sl*sl;
00277 }
00278 lds = std::sqrt(lds);
00279 if (lds > SmaLL) {
00280 for (i=kl-1; i>=k; i--) {
00281 for (j=0; j < n; j++) v[j][i+1] = v[j][i];
00282 d[i+1] = d[i];
00283 }
00284 d[k] = 0.0;
00285 for (i=0; i < n; i++) v[i][k] = y[i] / lds;
00286 min_first_dir(k,4,&d[k],&lds,f1,1,x,n,v,w,q01);
00287 if (lds <= 0.0) {
00288 lds = -lds;
00289 for (i=0; i<n; i++) v[i][k] = -v[i][k];
00290 }
00291 }
00292 ldt = ldfac * ldt;
00293 if (ldt < lds) ldt = lds;
00294 if (prtlevel > 1) print(n,x);
00295 t2 = 0.0;
00296 for (i=0; i<n; i++) t2 += x[i]*x[i];
00297 t2 = m2 * std::sqrt(t2) + t;
00298 if (ldt > (0.5 * t2))
00299 kt = 0;
00300 else
00301 kt++;
00302 if (kt > ktm) goto fret;
00303 }
00304 /////////////////////////////////////////////
00305 // try quadratic extrapolation in case //
00306 // we are stuck in a curved valley //
00307 /////////////////////////////////////////////
00308 quadratic(x,n,v,w,q01);
00309 dn = 0.0;
00310 for (i=0; i<n; i++) {
00311 d[i] = 1.0 / std::sqrt(d[i]);
00312 if (dn < d[i]) dn = d[i];
00313 }
00314 if (prtlevel > 2) {
00315 std::cout << "\n-- New Matrix of Directions:\n";
00316 v.print();
00317 }
00318
00319 for (j=0; j<n; j++) {
00320 s = d[j] / dn;
00321 for (i=0; i < n; i++) v[i][j] *= s;
00322 }
00323 if (scbd > 1.0) { // scale axis to reduce condition number
00324 s = vlarge;
00325 for (i=0; i<n; i++) {
00326 sl = 0.0;
00327 for (j=0; j < n; j++) sl += v[i][j]*v[i][j];
00328 z[i] = std::sqrt(sl);
00329 if (z[i] < m4) z[i] = m4;
00330 if (s > z[i]) s = z[i];
00331 }
00332 for (i=0; i<n; i++) {
00333 sl = s / z[i];
00334 z[i] = 1.0 / sl;
00335 if (z[i] > scbd) {
00336 sl = 1.0 / scbd;
00337 z[i] = scbd;
00338 }
00339 }
00340 }
00341 for (i=1; i<n; i++) {
00342 for (j=0; j<=i-1; j++) {
00343 s = v[i][j];
00344 v[i][j] = v[j][i];
00345 v[j][i] = s;
00346 }
00347 }
00348 svd(n, macheps, vsmall, v.begin(), d.begin());
00349 if (scbd > 1.0) {
00350 for (i=0; i<n; i++) {
00351 s = z[i];
00352 for (j=0; j<n; j++) v[i][j] *= s;
00353 }
00354 for (i=0; i<n; i++) {
00355 s = 0.0;
00356 for (j=0; j<n; j++) s += v[j][i]*v[j][i];
00357 s = std::sqrt(s);
00358 d[i] *= s;
00359 s = 1.0 / s;
00360 for (j=0; j<n; j++) v[j][i] *= s;
00361 }
00362 }
00363 for (i=0; i<n; i++) {
00364 if ((dn * d[i]) > large)
00365 d[i] = vsmall;
00366 else if ((dn * d[i]) < SmaLL)
00367 d[i] = vlarge;
00368 else
00369 d[i] = pow(dn * d[i],-2.0);
00370 }
00371 sort(n,d.begin(),v.begin()); // the new eigenvalues and eigenvectors
00372 dmin = d[n-1];
00373 if (dmin < SmaLL) dmin = SmaLL;
00374 illc = (m2 * d[0]) > dmin;
00375 if (prtlevel > 2 && scbd > 1.0) {
00376 std::cout << "\n-- Scale Factors:\n";
00377 z.print();
00378 }
00379 if (prtlevel > 2) {
00380 std::cout << "\n-- Eigenvalues of A:\n";
00381 d.print();
00382 std::cout << "\n-- Eigenvectors of A:\n";
00383 v.print();
00384 }
00385
00386 if (maxfun > 0 && nl > maxfun) {
00387 if (prtlevel) std::cout << "\n-- maxiimum # of function calls reached \n";
00388 goto fret;
00389 }
00390 if (nf > n*10) step = 1.0;
00391 goto mloop; // back to main loop
00392
00393 fret:
00394 if (prtlevel > 0) {
00395 std::cout << "\n--after " << nf << " function calls\n";
00396 std::cout << "-- your function f(x) has finally been minimized to "
00397 << std::setprecision(W) << fx << "\n";
00398 std::cout.precision(6);
00399 std::cout << "-- with final x:\n";
00400 x.print();
00401 }
00402 maxfun = nf;
00403 return fx;
00404 }
|
|
||||||||||||||||||||||||
|
Definition at line 485 of file min_praxis.cpp. References matvec::fx, min_first_dir(), matvec::nl, matvec::qd, and matvec::qf1.
00487 { // look for a minimum along the curve q0, q1 //
00488 int i;
00489 double l, s,qa,qb,qc;
00490 s = fx; fx = qf1; qf1 = s; // swap fx and qf1
00491 qd[1] = 0.0;
00492 for (i=0; i<n; i++) {
00493 s = x[i]; l = q01[1][i]; x[i] = l; q01[1][i] = s;
00494 qd[1] += (s-l)*(s-l);
00495 }
00496 s = 0.0; qd[1] = std::sqrt(qd[1]); l = qd[1];
00497 if ( qd[0]>0.0 && qd[1]>0.0 && nl>=3*n*n ) {
00498 min_first_dir(-1,2, &s,&l,qf1,1,x, n,v,w,q01);
00499 qa = l*(l-qd[1])/(qd[0]*(qd[0]+qd[1]));
00500 qb = (l+qd[0])*(qd[1]-l)/(qd[0]*qd[1]);
00501 qc = l*(l+qd[0])/(qd[1]*(qd[0]+qd[1]));
00502 }
00503 else { qa = qb = 0.0; qc = 1.0; }
00504 qd[0] = qd[1];
00505 for (i=0; i<n; i++) {
00506 s = q01[0][i]; q01[0][i] = x[i];
00507 x[i] = qa*s + qb*x[i] + qc*q01[1][i];
00508 }
00509 }
|
|
|
Definition at line 43 of file minimizer.h. Referenced by matvec::Model::minfun(), Minimizer(), matvec::Model::multipoint_ml_estimate(), and matvec::Model::vce_dfreml(). |
|
|
Definition at line 36 of file minimizer.h. Referenced by min_prtlevel(), Minimizer(), and praxis(). |
1.2.16