Skip to content

Commit

Permalink
Merge pull request #52 from arvigj/debug_fd
Browse files Browse the repository at this point in the history 
Add option to run finite difference on gradient of energy
  • Loading branch information
@arvigj
arvigj authored Dec 7, 2023
2 parents bd91fe5 + 035b97c commit 8a8c103
Show file tree
Hide file tree
Showing 6 changed files with 205 additions and 10 deletions.
5 changes: 5 additions & 0 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,6 +205,10 @@ target_link_libraries(polysolve_linear PUBLIC jse::jse)
include(spdlog)
target_link_libraries(polysolve_linear PUBLIC spdlog::spdlog)

# finite-diff
include(finite-diff)
target_link_libraries(polysolve PUBLIC finitediff::finitediff)

# Accelerate solver
if(POLYSOLVE_WITH_ACCELERATE)
include(CPM)
Expand Down Expand Up @@ -322,6 +326,7 @@ if(POLYSOLVE_WITH_CUSOLVER)
endif()
endif()


################################################################################
# Compiler options
################################################################################
Expand Down
11 changes: 11 additions & 0 deletions cmake/recipes/finite-diff.cmake
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
# finite-diff (https://github.com/zfergus/finite-diff)
# License: MIT

if(TARGET finitediff::finitediff)
return()
endif()

message(STATUS "Third-party: creating target 'finitediff::finitediff'")

include(CPM)
CPMAddPackage("gh:zfergus/finite-diff@1.0.2")
21 changes: 20 additions & 1 deletion non-linear-solver-spec.json
View file
Original file line number Diff line number Diff line change
Expand Up @@ -303,7 +303,9 @@
"type": "object",
"optional": [
"f_delta",
"f_delta_step_tol"
"f_delta_step_tol",
"apply_gradient_fd",
"gradient_fd_eps"
],
"doc": "Nonlinear solver advanced options"
},
Expand All @@ -319,5 +321,22 @@
"default": 100,
"type": "int",
"doc": "Dangerous Option: Quit the optimization if the solver reduces the energy by less than f_delta for consecutive f_delta_step_tol steps."
},
{
"pointer": "/advanced/apply_gradient_fd",
"default": "None",
"type": "string",
"options": [
"None",
"DirectionalDerivative",
"FullFiniteDiff"
],
"doc": "Expensive Option: For every iteration of the nonlinear solver, run finite difference to verify gradient of energy."
},
{
"pointer": "/advanced/gradient_fd_eps",
"default": 1e-7,
"type": "float",
"doc": "Expensive Option: Eps for finite difference to verify gradient of energy."
}
]
70 changes: 70 additions & 0 deletions src/polysolve/nonlinear/Solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,18 @@
#include <spdlog/fmt/bundled/color.h>
#include <spdlog/fmt/ostr.h>

#include <finitediff.hpp>

#include <iomanip>
#include <fstream>

namespace polysolve::nonlinear
{
NLOHMANN_JSON_SERIALIZE_ENUM(
FiniteDiffStrategy,
{{FiniteDiffStrategy::NONE, "None"},
{FiniteDiffStrategy::DIRECTIONAL_DERIVATIVE, "DirectionalDerivative"},
{FiniteDiffStrategy::FULL_FINITE_DIFF, "FullFiniteDiff"}})

// Static constructor
std::unique_ptr<Solver> Solver::create(
Expand Down Expand Up @@ -130,6 +137,9 @@ namespace polysolve::nonlinear
m_descent_strategy = 0;

set_line_search(solver_params);

gradient_fd_strategy = solver_params["advanced"]["apply_gradient_fd"];
gradient_fd_eps = solver_params["advanced"]["gradient_fd_eps"];
}

void Solver::set_strategies_iterations(const json &solver_params)
Expand Down Expand Up @@ -234,6 +244,11 @@ namespace polysolve::nonlinear
objFunc.gradient(x, grad);
}

{
POLYSOLVE_SCOPED_STOPWATCH("verify gradient", grad_time, m_logger);
this->verify_gradient(objFunc, x, grad);
}

