internal/ceres/inner_product_computer_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/inner_product_computer.h"
  32
  33 #include <numeric>
  34 #include "ceres/block_sparse_matrix.h"
  35 #include "ceres/internal/eigen.h"
  36 #include "ceres/internal/scoped_ptr.h"
  37 #include "ceres/random.h"
  38 #include "ceres/triplet_sparse_matrix.h"
  39 #include "glog/logging.h"
  40 #include "gtest/gtest.h"
  41
  42 #include "Eigen/SparseCore"
  43
  44 namespace ceres {
  45 namespace internal {
  46
  47 #define COMPUTE_AND_COMPARE                                                  \
  48   {                                                                          \
  49     inner_product_computer->Compute();                                       \
  50     CompressedRowSparseMatrix* actual_product_crsm =                         \
  51         inner_product_computer->mutable_result();                            \
  52     Matrix actual_inner_product =                                            \
  53         Eigen::MappedSparseMatrix<double, Eigen::ColMajor>(                  \
  54             actual_product_crsm->num_rows(),                                 \
  55             actual_product_crsm->num_rows(),                                 \
  56             actual_product_crsm->num_nonzeros(),                             \
  57             actual_product_crsm->mutable_rows(),                             \
  58             actual_product_crsm->mutable_cols(),                             \
  59             actual_product_crsm->mutable_values());                          \
  60     EXPECT_EQ(actual_inner_product.rows(), actual_inner_product.cols());     \
  61     EXPECT_EQ(expected_inner_product.rows(), expected_inner_product.cols()); \
  62     EXPECT_EQ(actual_inner_product.rows(), expected_inner_product.rows());   \
  63     Matrix expected_t, actual_t;                                             \
  64     if (actual_product_crsm->storage_type() ==                               \
  65         CompressedRowSparseMatrix::LOWER_TRIANGULAR) {                       \
  66       expected_t = expected_inner_product.triangularView<Eigen::Upper>();    \
  67       actual_t = actual_inner_product.triangularView<Eigen::Upper>();        \
  68     } else {                                                                 \
  69       expected_t = expected_inner_product.triangularView<Eigen::Lower>();    \
  70       actual_t = actual_inner_product.triangularView<Eigen::Lower>();        \
  71     }                                                                        \
  72     EXPECT_LE((expected_t - actual_t).norm() / actual_t.norm(),              \
  73               100 * std::numeric_limits<double>::epsilon())                  \
  74         << "expected: \n"                                                    \
  75         << expected_t << "\nactual: \n"                                      \
  76         << actual_t;                                                         \
  77   }
  78
  79 TEST(InnerProductComputer, NormalOperation) {
  80   // "Randomly generated seed."
  81   SetRandomState(29823);
  82   const int kMaxNumRowBlocks = 10;
  83   const int kMaxNumColBlocks = 10;
  84   const int kNumTrials = 10;
  85
  86   // Create a random matrix, compute its outer product using Eigen and
  87   // ComputeOuterProduct. Convert both matrices to dense matrices and
  88   // compare their upper triangular parts.
  89   for (int num_row_blocks = 1; num_row_blocks < kMaxNumRowBlocks;
  90        ++num_row_blocks) {
  91     for (int num_col_blocks = 1; num_col_blocks < kMaxNumColBlocks;
  92          ++num_col_blocks) {
  93       for (int trial = 0; trial < kNumTrials; ++trial) {
  94         BlockSparseMatrix::RandomMatrixOptions options;
  95         options.num_row_blocks = num_row_blocks;
  96         options.num_col_blocks = num_col_blocks;
  97         options.min_row_block_size = 1;
  98         options.max_row_block_size = 5;
  99         options.min_col_block_size = 1;
 100         options.max_col_block_size = 10;
 101         options.block_density = std::max(0.1, RandDouble());
 102
 103         VLOG(2) << "num row blocks: " << options.num_row_blocks;
 104         VLOG(2) << "num col blocks: " << options.num_col_blocks;
 105         VLOG(2) << "min row block size: " << options.min_row_block_size;
 106         VLOG(2) << "max row block size: " << options.max_row_block_size;
 107         VLOG(2) << "min col block size: " << options.min_col_block_size;
 108         VLOG(2) << "max col block size: " << options.max_col_block_size;
 109         VLOG(2) << "block density: " << options.block_density;
 110
 111         scoped_ptr<BlockSparseMatrix> random_matrix(
 112             BlockSparseMatrix::CreateRandomMatrix(options));
 113
 114         TripletSparseMatrix tsm(random_matrix->num_rows(),
 115                                 random_matrix->num_cols(),
 116                                 random_matrix->num_nonzeros());
 117         random_matrix->ToTripletSparseMatrix(&tsm);
 118         std::vector<Eigen::Triplet<double> > triplets;
 119         for (int i = 0; i < tsm.num_nonzeros(); ++i) {
 120           triplets.push_back(Eigen::Triplet<double>(
 121               tsm.rows()[i], tsm.cols()[i], tsm.values()[i]));
 122         }
 123         Eigen::SparseMatrix<double> eigen_random_matrix(
 124             random_matrix->num_rows(), random_matrix->num_cols());
 125         eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end());
 126         Matrix expected_inner_product =
 127             eigen_random_matrix.transpose() * eigen_random_matrix;
 128
 129         scoped_ptr<InnerProductComputer> inner_product_computer;
 130
 131         inner_product_computer.reset(InnerProductComputer::Create(
 132             *random_matrix, CompressedRowSparseMatrix::LOWER_TRIANGULAR));
 133         COMPUTE_AND_COMPARE;
 134         inner_product_computer.reset(InnerProductComputer::Create(
 135             *random_matrix, CompressedRowSparseMatrix::UPPER_TRIANGULAR));
 136         COMPUTE_AND_COMPARE;
 137
 138       }
 139     }
 140   }
 141 }
 142
 143
 144 TEST(InnerProductComputer, SubMatrix) {
 145   // "Randomly generated seed."
 146   SetRandomState(29823);
 147   const int kNumRowBlocks = 10;
 148   const int kNumColBlocks = 20;
 149   const int kNumTrials = 5;
 150
 151   // Create a random matrix, compute its outer product using Eigen and
 152   // ComputeInnerProductComputer. Convert both matrices to dense matrices and
 153   // compare their upper triangular parts.
 154   for (int trial = 0; trial < kNumTrials; ++trial) {
 155     BlockSparseMatrix::RandomMatrixOptions options;
 156     options.num_row_blocks = kNumRowBlocks;
 157     options.num_col_blocks = kNumColBlocks;
 158     options.min_row_block_size = 1;
 159     options.max_row_block_size = 5;
 160     options.min_col_block_size = 1;
 161     options.max_col_block_size = 10;
 162     options.block_density = std::max(0.1, RandDouble());
 163
 164     VLOG(2) << "num row blocks: " << options.num_row_blocks;
 165     VLOG(2) << "num col blocks: " << options.num_col_blocks;
 166     VLOG(2) << "min row block size: " << options.min_row_block_size;
 167     VLOG(2) << "max row block size: " << options.max_row_block_size;
 168     VLOG(2) << "min col block size: " << options.min_col_block_size;
 169     VLOG(2) << "max col block size: " << options.max_col_block_size;
 170     VLOG(2) << "block density: " << options.block_density;
 171
 172     scoped_ptr<BlockSparseMatrix> random_matrix(
 173         BlockSparseMatrix::CreateRandomMatrix(options));
 174
 175     const std::vector<CompressedRow>& row_blocks =
 176         random_matrix->block_structure()->rows;
 177     const int num_row_blocks = row_blocks.size();
 178
 179     for (int start_row_block = 0; start_row_block < num_row_blocks - 1;
 180          ++start_row_block) {
 181       for (int end_row_block = start_row_block + 1;
 182            end_row_block < num_row_blocks;
 183            ++end_row_block) {
 184         const int start_row = row_blocks[start_row_block].block.position;
 185         const int end_row = row_blocks[end_row_block].block.position;
 186
 187         TripletSparseMatrix tsm(random_matrix->num_rows(),
 188                                 random_matrix->num_cols(),
 189                                 random_matrix->num_nonzeros());
 190         random_matrix->ToTripletSparseMatrix(&tsm);
 191         std::vector<Eigen::Triplet<double> > triplets;
 192         for (int i = 0; i < tsm.num_nonzeros(); ++i) {
 193           if (tsm.rows()[i] >= start_row && tsm.rows()[i] < end_row) {
 194             triplets.push_back(Eigen::Triplet<double>(
 195                 tsm.rows()[i], tsm.cols()[i], tsm.values()[i]));
 196           }
 197         }
 198
 199         Eigen::SparseMatrix<double> eigen_random_matrix(
 200             random_matrix->num_rows(), random_matrix->num_cols());
 201         eigen_random_matrix.setFromTriplets(triplets.begin(), triplets.end());
 202
 203         Matrix expected_inner_product =
 204             eigen_random_matrix.transpose() * eigen_random_matrix;
 205
 206         scoped_ptr<InnerProductComputer> inner_product_computer;
 207         inner_product_computer.reset(InnerProductComputer::Create(
 208             *random_matrix,
 209             start_row_block,
 210             end_row_block,
 211             CompressedRowSparseMatrix::LOWER_TRIANGULAR));
 212         COMPUTE_AND_COMPARE;
 213         inner_product_computer.reset(InnerProductComputer::Create(
 214             *random_matrix,
 215             start_row_block,
 216             end_row_block,
 217             CompressedRowSparseMatrix::UPPER_TRIANGULAR));
 218         COMPUTE_AND_COMPARE;
 219
 220       }
 221     }
 222   }
 223 }
 224
 225 #undef COMPUTE_AND_COMPARE
 226
 227 }  // namespace internal
 228 }  // namespace ceres