eonc::ConjugateGradients Class Reference

Decleration of the Conjugate Gradients optimizer. More...

#include <ConjugateGradients.h>

Inheritance diagram for eonc::ConjugateGradients:

Public Member Functions
	ConjugateGradients (std::shared_ptr< ObjectiveFunction > a_objf, const Parameters &a_params)
	Conjugate Gradients optimizer constructor.
	~ConjugateGradients ()=default
	Conjugant Gradient deconstructor.
int	step (double a_maxMove) override
	Calls the next step in the algorithm.
int	run (size_t a_maxIterations, double a_maxMove) override
	Runs the conjugate gradient.
Eigen::VectorXd	getStep ()
	Gets the direction of the next step.
Public Member Functions inherited from eonc::Optimizer
	Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, const OptimizerConfig &a_config)
	Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, OptType a_optype, const OptimizerConfig &a_config)
	Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, const Parameters &a_params)
	Optimizer (std::shared_ptr< ObjectiveFunction > a_objf, OptType a_optype, const Parameters &a_params)
virtual	~Optimizer ()=default

Private Member Functions
int	single_step (double a_maxMove)
	Steps the conjugate gradient.
int	line_search (double a_maxMove)
	Steps the conjugate gradient.

Private Attributes
Eigen::VectorXd	m_direction
	Current step direction of the conjugate gradient.
Eigen::VectorXd	m_directionOld
	Algorithms previous step direction.
Eigen::VectorXd	m_directionNorm
	Normalised version of the current direction vector.
Eigen::VectorXd	m_force
	Current force vector.
Eigen::VectorXd	m_forceOld
	Previous force vector.
eonc::log::FileScoped	m_log {"cg", "_cg.log"}
size_t	m_cg_i {0}
	Counts the number of discrete steps until algorithm convergence.

Additional Inherited Members
Protected Attributes inherited from eonc::Optimizer
const OptimizerConfig	m_optConfig
std::shared_ptr< ObjectiveFunction >	m_objf

Detailed Description

Decleration of the Conjugate Gradients optimizer.

Definition at line 38 of file ConjugateGradients.h.

Constructor & Destructor Documentation

ConjugateGradients()

eonc::ConjugateGradients::ConjugateGradients ( std::shared_ptr< ObjectiveFunction > a_objf, const Parameters & a_params )

inline

Conjugate Gradients optimizer constructor.

Parameters

std::shared_ptr<ObjectiveFunction>	m_objf that tells the optimizer how to run
const	Parameters& m_params defined by the config.init file

Definition at line 46 of file ConjugateGradients.h.

      : Optimizer(a_objf, OptType::CG, OptimizerConfig::fromParams(a_params)),
        m_direction{Eigen::VectorXd::Zero(a_objf->degreesOfFreedom())},
        m_directionOld{Eigen::VectorXd::Zero(a_objf->degreesOfFreedom())},
        m_directionNorm{Eigen::VectorXd::Zero(a_objf->degreesOfFreedom())},
        m_force{Eigen::VectorXd::Zero(a_objf->degreesOfFreedom())},
        m_forceOld{Eigen::VectorXd::Zero(a_objf->degreesOfFreedom())} {}

~ConjugateGradients()

eonc::ConjugateGradients::~ConjugateGradients ( )

default

Conjugant Gradient deconstructor.

Member Function Documentation

getStep()

Eigen::VectorXd eonc::ConjugateGradients::getStep

(

)

Gets the direction of the next step.

Definition at line 19 of file ConjugateGradients.cpp.

                                          {
  double a = std::abs(m_force.dot(m_forceOld));
  double b = m_forceOld.squaredNorm();
  double gamma = 0.0;
  if (a < 0.5 * b) {
    // Polak-Ribiere way to determine how much to mix in of old direction
    gamma = eonc::safemath::safe_div(m_force.dot(m_force - m_forceOld), b, 0.0);
  } else {
    gamma = 0;
  }
  m_direction = m_force + gamma * m_directionOld;
  m_directionNorm = m_direction;
  eonc::safemath::safe_normalize_inplace(m_directionNorm);
  m_directionOld = m_direction;
  m_forceOld = m_force;
 
  // Only if value for max nr of iteration before reset
  if (m_optConfig.opts.cg.max_iter_before_reset > 0 &&
      m_optConfig.opts.cg.max_iter_before_reset <= m_cg_i) {
    m_cg_i = 0;
    m_forceOld.setZero();
    m_directionOld.setZero();
  }
  m_cg_i += 1;
 
  return m_direction;
}

