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
diff --git a/docs/source/nnls_modeling.rst b/docs/source/nnls_modeling.rst index 1aee32d..f451b48 100644 --- a/docs/source/nnls_modeling.rst +++ b/docs/source/nnls_modeling.rst
@@ -1661,13 +1661,27 @@ 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<double>::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<double>::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<double>::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<double>::max()``. .. function:: int Problem::NumParameterBlocks() const
diff --git a/include/ceres/problem.h b/include/ceres/problem.h index e7b60f0..43b27c8 100644 --- a/include/ceres/problem.h +++ b/include/ceres/problem.h
@@ -327,10 +327,17 @@ // 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<double>::max() and upper + // bound is std::numeric_limits<double>::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 2336ffe..6e2a82c 100644 --- a/internal/ceres/gradient_checking_cost_function.cc +++ b/internal/ceres/gradient_checking_cost_function.cc
@@ -211,6 +211,17 @@ 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 6ae77f6..fc8e071 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, 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<ProblemImpl> 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 9149a54..c07edef 100644 --- a/internal/ceres/parameter_block.h +++ b/internal/ceres/parameter_block.h
@@ -128,6 +128,19 @@ 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<double>::max()); + } + + double LowerBound(int index) const { + return (lower_bounds_ ? lower_bounds_[index] + : -std::numeric_limits<double>::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 @@ void SetUpperBound(int index, double upper_bound) { CHECK_LT(index, size_); - if (upper_bounds_.get() == NULL) { + if (upper_bound >= std::numeric_limits<double>::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 @@ void SetLowerBound(int index, double lower_bound) { CHECK_LT(index, size_); - if (lower_bounds_.get() == NULL) { + if (lower_bound <= -std::numeric_limits<double>::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 730ce64..956afa1 100644 --- a/internal/ceres/problem.cc +++ b/internal/ceres/problem.cc
@@ -201,6 +201,14 @@ 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<double>* residuals,
diff --git a/internal/ceres/problem_impl.cc b/internal/ceres/problem_impl.cc index 3fdbfed..f627659 100644 --- a/internal/ceres/problem_impl.cc +++ b/internal/ceres/problem_impl.cc
@@ -717,6 +717,28 @@ 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<double>* residuals,
diff --git a/internal/ceres/problem_impl.h b/internal/ceres/problem_impl.h index 88df82c..ff89e94 100644 --- a/internal/ceres/problem_impl.h +++ b/internal/ceres/problem_impl.h
@@ -140,6 +140,8 @@ 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 3473c59..8012583 100644 --- a/internal/ceres/problem_test.cc +++ b/internal/ceres/problem_test.cc
@@ -1527,5 +1527,59 @@ 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<double>::max()); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits<double>::max()); + + problem.SetParameterLowerBound(x, 0, -1.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), -1.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits<double>::max()); + + problem.SetParameterLowerBound(x, 0, -2.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), -2.0); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits<double>::max()); + + problem.SetParameterLowerBound(x, 0, -std::numeric_limits<double>::max()); + EXPECT_EQ(problem.GetParameterLowerBound(x, 0), + -std::numeric_limits<double>::max()); + EXPECT_EQ(problem.GetParameterLowerBound(x, 1), + -std::numeric_limits<double>::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<double>::max()); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits<double>::max()); + + problem.SetParameterUpperBound(x, 0, -1.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), -1.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits<double>::max()); + + problem.SetParameterUpperBound(x, 0, -2.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), -2.0); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits<double>::max()); + + problem.SetParameterUpperBound(x, 0, std::numeric_limits<double>::max()); + EXPECT_EQ(problem.GetParameterUpperBound(x, 0), + std::numeric_limits<double>::max()); + EXPECT_EQ(problem.GetParameterUpperBound(x, 1), + std::numeric_limits<double>::max()); +} + } // namespace internal } // namespace ceres
diff --git a/internal/ceres/solver.cc b/internal/ceres/solver.cc index d38982a..7ac60ac 100644 --- a/internal/ceres/solver.cc +++ b/internal/ceres/solver.cc
@@ -495,13 +495,6 @@ 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 @@ 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( Preprocessor::Create(modified_options.minimizer_type)); PreprocessedProblem pp;