Loading...
Searching...
No Matches
BFGSOpt.h
Go to the documentation of this file.
//
// Copyright (C) 2004-2008 Greg Landrum and Rational Discovery LLC
//
// @@ All Rights Reserved @@
// This file is part of the RDKit.
// The contents are covered by the terms of the BSD license
// which is included in the file license.txt, found at the root
// of the RDKit source tree.
//
#ifndef RD_BFGSOPT_H
#define RD_BFGSOPT_H
#include <RDGeneral/export.h>
#include <cmath>
#include <RDGeneral/Invariant.h>
#include <GraphMol/Trajectory/Snapshot.h>
#include <cstring>
#include <vector>
#include <algorithm>
#include "BFGSOpt_SVE.h"
namespace BFGSOpt {
RDKIT_OPTIMIZER_EXPORT extern int HEAD_ONLY_LIBRARY;
RDKIT_OPTIMIZER_EXPORT extern int REALLY_A_HEADER_ONLY_LIBRARY;
const double FUNCTOL =
1e-4; //!< Default tolerance for function convergence in the minimizer
const double MOVETOL =
1e-7; //!< Default tolerance for x changes in the minimizer
const int MAXITS = 200; //!< Default maximum number of iterations
const double EPS = 3e-8; //!< Default gradient tolerance in the minimizer
const double TOLX =
4. * EPS; //!< Default direction vector tolerance in the minimizer
const double MAXSTEP = 100.0; //!< Default maximum step size in the minimizer
/*!
See Numerical Recipes in C, Section 9.7 for a description of the algorithm.
\param dim the dimensionality of the space.
\param oldPt the current position, as an array.
\param oldVal the current function value.
\param grad the value of the function gradient at oldPt
\param dir the minimization direction
\param newPt used to return the final position
\param newVal used to return the final function value
\param func the function to minimize
\param maxStep the maximum allowable step size
\param resCode used to return the results of the search.
Possible values for resCode are on return are:
- 0: success
- 1: the stepsize got too small. This probably indicates success.
- -1: the direction is bad (orthogonal to the gradient)
*/
template <typename EnergyFunctor>
void linearSearch(unsigned int dim, double *oldPt, double oldVal, double *grad,
double *dir, double *newPt, double &newVal,
EnergyFunctor func, double maxStep, int &resCode) {
PRECONDITION(oldPt, "bad input array");
PRECONDITION(grad, "bad input array");
PRECONDITION(dir, "bad input array");
PRECONDITION(newPt, "bad input array");
const unsigned int MAX_ITER_LINEAR_SEARCH = 1000;
double sum = 0.0, slope = 0.0, test = 0.0, lambda = 0.0;
double lambda2 = 0.0, lambdaMin = 0.0, tmpLambda = 0.0, val2 = 0.0;
resCode = -1;
// get the length of the direction vector:
sum = 0.0;
for (unsigned int i = 0; i < dim; i++) {
sum += dir[i] * dir[i];
}
sum = sqrt(sum);
// Rescale if we're trying to move too far
if (sum > maxStep) {
for (unsigned int i = 0; i < dim; i++) {
dir[i] *= maxStep / sum;
}
}
// make sure our direction has at least some component along
// -grad
slope = 0.0;
for (unsigned int i = 0; i < dim; i++) {
slope += dir[i] * grad[i];
}
if (slope >= 0.0) {
return;
}
test = 0.0;
for (unsigned int i = 0; i < dim; i++) {
double temp = fabs(dir[i]) / std::max(fabs(oldPt[i]), 1.0);
if (temp > test) {
test = temp;
}
}
lambdaMin = MOVETOL / test;
lambda = 1.0;
unsigned int it = 0;
while (it < MAX_ITER_LINEAR_SEARCH) {
if (lambda < lambdaMin) {
// Step size is below the position-scaled threshold; treat as converged
resCode = 1;
break;
}
for (unsigned int i = 0; i < dim; i++) {
newPt[i] = oldPt[i] + lambda * dir[i];
}
newVal = func(newPt);
if (newVal - oldVal <= FUNCTOL * lambda * slope) {
// Armijo sufficient-decrease condition satisfied; accept the step
resCode = 0;
return;
}
// if we made it this far, we need to backtrack:
if (it == 0) {
// Quadratic model: only one prior function value available
tmpLambda = -slope / (2.0 * (newVal - oldVal - slope));
} else {
double rhs1 = newVal - oldVal - lambda * slope;
double rhs2 = val2 - oldVal - lambda2 * slope;
double a = (rhs1 / (lambda * lambda) - rhs2 / (lambda2 * lambda2)) /
(lambda - lambda2);
double b = (-lambda2 * rhs1 / (lambda * lambda) +
lambda * rhs2 / (lambda2 * lambda2)) /
(lambda - lambda2);
if (a == 0.0) {
tmpLambda = -slope / (2.0 * b);
} else {
double disc = b * b - 3 * a * slope;
if (disc < 0.0) {
tmpLambda = 0.5 * lambda;
} else if (b <= 0.0) {
tmpLambda = (-b + sqrt(disc)) / (3.0 * a);
} else {
tmpLambda = -slope / (b + sqrt(disc));
}
}
if (tmpLambda > 0.5 * lambda) {
tmpLambda = 0.5 * lambda;
}
}
lambda2 = lambda;
val2 = newVal;
lambda = std::max(tmpLambda, 0.1 * lambda);
++it;
}
// nothing was done
for (unsigned int i = 0; i < dim; i++) {
newPt[i] = oldPt[i];
}
}
//! Do a BFGS minimization of a function.
/*!
See Numerical Recipes in C, Section 10.7 for a description of the algorithm.
\param dim the dimensionality of the space.
\param pos the starting position, as an array.
\param gradTol tolerance for gradient convergence
\param numIters used to return the number of iterations required
\param funcVal used to return the final function value
\param func the function to minimize
\param gradFunc calculates the gradient of func
\param funcTol tolerance for changes in the function value for convergence.
\param maxIts maximum number of iterations allowed
\param snapshotFreq a snapshot of the minimization trajectory
will be stored after as many steps as indicated
through this parameter; defaults to 0 (no
snapshots stored)
\param snapshotVect pointer to a std::vector<Snapshot> object that will
receive the coordinates and energies every snapshotFreq steps; defaults to
NULL (no snapshots stored)
\return a flag indicating success (or type of failure). Possible values are:
- 0: success
- 1: too many iterations were required
*/
template <typename EnergyFunctor, typename GradientFunctor>
int minimize(unsigned int dim, double *pos, double gradTol,
unsigned int &numIters, double &funcVal, EnergyFunctor func,
GradientFunctor gradFunc, unsigned int snapshotFreq,
RDKit::SnapshotVect *snapshotVect, double funcTol = TOLX,
unsigned int maxIts = MAXITS) {
RDUNUSED_PARAM(funcTol);
PRECONDITION(pos, "bad input array");
PRECONDITION(gradTol > 0, "bad tolerance");
std::vector<double> grad(dim);
std::vector<double> dGrad(dim);
std::vector<double> hessDGrad(dim);
std::vector<double> xi(dim);
std::vector<double> invHessian(dim * dim, 0);
std::unique_ptr<double[]> newPos(new double[dim]);
snapshotFreq = std::min(snapshotFreq, maxIts);
double fp = func(pos);
gradFunc(pos, grad.data());
double sum = 0.0;
#ifdef RDK_SVE_AVAILABLE
if (cpuHasSVE()) {
// SVE path: initialise xi = -grad and compute ||pos||^2 in a single
// vectorised pass. The identity inverse Hessian is initialised separately
// (scalar, O(dim)) since it is a simple diagonal write and does not benefit
// from vectorisation over rows.
sveInitXiAndSum(dim, grad.data(), xi.data(), pos, &sum);
for (unsigned int i = 0; i < dim; i++) invHessian[i * dim + i] = 1.0;
} else
#endif
{
// Scalar path: initialise the inverse Hessian to the identity matrix,
// set the initial search direction xi = -grad (steepest descent step),
// and accumulate ||pos||^2 to set an appropriate maximum step size.
for (unsigned int i = 0; i < dim; i++) {
unsigned int itab = i * dim;
invHessian[itab + i] = 1.0;
xi[i] = -grad[i];
sum += pos[i] * pos[i];
}
}
double maxStep = MAXSTEP * std::max(sqrt(sum), static_cast<double>(dim));
for (unsigned int iter = 1; iter <= maxIts; ++iter) {
numIters = iter;
int status = -1;
linearSearch(dim, pos, fp, grad.data(), xi.data(), newPos.get(), funcVal,
func, maxStep, status);
CHECK_INVARIANT(status >= 0, "bad direction in linearSearch");
// save the function value for the next search:
fp = funcVal;
// set the direction of this line and save the gradient:
double test = 0.0;
for (unsigned int i = 0; i < dim; i++) {
xi[i] = newPos[i] - pos[i];
pos[i] = newPos[i];
double temp = fabs(xi[i]) / std::max(fabs(pos[i]), 1.0);
if (temp > test) {
test = temp;
}
dGrad[i] = grad[i];
}
if (test < TOLX) {
if (snapshotVect && snapshotFreq) {
RDKit::Snapshot s(boost::shared_array<double>(newPos.release()), fp);
snapshotVect->push_back(s);
}
return 0;
}
// update the gradient:
double gradScale = gradFunc(pos, grad.data());
test = 0.0;
// Use |funcVal| so that negative energies (which arise routinely
// mid-minimization in force fields containing stabilizing
// electrostatic or dispersion terms) do not drive
// funcVal * gradScale below zero and clamp the denominator to 1.0,
// which would artificially tighten the gradient convergence test.
double term = std::max(fabs(funcVal) * gradScale, 1.0);
for (unsigned int i = 0; i < dim; i++) {
double temp = fabs(grad[i]) * std::max(fabs(pos[i]), 1.0);
test = std::max(test, temp);
dGrad[i] = grad[i] - dGrad[i];
}
test /= term;
if (test < gradTol) {
if (snapshotVect && snapshotFreq) {
RDKit::Snapshot s(boost::shared_array<double>(newPos.release()), fp);
snapshotVect->push_back(s);
}
return 0;
}
// BFGS inverse Hessian update.
double fac = 0, fae = 0, sumDGrad = 0, sumXi = 0;
#ifdef RDK_SVE_AVAILABLE
if (cpuHasSVE()) {
// SVE path: matrix-vector multiply and all four dot products computed in
// one vectorised pass, saving two additional O(dim) traversals compared
// to separate scalar dot-product calls.
sveHessianVecMul(dim, invHessian.data(), dGrad.data(), hessDGrad.data(),
xi.data(), &fac, &fae, &sumDGrad, &sumXi);
} else
#endif
{
// Scalar path: fused matrix-vector multiply and dot-product accumulation.
// Pointer arithmetic (++ivh, ++dgj) avoids repeated index computations
// and helps the compiler generate efficient load sequences.
for (unsigned int i = 0; i < dim; i++) {
double *ivh = &(invHessian[i * dim]);
double &hdgradi = hessDGrad[i];
double *dgj = dGrad.data();
hdgradi = 0.0;
for (unsigned int j = 0; j < dim; ++j, ++ivh, ++dgj) {
hdgradi += *ivh * *dgj;
}
fac += dGrad[i] * xi[i];
fae += dGrad[i] * hessDGrad[i];
sumDGrad += dGrad[i] * dGrad[i];
sumXi += xi[i] * xi[i];
}
}
if (fac > sqrt(EPS * sumDGrad * sumXi)) {
fac = 1.0 / fac;
double fad = 1.0 / fae;
for (unsigned int i = 0; i < dim; i++) {
dGrad[i] = fac * xi[i] - fad * hessDGrad[i];
}
#ifdef RDK_SVE_AVAILABLE
if (cpuHasSVE()) {
// SVE path: symmetric rank-1 update with FMA, exploiting symmetry to
// halve memory writes and FLOPs versus a full-matrix update
sveHessianRank1Update(dim, invHessian.data(), xi.data(),
hessDGrad.data(), dGrad.data(), fac, fad, fae);
} else
#endif
{
// Scalar path: upper-triangle-only update (j >= i) followed by
// explicit symmetrisation. This halves the number of Hessian writes
// at the cost of one additional pass over a row to mirror elements.
for (unsigned int i = 0; i < dim; i++) {
unsigned int itab = i * dim;
double pxi = fac * xi[i], hdgi = fad * hessDGrad[i],
dgi = fae * dGrad[i];
double *pxj = &(xi[i]), *hdgj = &(hessDGrad[i]), *dgj = &(dGrad[i]);
for (unsigned int j = i; j < dim; ++j, ++pxj, ++hdgj, ++dgj) {
invHessian[itab + j] += pxi * *pxj - hdgi * *hdgj + dgi * *dgj;
invHessian[j * dim + i] = invHessian[itab + j];
}
}
}
}
#ifdef RDK_SVE_AVAILABLE
if (cpuHasSVE()) {
sveHessianVecMulNeg(dim, invHessian.data(), grad.data(), xi.data());
} else
#endif
{
for (unsigned int i = 0; i < dim; i++) {
unsigned int itab = i * dim;
xi[i] = 0.0;
double &pxi = xi[i];
double *ivh = &(invHessian[itab]);
double *gj = grad.data();
for (unsigned int j = 0; j < dim; ++j, ++ivh, ++gj) {
pxi -= *ivh * *gj;
}
}
}
if (snapshotVect && snapshotFreq && !(iter % snapshotFreq)) {
RDKit::Snapshot s(boost::shared_array<double>(newPos.release()), fp);
snapshotVect->push_back(s);
newPos.reset(new double[dim]);
}
}
return 1;
}
//! Do a BFGS minimization of a function.
/*!
\param dim the dimensionality of the space.
\param pos the starting position, as an array.
\param gradTol tolerance for gradient convergence
\param numIters used to return the number of iterations required
\param funcVal used to return the final function value
\param func the function to minimize
\param gradFunc calculates the gradient of func
\param funcTol tolerance for changes in the function value for convergence.
\param maxIts maximum number of iterations allowed
\return a flag indicating success (or type of failure). Possible values are:
- 0: success
- 1: too many iterations were required
*/
template <typename EnergyFunctor, typename GradientFunctor>
int minimize(unsigned int dim, double *pos, double gradTol,
unsigned int &numIters, double &funcVal, EnergyFunctor func,
GradientFunctor gradFunc, double funcTol = TOLX,
unsigned int maxIts = MAXITS) {
return minimize(dim, pos, gradTol, numIters, funcVal, func, gradFunc, 0,
nullptr, funcTol, maxIts);
}
} // namespace BFGSOpt
#endif // RD_BFGSOPT_H
BFGSOpt_SVE.h
Invariant.h
CHECK_INVARIANT
#define CHECK_INVARIANT(expr, mess)
Definition Invariant.h:100
RDUNUSED_PARAM
#define RDUNUSED_PARAM(x)
Definition Invariant.h:196
PRECONDITION
#define PRECONDITION(expr, mess)
Definition Invariant.h:108
Snapshot.h
RDKit::Snapshot
Definition Snapshot.h:26
export.h
RDKIT_OPTIMIZER_EXPORT
#define RDKIT_OPTIMIZER_EXPORT
Definition export.h:401
BFGSOpt
Definition BFGSOpt.h:22
BFGSOpt::EPS
const double EPS
Default gradient tolerance in the minimizer.
Definition BFGSOpt.h:30
BFGSOpt::TOLX
const double TOLX
Default direction vector tolerance in the minimizer.
Definition BFGSOpt.h:31
BFGSOpt::linearSearch
void linearSearch(unsigned int dim, double *oldPt, double oldVal, double *grad, double *dir, double *newPt, double &newVal, EnergyFunctor func, double maxStep, int &resCode)
Definition BFGSOpt.h:55
BFGSOpt::cpuHasSVE
static bool cpuHasSVE()
Definition BFGSOpt_SVE.h:16
BFGSOpt::FUNCTOL
const double FUNCTOL
Default tolerance for function convergence in the minimizer.
Definition BFGSOpt.h:25
BFGSOpt::minimize
int minimize(unsigned int dim, double *pos, double gradTol, unsigned int &numIters, double &funcVal, EnergyFunctor func, GradientFunctor gradFunc, unsigned int snapshotFreq, RDKit::SnapshotVect *snapshotVect, double funcTol=TOLX, unsigned int maxIts=MAXITS)
Do a BFGS minimization of a function.
Definition BFGSOpt.h:184
BFGSOpt::MAXSTEP
const double MAXSTEP
Default maximum step size in the minimizer.
Definition BFGSOpt.h:33
BFGSOpt::REALLY_A_HEADER_ONLY_LIBRARY
RDKIT_OPTIMIZER_EXPORT int REALLY_A_HEADER_ONLY_LIBRARY
BFGSOpt::MOVETOL
const double MOVETOL
Default tolerance for x changes in the minimizer.
Definition BFGSOpt.h:27
BFGSOpt::MAXITS
const int MAXITS
Default maximum number of iterations.
Definition BFGSOpt.h:29
BFGSOpt::HEAD_ONLY_LIBRARY
RDKIT_OPTIMIZER_EXPORT int HEAD_ONLY_LIBRARY
RDKit::SnapshotVect
std::vector< Snapshot > SnapshotVect
Definition Snapshot.h:21