line_search()

int eonc::ConjugateGradients::line_search ( double a_maxMove )

private

Steps the conjugate gradient.

Checks for convergence based on the ratio of the projected force along a line to the norm of the total force

Definition at line 59 of file ConjugateGradients.cpp.

                                                    {
  Eigen::VectorXd pos;
  Eigen::VectorXd posStep;
  Eigen::VectorXd forceBeforeStep;
  double stepSize;
  double projectedForce;
  double projectedForceBeforeStep;
  double curvature;
 
  forceBeforeStep = -m_objf->getGradient();
  m_force = forceBeforeStep;
  getStep();
  pos = m_objf->getPositions();
  projectedForceBeforeStep = m_force.dot(m_directionNorm);
 
  // move system an infinitesimal step to determine the optimal step size along
  // the search line
  posStep = pos + m_directionNorm * m_optConfig.finiteDifference;
  m_objf->setPositions(posStep);
  m_force = -m_objf->getGradient(true);
  projectedForce = m_force.dot(m_directionNorm);
  stepSize = m_optConfig.finiteDifference;
 
  int line_i = 0;
  do {
    // Determine curvature from last step (Secant method)
    curvature = std::abs(eonc::safemath::safe_div(
        projectedForceBeforeStep - projectedForce, stepSize, 0.0));
    // A negative max move is a length (bowl breakout), not a backward step.
    const double maxMove = std::fabs(a_maxMove);
    stepSize = eonc::safemath::safe_div(projectedForce, curvature, maxMove);
 
    if (maxMove < std::fabs(stepSize)) {
      // first part get the sign of stepSize
      stepSize = ((stepSize > 0) - (stepSize < 0)) * maxMove;
    }
 
    forceBeforeStep = m_force;
    projectedForceBeforeStep = projectedForce;
 
    pos += stepSize * m_directionNorm;
    m_objf->setPositions(pos);
    m_force = -m_objf->getGradient();
    projectedForce = m_force.dot(m_directionNorm);
 
    line_i += 1;
 
    // Line search considered converged based in the ratio between the projected
    // force and the norm of the true force
  } while (m_optConfig.opts.cg.line_converged <
               std::abs(projectedForce) /
                   (std::sqrt(m_force.dot(m_force) +
                              m_optConfig.opts.cg.line_converged)) &&
           line_i < m_optConfig.opts.cg.line_search_max_iter);
  if (m_objf->isConverged())
    return 1;
  return 0;
}

run()

int eonc::ConjugateGradients::run ( size_t a_maxIterations, double a_maxMove )

overridevirtual

Runs the conjugate gradient.

Todo

method should also return an error code and message if the algorithm errors out

Returns

algorithm convergence

Implements eonc::Optimizer.

Definition at line 196 of file ConjugateGradients.cpp.

                                                                    {
  size_t iterations = 0;
  while (!m_objf->isConverged() && iterations < a_maxIterations) {
    step(a_maxMove);
    iterations++;
  }
  return m_objf->isConverged() ? 1 : 0;
}

single_step()

int eonc::ConjugateGradients::single_step ( double a_maxMove )

private

Steps the conjugate gradient.

Checks for convergence based on the change in displacement or direction

