Covariance estimation using SuiteSparseQR. Change-Id: I70d1686e3288fdde5f9723e832e15ffb857d6d85
diff --git a/CMakeLists.txt b/CMakeLists.txt index 0ebd6ff..c4c89bc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt
@@ -215,6 +215,23 @@ SET(CHOLMOD_FOUND FALSE) ENDIF (EXISTS ${CHOLMOD_INCLUDE}) +SET(SUITESPARSEQR_FOUND TRUE) +FIND_LIBRARY(SUITESPARSEQR_LIB NAMES spqr PATHS ${SUITESPARSE_SEARCH_LIBS}) +IF (EXISTS ${SUITESPARSEQR_LIB}) + MESSAGE("-- Found SUITESPARSEQR library: ${SUITESPARSEQR_LIB}") +ELSE (EXISTS ${SUITESPARSEQR_LIB}) + MESSAGE("-- Did not find SUITESPARSEQR library") + SET(SUITESPARSEQR_FOUND FALSE) +ENDIF (EXISTS ${SUITESPARSEQR_LIB}) + +FIND_PATH(SUITESPARSEQR_INCLUDE NAMES SuiteSparseQR.hpp PATHS ${SUITESPARSE_SEARCH_HEADERS}) +IF (EXISTS ${SUITESPARSEQR_INCLUDE}) + MESSAGE("-- Found SUITESPARSEQR header in: ${SUITESPARSEQR_INCLUDE}") +ELSE (EXISTS ${SUITESPARSEQR_INCLUDE}) + MESSAGE("-- Did not find SUITESPARSEQR header") + SET(SUITESPARSEQR_FOUND FALSE) +ENDIF (EXISTS ${SUITESPARSEQR_INCLUDE}) + # If SuiteSparse version is >= 4 then SuiteSparse_config is required. # For SuiteSparse 3, UFconfig.h is required. SET(SUITESPARSE_CONFIG_FOUND TRUE) @@ -259,6 +276,24 @@ MESSAGE("-- Did not find METIS library") ENDIF (EXISTS ${METIS_LIB}) +# SuiteSparseQR may be compiled with Intel Threading Building Blocks. +SET(TBB_FOUND TRUE) +FIND_LIBRARY(TBB_LIB NAMES tbb PATHS ${SEARCH_LIBS}) +IF (EXISTS ${TBB_LIB}) + MESSAGE("-- Found TBB library: ${TBB_LIB}") +ELSE (EXISTS ${TBB_LIB}) + MESSAGE("-- Did not find TBB library") + SET(TBB_FOUND FALSE) +ENDIF (EXISTS ${TBB_LIB}) + +FIND_LIBRARY(TBB_MALLOC_LIB NAMES tbbmalloc PATHS ${SEARCH_LIBS}) +IF (EXISTS ${TBB_MALLOC_LIB}) + MESSAGE("-- Found TBB Malloc library: ${TBB_MALLOC_LIB}") +ELSE (EXISTS ${TBB_MALLOC_LIB}) + MESSAGE("-- Did not find TBB library") + SET(TBB_FOUND FALSE) +ENDIF (EXISTS ${TBB_MALLOC_LIB}) + SET(BLAS_AND_LAPACK_FOUND TRUE) IF (APPLE) # Mac OS X has LAPACK/BLAS bundled in a framework called @@ -557,6 +592,7 @@ INCLUDE_DIRECTORIES(${COLAMD_INCLUDE}) INCLUDE_DIRECTORIES(${CCOLAMD_INCLUDE}) INCLUDE_DIRECTORIES(${CHOLMOD_INCLUDE}) + INCLUDE_DIRECTORIES(${SUITESPARSEQR_INCLUDE}) IF (SUITESPARSE_CONFIG_FOUND) INCLUDE_DIRECTORIES(${SUITESPARSE_CONFIG_INCLUDE}) ENDIF (SUITESPARSE_CONFIG_FOUND)
diff --git a/include/ceres/covariance.h b/include/ceres/covariance.h index b70f98b..83126b5 100644 --- a/include/ceres/covariance.h +++ b/include/ceres/covariance.h
@@ -200,34 +200,76 @@ public: struct Options { Options() - : num_threads(1), #ifndef CERES_NO_SUITESPARSE - use_dense_linear_algebra(false), + : algorithm_type(SPARSE_QR), #else - use_dense_linear_algebra(true), + : algorithm_type(DENSE_SVD), #endif min_reciprocal_condition_number(1e-14), null_space_rank(0), + num_threads(1), apply_loss_function(true) { } - // Number of threads to be used for evaluating the Jacobian and - // estimation of covariance. - int num_threads; - - - // When use_dense_linear_algebra = true, Eigen's JacobiSVD - // algorithm is used to perform the computations. It is an - // accurate but slow method and should only be used for small to - // moderate sized problems. + // Ceres supports three different algorithms for covariance + // estimation, which represent different tradeoffs in speed, + // accuracy and reliability. // - // When use_dense_linear_algebra = false, SuiteSparse/CHOLMOD is - // used to perform the computation. 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 - // large scale covariance estimation, we strongly recommend using - // a recent version of SuiteSparse. - bool use_dense_linear_algebra; + // 1. DENSE_SVD uses Eigen's JacobiSVD to perform the + // computations. It computes the singular value decomposition + // + // U * S * V' = J + // + // and then uses it to compute the pseudo inverse of J'J as + // + // pseudoinverse[J'J]^ = V * pseudoinverse[S] * V' + // + // It is an accurate but slow method and should only be used + // 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 + // + // 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. + // + // Neither SPARSE_CHOLESKY or 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 // will result in unreliable results, e.g, if @@ -240,48 +282,46 @@ // inv(J'J) = [ 2.0471e+14 -2.0471e+14] // [-2.0471e+14 2.0471e+14] // - // This is not a useful result. + // This is not a useful result. Therefore, by default + // Covariance::Compute will return false if a rank deficient + // Jacobian is encountered. How rank deficiency is detected + // depends on the algorithm being used. // - // The reciprocal condition number of a matrix is a measure of - // ill-conditioning or how close the matrix is to being - // singular/rank deficient. It is defined as the ratio of the - // smallest eigenvalue of the matrix to the largest eigenvalue. In - // the above case the reciprocal condition number is about - // 1e-16. Which is close to machine precision and even though the - // inverse exists, it is meaningless, and care should be taken to - // interpet the results of such an inversion. + // 1. DENSE_SVD // - // Matrices with condition number lower than - // min_reciprocal_condition_number are considered rank deficient - // and by default Covariance::Compute will return false if it - // encounters such a matrix. + // min_sigma / max_sigma < sqrt(min_reciprocal_condition_number) // - // use_dense_linear_algebra = false - // -------------------------------- + // where min_sigma and max_sigma are the minimum and maxiumum + // singular values of J respectively. // - // When performing large scale sparse covariance estimation, - // computing the exact value of the reciprocal condition number is - // not possible as it would require computing the eigenvalues of - // J'J. + // 2. SPARSE_CHOLESKY // - // In this case we use cholmod_rcond, which uses the ratio of the - // smallest to the largest diagonal entries of the Cholesky - // factorization as an approximation to the reciprocal condition - // number. + // cholmod_rcond < min_reciprocal_conditioner_number // - // However, care must be taken as this is a heuristic and can - // sometimes be a very crude estimate. 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. + // 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. // - // use_dense_linear_algebra = true - // ------------------------------- + // 3. SPARSE_QR // - // When using dense linear algebra, the user has more control in - // dealing with singular and near singular covariance matrices. + // 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. + // + double min_reciprocal_condition_number; + + // When using DENSE_SVD, the user has more control in dealing with + // singular and near singular covariance matrices. // // As mentioned above, when the covariance matrix is near // singular, instead of computing the inverse of J'J, the @@ -311,19 +351,12 @@ // // lambda_i / lambda_max < min_reciprocal_condition_number. // - double min_reciprocal_condition_number; - - // Truncate the smallest "null_space_rank" eigenvectors when - // computing the pseudo inverse of J'J. - // - // If null_space_rank = -1, then all eigenvectors with eigenvalues s.t. - // - // lambda_i / lambda_max < min_reciprocal_condition_number. - // - // are dropped. See the documentation for - // min_reciprocal_condition_number for more details. + // This option has no effect on the SPARSE_CHOLESKY or SPARSE_QR + // algorithms. int null_space_rank; + int num_threads; + // Even though the residual blocks in the problem may contain loss // functions, setting apply_loss_function to false will turn off // the application of the loss function to the output of the cost
diff --git a/include/ceres/types.h b/include/ceres/types.h index 46c2fd2..a967541 100644 --- a/include/ceres/types.h +++ b/include/ceres/types.h
@@ -383,6 +383,12 @@ CUBIC }; +enum CovarianceAlgorithmType { + DENSE_SVD, + SPARSE_CHOLESKY, + SPARSE_QR +}; + const char* LinearSolverTypeToString(LinearSolverType type); bool StringToLinearSolverType(string value, LinearSolverType* type); @@ -424,6 +430,12 @@ string value, LineSearchInterpolationType* type); +const char* CovarianceAlgorithmTypeToString( + CovarianceAlgorithmType type); +bool StringToCovarianceAlgorithmType( + string value, + CovarianceAlgorithmType* type); + const char* LinearSolverTerminationTypeToString( LinearSolverTerminationType type);
diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index dfa567c..a7f1d97 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt
@@ -152,6 +152,7 @@ ENDIF (GFLAGS) IF (SUITESPARSE_FOUND) + LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${SUITESPARSEQR_LIB}) LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CHOLMOD_LIB}) LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CCOLAMD_LIB}) LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CAMD_LIB}) @@ -170,8 +171,14 @@ IF (EXISTS ${BLAS_LIB}) LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${BLAS_LIB}) ENDIF (EXISTS ${BLAS_LIB}) + + IF (TBB_FOUND) + LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${TBB_LIB}) + LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${TBB_MALLOC_LIB}) + ENDIF (TBB_FOUND) ENDIF (SUITESPARSE_FOUND) + IF (CXSPARSE_FOUND) LIST(APPEND CERES_LIBRARY_DEPENDENCIES ${CXSPARSE_LIB}) ENDIF (CXSPARSE_FOUND)
diff --git a/internal/ceres/covariance_impl.cc b/internal/ceres/covariance_impl.cc index e3401ea..bb04cd0 100644 --- a/internal/ceres/covariance_impl.cc +++ b/internal/ceres/covariance_impl.cc
@@ -48,6 +48,7 @@ #include "ceres/suitesparse.h" #include "ceres/wall_time.h" #include "glog/logging.h" +#include "SuiteSparseQR.hpp" namespace ceres { namespace internal { @@ -388,28 +389,34 @@ } bool CovarianceImpl::ComputeCovarianceValues() { - if (options_.use_dense_linear_algebra) { - return ComputeCovarianceValuesUsingEigen(); - } - + switch (options_.algorithm_type) { + case (DENSE_SVD): + return ComputeCovarianceValuesUsingDenseSVD(); #ifndef CERES_NO_SUITESPARSE - return ComputeCovarianceValuesUsingSuiteSparse(); -#else - LOG(ERROR) << "Ceres compiled without SuiteSparse. " - << "Large scale covariance computation is not possible."; - return false; + case (SPARSE_CHOLESKY): + return ComputeCovarianceValuesUsingSparseCholesky(); + case (SPARSE_QR): + return ComputeCovarianceValuesUsingSparseQR(); #endif + default: + LOG(ERROR) << "Unsupported covariance estimation algorithm type: " + << CovarianceAlgorithmTypeToString(options_.algorithm_type); + return false; + } + return false; } -bool CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse() { +bool CovarianceImpl::ComputeCovarianceValuesUsingSparseCholesky() { EventLogger event_logger( - "CovarianceImpl::ComputeCovarianceValuesUsingSuiteSparse"); + "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); @@ -431,12 +438,12 @@ cholmod_jacobian_view.sorted = 1; cholmod_jacobian_view.packed = 1; - cholmod_factor* factor = ss_.AnalyzeCholesky(&cholmod_jacobian_view); + cholmod_factor* factor = ss.AnalyzeCholesky(&cholmod_jacobian_view); event_logger.AddEvent("Symbolic Factorization"); - bool factorization_succeeded = ss_.Cholesky(&cholmod_jacobian_view, factor); + bool factorization_succeeded = ss.Cholesky(&cholmod_jacobian_view, factor); if (factorization_succeeded) { const double reciprocal_condition_number = - cholmod_rcond(factor, ss_.mutable_cc()); + cholmod_rcond(factor, ss.mutable_cc()); if (reciprocal_condition_number < options_.min_reciprocal_condition_number) { LOG(WARNING) << "Cholesky factorization of J'J is not reliable. " @@ -450,7 +457,7 @@ event_logger.AddEvent("Numeric Factorization"); if (!factorization_succeeded) { - ss_.Free(factor); + ss.Free(factor); LOG(WARNING) << "Cholesky factorization failed."; return false; } @@ -474,7 +481,7 @@ // 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); + cholmod_dense* rhs = ss.CreateDenseVector(NULL, num_rows, num_rows); double* rhs_x = reinterpret_cast<double*>(rhs->x); for (int r = 0; r < num_rows; ++r) { @@ -485,17 +492,17 @@ } rhs_x[r] = 1.0; - cholmod_dense* solution = ss_.Solve(factor, rhs); + cholmod_dense* solution = ss.Solve(factor, rhs); double* solution_x = reinterpret_cast<double*>(solution->x); for (int idx = row_begin; idx < row_end; ++idx) { const int c = cols[idx]; values[idx] = solution_x[c]; } - ss_.Free(solution); + ss.Free(solution); rhs_x[r] = 0.0; } - ss_.Free(rhs); + ss.Free(rhs); #else // SUITESPARSE_VERSION < 4002 const int num_threads = options_.num_threads; @@ -567,7 +574,7 @@ #endif // SUITESPARSE_VERSION < 4002 - ss_.Free(factor); + ss.Free(factor); event_logger.AddEvent("Inversion"); return true; @@ -578,9 +585,144 @@ #endif // CERES_NO_SUITESPARSE }; -bool CovarianceImpl::ComputeCovarianceValuesUsingEigen() { +bool CovarianceImpl::ComputeCovarianceValuesUsingSparseQR() { EventLogger event_logger( - "CovarianceImpl::ComputeCovarianceValuesUsingEigen"); + "CovarianceImpl::ComputeCovarianceValuesUsingSparseQR"); + +#ifndef CERES_NO_SUITESPARSE + if (covariance_matrix_.get() == NULL) { + // Nothing to do, all zeros covariance matrix. + return true; + } + + CRSMatrix jacobian; + problem_->Evaluate(evaluate_options_, NULL, NULL, NULL, &jacobian); + event_logger.AddEvent("Evaluate"); + + // Construct a compressed column form of the Jacobian. + const int num_rows = jacobian.num_rows; + const int num_cols = jacobian.num_cols; + const int num_nonzeros = jacobian.values.size(); + + vector<SuiteSparse_long> transpose_rows(num_cols + 1, 0); + vector<SuiteSparse_long> transpose_cols(num_nonzeros, 0); + vector<double> transpose_values(num_nonzeros, 0); + + for (int idx = 0; idx < num_nonzeros; ++idx) { + transpose_rows[jacobian.cols[idx] + 1] += 1; + } + + for (int i = 1; i < transpose_rows.size(); ++i) { + transpose_rows[i] += transpose_rows[i - 1]; + } + + for (int r = 0; r < num_rows; ++r) { + for (int idx = jacobian.rows[r]; idx < jacobian.rows[r + 1]; ++idx) { + const int c = jacobian.cols[idx]; + const int transpose_idx = transpose_rows[c]; + transpose_cols[transpose_idx] = r; + transpose_values[transpose_idx] = jacobian.values[idx]; + ++transpose_rows[c]; + } + } + + for (int i = transpose_rows.size() - 1; i > 0 ; --i) { + transpose_rows[i] = transpose_rows[i - 1]; + } + transpose_rows[0] = 0; + + cholmod_sparse cholmod_jacobian; + cholmod_jacobian.nrow = num_rows; + cholmod_jacobian.ncol = num_cols; + cholmod_jacobian.nzmax = num_nonzeros; + cholmod_jacobian.nz = NULL; + cholmod_jacobian.p = reinterpret_cast<void*>(&transpose_rows[0]); + cholmod_jacobian.i = reinterpret_cast<void*>(&transpose_cols[0]); + cholmod_jacobian.x = reinterpret_cast<void*>(&transpose_values[0]); + cholmod_jacobian.z = NULL; + cholmod_jacobian.stype = 0; // Matrix is not symmetric. + cholmod_jacobian.itype = CHOLMOD_LONG; + cholmod_jacobian.xtype = CHOLMOD_REAL; + cholmod_jacobian.dtype = CHOLMOD_DOUBLE; + cholmod_jacobian.sorted = 1; + cholmod_jacobian.packed = 1; + + cholmod_common cc; + cholmod_l_start(&cc); + + SuiteSparseQR_factorization<double>* factor = + SuiteSparseQR_factorize<double>(SPQR_ORDERING_BESTAMD, + SPQR_DEFAULT_TOL, + &cholmod_jacobian, + &cc); + event_logger.AddEvent("Numeric Factorization"); + + const int rank = cc.SPQR_istat[4]; + if (rank < cholmod_jacobian.ncol) { + LOG(WARNING) << "Jacobian matrix is rank deficient." + << "Number of columns: " << cholmod_jacobian.ncol + << " rank: " << rank; + SuiteSparseQR_free(&factor, &cc); + cholmod_l_finish(&cc); + return false; + } + + 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. + + cholmod_dense* rhs = cholmod_l_zeros(num_cols, 1, CHOLMOD_REAL, &cc); + double* rhs_x = reinterpret_cast<double*>(rhs->x); + + for (int r = 0; r < num_cols; ++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* y1 = SuiteSparseQR_solve<double>(SPQR_RTX_EQUALS_ETB, factor, rhs, &cc); + cholmod_dense* solution = SuiteSparseQR_solve<double>(SPQR_RETX_EQUALS_B, factor, y1, &cc); + + double* solution_x = reinterpret_cast<double*>(solution->x); + for (int idx = row_begin; idx < row_end; ++idx) { + const int c = cols[idx]; + values[idx] = solution_x[c]; + } + + cholmod_l_free_dense(&y1, &cc); + cholmod_l_free_dense(&solution, &cc); + rhs_x[r] = 0.0; + } + + cholmod_l_free_dense(&rhs, &cc); + SuiteSparseQR_free(&factor, &cc); + cholmod_l_finish(&cc); + event_logger.AddEvent("Inversion"); + return true; + +#else // CERES_NO_SUITESPARSE + + return false; + +#endif // CERES_NO_SUITESPARSE +}; + +bool CovarianceImpl::ComputeCovarianceValuesUsingDenseSVD() { + EventLogger event_logger( + "CovarianceImpl::ComputeCovarianceValuesUsingDenseSVD"); if (covariance_matrix_.get() == NULL) { // Nothing to do, all zeros covariance matrix. return true; @@ -602,6 +744,7 @@ Eigen::JacobiSVD<Matrix> svd(dense_jacobian, Eigen::ComputeThinU | Eigen::ComputeThinV); + event_logger.AddEvent("SingularValueDecomposition"); const Vector singular_values = svd.singularValues();
diff --git a/internal/ceres/covariance_impl.h b/internal/ceres/covariance_impl.h index 0be53e7..0e7e217 100644 --- a/internal/ceres/covariance_impl.h +++ b/internal/ceres/covariance_impl.h
@@ -64,8 +64,9 @@ ProblemImpl* problem); bool ComputeCovarianceValues(); - bool ComputeCovarianceValuesUsingSuiteSparse(); - bool ComputeCovarianceValuesUsingEigen(); + bool ComputeCovarianceValuesUsingSparseCholesky(); + bool ComputeCovarianceValuesUsingSparseQR(); + bool ComputeCovarianceValuesUsingDenseSVD(); const CompressedRowSparseMatrix* covariance_matrix() const { return covariance_matrix_.get(); @@ -80,9 +81,6 @@ map<const double*, int> parameter_block_to_row_index_; set<const double*> constant_parameter_blocks_; scoped_ptr<CompressedRowSparseMatrix> covariance_matrix_; -#ifndef CERES_NO_SUITESPARSE - SuiteSparse ss_; -#endif }; } // namespace internal
diff --git a/internal/ceres/covariance_test.cc b/internal/ceres/covariance_test.cc index f9c5cd4..e7d25a1 100644 --- a/internal/ceres/covariance_test.cc +++ b/internal/ceres/covariance_test.cc
@@ -400,11 +400,14 @@ Covariance::Options options; #ifndef CERES_NO_SUITESPARSE - options.use_dense_linear_algebra = false; + options.algorithm_type = SPARSE_CHOLESKY; + ComputeAndCompareCovarianceBlocks(options, expected_covariance); + + options.algorithm_type = SPARSE_QR; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; ComputeAndCompareCovarianceBlocks(options, expected_covariance); } @@ -445,11 +448,14 @@ options.num_threads = 4; #ifndef CERES_NO_SUITESPARSE - options.use_dense_linear_algebra = false; + options.algorithm_type = SPARSE_CHOLESKY; + ComputeAndCompareCovarianceBlocks(options, expected_covariance); + + options.algorithm_type = SPARSE_QR; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; ComputeAndCompareCovarianceBlocks(options, expected_covariance); } @@ -491,11 +497,11 @@ Covariance::Options options; #ifndef CERES_NO_SUITESPARSE - options.use_dense_linear_algebra = false; + options.algorithm_type = SPARSE_CHOLESKY; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; ComputeAndCompareCovarianceBlocks(options, expected_covariance); } @@ -544,11 +550,11 @@ Covariance::Options options; #ifndef CERES_NO_SUITESPARSE - options.use_dense_linear_algebra = false; + options.algorithm_type = SPARSE_CHOLESKY; ComputeAndCompareCovarianceBlocks(options, expected_covariance); #endif - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; ComputeAndCompareCovarianceBlocks(options, expected_covariance); } @@ -589,7 +595,7 @@ { Covariance::Options options; - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; // Force dropping of the smallest eigenvector. options.null_space_rank = 1; ComputeAndCompareCovarianceBlocks(options, expected_covariance); @@ -597,7 +603,7 @@ { Covariance::Options options; - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; // Force dropping of the smallest eigenvector via the ratio but // automatic truncation. options.min_reciprocal_condition_number = 0.044494; @@ -693,7 +699,7 @@ }; Covariance::Options options; - options.use_dense_linear_algebra = true; + options.algorithm_type = DENSE_SVD; options.null_space_rank = -1; ComputeAndCompareCovarianceBlocks(options, expected_covariance); } @@ -723,9 +729,10 @@ } } - void ComputeAndCompare(int num_threads) { + void ComputeAndCompare(CovarianceAlgorithmType algorithm_type, + int num_threads) { Covariance::Options options; - options.use_dense_linear_algebra = false; + options.algorithm_type = algorithm_type; options.num_threads = num_threads; Covariance covariance(options); EXPECT_TRUE(covariance.Compute(all_covariance_blocks_, &problem_)); @@ -768,7 +775,7 @@ #if !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP) TEST_F(LargeScaleCovarianceTest, Parallel) { - ComputeAndCompare(4); + ComputeAndCompare(SPARSE_CHOLESKY, 4); } #endif // !defined(CERES_NO_SUITESPARSE) && defined(CERES_USE_OPENMP)
diff --git a/internal/ceres/types.cc b/internal/ceres/types.cc index 42990e3..164185e 100644 --- a/internal/ceres/types.cc +++ b/internal/ceres/types.cc
@@ -239,6 +239,27 @@ return false; } +const char* CovarianceAlgorithmTypeToString( + CovarianceAlgorithmType type) { + switch (type) { + CASESTR(DENSE_SVD); + CASESTR(SPARSE_CHOLESKY); + CASESTR(SPARSE_QR); + default: + return "UNKNOWN"; + } +} + +bool StringToCovarianceAlgorithmType( + string value, + CovarianceAlgorithmType* type) { + UpperCase(&value); + STRENUM(DENSE_SVD); + STRENUM(SPARSE_CHOLESKY); + STRENUM(SPARSE_QR); + return false; +} + const char* SolverTerminationTypeToString(SolverTerminationType type) { switch (type) { CASESTR(NO_CONVERGENCE);