diff --git a/libsrc/linalg/bfgs.cpp b/libsrc/linalg/bfgs.cpp new file mode 100644 index 00000000..d106b6f1 --- /dev/null +++ b/libsrc/linalg/bfgs.cpp @@ -0,0 +1,410 @@ +/***************************************************************************/ +/* */ +/* Vorlesung Optimierung I, Gfrerer, WS94/95 */ +/* BFGS-Verfahren zur Lösung freier nichtlinearer Optimierungsprobleme */ +/* */ +/* Programmautor: Joachim Schöberl */ +/* Matrikelnummer: 9155284 */ +/* */ +/***************************************************************************/ + +#include +#include + +#include +#include "opti.hpp" + + +namespace netgen +{ + +void Cholesky (const DenseMatrix & a, + DenseMatrix & l, Vector & d) +{ + // Factors A = L D L^T + + double x; + + int i, j, k; + int n = a.Height(); + + // (*testout) << "a = " << a << endl; + + l = a; + + for (i = 1; i <= n; i++) + { + for (j = i; j <= n; j++) + { + x = l.Get(i, j); + + for (k = 1; k < i; k++) + x -= l.Get(i, k) * l.Get(j, k) * d.Get(k); + + if (i == j) + { + d.Elem(i) = x; + } + else + { + l.Elem(j, i) = x / d.Get(k); + } + } + } + + for (i = 1; i <= n; i++) + { + l.Elem(i, i) = 1; + for (j = i+1; j <= n; j++) + l.Elem(i, j) = 0; + } + + /* + // Multiply: + (*testout) << "multiplied factors: " << endl; + for (i = 1; i <= n; i++) + for (j = 1; j <= n; j++) + { + x = 0; + for (k = 1; k <= n; k++) + x += l.Get(i, k) * l.Get(j, k) * d.Get(k); + (*testout) << x << " "; + } + (*testout) << endl; + */ +} + + +void MultLDLt (const DenseMatrix & l, const Vector & d, const Vector & g, Vector & p) +{ + /* + int i, j, n; + double val; + + n = l.Height(); + p = g; + for (i = 1; i <= n; i++) + { + val = 0; + for (j = i; j <= n; j++) + val += p.Get(j) * l.Get(j, i); + p.Set(i, val); + } + for (i = 1; i <= n; i++) + p.Elem(i) *= d.Get(i); + + for (i = n; i >= 1; i--) + { + val = 0; + for (j = 1; j <= i; j++) + val += p.Get(j) * l.Get(i, j); + p.Set(i, val); + } + */ + + + + double val; + + int n = l.Height(); + p = g; + + for (int i = 0; i < n; i++) + { + val = 0; + for (int j = i; j < n; j++) + val += p(j) * l(j, i); + p(i) = val; + } + + for (int i = 0; i < n; i++) + p(i) *= d(i); + + for (int i = n-1; i >= 0; i--) + { + val = 0; + for (int j = 0; j <= i; j++) + val += p(j) * l(i, j); + p(i) = val; + } +} + +void SolveLDLt (const DenseMatrix & l, const Vector & d, const Vector & g, Vector & p) +{ + double val; + + int n = l.Height(); + p = g; + + for (int i = 0; i < n; i++) + { + val = 0; + for (int j = 0; j < i; j++) + val += p(j) * l(i,j); + p(i) -= val; + } + + for (int i = 0; i < n; i++) + p(i) /= d(i); + + for (int i = n-1; i >= 0; i--) + { + val = 0; + for (int j = i+1; j < n; j++) + val += p(j) * l(j, i); + p(i) -= val; + } +} + +int LDLtUpdate (DenseMatrix & l, Vector & d, double a, const Vector & u) +{ + // Bemerkung: Es wird a aus R erlaubt + // Rueckgabewert: 0 .. D bleibt positiv definit + // 1 .. sonst + + int i, j, n; + + n = l.Height(); + + Vector v(n); + double t, told, xi; + + told = 1; + v = u; + + for (j = 1; j <= n; j++) + { + t = told + a * sqr (v.Elem(j)) / d.Get(j); + + if (t <= 0) + { + (*testout) << "update err, t = " << t << endl; + return 1; + } + + xi = a * v.Elem(j) / (d.Get(j) * t); + + d.Elem(j) *= t / told; + + for (i = j + 1; i <= n; i++) + { + v.Elem(i) -= v.Elem(j) * l.Elem(i, j); + l.Elem(i, j) += xi * v.Elem(i); + } + + told = t; + } + + return 0; +} + + +double BFGS ( + Vector & x, // i: Startwert + // o: Loesung, falls IFAIL = 0 + const MinFunction & fun, + const OptiParameters & par, + double eps + ) + + +{ + int i, j, n = x.Size(); + long it; + char a1crit, a3acrit; + + + Vector d(n), g(n), p(n), temp(n), bs(n), xneu(n), y(n), s(n), x0(n); + DenseMatrix l(n); + DenseMatrix hesse(n); + + double /* normg, */ alphahat, hd, fold; + double a1, a2; + const double mu1 = 0.1, sigma = 0.1, xi1 = 1, xi2 = 10; + const double tau = 0.1, tau1 = 0.1, tau2 = 0.6; + + Vector typx(x.Size()); // i: typische Groessenordnung der Komponenten + double f, f0; + double typf; // i: typische Groessenordnung der Loesung + double fmin = -1e5; // i: untere Schranke fuer Funktionswert + // double eps = 1e-8; // i: Abbruchschranke fuer relativen Gradienten + double tauf = 0.1; // i: Abbruchschranke fuer die relative Aenderung der + // Funktionswerte + int ifail; // o: 0 .. Erfolg + // -1 .. Unterschreitung von fmin + // 1 .. kein Erfolg bei Liniensuche + // 2 .. Überschreitung von itmax + + typx = par.typx; + typf = par.typf; + + + l = 0; + for (i = 1; i <= n; i++) + l.Elem(i, i) = 1; + + f = fun.FuncGrad (x, g); + f0 = f; + x0 = x; + + it = 0; + do + { + // Restart + + if (it % (5 * n) == 0) + { + + for (i = 1; i <= n; i++) + d.Elem(i) = typf/ sqr (typx.Get(i)); // 1; + for (i = 2; i <= n; i++) + for (j = 1; j < i; j++) + l.Elem(i, j) = 0; + + /* + hesse = 0; + for (i = 1; i <= n; i++) + hesse.Elem(i, i) = typf / sqr (typx.Get(i)); + + fun.ApproximateHesse (x, hesse); + + Cholesky (hesse, l, d); + */ + } + + it++; + if (it > par.maxit_bfgs) + { + ifail = 2; + break; + } + + + // Solve with factorized B + + SolveLDLt (l, d, g, p); + + // (*testout) << "l " << l << endl +// << "d " << d << endl +// << "g " << g << endl +// << "p " << p << endl; + + + p *= -1; + y = g; + + fold = f; + + // line search + + alphahat = 1; + lines (x, xneu, p, f, g, fun, par, alphahat, fmin, + mu1, sigma, xi1, xi2, tau, tau1, tau2, ifail); + + if(ifail == 1) + (*testout) << "no success with linesearch" << endl; + + /* + // if (it > par.maxit_bfgs/2) + { + (*testout) << "x = " << x << endl; + (*testout) << "xneu = " << xneu << endl; + (*testout) << "f = " << f << endl; + (*testout) << "g = " << g << endl; + } + */ + + // (*testout) << "it = " << it << " f = " << f << endl; + // if (ifail != 0) break; + + s.Set2 (1, xneu, -1, x); + y *= -1; + y.Add (1,g); // y += g; + + x = xneu; + + // BFGS Update + + MultLDLt (l, d, s, bs); + + a1 = y * s; + a2 = s * bs; + + if (a1 > 0 && a2 > 0) + { + if (LDLtUpdate (l, d, 1 / a1, y) != 0) + { + cerr << "BFGS update error1" << endl; + (*testout) << "BFGS update error1" << endl; + (*testout) << "l " << endl << l << endl + << "d " << d << endl; + ifail = 1; + break; + } + + if (LDLtUpdate (l, d, -1 / a2, bs) != 0) + { + cerr << "BFGS update error2" << endl; + (*testout) << "BFGS update error2" << endl; + (*testout) << "l " << endl << l << endl + << "d " << d << endl; + ifail = 1; + break; + } + } + + // Calculate stop conditions + + hd = eps * max2 (typf, fabs (f)); + a1crit = 1; + for (i = 1; i <= n; i++) + if ( fabs (g.Elem(i)) * max2 (typx.Elem(i), fabs (x.Elem(i))) > hd) + a1crit = 0; + + + a3acrit = (fold - f <= tauf * max2 (typf, fabs (f))); + + // testout << "g = " << g << endl; + // testout << "a1crit, a3crit = " << int(a1crit) << ", " << int(a3acrit) << endl; + + /* + // Output for tests + + normg = sqrt (g * g); + + testout << "it =" << setw (5) << it + << " f =" << setw (12) << setprecision (5) << f + << " |g| =" << setw (12) << setprecision (5) << normg; + + testout << " x = (" << setw (12) << setprecision (5) << x.Elem(1); + for (i = 2; i <= n; i++) + testout << "," << setw (12) << setprecision (5) << x.Elem(i); + testout << ")" << endl; + */ + + //(*testout) << "it = " << it << " f = " << f << " x = " << x << endl + // << " g = " << g << " p = " << p << endl << endl; + + // (*testout) << "|g| = " << g.L2Norm() << endl; + + if (g.L2Norm() < fun.GradStopping (x)) break; + + } + while (!a1crit || !a3acrit); + + /* + (*testout) << "it = " << it << " g = " << g << " f = " << f + << " fail = " << ifail << endl; + */ + if (f0 < f || (ifail == 1)) + { + (*testout) << "fail, f = " << f << " f0 = " << f0 << endl; + f = f0; + x = x0; + } + + // (*testout) << "x = " << x << ", x0 = " << x0 << endl; + return f; +} + +} diff --git a/libsrc/linalg/linopt.cpp b/libsrc/linalg/linopt.cpp new file mode 100644 index 00000000..a5d381e6 --- /dev/null +++ b/libsrc/linalg/linopt.cpp @@ -0,0 +1,73 @@ +#include +#include + +#include +#include "opti.hpp" + +namespace netgen +{ + +void LinearOptimize (const DenseMatrix & a, const Vector & b, + const Vector & c, Vector & x) + + { + int i1, i2, i3, j; + DenseMatrix m(3), inv(3); + Vector rs(3), hx(3), res(a.Height()), res2(3); + double f, fmin; + int nrest; + + if (a.Width() != 3) + { + cerr << "LinearOptimize only implemented for 3 unknowns" << endl; + return; + } + + fmin = 1e10; + x = 0; + nrest = a.Height(); + for (i1 = 1; i1 <= nrest; i1++) + for (i2 = i1 + 1; i2 <= nrest; i2++) + for (i3 = i2 + 1; i3 <= nrest; i3++) + { + for (j = 1; j <= 3; j++) + { + m.Elem(1, j) = a.Get(i1, j); + m.Elem(2, j) = a.Get(i2, j); + m.Elem(3, j) = a.Get(i3, j); + } + + rs.Elem(1) = b.Get(i1); + rs.Elem(2) = b.Get(i2); + rs.Elem(3) = b.Get(i3); + + if (fabs (m.Det()) < 1e-12) continue; + + CalcInverse (m, inv); + inv.Mult (rs, hx); + + a.Residuum (hx, b, res); +// m.Residuum (hx, rs, res2); + f = c * hx; + +/* + testout -> precision(12); + (*testout) << "i = (" << i1 << "," << i2 << "," << i3 + << "), f = " << f << " x = " << x << " res = " << res + << " resmin = " << res.Min() + << " res2 = " << res2 << " prod = " << prod << endl; +*/ + + + double rmin = res.Elem(1); + for (int hi = 2; hi <= res.Size(); hi++) + if (res.Elem(hi) < rmin) rmin = res.Elem(hi); + + if ( (f < fmin) && rmin >= -1e-8) + { + fmin = f; + x = hx; + } + } + } +} diff --git a/libsrc/linalg/linsearch.cpp b/libsrc/linalg/linsearch.cpp new file mode 100644 index 00000000..4c297fd9 --- /dev/null +++ b/libsrc/linalg/linsearch.cpp @@ -0,0 +1,351 @@ +/***************************************************************************/ +/* */ +/* Problem: Liniensuche */ +/* */ +/* Programmautor: Joachim Schöberl */ +/* Matrikelnummer: 9155284 */ +/* */ +/* Algorithmus nach: */ +/* */ +/* Optimierung I, Gfrerer, WS94/95 */ +/* Algorithmus 2.1: Liniensuche Problem (ii) */ +/* */ +/***************************************************************************/ + + + +#include + +#include // min, max, sqr + +#include +#include "opti.hpp" + + +namespace netgen +{ +const double eps0 = 1E-15; + +// Liniensuche + + +double MinFunction :: Func (const Vector & /* x */) const +{ + cerr << "Func of MinFunction called" << endl; + return 0; +} + +void MinFunction :: Grad (const Vector & /* x */, Vector & /* g */) const +{ + cerr << "Grad of MinFunction called" << endl; +} + +double MinFunction :: FuncGrad (const Vector & x, Vector & g) const +{ + cerr << "Grad of MinFunction called" << endl; + return 0; + /* + int n = x.Size(); + + static Vector xr; + static Vector xl; + xr.SetSize(n); + xl.SetSize(n); + + double eps = 1e-6; + double fl, fr; + + for (int i = 1; i <= n; i++) + { + xr.Set (1, x); + xl.Set (1, x); + + xr.Elem(i) += eps; + fr = Func (xr); + + xl.Elem(i) -= eps; + fl = Func (xl); + + g.Elem(i) = (fr - fl) / (2 * eps); + } + + double f = Func(x); + // (*testout) << "f = " << f << " grad = " << g << endl; + return f; + */ +} + + +double MinFunction :: FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const +{ + Vector g(x.Size()); + double f = FuncGrad (x, g); + deriv = (g * dir); + + // (*testout) << "g = " << g << ", dir = " << dir << ", deriv = " << deriv << endl; + return f; +} + +void MinFunction :: ApproximateHesse (const Vector & x, + DenseMatrix & hesse) const +{ + int n = x.Size(); + int i, j; + + static Vector hx; + hx.SetSize(n); + + double eps = 1e-6; + double f, f11, f12, f21, f22; + + for (i = 1; i <= n; i++) + { + for (j = 1; j < i; j++) + { + hx = x; + hx.Elem(i) = x.Get(i) + eps; + hx.Elem(j) = x.Get(j) + eps; + f11 = Func(hx); + hx.Elem(i) = x.Get(i) + eps; + hx.Elem(j) = x.Get(j) - eps; + f12 = Func(hx); + hx.Elem(i) = x.Get(i) - eps; + hx.Elem(j) = x.Get(j) + eps; + f21 = Func(hx); + hx.Elem(i) = x.Get(i) - eps; + hx.Elem(j) = x.Get(j) - eps; + f22 = Func(hx); + + hesse.Elem(i, j) = hesse.Elem(j, i) = + (f11 + f22 - f12 - f21) / (2 * eps * eps); + } + + hx = x; + f = Func(x); + hx.Elem(i) = x.Get(i) + eps; + f11 = Func(hx); + hx.Elem(i) = x.Get(i) - eps; + f22 = Func(hx); + + hesse.Elem(i, i) = (f11 + f22 - 2 * f) / (eps * eps); + } + // (*testout) << "hesse = " << hesse << endl; +} + + + + + + + +/// Line search, modified Mangasarien conditions +void lines (Vector & x, // i: initial point of line-search + Vector & xneu, // o: solution, if successful + Vector & p, // i: search direction + double & f, // i: function-value at x + // o: function-value at xneu, iff ifail = 0 + Vector & g, // i: gradient at x + // o: gradient at xneu, iff ifail = 0 + const MinFunction & fun, // function to minimize + const OptiParameters & par, + double & alphahat, // i: initial value for alpha_hat + // o: solution alpha iff ifail = 0 + double fmin, // i: lower bound for f + double mu1, // i: Parameter mu_1 of Alg.2.1 + double sigma, // i: Parameter sigma of Alg.2.1 + double xi1, // i: Parameter xi_1 of Alg.2.1 + double xi2, // i: Parameter xi_1 of Alg.2.1 + double tau, // i: Parameter tau of Alg.2.1 + double tau1, // i: Parameter tau_1 of Alg.2.1 + double tau2, // i: Parameter tau_2 of Alg.2.1 + int & ifail) // o: 0 on success + // -1 bei termination because lower limit fmin + // 1 bei illegal termination due to different reasons + +{ + double phi0, phi0prime, phi1, phi1prime, phihatprime; + double alpha1, alpha2, alphaincr, c; + char flag = 1; + long it; + + alpha1 = 0; + alpha2 = 1e50; + phi0 = phi1 = f; + + phi0prime = g * p; + + if (phi0prime > 0) + { + ifail = 1; + return; + } + + ifail = 1; // Markus + + phi1prime = phi0prime; + + // (*testout) << "phi0prime = " << phi0prime << endl; + + // it = 100000l; + it = 0; + + while (it++ <= par.maxit_linsearch) + { + + xneu.Set2 (1, x, alphahat, p); + + + // f = fun.FuncGrad (xneu, g); + // f = fun.Func (xneu); + f = fun.FuncDeriv (xneu, p, phihatprime); + + // (*testout) << "lines, f = " << f << " phip = " << phihatprime << endl; + + if (f < fmin) + { + ifail = -1; + break; + } + + + if (alpha2 - alpha1 < eps0 * alpha2) + { + ifail = 0; + break; + } + + // (*testout) << "i = " << it << " al = " << alphahat << " f = " << f << " fprime " << phihatprime << endl;; + + if (f - phi0 > mu1 * alphahat * phi1prime + eps0 * fabs (phi0)) + + { + + flag = 0; + alpha2 = alphahat; + + c = + (f - phi1 - phi1prime * (alphahat-alpha1)) / + sqr (alphahat-alpha1); + + alphahat = alpha1 - 0.5 * phi1prime / c; + + if (alphahat > alpha2) + alphahat = alpha1 + 1/(4*c) * + ( (sigma+mu1) * phi0prime - 2*phi1prime + + sqrt (sqr(phi1prime - mu1 * phi0prime) - + 4 * (phi1 - phi0 - mu1 * alpha1 * phi0prime) * c)); + + alphahat = max2 (alphahat, alpha1 + tau * (alpha2 - alpha1)); + alphahat = min2 (alphahat, alpha2 - tau * (alpha2 - alpha1)); + + // (*testout) << " if-branch" << endl; + + } + + else + + { + /* + f = fun.FuncGrad (xneu, g); + phihatprime = g * p; + */ + f = fun.FuncDeriv (xneu, p, phihatprime); + + if (phihatprime < sigma * phi0prime * (1 + eps0)) + + { + if (phi1prime < phihatprime) + // Approximationsfunktion ist konvex + + alphaincr = (alphahat - alpha1) * phihatprime / + (phi1prime - phihatprime); + + else + alphaincr = 1e99; // MAXDOUBLE; + + if (flag) + { + alphaincr = max2 (alphaincr, xi1 * (alphahat-alpha1)); + alphaincr = min2 (alphaincr, xi2 * (alphahat-alpha1)); + } + else + { + alphaincr = max2 (alphaincr, tau1 * (alpha2 - alphahat)); + alphaincr = min2 (alphaincr, tau2 * (alpha2 - alphahat)); + } + + alpha1 = alphahat; + alphahat += alphaincr; + phi1 = f; + phi1prime = phihatprime; + } + + else + + { + ifail = 0; // Erfolg !! + break; + } + + // (*testout) << " else, " << endl; + + } + + } + + // (*testout) << "linsearch: it = " << it << " ifail = " << ifail << endl; + + fun.FuncGrad (xneu, g); + + + if (it < 0) + ifail = 1; + + // (*testout) << "fail = " << ifail << endl; +} + + + + + + + + + + + + + + + + + + + +void SteepestDescent (Vector & x, const MinFunction & fun, + const OptiParameters & par) +{ + int it, n = x.Size(); + Vector xnew(n), p(n), g(n), g2(n); + double val, alphahat; + int fail; + + val = fun.FuncGrad(x, g); + + alphahat = 1; + // testout << "f = "; + for (it = 0; it < 10; it++) + { + // testout << val << " "; + + // p = -g; + p.Set (-1, g); + + lines (x, xnew, p, val, g, fun, par, alphahat, -1e5, + 0.1, 0.1, 1, 10, 0.1, 0.1, 0.6, fail); + + x = xnew; + } + // testout << endl; +} +} diff --git a/libsrc/linalg/opti.hpp b/libsrc/linalg/opti.hpp new file mode 100644 index 00000000..98757869 --- /dev/null +++ b/libsrc/linalg/opti.hpp @@ -0,0 +1,142 @@ +#ifndef FILE_OPTI +#define FILE_OPTI + +/**************************************************************************/ +/* File: opti.hpp */ +/* Author: Joachim Schoeberl */ +/* Date: 01. Jun. 95 */ +/**************************************************************************/ + + + +namespace netgen +{ + + /** + Function to be minimized. + */ + class MinFunction + { + public: + /// + virtual double Func (const Vector & x) const; + /// + virtual void Grad (const Vector & x, Vector & g) const; + /// function and gradient + virtual double FuncGrad (const Vector & x, Vector & g) const; + /// directional derivative + virtual double FuncDeriv (const Vector & x, const Vector & dir, double & deriv) const; + /// if |g| < gradaccuray, then stop bfgs + virtual double GradStopping (const Vector & /* x */) const { return 0; } + + /// + virtual void ApproximateHesse (const Vector & /* x */, + DenseMatrix & /* hesse */) const; + }; + + + class OptiParameters + { + public: + int maxit_linsearch; + int maxit_bfgs; + double typf; + double typx; + + OptiParameters () + { + maxit_linsearch = 100; + maxit_bfgs = 100; + typf = 1; + typx = 1; + } + }; + + + /** Implementation of BFGS method. + Efficient method for non-linear minimiztion problems. + @param x initial value and solution + @param fun function to be minimized + */ + extern double BFGS (Vector & x, const MinFunction & fun, + const OptiParameters & par, + double eps = 1e-8); + + /** Steepest descent method. + Simple method for non-linear minimization problems. + @param x initial value and solution + @param fun function to be minimized + */ + void SteepestDescent (Vector & x, const MinFunction & fun, + const OptiParameters & par); + + + extern void lines ( + Vector & x, // i: Ausgangspunkt der Liniensuche + Vector & xneu, // o: Loesung der Liniensuche bei Erfolg + Vector & p, // i: Suchrichtung + double & f, // i: Funktionswert an der Stelle x + // o: Funktionswert an der Stelle xneu, falls ifail = 0 + Vector & g, // i: Gradient an der Stelle x + // o: Gradient an der Stelle xneu, falls ifail = 0 + + const MinFunction & fun, // function to minmize + const OptiParameters & par, // parameters + double & alphahat, // i: Startwert für alpha_hat + // o: Loesung falls ifail = 0 + double fmin, // i: untere Schranke für f + double mu1, // i: Parameter mu_1 aus Alg.2.1 + double sigma, // i: Parameter sigma aus Alg.2.1 + double xi1, // i: Parameter xi_1 aus Alg.2.1 + double xi2, // i: Parameter xi_1 aus Alg.2.1 + double tau, // i: Parameter tau aus Alg.2.1 + double tau1, // i: Parameter tau_1 aus Alg.2.1 + double tau2, // i: Parameter tau_2 aus Alg.2.1 + int & ifail); // o: 0 bei erfolgreicher Liniensuche + // -1 bei Abbruch wegen Unterschreiten von fmin + // 1 bei Abbruch, aus sonstigen Gründen + + + + + /** + Solver for linear programming problem. + + \begin{verbatim} + min c^t x + A x <= b + \end{verbatim} + */ + extern void LinearOptimize (const DenseMatrix & a, const Vector & b, + const Vector & c, Vector & x); + + +#ifdef NONE + + /** + Simple projection iteration. + + find $u = argmin_{v >= 0} 0.5 u A u - f u$ + */ + extern void ApproxProject (const BaseMatrix & a, Vector & u, + const Vector & f, + double tau, int its); + + + /** + CG Algorithm for quadratic programming problem. + See: Dostal ... + + d ... diag(A) ^{-1} + */ + extern void ApproxProjectCG (const BaseMatrix & a, Vector & x, + const Vector & b, const class DiagMatrix & d, + double gamma, int & steps, int & changes); + +#endif + + +} + +#endif +