Definition at line 118 of file ConjugateGradients.cpp.

                                                    {
  Eigen::VectorXd pos;
  Eigen::VectorXd posStep;
  Eigen::VectorXd forceAfterStep;
 
  m_force = -m_objf->getGradient();
  pos = m_objf->getPositions();
  getStep();
 
  // move system an infinitesimal step to determine the optimal step size along
  // the search line
  posStep = pos + m_directionNorm * m_optConfig.finiteDifference;
  m_objf->setPositions(posStep);
  forceAfterStep = -m_objf->getGradient(true);
 
  // Determine curvature
  double projectedForce1 = m_force.dot(m_directionNorm);
  double projectedForce2 = forceAfterStep.dot(m_directionNorm);
  double curvature = eonc::safemath::safe_div(
      projectedForce1 - projectedForce2, m_optConfig.finiteDifference, 0.0);
 
  double stepSize = a_maxMove;
 
  if (curvature > 0.0) {
    stepSize = eonc::safemath::safe_div(projectedForce1, curvature, a_maxMove);
  }
 
  if (m_optConfig.bowlBreakout && a_maxMove < 0.0) {
    stepSize = -a_maxMove;
    a_maxMove = -a_maxMove;
  }
 
  if (!m_optConfig.opts.cg.no_overshooting) {
    if (m_optConfig.bowlBreakout) {
      // max displacement is based on system not single atom
      pos += eonc::geometry::maxMotionAppliedV(stepSize * m_directionNorm,
                                               a_maxMove);
    } else {
      pos += eonc::geometry::maxAtomMotionAppliedV(stepSize * m_directionNorm,
                                                   a_maxMove);
    }
    m_objf->setPositions(pos);
  } else {
    // negative if product of the projected forces before and after the step are
    // in opposite directions
    double passedMinimum = -1.;
    double forceChange = 0.;
    while (passedMinimum < 0.0 &&
           0.1 * std::abs(projectedForce1) < std::abs(projectedForce2)) {
      posStep = pos + eonc::geometry::maxAtomMotionAppliedV(
                          stepSize * m_directionNorm, a_maxMove);
      m_objf->setPositions(posStep);
      forceAfterStep = -m_objf->getGradient(true);
      projectedForce2 = forceAfterStep.dot(m_directionNorm);
 
      passedMinimum = projectedForce1 * projectedForce2;
      if (passedMinimum < 0.0 &&
          0.1 * std::abs(projectedForce1) < std::abs(projectedForce2)) {
        forceChange = (projectedForce1 - projectedForce2);
        stepSize = eonc::safemath::safe_div(projectedForce1, forceChange, 0.0) *
                   stepSize;
        QUILL_LOG_DEBUG(m_log, "Force changed {}, step size adjusted to {}",
                        forceChange, stepSize);
      }
    }
  }
  if (m_optConfig.opts.cg.knock_out_max_move) {
    if (stepSize >= a_maxMove) {
      // knockout old search direction
      m_directionOld.setZero();
      m_forceOld.setZero();
      QUILL_LOG_DEBUG(m_log, "Resetting the old search direction");
    }
  }
 
  return m_objf->isConverged() ? 1 : 0;
}

step()

int eonc::ConjugateGradients::step ( double a_maxMove )

overridevirtual

Calls the next step in the algorithm.

Either calls the single_step or line_search method depending on the parameters

Returns

whether or not the algorithm has converged

Implements eonc::Optimizer.

Definition at line 47 of file ConjugateGradients.cpp.

                                             {
  bool converged;
  if (m_optConfig.opts.cg.line_search) {
    converged = line_search(a_maxMove);
  } else {
    converged = single_step(a_maxMove);
  }
  if (converged)
    return 1;
  return 0;
}

Member Data Documentation

m_cg_i

size_t eonc::ConjugateGradients::m_cg_i {0}

private

Counts the number of discrete steps until algorithm convergence.

Definition at line 86 of file ConjugateGradients.h.

{0};

m_direction

Eigen::VectorXd eonc::ConjugateGradients::m_direction

private

Current step direction of the conjugate gradient.

Definition at line 74 of file ConjugateGradients.h.

m_directionNorm

Eigen::VectorXd eonc::ConjugateGradients::m_directionNorm

private

Normalised version of the current direction vector.

Definition at line 78 of file ConjugateGradients.h.

m_directionOld

Eigen::VectorXd eonc::ConjugateGradients::m_directionOld

private

Algorithms previous step direction.

Definition at line 76 of file ConjugateGradients.h.

m_force

Eigen::VectorXd eonc::ConjugateGradients::m_force

private

Current force vector.

Definition at line 80 of file ConjugateGradients.h.

m_forceOld

Eigen::VectorXd eonc::ConjugateGradients::m_forceOld

private

Previous force vector.

Definition at line 82 of file ConjugateGradients.h.

m_log

eonc::log::FileScoped eonc::ConjugateGradients::m_log {"cg", "_cg.log"}

private

Definition at line 83 of file ConjugateGradients.h.

{"cg", "_cg.log"};

---

The documentation for this class was generated from the following files:

• /home/runner/work/eOn/eOn/include/eon/ConjugateGradients.h
• /home/runner/work/eOn/eOn/client/ConjugateGradients.cpp