saberbud/ceres-solver
1
0
Fork 0
mirror of https://github.com/ceres-solver/ceres-solver.git synced 2026-08-30 17:10:38 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
8fd7828e3dd2ac6dccb57985249d40fdb2ff9428
ceres-solver/internal/ceres/solver.cc
T
Dmitriy Korchemkin b7116824b7 Evaluation benchmark 
Benchmark for evaluation of residuals and evaluation of both residuals
and jacobian.

For each input file (in BAL format) specified on the command line two
sets of benchmarks are created for different number of threads.

BAL data is preloaded before starting benchmark.

Usage:
./bin/evaluation_benchmark [flags] input_1.txt ... input_N.txt

Change-Id: I543df65f483e3b186f52785b68bd5e2b3e3552a3
2022-09-21 13:40:03 +03:00

1138 lines
43 KiB
C++
Raw Blame History

// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2022 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
// this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
// this list of conditions and the following disclaimer in the documentation
// and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
// used to endorse or promote products derived from this software without
// specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: keir@google.com (Keir Mierle)
// sameeragarwal@google.com (Sameer Agarwal)
#include "ceres/solver.h"
#include <algorithm>
#include <memory>
#include <sstream> // NOLINT
#include <string>
#include <vector>
#include "ceres/casts.h"
#include "ceres/context.h"
#include "ceres/context_impl.h"
#include "ceres/detect_structure.h"
#include "ceres/eigensparse.h"
#include "ceres/gradient_checking_cost_function.h"
#include "ceres/internal/export.h"
#include "ceres/parameter_block_ordering.h"
#include "ceres/preprocessor.h"
#include "ceres/problem.h"
#include "ceres/problem_impl.h"
#include "ceres/program.h"
#include "ceres/schur_templates.h"
#include "ceres/solver_utils.h"
#include "ceres/stringprintf.h"
#include "ceres/suitesparse.h"
#include "ceres/types.h"
#include "ceres/wall_time.h"
namespace ceres {
namespace {
using internal::StringAppendF;
using internal::StringPrintf;
using std::map;
using std::string;
using std::vector;
#define OPTION_OP(x, y, OP) \
if (!(options.x OP y)) { \
std::stringstream ss; \
ss << "Invalid configuration. "; \
ss << string("Solver::Options::" #x " = ") << options.x << ". "; \
ss << "Violated constraint: "; \
ss << string("Solver::Options::" #x " " #OP " " #y); \
*error = ss.str(); \
return false; \
}
#define OPTION_OP_OPTION(x, y, OP) \
if (!(options.x OP options.y)) { \
std::stringstream ss; \
ss << "Invalid configuration. "; \
ss << string("Solver::Options::" #x " = ") << options.x << ". "; \
ss << string("Solver::Options::" #y " = ") << options.y << ". "; \
ss << "Violated constraint: "; \
ss << string("Solver::Options::" #x); \
ss << string(#OP " Solver::Options::" #y "."); \
*error = ss.str(); \
return false; \
}
#define OPTION_GE(x, y) OPTION_OP(x, y, >=);
#define OPTION_GT(x, y) OPTION_OP(x, y, >);
#define OPTION_LE(x, y) OPTION_OP(x, y, <=);
#define OPTION_LT(x, y) OPTION_OP(x, y, <);
#define OPTION_LE_OPTION(x, y) OPTION_OP_OPTION(x, y, <=)
#define OPTION_LT_OPTION(x, y) OPTION_OP_OPTION(x, y, <)
bool CommonOptionsAreValid(const Solver::Options& options, string* error) {
OPTION_GE(max_num_iterations, 0);
OPTION_GE(max_solver_time_in_seconds, 0.0);
OPTION_GE(function_tolerance, 0.0);
OPTION_GE(gradient_tolerance, 0.0);
OPTION_GE(parameter_tolerance, 0.0);
OPTION_GT(num_threads, 0);
if (options.check_gradients) {
OPTION_GT(gradient_check_relative_precision, 0.0);
OPTION_GT(gradient_check_numeric_derivative_relative_step_size, 0.0);
}
return true;
}
bool IsNestedDissectionAvailable(SparseLinearAlgebraLibraryType type) {
return (((type == SUITE_SPARSE) &&
internal::SuiteSparse::IsNestedDissectionAvailable()) ||
(type == ACCELERATE_SPARSE) ||
((type == EIGEN_SPARSE) &&
internal::EigenSparse::IsNestedDissectionAvailable()));
}
bool IsIterativeSolver(LinearSolverType type) {
return (type == CGNR || type == ITERATIVE_SCHUR);
}
bool OptionsAreValidForDenseSolver(const Solver::Options& options,
string* error) {
const char* library_name = DenseLinearAlgebraLibraryTypeToString(
options.dense_linear_algebra_library_type);
const char* solver_name =
LinearSolverTypeToString(options.linear_solver_type);
constexpr char kFormat[] =
"Can't use %s with dense_linear_algebra_library_type = %s "
"because support not enabled when Ceres was built.";
if (!IsDenseLinearAlgebraLibraryTypeAvailable(
options.dense_linear_algebra_library_type)) {
*error = StringPrintf(kFormat, solver_name, library_name);
return false;
}
return true;
}
bool OptionsAreValidForSparseCholeskyBasedSolver(const Solver::Options& options,
string* error) {
const char* library_name = SparseLinearAlgebraLibraryTypeToString(
options.sparse_linear_algebra_library_type);
// Sparse factorization based solvers and some preconditioners require a
// sparse Cholesky factorization.
const char* solver_name =
IsIterativeSolver(options.linear_solver_type)
? PreconditionerTypeToString(options.preconditioner_type)
: LinearSolverTypeToString(options.linear_solver_type);
constexpr char kNoSparseFormat[] =
"Can't use %s with sparse_linear_algebra_library_type = %s.";
constexpr char kNoLibraryFormat[] =
"Can't use %s sparse_linear_algebra_library_type = %s, because support "
"was not enabled when Ceres Solver was built.";
constexpr char kNoNesdisFormat[] =
"NESDIS is not available with sparse_linear_algebra_library_type = %s.";
constexpr char kMixedFormat[] =
"use_mixed_precision_solves with %s is not supported with "
"sparse_linear_algebra_library_type = %s";
constexpr char kDynamicSparsityFormat[] =
"dynamic sparsity is not supported with "
"sparse_linear_algebra_library_type = %s";
if (options.sparse_linear_algebra_library_type == NO_SPARSE) {
*error = StringPrintf(kNoSparseFormat, solver_name, library_name);
return false;
}
if (!IsSparseLinearAlgebraLibraryTypeAvailable(
options.sparse_linear_algebra_library_type)) {
*error = StringPrintf(kNoLibraryFormat, solver_name, library_name);
return false;
}
if (options.linear_solver_ordering_type == ceres::NESDIS &&
!IsNestedDissectionAvailable(
options.sparse_linear_algebra_library_type)) {
*error = StringPrintf(kNoNesdisFormat, library_name);
return false;
}
if (options.use_mixed_precision_solves &&
options.sparse_linear_algebra_library_type == SUITE_SPARSE) {
*error = StringPrintf(kMixedFormat, solver_name, library_name);
return false;
}
if (options.dynamic_sparsity &&
options.sparse_linear_algebra_library_type == ACCELERATE_SPARSE) {
*error = StringPrintf(kDynamicSparsityFormat, library_name);
return false;
}
return true;
}
bool OptionsAreValidForDenseNormalCholesky(const Solver::Options& options,
string* error) {
CHECK_EQ(options.linear_solver_type, DENSE_NORMAL_CHOLESKY);
return OptionsAreValidForDenseSolver(options, error);
}
bool OptionsAreValidForDenseQr(const Solver::Options& options, string* error) {
CHECK_EQ(options.linear_solver_type, DENSE_QR);
if (!OptionsAreValidForDenseSolver(options, error)) {
return false;
}
if (options.use_mixed_precision_solves) {
*error = "Can't use use_mixed_precision_solves with DENSE_QR.";
return false;
}
return true;
}
bool OptionsAreValidForSparseNormalCholesky(const Solver::Options& options,
string* error) {
CHECK_EQ(options.linear_solver_type, SPARSE_NORMAL_CHOLESKY);
return OptionsAreValidForSparseCholeskyBasedSolver(options, error);
}
bool OptionsAreValidForDenseSchur(const Solver::Options& options,
string* error) {
CHECK_EQ(options.linear_solver_type, DENSE_SCHUR);
if (options.dynamic_sparsity) {
*error = "dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY";
return false;
}
if (!OptionsAreValidForDenseSolver(options, error)) {
return false;
}
return true;
}
bool OptionsAreValidForSparseSchur(const Solver::Options& options,
string* error) {
CHECK_EQ(options.linear_solver_type, SPARSE_SCHUR);
if (options.dynamic_sparsity) {
*error = "Dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY.";
return false;
}
return OptionsAreValidForSparseCholeskyBasedSolver(options, error);
}
bool OptionsAreValidForIterativeSchur(const Solver::Options& options,
string* error) {
CHECK_EQ(options.linear_solver_type, ITERATIVE_SCHUR);
if (options.dynamic_sparsity) {
*error = "Dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY.";
return false;
}
if (options.use_explicit_schur_complement) {
if (options.preconditioner_type != SCHUR_JACOBI) {
*error =
"use_explicit_schur_complement only supports "
"SCHUR_JACOBI as the preconditioner.";
return false;
}
if (options.use_spse_initialization) {
*error =
"use_explicit_schur_complement does not support "
"use_spse_initialization.";
return false;
}
}
if (options.use_spse_initialization ||
options.preconditioner_type == SCHUR_POWER_SERIES_EXPANSION) {
OPTION_GE(max_num_spse_iterations, 1)
OPTION_GE(spse_tolerance, 0.0)
}
if (options.use_mixed_precision_solves) {
*error = "Can't use use_mixed_precision_solves with ITERATIVE_SCHUR";
return false;
}
if (options.dynamic_sparsity) {
*error = "Dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY.";
return false;
}
if (options.preconditioner_type == SUBSET) {
*error = "Can't use SUBSET preconditioner with ITERATIVE_SCHUR";
return false;
}
// CLUSTER_JACOBI and CLUSTER_TRIDIAGONAL require sparse Cholesky
// factorization.
if (options.preconditioner_type == CLUSTER_JACOBI ||
options.preconditioner_type == CLUSTER_TRIDIAGONAL) {
return OptionsAreValidForSparseCholeskyBasedSolver(options, error);
}
return true;
}
bool OptionsAreValidForCgnr(const Solver::Options& options, string* error) {
CHECK_EQ(options.linear_solver_type, CGNR);
if (options.preconditioner_type != IDENTITY &&
options.preconditioner_type != JACOBI &&
options.preconditioner_type != SUBSET) {
*error =
StringPrintf("Can't use CGNR with preconditioner_type = %s.",
PreconditionerTypeToString(options.preconditioner_type));
return false;
}
if (options.use_mixed_precision_solves) {
*error = "use_mixed_precision_solves cannot be used with CGNR";
return false;
}
if (options.dynamic_sparsity) {
*error = "Dynamic sparsity is only supported with SPARSE_NORMAL_CHOLESKY.";
return false;
}
if (options.preconditioner_type == SUBSET) {
if (options.sparse_linear_algebra_library_type == CUDA_SPARSE) {
*error =
"Can't use CGNR with preconditioner_type = SUBSET when "
"sparse_linear_algebra_library_type = CUDA_SPARSE.";
return false;
}
if (options.residual_blocks_for_subset_preconditioner.empty()) {
*error =
"When using SUBSET preconditioner, "
"residual_blocks_for_subset_preconditioner cannot be empty";
return false;
}
// SUBSET preconditioner requires sparse Cholesky factorization.
if (!OptionsAreValidForSparseCholeskyBasedSolver(options, error)) {
return false;
}
}
// Check options for CGNR with CUDA_SPARSE.
if (options.sparse_linear_algebra_library_type == CUDA_SPARSE) {
if (!IsSparseLinearAlgebraLibraryTypeAvailable(CUDA_SPARSE)) {
*error =
"Can't use CGNR with sparse_linear_algebra_library_type = "
"CUDA_SPARSE because support was not enabled when Ceres was built.";
return false;
}
}
return true;
}
bool OptionsAreValidForLinearSolver(const Solver::Options& options,
string* error) {
switch (options.linear_solver_type) {
case DENSE_NORMAL_CHOLESKY:
return OptionsAreValidForDenseNormalCholesky(options, error);
case DENSE_QR:
return OptionsAreValidForDenseQr(options, error);
case SPARSE_NORMAL_CHOLESKY:
return OptionsAreValidForSparseNormalCholesky(options, error);
case DENSE_SCHUR:
return OptionsAreValidForDenseSchur(options, error);
case SPARSE_SCHUR:
return OptionsAreValidForSparseSchur(options, error);
case ITERATIVE_SCHUR:
return OptionsAreValidForIterativeSchur(options, error);
case CGNR:
return OptionsAreValidForCgnr(options, error);
default:
LOG(FATAL) << "Congratulations you have found a bug. Please report "
"this to the "
"Ceres Solver developers. Unknown linear solver type: "
<< LinearSolverTypeToString(options.linear_solver_type);
}
return false;
}
bool TrustRegionOptionsAreValid(const Solver::Options& options, string* error) {
OPTION_GT(initial_trust_region_radius, 0.0);
OPTION_GT(min_trust_region_radius, 0.0);
OPTION_GT(max_trust_region_radius, 0.0);
OPTION_LE_OPTION(min_trust_region_radius, max_trust_region_radius);
OPTION_LE_OPTION(min_trust_region_radius, initial_trust_region_radius);
OPTION_LE_OPTION(initial_trust_region_radius, max_trust_region_radius);
OPTION_GE(min_relative_decrease, 0.0);
OPTION_GE(min_lm_diagonal, 0.0);
OPTION_GE(max_lm_diagonal, 0.0);
OPTION_LE_OPTION(min_lm_diagonal, max_lm_diagonal);
OPTION_GE(max_num_consecutive_invalid_steps, 0);
OPTION_GT(eta, 0.0);
OPTION_GE(min_linear_solver_iterations, 0);
OPTION_GE(max_linear_solver_iterations, 0);
OPTION_LE_OPTION(min_linear_solver_iterations, max_linear_solver_iterations);
if (options.use_inner_iterations) {
OPTION_GE(inner_iteration_tolerance, 0.0);
}
if (options.use_nonmonotonic_steps) {
OPTION_GT(max_consecutive_nonmonotonic_steps, 0);
}
if ((options.trust_region_strategy_type == DOGLEG) &&
IsIterativeSolver(options.linear_solver_type)) {
*error =
"DOGLEG only supports exact factorization based linear "
"solvers. If you want to use an iterative solver please "
"use LEVENBERG_MARQUARDT as the trust_region_strategy_type";
return false;
}
if (!OptionsAreValidForLinearSolver(options, error)) {
return false;
}
if (!options.trust_region_minimizer_iterations_to_dump.empty() &&
options.trust_region_problem_dump_format_type != CONSOLE &&
options.trust_region_problem_dump_directory.empty()) {
*error = "Solver::Options::trust_region_problem_dump_directory is empty.";
return false;
}
return true;
}
bool LineSearchOptionsAreValid(const Solver::Options& options, string* error) {
OPTION_GT(max_lbfgs_rank, 0);
OPTION_GT(min_line_search_step_size, 0.0);
OPTION_GT(max_line_search_step_contraction, 0.0);
OPTION_LT(max_line_search_step_contraction, 1.0);
OPTION_LT_OPTION(max_line_search_step_contraction,
min_line_search_step_contraction);
OPTION_LE(min_line_search_step_contraction, 1.0);
OPTION_GE(max_num_line_search_step_size_iterations,
(options.minimizer_type == ceres::TRUST_REGION ? 0 : 1));
OPTION_GT(line_search_sufficient_function_decrease, 0.0);
OPTION_LT_OPTION(line_search_sufficient_function_decrease,
line_search_sufficient_curvature_decrease);
OPTION_LT(line_search_sufficient_curvature_decrease, 1.0);
OPTION_GT(max_line_search_step_expansion, 1.0);
if ((options.line_search_direction_type == ceres::BFGS ||
options.line_search_direction_type == ceres::LBFGS) &&
options.line_search_type != ceres::WOLFE) {
*error =
string("Invalid configuration: Solver::Options::line_search_type = ") +
string(LineSearchTypeToString(options.line_search_type)) +
string(
". When using (L)BFGS, "
"Solver::Options::line_search_type must be set to WOLFE.");
return false;
}
// Warn user if they have requested BISECTION interpolation, but constraints
// on max/min step size change during line search prevent bisection scaling
// from occurring. Warn only, as this is likely a user mistake, but one
// which does not prevent us from continuing.
if (options.line_search_interpolation_type == ceres::BISECTION &&
(options.max_line_search_step_contraction > 0.5 ||
options.min_line_search_step_contraction < 0.5)) {
LOG(WARNING)
<< "Line search interpolation type is BISECTION, but specified "
<< "max_line_search_step_contraction: "
<< options.max_line_search_step_contraction << ", and "
<< "min_line_search_step_contraction: "
<< options.min_line_search_step_contraction
<< ", prevent bisection (0.5) scaling, continuing with solve "
"regardless.";
}
return true;
}
#undef OPTION_OP
#undef OPTION_OP_OPTION
#undef OPTION_GT
#undef OPTION_GE
#undef OPTION_LE
#undef OPTION_LT
#undef OPTION_LE_OPTION
#undef OPTION_LT_OPTION
void StringifyOrdering(const vector<int>& ordering, string* report) {
if (ordering.empty()) {
internal::StringAppendF(report, "AUTOMATIC");
return;
}
for (int i = 0; i < ordering.size() - 1; ++i) {
internal::StringAppendF(report, "%d,", ordering[i]);
}
internal::StringAppendF(report, "%d", ordering.back());
}
void SummarizeGivenProgram(const internal::Program& program,
Solver::Summary* summary) {
// clang-format off
summary->num_parameter_blocks = program.NumParameterBlocks();
summary->num_parameters = program.NumParameters();
summary->num_effective_parameters = program.NumEffectiveParameters();
summary->num_residual_blocks = program.NumResidualBlocks();
summary->num_residuals = program.NumResiduals();
// clang-format on
}
void SummarizeReducedProgram(const internal::Program& program,
Solver::Summary* summary) {
// clang-format off
summary->num_parameter_blocks_reduced = program.NumParameterBlocks();
summary->num_parameters_reduced = program.NumParameters();
summary->num_effective_parameters_reduced = program.NumEffectiveParameters();
summary->num_residual_blocks_reduced = program.NumResidualBlocks();
summary->num_residuals_reduced = program.NumResiduals();
// clang-format on
}
void PreSolveSummarize(const Solver::Options& options,
const internal::ProblemImpl* problem,
Solver::Summary* summary) {
SummarizeGivenProgram(problem->program(), summary);
internal::OrderingToGroupSizes(options.linear_solver_ordering.get(),
&(summary->linear_solver_ordering_given));
internal::OrderingToGroupSizes(options.inner_iteration_ordering.get(),
&(summary->inner_iteration_ordering_given));
// clang-format off
summary->dense_linear_algebra_library_type = options.dense_linear_algebra_library_type;
summary->dogleg_type = options.dogleg_type;
summary->inner_iteration_time_in_seconds = 0.0;
summary->num_line_search_steps = 0;
summary->line_search_cost_evaluation_time_in_seconds = 0.0;
summary->line_search_gradient_evaluation_time_in_seconds = 0.0;
summary->line_search_polynomial_minimization_time_in_seconds = 0.0;
summary->line_search_total_time_in_seconds = 0.0;
summary->inner_iterations_given = options.use_inner_iterations;
summary->line_search_direction_type = options.line_search_direction_type;
summary->line_search_interpolation_type = options.line_search_interpolation_type;
summary->line_search_type = options.line_search_type;
summary->linear_solver_type_given = options.linear_solver_type;
summary->max_lbfgs_rank = options.max_lbfgs_rank;
summary->minimizer_type = options.minimizer_type;
summary->nonlinear_conjugate_gradient_type = options.nonlinear_conjugate_gradient_type;
summary->num_threads_given = options.num_threads;
summary->preconditioner_type_given = options.preconditioner_type;
summary->sparse_linear_algebra_library_type = options.sparse_linear_algebra_library_type;
summary->linear_solver_ordering_type = options.linear_solver_ordering_type;
summary->trust_region_strategy_type = options.trust_region_strategy_type;
summary->visibility_clustering_type = options.visibility_clustering_type;
// clang-format on
}
void PostSolveSummarize(const internal::PreprocessedProblem& pp,
Solver::Summary* summary) {
internal::OrderingToGroupSizes(pp.options.linear_solver_ordering.get(),
&(summary->linear_solver_ordering_used));
// TODO(sameeragarwal): Update the preprocessor to collapse the
// second and higher groups into one group when nested dissection is
// used.
internal::OrderingToGroupSizes(pp.options.inner_iteration_ordering.get(),
&(summary->inner_iteration_ordering_used));
// clang-format off
summary->inner_iterations_used = pp.inner_iteration_minimizer != nullptr;
summary->linear_solver_type_used = pp.linear_solver_options.type;
summary->mixed_precision_solves_used = pp.options.use_mixed_precision_solves;
summary->num_threads_used = pp.options.num_threads;
summary->preconditioner_type_used = pp.options.preconditioner_type;
// clang-format on
internal::SetSummaryFinalCost(summary);
if (pp.reduced_program != nullptr) {
SummarizeReducedProgram(*pp.reduced_program, summary);
}
using internal::CallStatistics;
// It is possible that no evaluator was created. This would be the
// case if the preprocessor failed, or if the reduced problem did
// not contain any parameter blocks. Thus, only extract the
// evaluator statistics if one exists.
if (pp.evaluator != nullptr) {
const map<string, CallStatistics>& evaluator_statistics =
pp.evaluator->Statistics();
{
const CallStatistics& call_stats = FindWithDefault(
evaluator_statistics, "Evaluator::Residual", CallStatistics());
summary->residual_evaluation_time_in_seconds = call_stats.time;
summary->num_residual_evaluations = call_stats.calls;
}
{
const CallStatistics& call_stats = FindWithDefault(
evaluator_statistics, "Evaluator::Jacobian", CallStatistics());
summary->jacobian_evaluation_time_in_seconds = call_stats.time;
summary->num_jacobian_evaluations = call_stats.calls;
}
}
// Again, like the evaluator, there may or may not be a linear
// solver from which we can extract run time statistics. In
// particular the line search solver does not use a linear solver.
if (pp.linear_solver != nullptr) {
const map<string, CallStatistics>& linear_solver_statistics =
pp.linear_solver->Statistics();
const CallStatistics& call_stats = FindWithDefault(
linear_solver_statistics, "LinearSolver::Solve", CallStatistics());
summary->num_linear_solves = call_stats.calls;
summary->linear_solver_time_in_seconds = call_stats.time;
}
}
void Minimize(internal::PreprocessedProblem* pp, Solver::Summary* summary) {
using internal::Minimizer;
using internal::Program;
Program* program = pp->reduced_program.get();
if (pp->reduced_program->NumParameterBlocks() == 0) {
summary->message =
"Function tolerance reached. "
"No non-constant parameter blocks found.";
summary->termination_type = CONVERGENCE;
if (pp->options.logging_type != SILENT) {
VLOG(1) << summary->message;
}
summary->initial_cost = summary->fixed_cost;
summary->final_cost = summary->fixed_cost;
return;
}
const Vector original_reduced_parameters = pp->reduced_parameters;
auto minimizer = Minimizer::Create(pp->options.minimizer_type);
minimizer->Minimize(
pp->minimizer_options, pp->reduced_parameters.data(), summary);
program->StateVectorToParameterBlocks(
summary->IsSolutionUsable() ? pp->reduced_parameters.data()
: original_reduced_parameters.data());
program->CopyParameterBlockStateToUserState();
}
std::string SchurStructureToString(const int row_block_size,
const int e_block_size,
const int f_block_size) {
const std::string row = (row_block_size == Eigen::Dynamic)
? "d"
: internal::StringPrintf("%d", row_block_size);
const std::string e = (e_block_size == Eigen::Dynamic)
? "d"
: internal::StringPrintf("%d", e_block_size);
const std::string f = (f_block_size == Eigen::Dynamic)
? "d"
: internal::StringPrintf("%d", f_block_size);
return internal::StringPrintf("%s,%s,%s", row.c_str(), e.c_str(), f.c_str());
}
#ifndef CERES_NO_CUDA
bool IsCudaRequired(const Solver::Options& options) {
if (options.linear_solver_type == DENSE_NORMAL_CHOLESKY ||
options.linear_solver_type == DENSE_SCHUR ||
options.linear_solver_type == DENSE_QR) {
return (options.dense_linear_algebra_library_type == CUDA);
}
if (options.linear_solver_type == CGNR) {
return (options.sparse_linear_algebra_library_type == CUDA_SPARSE);
}
return false;
}
#endif
} // namespace
bool Solver::Options::IsValid(string* error) const {
if (!CommonOptionsAreValid(*this, error)) {
return false;
}
if (minimizer_type == TRUST_REGION &&
!TrustRegionOptionsAreValid(*this, error)) {
return false;
}
// We do not know if the problem is bounds constrained or not, if it
// is then the trust region solver will also use the line search
// solver to do a projection onto the box constraints, so make sure
// that the line search options are checked independent of what
// minimizer algorithm is being used.
return LineSearchOptionsAreValid(*this, error);
}
Solver::~Solver() = default;
void Solver::Solve(const Solver::Options& options,
Problem* problem,
Solver::Summary* summary) {
using internal::PreprocessedProblem;
using internal::Preprocessor;
using internal::ProblemImpl;
using internal::Program;
using internal::WallTimeInSeconds;
CHECK(problem != nullptr);
CHECK(summary != nullptr);
double start_time = WallTimeInSeconds();
*summary = Summary();
if (!options.IsValid(&summary->message)) {
LOG(ERROR) << "Terminating: " << summary->message;
return;
}
ProblemImpl* problem_impl = problem->mutable_impl();
Program* program = problem_impl->mutable_program();
PreSolveSummarize(options, problem_impl, summary);
#ifndef CERES_NO_CUDA
if (IsCudaRequired(options)) {
if (!problem_impl->context()->InitCuda(&summary->message)) {
LOG(ERROR) << "Terminating: " << summary->message;
return;
}
}
#endif // CERES_NO_CUDA
// 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.
std::unique_ptr<internal::ProblemImpl> gradient_checking_problem;
internal::GradientCheckingIterationCallback gradient_checking_callback;
Solver::Options modified_options = options;
if (options.check_gradients) {
modified_options.callbacks.push_back(&gradient_checking_callback);
gradient_checking_problem = CreateGradientCheckingProblemImpl(
problem_impl,
options.gradient_check_numeric_derivative_relative_step_size,
options.gradient_check_relative_precision,
&gradient_checking_callback);
problem_impl = gradient_checking_problem.get();
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);
auto preprocessor = Preprocessor::Create(modified_options.minimizer_type);
PreprocessedProblem pp;
const bool status =
preprocessor->Preprocess(modified_options, problem_impl, &pp);
// We check the linear_solver_options.type rather than
// modified_options.linear_solver_type because, depending on the
// lack of a Schur structure, the preprocessor may change the linear
// solver type.
if (IsSchurType(pp.linear_solver_options.type)) {
// TODO(sameeragarwal): We can likely eliminate the duplicate call
// to DetectStructure here and inside the linear solver, by
// calling this in the preprocessor.
int row_block_size;
int e_block_size;
int f_block_size;
DetectStructure(*static_cast<internal::BlockSparseMatrix*>(
pp.minimizer_options.jacobian.get())
->block_structure(),
pp.linear_solver_options.elimination_groups[0],
&row_block_size,
&e_block_size,
&f_block_size);
summary->schur_structure_given =
SchurStructureToString(row_block_size, e_block_size, f_block_size);
internal::GetBestSchurTemplateSpecialization(
&row_block_size, &e_block_size, &f_block_size);
summary->schur_structure_used =
SchurStructureToString(row_block_size, e_block_size, f_block_size);
}
summary->fixed_cost = pp.fixed_cost;
summary->preprocessor_time_in_seconds = WallTimeInSeconds() - start_time;
if (status) {
const double minimizer_start_time = WallTimeInSeconds();
Minimize(&pp, summary);
summary->minimizer_time_in_seconds =
WallTimeInSeconds() - minimizer_start_time;
} else {
summary->message = pp.error;
}
const double postprocessor_start_time = WallTimeInSeconds();
problem_impl = problem->mutable_impl();
program = problem_impl->mutable_program();
// On exit, ensure that the parameter blocks again point at the user
// provided values and the parameter blocks are numbered according
// to their position in the original user provided program.
program->SetParameterBlockStatePtrsToUserStatePtrs();
program->SetParameterOffsetsAndIndex();
PostSolveSummarize(pp, summary);
summary->postprocessor_time_in_seconds =
WallTimeInSeconds() - postprocessor_start_time;
// If the gradient checker reported an error, we want to report FAILURE
// instead of USER_FAILURE and provide the error log.
if (gradient_checking_callback.gradient_error_detected()) {
summary->termination_type = FAILURE;
summary->message = gradient_checking_callback.error_log();
}
summary->total_time_in_seconds = WallTimeInSeconds() - start_time;
}
void Solve(const Solver::Options& options,
Problem* problem,
Solver::Summary* summary) {
Solver solver;
solver.Solve(options, problem, summary);
}
string Solver::Summary::BriefReport() const {
return StringPrintf(
"Ceres Solver Report: "
"Iterations: %d, "
"Initial cost: %e, "
"Final cost: %e, "
"Termination: %s",
num_successful_steps + num_unsuccessful_steps,
initial_cost,
final_cost,
TerminationTypeToString(termination_type));
}
string Solver::Summary::FullReport() const {
using internal::VersionString;
// NOTE operator+ is not usable for concatenating a string and a string_view.
string report =
string{"\nSolver Summary (v "}.append(VersionString()) + ")\n\n";
StringAppendF(&report, "%45s %21s\n", "Original", "Reduced");
StringAppendF(&report,
"Parameter blocks % 25d% 25d\n",
num_parameter_blocks,
num_parameter_blocks_reduced);
StringAppendF(&report,
"Parameters % 25d% 25d\n",
num_parameters,
num_parameters_reduced);
if (num_effective_parameters_reduced != num_parameters_reduced) {
StringAppendF(&report,
"Effective parameters% 25d% 25d\n",
num_effective_parameters,
num_effective_parameters_reduced);
}
StringAppendF(&report,
"Residual blocks % 25d% 25d\n",
num_residual_blocks,
num_residual_blocks_reduced);
StringAppendF(&report,
"Residuals % 25d% 25d\n",
num_residuals,
num_residuals_reduced);
if (minimizer_type == TRUST_REGION) {
// TRUST_SEARCH HEADER
StringAppendF(
&report, "\nMinimizer %19s\n", "TRUST_REGION");
if (linear_solver_type_used == DENSE_NORMAL_CHOLESKY ||
linear_solver_type_used == DENSE_SCHUR ||
linear_solver_type_used == DENSE_QR) {
const char* mixed_precision_suffix =
(mixed_precision_solves_used ? "(Mixed Precision)" : "");
StringAppendF(&report,
"\nDense linear algebra library %15s %s\n",
DenseLinearAlgebraLibraryTypeToString(
dense_linear_algebra_library_type),
mixed_precision_suffix);
}
StringAppendF(&report,
"Trust region strategy %19s",
TrustRegionStrategyTypeToString(trust_region_strategy_type));
if (trust_region_strategy_type == DOGLEG) {
if (dogleg_type == TRADITIONAL_DOGLEG) {
StringAppendF(&report, " (TRADITIONAL)");
} else {
StringAppendF(&report, " (SUBSPACE)");
}
}
const bool used_sparse_linear_algebra_library =
linear_solver_type_used == SPARSE_NORMAL_CHOLESKY ||
linear_solver_type_used == SPARSE_SCHUR ||
linear_solver_type_used == CGNR ||
(linear_solver_type_used == ITERATIVE_SCHUR &&
(preconditioner_type_used == CLUSTER_JACOBI ||
preconditioner_type_used == CLUSTER_TRIDIAGONAL));
const bool linear_solver_ordering_required =
linear_solver_type_used == SPARSE_SCHUR ||
(linear_solver_type_used == ITERATIVE_SCHUR &&
(preconditioner_type_used == CLUSTER_JACOBI ||
preconditioner_type_used == CLUSTER_TRIDIAGONAL)) ||
(linear_solver_type_used == CGNR && preconditioner_type_used == SUBSET);
if (used_sparse_linear_algebra_library) {
const char* mixed_precision_suffix =
(mixed_precision_solves_used ? "(Mixed Precision)" : "");
if (linear_solver_ordering_required) {
StringAppendF(
&report,
"\nSparse linear algebra library %15s + %s %s\n",
SparseLinearAlgebraLibraryTypeToString(
sparse_linear_algebra_library_type),
LinearSolverOrderingTypeToString(linear_solver_ordering_type),
mixed_precision_suffix);
} else {
StringAppendF(&report,
"\nSparse linear algebra library %15s %s\n",
SparseLinearAlgebraLibraryTypeToString(
sparse_linear_algebra_library_type),
mixed_precision_suffix);
}
}
StringAppendF(&report, "\n");
StringAppendF(&report, "%45s %21s\n", "Given", "Used");
StringAppendF(&report,
"Linear solver %25s%25s\n",
LinearSolverTypeToString(linear_solver_type_given),
LinearSolverTypeToString(linear_solver_type_used));
if (IsIterativeSolver(linear_solver_type_given)) {
StringAppendF(&report,
"Preconditioner %25s%25s\n",
PreconditionerTypeToString(preconditioner_type_given),
PreconditionerTypeToString(preconditioner_type_used));
}
if (preconditioner_type_used == CLUSTER_JACOBI ||
preconditioner_type_used == CLUSTER_TRIDIAGONAL) {
StringAppendF(
&report,
"Visibility clustering%24s%25s\n",
VisibilityClusteringTypeToString(visibility_clustering_type),
VisibilityClusteringTypeToString(visibility_clustering_type));
}
StringAppendF(&report,
"Threads % 25d% 25d\n",
num_threads_given,
num_threads_used);
string given;
StringifyOrdering(linear_solver_ordering_given, &given);
string used;
StringifyOrdering(linear_solver_ordering_used, &used);
StringAppendF(&report,
"Linear solver ordering %22s %24s\n",
given.c_str(),
used.c_str());
if (IsSchurType(linear_solver_type_used)) {
StringAppendF(&report,
"Schur structure %22s %24s\n",
schur_structure_given.c_str(),
schur_structure_used.c_str());
}
if (inner_iterations_given) {
StringAppendF(&report,
"Use inner iterations %20s %20s\n",
inner_iterations_given ? "True" : "False",
inner_iterations_used ? "True" : "False");
}
if (inner_iterations_used) {
string given;
StringifyOrdering(inner_iteration_ordering_given, &given);
string used;
StringifyOrdering(inner_iteration_ordering_used, &used);
StringAppendF(&report,
"Inner iteration ordering %20s %24s\n",
given.c_str(),
used.c_str());
}
} else {
// LINE_SEARCH HEADER
StringAppendF(&report, "\nMinimizer %19s\n", "LINE_SEARCH");
string line_search_direction_string;
if (line_search_direction_type == LBFGS) {
line_search_direction_string = StringPrintf("LBFGS (%d)", max_lbfgs_rank);
} else if (line_search_direction_type == NONLINEAR_CONJUGATE_GRADIENT) {
line_search_direction_string = NonlinearConjugateGradientTypeToString(
nonlinear_conjugate_gradient_type);
} else {
line_search_direction_string =
LineSearchDirectionTypeToString(line_search_direction_type);
}
StringAppendF(&report,
"Line search direction %19s\n",
line_search_direction_string.c_str());
const string line_search_type_string = StringPrintf(
"%s %s",
LineSearchInterpolationTypeToString(line_search_interpolation_type),
LineSearchTypeToString(line_search_type));
StringAppendF(&report,
"Line search type %19s\n",
line_search_type_string.c_str());
StringAppendF(&report, "\n");
StringAppendF(&report, "%45s %21s\n", "Given", "Used");
StringAppendF(&report,
"Threads % 25d% 25d\n",
num_threads_given,
num_threads_used);
}
StringAppendF(&report, "\nCost:\n");
StringAppendF(&report, "Initial % 30e\n", initial_cost);
if (termination_type != FAILURE && termination_type != USER_FAILURE) {
StringAppendF(&report, "Final % 30e\n", final_cost);
StringAppendF(&report, "Change % 30e\n", initial_cost - final_cost);
}
StringAppendF(&report,
"\nMinimizer iterations % 16d\n",
num_successful_steps + num_unsuccessful_steps);
// Successful/Unsuccessful steps only matter in the case of the
// trust region solver. Line search terminates when it encounters
// the first unsuccessful step.
if (minimizer_type == TRUST_REGION) {
StringAppendF(&report,
"Successful steps % 14d\n",
num_successful_steps);
StringAppendF(&report,
"Unsuccessful steps % 14d\n",
num_unsuccessful_steps);
}
if (inner_iterations_used) {
StringAppendF(&report,
"Steps with inner iterations % 14d\n",
num_inner_iteration_steps);
}
const bool line_search_used =
(minimizer_type == LINE_SEARCH ||
(minimizer_type == TRUST_REGION && is_constrained));
if (line_search_used) {
StringAppendF(&report,
"Line search steps % 14d\n",
num_line_search_steps);
}
StringAppendF(&report, "\nTime (in seconds):\n");
StringAppendF(
&report, "Preprocessor %25.6f\n", preprocessor_time_in_seconds);
StringAppendF(&report,
"\n Residual only evaluation %18.6f (%d)\n",
residual_evaluation_time_in_seconds,
num_residual_evaluations);
if (line_search_used) {
StringAppendF(&report,
" Line search cost evaluation %10.6f\n",
line_search_cost_evaluation_time_in_seconds);
}
StringAppendF(&report,
" Jacobian & residual evaluation %12.6f (%d)\n",
jacobian_evaluation_time_in_seconds,
num_jacobian_evaluations);
if (line_search_used) {
StringAppendF(&report,
" Line search gradient evaluation %6.6f\n",
line_search_gradient_evaluation_time_in_seconds);
}
if (minimizer_type == TRUST_REGION) {
StringAppendF(&report,
" Linear solver %23.6f (%d)\n",
linear_solver_time_in_seconds,
num_linear_solves);
}
if (inner_iterations_used) {
StringAppendF(&report,
" Inner iterations %23.6f\n",
inner_iteration_time_in_seconds);
}
if (line_search_used) {
StringAppendF(&report,
" Line search polynomial minimization %.6f\n",
line_search_polynomial_minimization_time_in_seconds);
}
StringAppendF(
&report, "Minimizer %25.6f\n\n", minimizer_time_in_seconds);
StringAppendF(
&report, "Postprocessor %24.6f\n", postprocessor_time_in_seconds);
StringAppendF(
&report, "Total %25.6f\n\n", total_time_in_seconds);
StringAppendF(&report,
"Termination: %25s (%s)\n",
TerminationTypeToString(termination_type),
message.c_str());
return report;
}
bool Solver::Summary::IsSolutionUsable() const {
return internal::IsSolutionUsable(*this);
}
} // namespace ceres
Reference in New Issue View Git Blame Copy Permalink