From c0de4de661c8763d25ce6f5a013ad0a1a1fe0f87 Mon Sep 17 00:00:00 2001 From: Bogdan Burlacu Date: Wed, 27 Nov 2024 13:05:00 +0100 Subject: [PATCH] fix #41 (issue with ceres optimizer) --- .../optimizer/dynamic_cost_function.hpp | 4 +- include/operon/optimizer/optimizer.hpp | 81 ++++++++----- .../operon/optimizer/tiny_cost_function.hpp | 110 ------------------ 3 files changed, 54 insertions(+), 141 deletions(-) delete mode 100644 include/operon/optimizer/tiny_cost_function.hpp diff --git a/include/operon/optimizer/dynamic_cost_function.hpp b/include/operon/optimizer/dynamic_cost_function.hpp index 2b449336..fbe2b979 100644 --- a/include/operon/optimizer/dynamic_cost_function.hpp +++ b/include/operon/optimizer/dynamic_cost_function.hpp @@ -5,11 +5,13 @@ #define OPERON_OPTIMIZER_COST_FUNCTION_HPP #ifdef HAVE_CERES +#include +#include #ifndef CERES_EXPORT #define CERES_EXPORT OPERON_EXPORT #endif -#include + #include "operon/core/contracts.hpp" namespace Operon { diff --git a/include/operon/optimizer/optimizer.hpp b/include/operon/optimizer/optimizer.hpp index cbe5e768..8f4c856d 100644 --- a/include/operon/optimizer/optimizer.hpp +++ b/include/operon/optimizer/optimizer.hpp @@ -151,7 +151,7 @@ struct LevenbergMarquardtOptimizer final : public Operon::Interpreter interpreter{dtable, dataset, tree}; Operon::LMCostFunction cf{interpreter, target, range}; Eigen::LevenbergMarquardt lm(cf); - lm.setMaxfev(static_cast(iterations+1)); + lm.setMaxfev(static_cast(iterations)); auto x0 = tree.GetCoefficients(); OptimizerSummary summary; @@ -195,35 +195,39 @@ struct LevenbergMarquardtOptimizer final : public }; #if defined(HAVE_CERES) -template -struct NonlinearLeastSquaresOptimizer : public OptimizerBase { - explicit NonlinearLeastSquaresOptimizer(InterpreterBase& interpreter) - : OptimizerBase{interpreter} +template +struct LevenbergMarquardtOptimizer final : public OptimizerBase { + explicit LevenbergMarquardtOptimizer(DTable const& dtable, Problem const& problem) + : OptimizerBase{problem}, dtable_{dtable} { } - auto Optimize(Operon::Span target, Range range, size_t iterations, OptimizerSummary& summary) -> std::vector final + [[nodiscard]] auto Optimize(Operon::RandomGenerator& /*unused*/, Operon::Tree const& tree) const -> OptimizerSummary final { - auto const& tree = this->GetTree(); - auto const& ds = this->GetDataset(); - auto const& dt = this->GetDispatchTable(); - - auto x0 = tree.GetCoefficients(); + auto const& dtable = this->GetDispatchTable(); + auto const& problem = this->GetProblem(); + auto const& dataset = problem.GetDataset(); + auto range = problem.TrainingRange(); + auto target = problem.TargetValues(range); + auto iterations = this->Iterations();; - Operon::CostFunction cf(tree, ds, target, range, dt); - auto costFunction = new Operon::DynamicCostFunction(cf); // NOLINT + auto initialParameters = tree.GetCoefficients(); + auto finalParameters = initialParameters; + Operon::Interpreter interpreter{dtable, dataset, tree}; + Operon::LMCostFunction cf{interpreter, target, range}; + auto* dynamicCostFunction = new Operon::DynamicCostFunction{cf}; ceres::Solver::Summary s; - if (!x0.empty()) { - Eigen::Map> m0(x0.data(), std::ssize(x0)); - auto sz = static_cast(x0.size()); - Eigen::VectorXd params = Eigen::Map>(x0.data(), sz).template cast(); + if (!initialParameters.empty()) { + auto sz = std::ssize(finalParameters); + Eigen::Map> m0(finalParameters.data(), sz); + Eigen::VectorXd params = m0.template cast(); ceres::Problem problem; - problem.AddResidualBlock(costFunction, nullptr, params.data()); + problem.AddResidualBlock(dynamicCostFunction, nullptr, params.data()); ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; options.logging_type = ceres::LoggingType::SILENT; - options.max_num_iterations = static_cast(iterations - 1); // workaround since for some reason ceres sometimes does 1 more iteration + options.max_num_iterations = static_cast(iterations); options.minimizer_progress_to_stdout = false; options.num_threads = 1; options.trust_region_strategy_type = ceres::LEVENBERG_MARQUARDT; @@ -231,14 +235,31 @@ struct NonlinearLeastSquaresOptimizer : public Optimize Solve(options, &problem, &s); m0 = params.cast(); } - summary.InitialCost = s.initial_cost; - summary.FinalCost = s.final_cost; - summary.Iterations = static_cast(s.iterations.size()); - summary.FunctionEvaluations = s.num_residual_evaluations; - summary.JacobianEvaluations = s.num_jacobian_evaluations; - summary.Success = detail::CheckSuccess(summary.InitialCost, summary.FinalCost); - return x0; + return Operon::OptimizerSummary { + .InitialParameters = initialParameters, + .FinalParameters = finalParameters, + .InitialCost = static_cast(s.initial_cost), + .FinalCost = static_cast(s.final_cost), + .Iterations = static_cast(s.iterations.size()), + .FunctionEvaluations = s.num_residual_evaluations, + .JacobianEvaluations = s.num_jacobian_evaluations, + .Success = detail::CheckSuccess(s.initial_cost, s.final_cost) + }; } + + auto GetDispatchTable() const -> DTable const& { return dtable_.get(); } + + [[nodiscard]] auto ComputeLikelihood(Operon::Span x, Operon::Span y, Operon::Span w) const -> Operon::Scalar final + { + return GaussianLikelihood::ComputeLikelihood(x, y, w); + } + + [[nodiscard]] auto ComputeFisherMatrix(Operon::Span pred, Operon::Span jac, Operon::Span sigma) const -> Eigen::Matrix final { + return GaussianLikelihood::ComputeFisherMatrix(pred, jac, sigma); + } + + private: + std::reference_wrapper dtable_; }; #endif @@ -298,12 +319,12 @@ struct LBFGSOptimizer final : public OptimizerBase { auto GetDispatchTable() const -> DTable const& { return dtable_.get(); } - [[nodiscard]] virtual auto ComputeLikelihood(Operon::Span x, Operon::Span y, Operon::Span w) const -> Operon::Scalar + [[nodiscard]] auto ComputeLikelihood(Operon::Span x, Operon::Span y, Operon::Span w) const -> Operon::Scalar override { return LossFunction::ComputeLikelihood(x, y, w); } - [[nodiscard]] virtual auto ComputeFisherMatrix(Operon::Span pred, Operon::Span jac, Operon::Span sigma) const -> Eigen::Matrix final { + [[nodiscard]] auto ComputeFisherMatrix(Operon::Span pred, Operon::Span jac, Operon::Span sigma) const -> Eigen::Matrix final { return LossFunction::ComputeFisherMatrix(pred, jac, sigma); } @@ -371,12 +392,12 @@ struct SGDOptimizer final : public OptimizerBase { return summary; } - [[nodiscard]] virtual auto ComputeLikelihood(Operon::Span x, Operon::Span y, Operon::Span w) const -> Operon::Scalar + [[nodiscard]] auto ComputeLikelihood(Operon::Span x, Operon::Span y, Operon::Span w) const -> Operon::Scalar override { return LossFunction::ComputeLikelihood(x, y, w); } - [[nodiscard]] virtual auto ComputeFisherMatrix(Operon::Span pred, Operon::Span jac, Operon::Span sigma) const -> Eigen::Matrix final { + [[nodiscard]] auto ComputeFisherMatrix(Operon::Span pred, Operon::Span jac, Operon::Span sigma) const -> Eigen::Matrix final { return LossFunction::ComputeFisherMatrix(pred, jac, sigma); } diff --git a/include/operon/optimizer/tiny_cost_function.hpp b/include/operon/optimizer/tiny_cost_function.hpp deleted file mode 100644 index 283a9a1f..00000000 --- a/include/operon/optimizer/tiny_cost_function.hpp +++ /dev/null @@ -1,110 +0,0 @@ -// SPDX-License-Identifier: MIT -// SPDX-FileCopyrightText: Copyright 2019-2023 Heal Research - -#ifndef OPERON_OPTIMIZER_TINY_HPP -#define OPERON_OPTIMIZER_TINY_HPP - -#include -#include "operon/interpreter/interpreter.hpp" - -namespace Operon { - -// this cost function is adapted to work with both solvers from Ceres: the normal one and the tiny solver -// for this, a number of template parameters are necessary: -// - the AutodiffCalculator will compute and return the Jacobian matrix -// - the StorageOrder specifies the format of the jacobian (row-major for the big Ceres solver, column-major for the tiny solver) - -template -struct CostFunction { - static auto constexpr Storage{ StorageOrder }; - using Scalar = Operon::Scalar; - - enum { - NUM_RESIDUALS = Eigen::Dynamic, // NOLINT - NUM_PARAMETERS = Eigen::Dynamic, // NOLINT - }; - - explicit CostFunction(Operon::Tree const& tree, Operon::Dataset const& dataset, Operon::Span target, Operon::Range const range, DTable const& dtable) - : tree_{tree} - , dataset_{dataset} - , target_{target} - , range_{range} - , dtable_{dtable} - , numResiduals_{range.Size()} - , numParameters_{ParameterCount(tree)} - { } - - inline auto Evaluate(Scalar const* parameters, Scalar* residuals, Scalar* jacobian) const -> bool // NOLINT - { - EXPECT(target_.size() == numResiduals_); - EXPECT(parameters != nullptr); - Operon::Span params{ parameters, numParameters_ }; - - Interpreter interpreter{dtable_.get(), dataset_, tree_}; - - if (jacobian != nullptr) { - Operon::Span jac{jacobian, static_cast(numResiduals_ * numParameters_)}; - interpreter.JacRev(params, range_, jac); - } - - if (residuals != nullptr) { - Operon::Span res{ residuals, static_cast(numResiduals_) }; - interpreter.Evaluate(params, range_, res); - Eigen::Map> x(residuals, numResiduals_); - Eigen::Map const> y(target_.data(), numResiduals_); - x -= y; - } - return true; - } - - // ceres solver - jacobian must be in row-major format - // tiny solver - jacobian must be in col-major format - auto operator()(Scalar const* parameters, Scalar* residuals, Scalar* jacobian) const -> bool - { - return Evaluate(parameters, residuals, jacobian); - } - - [[nodiscard]] auto NumResiduals() const -> int { return numResiduals_; } - [[nodiscard]] auto NumParameters() const -> int { return numParameters_; } - - // required by Eigen::LevenbergMarquardt - using JacobianType = Eigen::Matrix; - using QRSolver = Eigen::ColPivHouseholderQR; - - // there is no real documentation but looking at Eigen unit tests, these functions should return zero - // see: https://gitlab.com/libeigen/eigen/-/blob/master/unsupported/test/NonLinearOptimization.cpp - auto operator()(Eigen::Matrix const& input, Eigen::Matrix& residual) const -> int - { - Evaluate(input.data(), residual.data(), nullptr); - return 0; - } - - auto df(Eigen::Matrix const& input, Eigen::Matrix& jacobian) const -> int // NOLINT - { - static_assert(StorageOrder == Eigen::ColMajor, "Eigen::LevenbergMarquardt requires the Jacobian to be stored in column-major format."); - Evaluate(input.data(), nullptr, jacobian.data()); - //std::cout << "jacobian:\n" << jacobian << "\n"; - return 0; - } - - [[nodiscard]] auto values() const -> int { return NumResiduals(); } // NOLINT - [[nodiscard]] auto inputs() const -> int { return NumParameters(); } // NOLINT - -private: - std::reference_wrapper tree_; - std::reference_wrapper dataset_; - std::reference_wrapper dtable_; - - Operon::Range const range_; // NOLINT - Operon::Span target_; - std::size_t numResiduals_; - std::size_t numParameters_; - - inline auto ParameterCount(auto const& tree) const -> std::size_t { - auto const& nodes = tree.Nodes(); - return std::count_if(nodes.cbegin(), nodes.cend(), [](auto const& n) { return n.Optimize; }); - } -}; -} // namespace Operon - -#endif