diff --git a/docs/source/solving.rst b/docs/source/solving.rst index 89c837e9a..fc8219ea7 100644 --- a/docs/source/solving.rst +++ b/docs/source/solving.rst @@ -2262,7 +2262,8 @@ cases. .. member:: CovarianceAlgorithmType Covariance::Options::algorithm_type - Default: ``SPARSE_QR`` or ``DENSE_SVD`` + Default: ``SUITE_SPARSE_QR`` if ``SuiteSparseQR`` is installed and + ``EIGEN_SPARSE_QR`` otherwise. Ceres supports three different algorithms for covariance estimation, which represent different tradeoffs in speed, accuracy @@ -2281,47 +2282,23 @@ cases. small to moderate sized problems. It can handle full-rank as well as rank deficient Jacobians. - 2. ``SPARSE_CHOLESKY`` uses the ``CHOLMOD`` sparse Cholesky - factorization library to compute the decomposition : - - .. math:: R^\top R = J^\top J - - and then - - .. math:: \left(J^\top J\right)^{-1} = \left(R^\top R\right)^{-1} - - It a fast algorithm for sparse matrices that should be used when - the Jacobian matrix J is well conditioned. For ill-conditioned - matrices, this algorithm can fail unpredictabily. This is - because Cholesky factorization is not a rank-revealing - factorization, i.e., it cannot reliably detect when the matrix - being factorized is not of full - rank. ``SuiteSparse``/``CHOLMOD`` supplies a heuristic for - checking if the matrix is rank deficient (cholmod_rcond), but it - is only a heuristic and can have both false positive and false - negatives. - - Recent versions of ``SuiteSparse`` (>= 4.2.0) provide a much more - efficient method for solving for rows of the covariance - matrix. Therefore, if you are doing ``SPARSE_CHOLESKY``, we strongly - recommend using a recent version of ``SuiteSparse``. - - 3. ``SPARSE_QR`` uses the ``SuiteSparseQR`` sparse QR factorization - library to compute the decomposition + 2. ``EIGEN_SPARSE_QR`` uses the sparse QR factorization algorithm + in ``Eigen`` to compute the decomposition .. math:: QR &= J\\ \left(J^\top J\right)^{-1} &= \left(R^\top R\right)^{-1} - It is a moderately fast algorithm for sparse matrices, which at - the price of more time and memory than the ``SPARSE_CHOLESKY`` - algorithm is numerically better behaved and is rank revealing, - i.e., it can reliably detect when the Jacobian matrix is rank - deficient. + It is a moderately fast algorithm for sparse matrices. - Neither ``SPARSE_CHOLESKY`` or ``SPARSE_QR`` are capable of computing - the covariance if the Jacobian is rank deficient. + 3. ``SUITE_SPARSE_QR`` uses the sparse QR factorization algorithm + in ``SuiteSparse``. It uses dense linear algebra and is multi + threaded, so for large sparse sparse matrices it is + significantly faster than ``EIGEN_SPARSE_QR``. + + Neither ``EIGEN_SPARSE_QR`` nor ``SUITE_SPARSE_QR`` are capable of + computing the covariance if the Jacobian is rank deficient. .. member:: int Covariance::Options::min_reciprocal_condition_number @@ -2360,29 +2337,14 @@ cases. :math:`\sigma_{\text{max}}` are the minimum and maxiumum singular values of :math:`J` respectively. - 2. ``SPARSE_CHOLESKY`` - - .. math:: \text{cholmod_rcond} < \text{min_reciprocal_conditioner_number} - - Here cholmod_rcond is a crude estimate of the reciprocal - condition number of :math:`J^\top J` by using the maximum and - minimum diagonal entries of the Cholesky factor :math:`R`. There - are no theoretical guarantees associated with this test. It can - give false positives and negatives. Use at your own risk. The - default value of ``min_reciprocal_condition_number`` has been - set to a conservative value, and sometimes the - :func:`Covariance::Compute` may return false even if it is - possible to estimate the covariance reliably. In such cases, the - user should exercise their judgement before lowering the value - of ``min_reciprocal_condition_number``. - - 3. ``SPARSE_QR`` + 2. ``EIGEN_SPARSE_QR`` and ``SUITE_SPARSE_QR`` .. math:: \operatorname{rank}(J) < \operatorname{num\_col}(J) Here :\math:`\operatorname{rank}(J)` is the estimate of the - rank of `J` returned by the ``SuiteSparseQR`` algorithm. It is - a fairly reliable indication of rank deficiency. + rank of `J` returned by the sparse QR factorization + algorithm. It is a fairly reliable indication of rank + deficiency. .. member:: int Covariance::Options::null_space_rank @@ -2417,8 +2379,8 @@ cases. .. math:: \frac{\lambda_i}{\lambda_{\textrm{max}}} < \textrm{min_reciprocal_condition_number} - This option has no effect on ``SPARSE_QR`` and ``SPARSE_CHOLESKY`` - algorithms. + This option has no effect on ``EIGEN_SPARSE_QR`` and + ``SUITE_SPARSE_QR``. .. member:: bool Covariance::Options::apply_loss_function diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst index 1dcc002b1..a52ab30e0 100644 --- a/docs/source/version_history.rst +++ b/docs/source/version_history.rst @@ -10,12 +10,23 @@ HEAD #. Added ``Solver::Options::IsValid`` which allows users to validate their solver configuration before calling ``Solve``. +#. Added ``EIGEN_SPARSE_QR`` algorithm for covariance estimation using + ``Eigen``'s sparse QR factorization. (Michael Vitus) + Backward Incompatible API Changes --------------------------------- #. ``Solver::Options::solver_log`` has been removed. If needed this iteration callback can easily be implemented in user code. +#. The ``SPARSE_CHOLESKY`` algorithm for covariance estimation has + been removed. It is not rank revealing and numerically poorly + behaved. Sparse QR factorization is a much better way to do this. + +#. The ``SPARSE_QR`` algorithm for covariance estimation has been + renamed to ``SUITE_SPARSE_QR`` to be consistent with + ``EIGEN_SPARSE_QR``. + 1.9.0 ===== diff --git a/include/ceres/covariance.h b/include/ceres/covariance.h index 245381a2d..35fde4de0 100644 --- a/include/ceres/covariance.h +++ b/include/ceres/covariance.h @@ -202,9 +202,9 @@ class CERES_EXPORT Covariance { struct CERES_EXPORT Options { Options() #ifndef CERES_NO_SUITESPARSE - : algorithm_type(SPARSE_QR), + : algorithm_type(SUITE_SPARSE_QR), #else - : algorithm_type(DENSE_SVD), + : algorithm_type(EIGEN_SPARSE_QR), #endif min_reciprocal_condition_number(1e-14), null_space_rank(0), @@ -229,47 +229,22 @@ class CERES_EXPORT Covariance { // for small to moderate sized problems. It can handle // full-rank as well as rank deficient Jacobians. // - // 2. SPARSE_CHOLESKY uses the CHOLMOD sparse Cholesky - // factorization library to compute the decomposition : - // - // R'R = J'J - // - // and then - // - // [J'J]^-1 = [R'R]^-1 - // - // It a fast algorithm for sparse matrices that should be used - // when the Jacobian matrix J is well conditioned. For - // ill-conditioned matrices, this algorithm can fail - // unpredictabily. This is because Cholesky factorization is - // not a rank-revealing factorization, i.e., it cannot reliably - // detect when the matrix being factorized is not of full - // rank. SuiteSparse/CHOLMOD supplies a heuristic for checking - // if the matrix is rank deficient (cholmod_rcond), but it is - // only a heuristic and can have both false positive and false - // negatives. - // - // Recent versions of SuiteSparse (>= 4.2.0) provide a much - // more efficient method for solving for rows of the covariance - // matrix. Therefore, if you are doing SPARSE_CHOLESKY, we - // strongly recommend using a recent version of SuiteSparse. - // - // 3. SPARSE_QR uses the SuiteSparseQR sparse QR factorization - // library to compute the decomposition + // 2. EIGEN_SPARSE_QR uses the sparse QR factorization algorithm + // in Eigen to compute the decomposition // // Q * R = J // // [J'J]^-1 = [R*R']^-1 // - // It is a moderately fast algorithm for sparse matrices, which - // at the price of more time and memory than the - // SPARSE_CHOLESKY algorithm is numerically better behaved and - // is rank revealing, i.e., it can reliably detect when the - // Jacobian matrix is rank deficient. + // It is a moderately fast algorithm for sparse matrices. // - // Neither SPARSE_CHOLESKY or SPARSE_QR are capable of computing - // the covariance if the Jacobian is rank deficient. - + // 3. SUITE_SPARSE_QR uses the SuiteSparseQR sparse QR + // factorization algorithm. It uses dense linear algebra and is + // multi threaded, so for large sparse sparse matrices it is + // significantly faster than EIGEN_SPARSE_QR. + // + // Neither EIGEN_SPARSE_QR not SUITE_SPARSE_QR are capable of + // computing the covariance if the Jacobian is rank deficient. CovarianceAlgorithmType algorithm_type; // If the Jacobian matrix is near singular, then inverting J'J @@ -295,29 +270,13 @@ class CERES_EXPORT Covariance { // where min_sigma and max_sigma are the minimum and maxiumum // singular values of J respectively. // - // 2. SPARSE_CHOLESKY - // - // cholmod_rcond < min_reciprocal_conditioner_number - // - // Here cholmod_rcond is a crude estimate of the reciprocal - // condition number of J'J by using the maximum and minimum - // diagonal entries of the Cholesky factor R. There are no - // theoretical guarantees associated with this test. It can - // give false positives and negatives. Use at your own - // risk. The default value of min_reciprocal_condition_number - // has been set to a conservative value, and sometimes the - // Covariance::Compute may return false even if it is possible - // to estimate the covariance reliably. In such cases, the user - // should exercise their judgement before lowering the value of - // min_reciprocal_condition_number. - // - // 3. SPARSE_QR + // 2. SUITE_SPARSE_QR and EIGEN_SPARSE_QR // // rank(J) < num_col(J) // // Here rank(J) is the estimate of the rank of J returned by the - // SuiteSparseQR algorithm. It is a fairly reliable indication - // of rank deficiency. + // sparse QR factorization algorithm. It is a fairly reliable + // indication of rank deficiency. // double min_reciprocal_condition_number; @@ -352,8 +311,8 @@ class CERES_EXPORT Covariance { // // lambda_i / lambda_max < min_reciprocal_condition_number. // - // This option has no effect on the SPARSE_CHOLESKY or SPARSE_QR - // algorithms. + // This option has no effect on the SUITE_SPARSE_QR and + // EIGEN_SPARSE_QR algorithms. int null_space_rank; int num_threads; diff --git a/include/ceres/types.h b/include/ceres/types.h index ba31236a4..5fccd0363 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h @@ -399,8 +399,7 @@ enum LineSearchInterpolationType { enum CovarianceAlgorithmType { DENSE_SVD, - SPARSE_CHOLESKY, - SPARSE_QR, + SUITE_SPARSE_QR, EIGEN_SPARSE_QR }; diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc index 437a03db5..0fe481d43 100644 --- a/internal/ceres/covariance_impl.cc +++ b/internal/ceres/covariance_impl.cc @@ -55,40 +55,6 @@ namespace ceres { namespace internal { -namespace { - -// Per thread storage for SuiteSparse. -#ifndef CERES_NO_SUITESPARSE - -struct PerThreadContext { - explicit PerThreadContext(int num_rows) - : solution(NULL), - solution_set(NULL), - y_workspace(NULL), - e_workspace(NULL), - rhs(NULL) { - rhs = ss.CreateDenseVector(NULL, num_rows, num_rows); - } - - ~PerThreadContext() { - ss.Free(solution); - ss.Free(solution_set); - ss.Free(y_workspace); - ss.Free(e_workspace); - ss.Free(rhs); - } - - cholmod_dense* solution; - cholmod_sparse* solution_set; - cholmod_dense* y_workspace; - cholmod_dense* e_workspace; - cholmod_dense* rhs; - SuiteSparse ss; -}; - -#endif - -} // namespace typedef vector > CovarianceBlocks; @@ -398,10 +364,12 @@ bool CovarianceImpl::ComputeCovarianceValues() { case DENSE_SVD: return ComputeCovarianceValuesUsingDenseSVD(); #ifndef CERES_NO_SUITESPARSE - case SPARSE_CHOLESKY: - return ComputeCovarianceValuesUsingSparseCholesky(); - case SPARSE_QR: - return ComputeCovarianceValuesUsingSparseQR(); + case SUITE_SPARSE_QR: + return ComputeCovarianceValuesUsingSuiteSparseQR(); +#else + LOG(ERROR) << "SuiteSparse is required to use the " + << "SUITE_SPARSE_QR algorithm."; + return false; #endif case EIGEN_SPARSE_QR: return ComputeCovarianceValuesUsingEigenSparseQR(); @@ -413,197 +381,7 @@ bool CovarianceImpl::ComputeCovarianceValues() { return false; } -bool CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky() { - EventLogger event_logger( - "CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky"); -#ifndef CERES_NO_SUITESPARSE - if (covariance_matrix_.get() == NULL) { - // Nothing to do, all zeros covariance matrix. - return true; - } - - SuiteSparse ss; - - CRSMatrix jacobian; - problem_->Evaluate(evaluate_options_, NULL, NULL, NULL, &jacobian); - - event_logger.AddEvent("Evaluate"); - // m is a transposed view of the Jacobian. - cholmod_sparse cholmod_jacobian_view; - cholmod_jacobian_view.nrow = jacobian.num_cols; - cholmod_jacobian_view.ncol = jacobian.num_rows; - cholmod_jacobian_view.nzmax = jacobian.values.size(); - cholmod_jacobian_view.nz = NULL; - cholmod_jacobian_view.p = reinterpret_cast(&jacobian.rows[0]); - cholmod_jacobian_view.i = reinterpret_cast(&jacobian.cols[0]); - cholmod_jacobian_view.x = reinterpret_cast(&jacobian.values[0]); - cholmod_jacobian_view.z = NULL; - cholmod_jacobian_view.stype = 0; // Matrix is not symmetric. - cholmod_jacobian_view.itype = CHOLMOD_INT; - cholmod_jacobian_view.xtype = CHOLMOD_REAL; - cholmod_jacobian_view.dtype = CHOLMOD_DOUBLE; - cholmod_jacobian_view.sorted = 1; - cholmod_jacobian_view.packed = 1; - - string message; - cholmod_factor* factor = ss.AnalyzeCholesky(&cholmod_jacobian_view, &message); - event_logger.AddEvent("Symbolic Factorization"); - if (factor == NULL) { - LOG(ERROR) << "Covariance estimation failed. " - << "CHOLMOD symbolic cholesky factorization returned with: " - << message; - return false; - } - - LinearSolverTerminationType termination_type = - ss.Cholesky(&cholmod_jacobian_view, factor, &message); - event_logger.AddEvent("Numeric Factorization"); - if (termination_type != LINEAR_SOLVER_SUCCESS) { - LOG(ERROR) << "Covariance estimation failed. " - << "CHOLMOD numeric cholesky factorization returned with: " - << message; - ss.Free(factor); - return false; - } - - const double reciprocal_condition_number = - cholmod_rcond(factor, ss.mutable_cc()); - - if (reciprocal_condition_number < - options_.min_reciprocal_condition_number) { - LOG(ERROR) << "Cholesky factorization of J'J is not reliable. " - << "Reciprocal condition number: " - << reciprocal_condition_number << " " - << "min_reciprocal_condition_number: " - << options_.min_reciprocal_condition_number; - ss.Free(factor); - return false; - } - - const int num_rows = covariance_matrix_->num_rows(); - const int* rows = covariance_matrix_->rows(); - const int* cols = covariance_matrix_->cols(); - double* values = covariance_matrix_->mutable_values(); - - // The following loop exploits the fact that the i^th column of A^{-1} - // is given by the solution to the linear system - // - // A x = e_i - // - // where e_i is a vector with e(i) = 1 and all other entries zero. - // - // Since the covariance matrix is symmetric, the i^th row and column - // are equal. - // - // The ifdef separates two different version of SuiteSparse. Newer - // versions of SuiteSparse have the cholmod_solve2 function which - // re-uses memory across calls. -#if (SUITESPARSE_VERSION < 4002) - cholmod_dense* rhs = ss.CreateDenseVector(NULL, num_rows, num_rows); - double* rhs_x = reinterpret_cast(rhs->x); - - for (int r = 0; r < num_rows; ++r) { - int row_begin = rows[r]; - int row_end = rows[r + 1]; - if (row_end == row_begin) { - continue; - } - - rhs_x[r] = 1.0; - cholmod_dense* solution = ss.Solve(factor, rhs, &message); - double* solution_x = reinterpret_cast(solution->x); - for (int idx = row_begin; idx < row_end; ++idx) { - const int c = cols[idx]; - values[idx] = solution_x[c]; - } - ss.Free(solution); - rhs_x[r] = 0.0; - } - - ss.Free(rhs); -#else // SUITESPARSE_VERSION < 4002 - - const int num_threads = options_.num_threads; - vector contexts(num_threads); - for (int i = 0; i < num_threads; ++i) { - contexts[i] = new PerThreadContext(num_rows); - } - - // The first call to cholmod_solve2 is not thread safe, since it - // changes the factorization from supernodal to simplicial etc. - { - PerThreadContext* context = contexts[0]; - double* context_rhs_x = reinterpret_cast(context->rhs->x); - context_rhs_x[0] = 1.0; - cholmod_solve2(CHOLMOD_A, - factor, - context->rhs, - NULL, - &context->solution, - &context->solution_set, - &context->y_workspace, - &context->e_workspace, - context->ss.mutable_cc()); - context_rhs_x[0] = 0.0; - } - -#pragma omp parallel for num_threads(num_threads) schedule(dynamic) - for (int r = 0; r < num_rows; ++r) { - int row_begin = rows[r]; - int row_end = rows[r + 1]; - if (row_end == row_begin) { - continue; - } - -# ifdef CERES_USE_OPENMP - int thread_id = omp_get_thread_num(); -# else - int thread_id = 0; -# endif - - PerThreadContext* context = contexts[thread_id]; - double* context_rhs_x = reinterpret_cast(context->rhs->x); - context_rhs_x[r] = 1.0; - - // TODO(sameeragarwal) There should be a more efficient way - // involving the use of Bset but I am unable to make it work right - // now. - cholmod_solve2(CHOLMOD_A, - factor, - context->rhs, - NULL, - &context->solution, - &context->solution_set, - &context->y_workspace, - &context->e_workspace, - context->ss.mutable_cc()); - - double* solution_x = reinterpret_cast(context->solution->x); - for (int idx = row_begin; idx < row_end; ++idx) { - const int c = cols[idx]; - values[idx] = solution_x[c]; - } - context_rhs_x[r] = 0.0; - } - - for (int i = 0; i < num_threads; ++i) { - delete contexts[i]; - } - -#endif // SUITESPARSE_VERSION < 4002 - - ss.Free(factor); - event_logger.AddEvent("Inversion"); - return true; - -#else // CERES_NO_SUITESPARSE - - return false; - -#endif // CERES_NO_SUITESPARSE -} - -bool CovarianceImpl::ComputeCovarianceValuesUsingSparseQR() { +bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparseQR() { EventLogger event_logger( "CovarianceImpl::ComputeCovarianceValuesUsingSparseQR"); diff --git a/internal/ceres/covariance_impl.h b/internal/ceres/covariance_impl.h index ff640b9d6..135f4a1d6 100644 --- a/internal/ceres/covariance_impl.h +++ b/internal/ceres/covariance_impl.h @@ -64,9 +64,8 @@ class CovarianceImpl { ProblemImpl* problem); bool ComputeCovarianceValues(); - bool ComputeCovarianceValuesUsingSparseCholesky(); - bool ComputeCovarianceValuesUsingSparseQR(); bool ComputeCovarianceValuesUsingDenseSVD(); + bool ComputeCovarianceValuesUsingSuiteSparseQR(); bool ComputeCovarianceValuesUsingEigenSparseQR(); const CompressedRowSparseMatrix* covariance_matrix() const { diff --git a/internal/ceres/covariance_test.cc b/internal/ceres/covariance_test.cc index 6de3fecd4..6c506b72d 100644 --- a/internal/ceres/covariance_test.cc +++ b/internal/ceres/covariance_test.cc @@ -400,10 +400,7 @@ TEST_F(CovarianceTest, NormalBehavior) { Covariance::Options options; #ifndef CERES_NO_SUITESPARSE - options.algorithm_type = SPARSE_CHOLESKY; - ComputeAndCompareCovarianceBlocks(options, expected_covariance); - - options.algorithm_type = SPARSE_QR; + options.algorithm_type = SUITE_SPARSE_QR; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif @@ -451,10 +448,7 @@ TEST_F(CovarianceTest, ThreadedNormalBehavior) { options.num_threads = 4; #ifndef CERES_NO_SUITESPARSE - options.algorithm_type = SPARSE_CHOLESKY; - ComputeAndCompareCovarianceBlocks(options, expected_covariance); - - options.algorithm_type = SPARSE_QR; + options.algorithm_type = SUITE_SPARSE_QR; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif @@ -503,10 +497,7 @@ TEST_F(CovarianceTest, ConstantParameterBlock) { Covariance::Options options; #ifndef CERES_NO_SUITESPARSE - options.algorithm_type = SPARSE_CHOLESKY; - ComputeAndCompareCovarianceBlocks(options, expected_covariance); - - options.algorithm_type = SPARSE_QR; + options.algorithm_type = SUITE_SPARSE_QR; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif @@ -562,10 +553,7 @@ TEST_F(CovarianceTest, LocalParameterization) { Covariance::Options options; #ifndef CERES_NO_SUITESPARSE - options.algorithm_type = SPARSE_CHOLESKY; - ComputeAndCompareCovarianceBlocks(options, expected_covariance); - - options.algorithm_type = SPARSE_QR; + options.algorithm_type = SUITE_SPARSE_QR; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif @@ -793,8 +781,7 @@ class LargeScaleCovarianceTest : public ::testing::Test { #if !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP) TEST_F(LargeScaleCovarianceTest, Parallel) { - ComputeAndCompare(SPARSE_CHOLESKY, 4); - ComputeAndCompare(SPARSE_QR, 4); + ComputeAndCompare(SUITE_SPARSE_QR, 4); } #endif // !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP) diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 8eaadc3e1..adf109089 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc @@ -261,8 +261,8 @@ const char* CovarianceAlgorithmTypeToString( CovarianceAlgorithmType type) { switch (type) { CASESTR(DENSE_SVD); - CASESTR(SPARSE_CHOLESKY); - CASESTR(SPARSE_QR); + CASESTR(EIGEN_SPARSE_QR); + CASESTR(SUITE_SPARSE_QR); default: return "UNKNOWN"; } @@ -273,8 +273,8 @@ bool StringToCovarianceAlgorithmType( CovarianceAlgorithmType* type) { UpperCase(&value); STRENUM(DENSE_SVD); - STRENUM(SPARSE_CHOLESKY); - STRENUM(SPARSE_QR); + STRENUM(EIGEN_SPARSE_QR); + STRENUM(SUITE_SPARSE_QR); return false; }