const double grad_norm = compute_grad_norm(x, grad);
if (std::isnan(grad_norm))
{
Expand Down Expand Up @@ -495,4 +510,59 @@ namespace polysolve::nonlinear
m_line_search ? m_line_search->broad_phase_ccd_time : 0, m_line_search ? m_line_search->ccd_time : 0,
m_line_search ? m_line_search->classical_line_search_time : 0);
}

void Solver::verify_gradient(Problem &objFunc, const TVector &x, const TVector &grad)
{
bool match = false;

switch (gradient_fd_strategy)
{
case FiniteDiffStrategy::NONE:
return;
case FiniteDiffStrategy::DIRECTIONAL_DERIVATIVE:
{
Eigen::VectorXd direc = grad.normalized();
Eigen::VectorXd x2 = x + direc * gradient_fd_eps;
Eigen::VectorXd x1 = x - direc * gradient_fd_eps;

objFunc.solution_changed(x2);
double J2 = objFunc.value(x2);

objFunc.solution_changed(x1);
double J1 = objFunc.value(x1);

double fd = (J2 - J1) / 2 / gradient_fd_eps;
double analytic = direc.dot(grad);

match = abs(fd - analytic) < 1e-8 || abs(fd - analytic) < 1e-4 * abs(analytic);

// Log error in either case to make it more visible in the logs.
if (match)
m_logger.debug("step size: {}, finite difference: {}, derivative: {}", gradient_fd_eps, fd, analytic);
else
m_logger.error("step size: {}, finite difference: {}, derivative: {}", gradient_fd_eps, fd, analytic);
}
break;
case FiniteDiffStrategy::FULL_FINITE_DIFF:
{
Eigen::VectorXd grad_fd;
fd::finite_gradient(
x, [&](const Eigen::VectorXd &x_) {
objFunc.solution_changed(x_);
return objFunc.value(x_);
},
grad_fd, fd::AccuracyOrder::SECOND, gradient_fd_eps);

match = (grad_fd - grad).norm() < 1e-8 || (grad_fd - grad).norm() < 1e-4 * (grad).norm();

if (match)
m_logger.debug("step size: {}, all gradient components match finite difference", gradient_fd_eps);
else
m_logger.error("step size: {}, all gradient components do not match finite difference", gradient_fd_eps);
}
break;
}

objFunc.solution_changed(x);
}
} // namespace polysolve::nonlinear
15 changes: 15 additions & 0 deletions src/polysolve/nonlinear/Solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,13 @@ namespace polysolve::nonlinear
SUCCESS = 0,
};

enum class FiniteDiffStrategy
{
NONE = 0,
DIRECTIONAL_DERIVATIVE = 1,
FULL_FINITE_DIFF = 2
};

class Solver : public cppoptlib::ISolver<Problem, /*Ord=*/-1>
{
public:
Expand Down Expand Up @@ -99,6 +106,14 @@ namespace polysolve::nonlinear

spdlog::logger &m_logger;

// ====================================================================
// Finite Difference Utilities
// ====================================================================

FiniteDiffStrategy gradient_fd_strategy = FiniteDiffStrategy::NONE;
double gradient_fd_eps = 1e-7;
virtual void verify_gradient(Problem &objFunc, const TVector &x, const TVector &grad) final;

private:
// ====================================================================
// Solver parameters
Expand Down
93 changes: 84 additions & 9 deletions tests/test_nonlinear_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -240,14 +240,18 @@ class Beale : public AnalyticTestProblem
}
};


class InequalityConstraint : public Problem
{
public:
InequalityConstraint(const double upper_bound): upper_bound_(upper_bound) {}
InequalityConstraint(const double upper_bound) : upper_bound_(upper_bound) {}
double value(const TVector &x) override { return x(0) - upper_bound_; }
void gradient(const TVector &x, TVector &gradv) override { gradv.setZero(x.size()); gradv(0) = 1; }
void gradient(const TVector &x, TVector &gradv) override
{
gradv.setZero(x.size());
gradv(0) = 1;
}
void hessian(const TVector &x, THessian &hessian) override {}

private:
const double upper_bound_;
};
Expand Down Expand Up @@ -351,12 +355,83 @@ void test_solvers(const std::vector<std::string> &solvers, const int iters, cons
}
}

