internal/ceres/sparse_cholesky_test.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2017 Google Inc. All rights reserved.
   3 // http://ceres-solver.org/
   4 //
   5 // Redistribution and use in source and binary forms, with or without
   6 // modification, are permitted provided that the following conditions are met:
   7 //
   8 // * Redistributions of source code must retain the above copyright notice,
   9 //   this list of conditions and the following disclaimer.
  10 // * Redistributions in binary form must reproduce the above copyright notice,
  11 //   this list of conditions and the following disclaimer in the documentation
  12 //   and/or other materials provided with the distribution.
  13 // * Neither the name of Google Inc. nor the names of its contributors may be
  14 //   used to endorse or promote products derived from this software without
  15 //   specific prior written permission.
  16 //
  17 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
  18 // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  19 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  20 // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
  21 // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
  22 // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
  23 // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
  24 // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
  25 // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  26 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
  27 // POSSIBILITY OF SUCH DAMAGE.
  28 //
  29 // Author: sameeragarwal@google.com (Sameer Agarwal)
  30
  31 #include "ceres/sparse_cholesky.h"
  32
  33 #include <numeric>
  34 #include <vector>
  35
  36 #include "Eigen/Dense"
  37 #include "Eigen/SparseCore"
  38 #include "ceres/block_sparse_matrix.h"
  39 #include "ceres/compressed_row_sparse_matrix.h"
  40 #include "ceres/inner_product_computer.h"
  41 #include "ceres/internal/eigen.h"
  42 #include "ceres/internal/scoped_ptr.h"
  43 #include "ceres/random.h"
  44 #include "glog/logging.h"
  45 #include "gtest/gtest.h"
  46
  47 namespace ceres {
  48 namespace internal {
  49
  50 BlockSparseMatrix* CreateRandomFullRankMatrix(const int num_col_blocks,
  51                                               const int min_col_block_size,
  52                                               const int max_col_block_size,
  53                                               const double block_density) {
  54   // Create a random matrix
  55   BlockSparseMatrix::RandomMatrixOptions options;
  56   options.num_col_blocks = num_col_blocks;
  57   options.min_col_block_size = min_col_block_size;
  58   options.max_col_block_size = max_col_block_size;
  59
  60   options.num_row_blocks = 2 * num_col_blocks;
  61   options.min_row_block_size = 1;
  62   options.max_row_block_size = max_col_block_size;
  63   options.block_density = block_density;
  64   scoped_ptr<BlockSparseMatrix> random_matrix(
  65       BlockSparseMatrix::CreateRandomMatrix(options));
  66
  67   // Add a diagonal block sparse matrix to make it full rank.
  68   Vector diagonal = Vector::Ones(random_matrix->num_cols());
  69   scoped_ptr<BlockSparseMatrix> block_diagonal(
  70       BlockSparseMatrix::CreateDiagonalMatrix(
  71           diagonal.data(), random_matrix->block_structure()->cols));
  72   random_matrix->AppendRows(*block_diagonal);
  73   return random_matrix.release();
  74 }
  75
  76 bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs,
  77                              const Vector& rhs,
  78                              Vector* solution) {
  79   Matrix eigen_lhs;
  80   lhs.ToDenseMatrix(&eigen_lhs);
  81   if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) {
  82     Matrix full_lhs = eigen_lhs.selfadjointView<Eigen::Upper>();
  83     Eigen::LLT<Matrix, Eigen::Upper> llt =
  84         eigen_lhs.selfadjointView<Eigen::Upper>().llt();
  85     if (llt.info() != Eigen::Success) {
  86       return false;
  87     }
  88     *solution = llt.solve(rhs);
  89     return (llt.info() == Eigen::Success);
  90   }
  91
  92   Matrix full_lhs = eigen_lhs.selfadjointView<Eigen::Lower>();
  93   Eigen::LLT<Matrix, Eigen::Lower> llt =
  94       eigen_lhs.selfadjointView<Eigen::Lower>().llt();
  95   if (llt.info() != Eigen::Success) {
  96     return false;
  97   }
  98   *solution = llt.solve(rhs);
  99   return (llt.info() == Eigen::Success);
 100 }
 101
 102 void SparseCholeskySolverUnitTest(
 103     const SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
 104     const OrderingType ordering_type,
 105     const bool use_block_structure,
 106     const int num_blocks,
 107     const int min_block_size,
 108     const int max_block_size,
 109     const double block_density) {
 110   scoped_ptr<SparseCholesky> sparse_cholesky(SparseCholesky::Create(
 111       sparse_linear_algebra_library_type, ordering_type));
 112   const CompressedRowSparseMatrix::StorageType storage_type =
 113       sparse_cholesky->StorageType();
 114
 115   scoped_ptr<BlockSparseMatrix> m(CreateRandomFullRankMatrix(
 116       num_blocks, min_block_size, max_block_size, block_density));
 117   scoped_ptr<InnerProductComputer> inner_product_computer(
 118       InnerProductComputer::Create(*m, storage_type));
 119   inner_product_computer->Compute();
 120   CompressedRowSparseMatrix* lhs = inner_product_computer->mutable_result();
 121
 122   if (!use_block_structure) {
 123     lhs->mutable_row_blocks()->clear();
 124     lhs->mutable_col_blocks()->clear();
 125   }
 126
 127   Vector rhs = Vector::Random(lhs->num_rows());
 128   Vector expected(lhs->num_rows());
 129   Vector actual(lhs->num_rows());
 130
 131   EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected));
 132   std::string message;
 133   EXPECT_EQ(sparse_cholesky->FactorAndSolve(
 134                 lhs, rhs.data(), actual.data(), &message),
 135             LINEAR_SOLVER_SUCCESS);
 136   Matrix eigen_lhs;
 137   lhs->ToDenseMatrix(&eigen_lhs);
 138   EXPECT_NEAR((actual - expected).norm() / actual.norm(),
 139               0.0,
 140               std::numeric_limits<double>::epsilon() * 10)
 141       << "\n"
 142       << eigen_lhs;
 143 }
 144
 145 typedef ::testing::tuple<SparseLinearAlgebraLibraryType, OrderingType, bool>
 146     Param;
 147
 148 std::string ParamInfoToString(testing::TestParamInfo<Param> info) {
 149   Param param = info.param;
 150   std::stringstream ss;
 151   ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_"
 152      << (::testing::get<1>(param) == AMD ? "AMD" : "NATURAL") << "_"
 153      << (::testing::get<2>(param) ? "UseBlockStructure" : "NoBlockStructure");
 154   return ss.str();
 155 }
 156
 157 class SparseCholeskyTest : public ::testing::TestWithParam<Param> {};
 158
 159 TEST_P(SparseCholeskyTest, FactorAndSolve) {
 160   SetRandomState(2982);
 161   const int kMinNumBlocks = 1;
 162   const int kMaxNumBlocks = 10;
 163   const int kNumTrials = 10;
 164   const int kMinBlockSize = 1;
 165   const int kMaxBlockSize = 5;
 166
 167   for (int num_blocks = kMinNumBlocks; num_blocks < kMaxNumBlocks;
 168        ++num_blocks) {
 169     for (int trial = 0; trial < kNumTrials; ++trial) {
 170       const double block_density = std::max(0.1, RandDouble());
 171       Param param = GetParam();
 172       SparseCholeskySolverUnitTest(::testing::get<0>(param),
 173                                    ::testing::get<1>(param),
 174                                    ::testing::get<2>(param),
 175                                    num_blocks,
 176                                    kMinBlockSize,
 177                                    kMaxBlockSize,
 178                                    block_density);
 179     }
 180   }
 181 }
 182
 183 #ifndef CERES_NO_SUITESPARSE
 184 INSTANTIATE_TEST_CASE_P(SuiteSparseCholesky,
 185                         SparseCholeskyTest,
 186                         ::testing::Combine(::testing::Values(SUITE_SPARSE),
 187                                            ::testing::Values(AMD, NATURAL),
 188                                            ::testing::Values(true, false)),
 189                         ParamInfoToString);
 190 #endif
 191
 192 #ifndef CERES_NO_CXSPARSE
 193 INSTANTIATE_TEST_CASE_P(CXSparseCholesky,
 194                         SparseCholeskyTest,
 195                         ::testing::Combine(::testing::Values(CX_SPARSE),
 196                                            ::testing::Values(AMD, NATURAL),
 197                                            ::testing::Values(true, false)),
 198                         ParamInfoToString);
 199 #endif
 200
 201 #ifdef CERES_USE_EIGEN_SPARSE
 202 INSTANTIATE_TEST_CASE_P(EigenSparseCholesky,
 203                         SparseCholeskyTest,
 204                         ::testing::Combine(::testing::Values(EIGEN_SPARSE),
 205                                            ::testing::Values(AMD, NATURAL),
 206                                            ::testing::Values(true, false)),
 207                         ParamInfoToString);
 208 #endif
 209
 210 }  // namespace internal
 211 }  // namespace ceres