internal/ceres/schur_complement_solver_test.cc

   1 // Ceres Solver - A fast non-linear least squares minimizer
   2 // Copyright 2015 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/schur_complement_solver.h"
  32
  33 #include <cstddef>
  34
  35 #include "ceres/block_sparse_matrix.h"
  36 #include "ceres/block_structure.h"
  37 #include "ceres/casts.h"
  38 #include "ceres/detect_structure.h"
  39 #include "ceres/internal/scoped_ptr.h"
  40 #include "ceres/linear_least_squares_problems.h"
  41 #include "ceres/linear_solver.h"
  42 #include "ceres/triplet_sparse_matrix.h"
  43 #include "ceres/types.h"
  44 #include "glog/logging.h"
  45 #include "gtest/gtest.h"
  46
  47 namespace ceres {
  48 namespace internal {
  49
  50 class SchurComplementSolverTest : public ::testing::Test {
  51  protected:
  52   void SetUpFromProblemId(int problem_id) {
  53     scoped_ptr<LinearLeastSquaresProblem> problem(
  54         CreateLinearLeastSquaresProblemFromId(problem_id));
  55
  56     CHECK_NOTNULL(problem.get());
  57     A.reset(down_cast<BlockSparseMatrix*>(problem->A.release()));
  58     b.reset(problem->b.release());
  59     D.reset(problem->D.release());
  60
  61     num_cols = A->num_cols();
  62     num_rows = A->num_rows();
  63     num_eliminate_blocks = problem->num_eliminate_blocks;
  64
  65     x.resize(num_cols);
  66     sol.resize(num_cols);
  67     sol_d.resize(num_cols);
  68
  69     LinearSolver::Options options;
  70     options.type = DENSE_QR;
  71
  72     scoped_ptr<LinearSolver> qr(LinearSolver::Create(options));
  73
  74     TripletSparseMatrix triplet_A(A->num_rows(),
  75                                   A->num_cols(),
  76                                   A->num_nonzeros());
  77     A->ToTripletSparseMatrix(&triplet_A);
  78
  79     // Gold standard solutions using dense QR factorization.
  80     DenseSparseMatrix dense_A(triplet_A);
  81     qr->Solve(&dense_A, b.get(), LinearSolver::PerSolveOptions(), sol.data());
  82
  83     // Gold standard solution with appended diagonal.
  84     LinearSolver::PerSolveOptions per_solve_options;
  85     per_solve_options.D = D.get();
  86     qr->Solve(&dense_A, b.get(), per_solve_options, sol_d.data());
  87   }
  88
  89   void ComputeAndCompareSolutions(
  90       int problem_id,
  91       bool regularization,
  92       ceres::LinearSolverType linear_solver_type,
  93       ceres::DenseLinearAlgebraLibraryType dense_linear_algebra_library_type,
  94       ceres::SparseLinearAlgebraLibraryType sparse_linear_algebra_library_type,
  95       bool use_postordering) {
  96     SetUpFromProblemId(problem_id);
  97     LinearSolver::Options options;
  98     options.elimination_groups.push_back(num_eliminate_blocks);
  99     options.elimination_groups.push_back(
 100         A->block_structure()->cols.size() - num_eliminate_blocks);
 101     options.type = linear_solver_type;
 102     options.dense_linear_algebra_library_type =
 103         dense_linear_algebra_library_type;
 104     options.sparse_linear_algebra_library_type =
 105         sparse_linear_algebra_library_type;
 106     options.use_postordering = use_postordering;
 107     DetectStructure(*A->block_structure(),
 108                     num_eliminate_blocks,
 109                     &options.row_block_size,
 110                     &options.e_block_size,
 111                     &options.f_block_size);
 112
 113     scoped_ptr<LinearSolver> solver(LinearSolver::Create(options));
 114
 115     LinearSolver::PerSolveOptions per_solve_options;
 116     LinearSolver::Summary summary;
 117     if (regularization) {
 118       per_solve_options.D = D.get();
 119     }
 120
 121     summary = solver->Solve(A.get(), b.get(), per_solve_options, x.data());
 122     EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS);
 123
 124     if (regularization) {
 125
 126       ASSERT_NEAR((sol_d - x).norm() / num_cols, 0, 1e-10)
 127           << "Regularized Expected solution: " << sol_d.transpose()
 128           << " Actual solution: " << x.transpose();
 129     } else {
 130       ASSERT_NEAR((sol - x).norm() / num_cols, 0, 1e-10)
 131           << "Unregularized Expected solution: " << sol.transpose()
 132           << " Actual solution: " << x.transpose();
 133     }
 134   }
 135
 136   int num_rows;
 137   int num_cols;
 138   int num_eliminate_blocks;
 139
 140   scoped_ptr<BlockSparseMatrix> A;
 141   scoped_array<double> b;
 142   scoped_array<double> D;
 143   Vector x;
 144   Vector sol;
 145   Vector sol_d;
 146 };
 147
 148 // TODO(sameeragarwal): Refactor these using value parameterized tests.
 149 // TODO(sameeragarwal): More extensive tests using random matrices.
 150 TEST_F(SchurComplementSolverTest, DenseSchurWithEigenSmallProblem) {
 151   ComputeAndCompareSolutions(2, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 152   ComputeAndCompareSolutions(2, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 153 }
 154
 155 TEST_F(SchurComplementSolverTest, DenseSchurWithEigenLargeProblem) {
 156   ComputeAndCompareSolutions(3, false, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 157   ComputeAndCompareSolutions(3, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 158 }
 159
 160 TEST_F(SchurComplementSolverTest, DenseSchurWithEigenVaryingFBlockSize) {
 161   ComputeAndCompareSolutions(4, true, DENSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 162 }
 163
 164 #ifndef CERES_NO_LAPACK
 165 TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKSmallProblem) {
 166   ComputeAndCompareSolutions(2, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
 167   ComputeAndCompareSolutions(2, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
 168 }
 169
 170 TEST_F(SchurComplementSolverTest, DenseSchurWithLAPACKLargeProblem) {
 171   ComputeAndCompareSolutions(3, false, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
 172   ComputeAndCompareSolutions(3, true, DENSE_SCHUR, LAPACK, SUITE_SPARSE, true);
 173 }
 174 #endif
 175
 176 #ifndef CERES_NO_SUITESPARSE
 177 TEST_F(SchurComplementSolverTest,
 178        SparseSchurWithSuiteSparseSmallProblemNoPostOrdering) {
 179   ComputeAndCompareSolutions(
 180       2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
 181   ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
 182 }
 183
 184 TEST_F(SchurComplementSolverTest,
 185        SparseSchurWithSuiteSparseSmallProblemPostOrdering) {
 186   ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 187   ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 188 }
 189
 190 TEST_F(SchurComplementSolverTest,
 191        SparseSchurWithSuiteSparseLargeProblemNoPostOrdering) {
 192   ComputeAndCompareSolutions(
 193       3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
 194   ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, false);
 195 }
 196
 197 TEST_F(SchurComplementSolverTest,
 198        SparseSchurWithSuiteSparseLargeProblemPostOrdering) {
 199   ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 200   ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, SUITE_SPARSE, true);
 201 }
 202 #endif  // CERES_NO_SUITESPARSE
 203
 204 #ifndef CERES_NO_CXSPARSE
 205 TEST_F(SchurComplementSolverTest,
 206        SparseSchurWithCXSparseSmallProblem) {
 207   ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
 208   ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
 209 }
 210
 211 TEST_F(SchurComplementSolverTest,
 212        SparseSchurWithCXSparseLargeProblem) {
 213   ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
 214   ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, CX_SPARSE, true);
 215 }
 216 #endif  // CERES_NO_CXSPARSE
 217
 218 #ifdef CERES_USE_EIGEN_SPARSE
 219 TEST_F(SchurComplementSolverTest,
 220        SparseSchurWithEigenSparseSmallProblem) {
 221   ComputeAndCompareSolutions(2, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
 222   ComputeAndCompareSolutions(2, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
 223 }
 224
 225 TEST_F(SchurComplementSolverTest,
 226        SparseSchurWithEigenSparseLargeProblem) {
 227   ComputeAndCompareSolutions(3, false, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
 228   ComputeAndCompareSolutions(3, true, SPARSE_SCHUR, EIGEN, EIGEN_SPARSE, true);
 229 }
 230 #endif  // CERES_USE_EIGEN_SPARSE
 231
 232 }  // namespace internal
 233 }  // namespace ceres