From 32cb9e4a12f19a10c6987e6f9769030f0f2c824c Mon Sep 17 00:00:00 2001 From: Sameer Agarwal Date: Sat, 7 Jul 2018 17:11:38 -0700 Subject: [PATCH] Respect bounds when using Solver::Options::check_gradients When Solver::Options::check_gradients is true, Ceres internally creates a new ProblemImpl object which wraps each CostFunction in the user's problem with a GradientCheckingCostFunction. Doing this also requires creating new ParameterBlock objects, and when support for upper and lower bounds was added to Ceres, CreateGradientCheckingProblemImpl should also have been updated to create a problem with the same parameter bounds. As a result, if check_gradients is enabled for a bounded problem, it constructs an unconstrained problem and solves it. This CL fixes this, by introducing Problem::GetParameterLowerBound, and Problem::GetParameterUpperBound and using them to create a bounded problem when checking gradients. Thanks to @pbeeson for not only reporting this problem, but also providing a small standalone reproduction which made debugging this possible. https://github.com/ceres-solver/ceres-solver/issues/379 Change-Id: Id18eb858a7009bf4fa452a21b925922d13f3249f --- docs/source/nnls_modeling.rst | 18 ++++++- include/ceres/problem.h | 9 +++- .../ceres/gradient_checking_cost_function.cc | 11 ++++ .../gradient_checking_cost_function_test.cc | 33 ++++++++++++ internal/ceres/parameter_block.h | 25 ++++++++- internal/ceres/problem.cc | 8 +++ internal/ceres/problem_impl.cc | 22 ++++++++ internal/ceres/problem_impl.h | 2 + internal/ceres/problem_test.cc | 54 +++++++++++++++++++ internal/ceres/solver.cc | 14 ++--- 10 files changed, 184 insertions(+), 12 deletions(-) diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst index 1aee32d3c..f451b4856 100644 --- a/docs/source/nnls_modeling.rst +++ b/docs/source/nnls_modeling.rst @@ -1661,13 +1661,27 @@ Instances Set the lower bound for the parameter at position `index` in the parameter block corresponding to `values`. By default the lower - bound is :math:`-\infty`. + bound is ``-std::numeric_limits::max()``, which is treated + by the solver as the same as :math:`-\infty`. .. function:: void Problem::SetParameterUpperBound(double* values, int index, double upper_bound) Set the upper bound for the parameter at position `index` in the parameter block corresponding to `values`. By default the value is - :math:`\infty`. + ``std::numeric_limits::max()``, which is treated by the + solver as the same as :math:`\infty`. + +.. function:: double Problem::GetParameterLowerBound(double* values, int index) + + Get the lower bound for the parameter with position `index`. If the + parameter is not bounded by the user, then its lower bound is + ``-std::numeric_limits::max()``. + +.. function:: double Problem::GetParameterUpperBound(double* values, int index) + + Get the upper bound for the parameter with position `index`. If the + parameter is not bounded by the user, then its upper bound is + ``std::numeric_limits::max()``. .. function:: int Problem::NumParameterBlocks() const diff --git a/include/ceres/problem.h b/include/ceres/problem.h index e7b60f02b..43b27c89d 100644 --- a/include/ceres/problem.h +++ b/include/ceres/problem.h @@ -327,10 +327,17 @@ class CERES_EXPORT Problem { // associated then NULL is returned. const LocalParameterization* GetParameterization(double* values) const; - // Set the lower/upper bound for the parameter with position "index". + // Set the lower/upper bound for the parameter at position "index". void SetParameterLowerBound(double* values, int index, double lower_bound); void SetParameterUpperBound(double* values, int index, double upper_bound); + // Get the lower/upper bound for the parameter at position + // "index". If the parameter is not bounded by the user, then its + // lower bound is -std::numeric_limits::max() and upper + // bound is std::numeric_limits::max(). + double GetParameterLowerBound(double* values, int index) const; + double GetParameterUpperBound(double* values, int index) const; + // Number of parameter blocks in the problem. Always equals // parameter_blocks().size() and parameter_block_sizes().size(). int NumParameterBlocks() const; diff --git a/internal/ceres/gradient_checking_cost_function.cc b/internal/ceres/gradient_checking_cost_function.cc index 2336ffe9e..6e2a82cdb 100644 --- a/internal/ceres/gradient_checking_cost_function.cc +++ b/internal/ceres/gradient_checking_cost_function.cc @@ -211,6 +211,17 @@ ProblemImpl* CreateGradientCheckingProblemImpl( gradient_checking_problem_impl->SetParameterBlockConstant( parameter_block->mutable_user_state()); } + + for (int i = 0; i < parameter_block->Size(); ++i) { + gradient_checking_problem_impl->SetParameterUpperBound( + parameter_block->mutable_user_state(), + i, + parameter_block->UpperBound(i)); + gradient_checking_problem_impl->SetParameterLowerBound( + parameter_block->mutable_user_state(), + i, + parameter_block->LowerBound(i)); + } } // For every ResidualBlock in problem_impl, create a new diff --git a/internal/ceres/gradient_checking_cost_function_test.cc b/internal/ceres/gradient_checking_cost_function_test.cc index 6ae77f62f..fc8e07197 100644 --- a/internal/ceres/gradient_checking_cost_function_test.cc +++ b/internal/ceres/gradient_checking_cost_function_test.cc @@ -411,5 +411,38 @@ TEST(GradientCheckingProblemImpl, ProblemDimensionsMatch) { } } + +TEST(GradientCheckingProblemImpl, ConstrainedProblemBoundsArePropagated) { + // Parameter blocks with arbitrarily chosen initial values. + double x[] = {1.0, 2.0, 3.0}; + ProblemImpl problem_impl; + problem_impl.AddParameterBlock(x, 3); + problem_impl.AddResidualBlock(new UnaryCostFunction(2, 3), NULL, x); + problem_impl.SetParameterLowerBound(x,0,0.9); + problem_impl.SetParameterUpperBound(x,1,2.5); + + GradientCheckingIterationCallback callback; + std::unique_ptr gradient_checking_problem_impl( + CreateGradientCheckingProblemImpl(&problem_impl, 1.0, 1.0, &callback)); + + // The dimensions of the two problems match. + EXPECT_EQ(problem_impl.NumParameterBlocks(), + gradient_checking_problem_impl->NumParameterBlocks()); + EXPECT_EQ(problem_impl.NumResidualBlocks(), + gradient_checking_problem_impl->NumResidualBlocks()); + + EXPECT_EQ(problem_impl.NumParameters(), + gradient_checking_problem_impl->NumParameters()); + EXPECT_EQ(problem_impl.NumResiduals(), + gradient_checking_problem_impl->NumResiduals()); + + for (int i = 0; i < 3; ++i) { + EXPECT_EQ(problem_impl.GetParameterLowerBound(x, i), + gradient_checking_problem_impl->GetParameterLowerBound(x, i)); + EXPECT_EQ(problem_impl.GetParameterUpperBound(x, i), + gradient_checking_problem_impl->GetParameterUpperBound(x, i)); + } +} + } // namespace internal } // namespace ceres diff --git a/internal/ceres/parameter_block.h b/internal/ceres/parameter_block.h index 9149a5488..c07edef0c 100644 --- a/internal/ceres/parameter_block.h +++ b/internal/ceres/parameter_block.h @@ -128,6 +128,19 @@ class ParameterBlock { void SetVarying() { is_constant_ = false; } bool IsConstant() const { return is_constant_; } + double UpperBound(int index) const { + return (upper_bounds_ ? upper_bounds_[index] + : std::numeric_limits::max()); + } + + double LowerBound(int index) const { + return (lower_bounds_ ? lower_bounds_[index] + : -std::numeric_limits::max()); + } + + bool IsUpperBounded() const { return (upper_bounds_ == nullptr); } + bool IsLowerBounded() const { return (lower_bounds_ == nullptr); } + // This parameter block's index in an array. int index() const { return index_; } void set_index(int index) { index_ = index; } @@ -194,7 +207,11 @@ class ParameterBlock { void SetUpperBound(int index, double upper_bound) { CHECK_LT(index, size_); - if (upper_bounds_.get() == NULL) { + if (upper_bound >= std::numeric_limits::max() && !upper_bounds_) { + return; + } + + if (!upper_bounds_) { upper_bounds_.reset(new double[size_]); std::fill(upper_bounds_.get(), upper_bounds_.get() + size_, @@ -207,7 +224,11 @@ class ParameterBlock { void SetLowerBound(int index, double lower_bound) { CHECK_LT(index, size_); - if (lower_bounds_.get() == NULL) { + if (lower_bound <= -std::numeric_limits::max() && !lower_bounds_) { + return; + } + + if (!lower_bounds_) { lower_bounds_.reset(new double[size_]); std::fill(lower_bounds_.get(), lower_bounds_.get() + size_, diff --git a/internal/ceres/problem.cc b/internal/ceres/problem.cc index 730ce6420..956afa1ed 100644 --- a/internal/ceres/problem.cc +++ b/internal/ceres/problem.cc @@ -201,6 +201,14 @@ void Problem::SetParameterUpperBound(double* values, problem_impl_->SetParameterUpperBound(values, index, upper_bound); } +double Problem::GetParameterUpperBound(double* values, int index) const { + return problem_impl_->GetParameterUpperBound(values, index); +} + +double Problem::GetParameterLowerBound(double* values, int index) const { + return problem_impl_->GetParameterLowerBound(values, index); +} + bool Problem::Evaluate(const EvaluateOptions& evaluate_options, double* cost, vector* residuals, diff --git a/internal/ceres/problem_impl.cc b/internal/ceres/problem_impl.cc index 3fdbfedf2..f627659f2 100644 --- a/internal/ceres/problem_impl.cc +++ b/internal/ceres/problem_impl.cc @@ -717,6 +717,28 @@ void ProblemImpl::SetParameterUpperBound(double* values, parameter_block->SetUpperBound(index, upper_bound); } +double ProblemImpl::GetParameterLowerBound(double* values, int index) const { + ParameterBlock* parameter_block = + FindWithDefault(parameter_block_map_, values, NULL); + if (parameter_block == NULL) { + LOG(FATAL) << "Parameter block not found: " << values + << ". You must add the parameter block to the problem before " + << "you can get the lower bound on one of its components."; + } + return parameter_block->LowerBound(index); +} + +double ProblemImpl::GetParameterUpperBound(double* values, int index) const { + ParameterBlock* parameter_block = + FindWithDefault(parameter_block_map_, values, NULL); + if (parameter_block == NULL) { + LOG(FATAL) << "Parameter block not found: " << values + << ". You must add the parameter block to the problem before " + << "you can set an upper bound on one of its components."; + } + return parameter_block->UpperBound(index); +} + bool ProblemImpl::Evaluate(const Problem::EvaluateOptions& evaluate_options, double* cost, vector* residuals, diff --git a/internal/ceres/problem_impl.h b/internal/ceres/problem_impl.h index 88df82cc9..ff89e947c 100644 --- a/internal/ceres/problem_impl.h +++ b/internal/ceres/problem_impl.h @@ -140,6 +140,8 @@ class ProblemImpl { void SetParameterLowerBound(double* values, int index, double lower_bound); void SetParameterUpperBound(double* values, int index, double upper_bound); + double GetParameterLowerBound(double* values, int index) const; + double GetParameterUpperBound(double* values, int index) const; bool Evaluate(const Problem::EvaluateOptions& options, double* cost, diff --git a/internal/ceres/problem_test.cc b/internal/ceres/problem_test.cc index 3473c59a8..80125834f 100644 --- a/internal/ceres/problem_test.cc +++ b/internal/ceres/problem_test.cc @@ -1527,5 +1527,59 @@ TEST_F(ProblemEvaluateTest, LocalParameterization) { CheckAllEvaluationCombinations(Problem::EvaluateOptions(), expected); } +TEST(Problem, SetAndGetParameterLowerBound) { + Problem problem; + double x[] = {1.0, 2.0}; + problem.AddParameterBlock(x, 2); + + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), + -std::numeric_limits::max()); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits::max()); + + problem.SetParameterLowerBound(x, 0, -1.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), -1.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits::max()); + + problem.SetParameterLowerBound(x, 0, -2.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), -2.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits::max()); + + problem.SetParameterLowerBound(x, 0, -std::numeric_limits::max()); + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), + -std::numeric_limits::max()); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits::max()); +} + +TEST(Problem, SetAndGetParameterUpperBound) { + Problem problem; + double x[] = {1.0, 2.0}; + problem.AddParameterBlock(x, 2); + + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), + std::numeric_limits::max()); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits::max()); + + problem.SetParameterUpperBound(x, 0, -1.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), -1.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits::max()); + + problem.SetParameterUpperBound(x, 0, -2.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), -2.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits::max()); + + problem.SetParameterUpperBound(x, 0, std::numeric_limits::max()); + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), + std::numeric_limits::max()); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits::max()); +} + } // namespace internal } // namespace ceres diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index d38982a0d..7ac60ac99 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc @@ -495,13 +495,6 @@ void Solver::Solve(const Solver::Options& options, Program* program = problem_impl->mutable_program(); PreSolveSummarize(options, problem_impl, summary); - // The main thread also does work so we only need to launch num_threads - 1. - problem_impl->context()->EnsureMinimumThreads(options.num_threads - 1); - - // Make sure that all the parameter blocks states are set to the - // values provided by the user. - program->SetParameterBlockStatePtrsToUserStatePtrs(); - // If gradient_checking is enabled, wrap all cost functions in a // gradient checker and install a callback that terminates if any gradient // error is detected. @@ -520,6 +513,13 @@ void Solver::Solve(const Solver::Options& options, program = problem_impl->mutable_program(); } + // Make sure that all the parameter blocks states are set to the + // values provided by the user. + program->SetParameterBlockStatePtrsToUserStatePtrs(); + + // The main thread also does work so we only need to launch num_threads - 1. + problem_impl->context()->EnsureMinimumThreads(options.num_threads - 1); + std::unique_ptr preprocessor( Preprocessor::Create(modified_options.minimizer_type)); PreprocessedProblem pp;