void test_solvers_gradient_fd(const bool full_fd)
{
std::unique_ptr<TestProblem> problem = std::make_unique<QuadraticProblem>();
std::string solver_name = "L-BFGS";

json solver_params, linear_solver_params;
solver_params["line_search"] = {};
solver_params["max_iterations"] = 100;
if (full_fd)
solver_params["advanced"]["apply_gradient_fd"] = "FullFiniteDiff";
else
solver_params["advanced"]["apply_gradient_fd"] = "DirectionalDerivative";
solver_params["solver"] = solver_name;

const double characteristic_length = 1;

static std::shared_ptr<spdlog::logger> logger = spdlog::stdout_color_mt("test_logger");
logger->set_level(spdlog::level::info);
TestProblem::TVector g;
linear_solver_params["solver"] = "Eigen::LDLT";

for (const auto &ls : line_search::LineSearch::available_methods())
{
solver_params["line_search"]["method"] = ls;

TestProblem::TVector x(problem->size());
x.setZero();

for (int i = 0; i < N_RANDOM; ++i)
{
auto solver = Solver::create(solver_params,
linear_solver_params,
characteristic_length,
*logger);
REQUIRE(solver->name() == solver_name);

try
{
solver->minimize(*problem, x);

double err = std::numeric_limits<double>::max();
for (auto sol : problem->solutions())
err = std::min(err, (x - sol).norm());
if (err >= 1e-7)
{
problem->gradient(x, g);
err = g.norm();
}
INFO("solver: " + solver_name + " LS: " + ls + " problem " + problem->name());
CHECK(err < 1e-7);
if (err >= 1e-7)
break;
}
catch (const std::exception &)
{
break;
}

x.setRandom();
x += problem->min();
x.array() *= (problem->max() - problem->min()).array();
}
}
}

TEST_CASE("non-linear", "[solver]")
{
test_solvers(Solver::available_solvers(), 1000, false);
// test_solvers({"L-BFGS"}, 1000, false);
}

TEST_CASE("non-linear-gradient-fd", "[solver]")
{
test_solvers_gradient_fd(false);
test_solvers_gradient_fd(true);
}

TEST_CASE("non-linear-easier", "[solver]")
{
test_solvers(Solver::available_solvers(), 5000, true);
Expand Down Expand Up @@ -395,7 +470,7 @@ TEST_CASE("non-linear-box-constraint", "[solver]")
if (ls == "None" && solver_name != "MMA")
continue;
if (solver_name == "MMA" && ls != "None")
continue;
continue;
solver_params["line_search"]["method"] = ls;

auto solver = BoxConstraintSolver::create(solver_params,
Expand Down Expand Up @@ -461,7 +536,7 @@ TEST_CASE("non-linear-box-constraint-input", "[solver]")
if (ls == "None" && solver_name != "MMA")
continue;
if (solver_name == "MMA" && ls != "None")
continue;
continue;
solver_params["line_search"]["method"] = ls;

QuadraticProblem::TVector x(prob->size());
Expand Down Expand Up @@ -529,12 +604,12 @@ TEST_CASE("MMA", "[solver]")
solver_params["line_search"]["method"] = "None";

auto solver = BoxConstraintSolver::create(solver_params,
linear_solver_params,
characteristic_length,
*logger);
linear_solver_params,
characteristic_length,
*logger);

auto c = std::make_shared<InequalityConstraint>(solver_params["box_constraints"]["bounds"][1]);
dynamic_cast<BoxConstraintSolver&>(*solver).add_constraint(c);
dynamic_cast<BoxConstraintSolver &>(*solver).add_constraint(c);

QuadraticProblem::TVector x(prob->size());
x.setConstant(3);
Expand Down

0 comments on commit 8a8c103

Please sign in to comment.