internal/ceres/compressed_row_sparse_matrix_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/compressed_row_sparse_matrix.h"
  32
  33 #include <numeric>
  34 #include "ceres/casts.h"
  35 #include "ceres/crs_matrix.h"
  36 #include "ceres/internal/eigen.h"
  37 #include "ceres/internal/scoped_ptr.h"
  38 #include "ceres/linear_least_squares_problems.h"
  39 #include "ceres/random.h"
  40 #include "ceres/triplet_sparse_matrix.h"
  41 #include "glog/logging.h"
  42 #include "gtest/gtest.h"
  43
  44 #include "Eigen/SparseCore"
  45
  46 namespace ceres {
  47 namespace internal {
  48
  49 using std::vector;
  50
  51 void CompareMatrices(const SparseMatrix* a, const SparseMatrix* b) {
  52   EXPECT_EQ(a->num_rows(), b->num_rows());
  53   EXPECT_EQ(a->num_cols(), b->num_cols());
  54
  55   int num_rows = a->num_rows();
  56   int num_cols = a->num_cols();
  57
  58   for (int i = 0; i < num_cols; ++i) {
  59     Vector x = Vector::Zero(num_cols);
  60     x(i) = 1.0;
  61
  62     Vector y_a = Vector::Zero(num_rows);
  63     Vector y_b = Vector::Zero(num_rows);
  64
  65     a->RightMultiply(x.data(), y_a.data());
  66     b->RightMultiply(x.data(), y_b.data());
  67
  68     EXPECT_EQ((y_a - y_b).norm(), 0);
  69   }
  70 }
  71
  72 class CompressedRowSparseMatrixTest : public ::testing::Test {
  73  protected:
  74   virtual void SetUp() {
  75     scoped_ptr<LinearLeastSquaresProblem> problem(
  76         CreateLinearLeastSquaresProblemFromId(1));
  77
  78     CHECK_NOTNULL(problem.get());
  79
  80     tsm.reset(down_cast<TripletSparseMatrix*>(problem->A.release()));
  81     crsm.reset(CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
  82
  83     num_rows = tsm->num_rows();
  84     num_cols = tsm->num_cols();
  85
  86     vector<int>* row_blocks = crsm->mutable_row_blocks();
  87     row_blocks->resize(num_rows);
  88     std::fill(row_blocks->begin(), row_blocks->end(), 1);
  89
  90     vector<int>* col_blocks = crsm->mutable_col_blocks();
  91     col_blocks->resize(num_cols);
  92     std::fill(col_blocks->begin(), col_blocks->end(), 1);
  93   }
  94
  95   int num_rows;
  96   int num_cols;
  97
  98   scoped_ptr<TripletSparseMatrix> tsm;
  99   scoped_ptr<CompressedRowSparseMatrix> crsm;
 100 };
 101
 102 TEST_F(CompressedRowSparseMatrixTest, RightMultiply) {
 103   CompareMatrices(tsm.get(), crsm.get());
 104 }
 105
 106 TEST_F(CompressedRowSparseMatrixTest, LeftMultiply) {
 107   for (int i = 0; i < num_rows; ++i) {
 108     Vector a = Vector::Zero(num_rows);
 109     a(i) = 1.0;
 110
 111     Vector b1 = Vector::Zero(num_cols);
 112     Vector b2 = Vector::Zero(num_cols);
 113
 114     tsm->LeftMultiply(a.data(), b1.data());
 115     crsm->LeftMultiply(a.data(), b2.data());
 116
 117     EXPECT_EQ((b1 - b2).norm(), 0);
 118   }
 119 }
 120
 121 TEST_F(CompressedRowSparseMatrixTest, ColumnNorm) {
 122   Vector b1 = Vector::Zero(num_cols);
 123   Vector b2 = Vector::Zero(num_cols);
 124
 125   tsm->SquaredColumnNorm(b1.data());
 126   crsm->SquaredColumnNorm(b2.data());
 127
 128   EXPECT_EQ((b1 - b2).norm(), 0);
 129 }
 130
 131 TEST_F(CompressedRowSparseMatrixTest, Scale) {
 132   Vector scale(num_cols);
 133   for (int i = 0; i < num_cols; ++i) {
 134     scale(i) = i + 1;
 135   }
 136
 137   tsm->ScaleColumns(scale.data());
 138   crsm->ScaleColumns(scale.data());
 139   CompareMatrices(tsm.get(), crsm.get());
 140 }
 141
 142 TEST_F(CompressedRowSparseMatrixTest, DeleteRows) {
 143   // Clear the row and column blocks as these are purely scalar tests.
 144   crsm->mutable_row_blocks()->clear();
 145   crsm->mutable_col_blocks()->clear();
 146
 147   for (int i = 0; i < num_rows; ++i) {
 148     tsm->Resize(num_rows - i, num_cols);
 149     crsm->DeleteRows(crsm->num_rows() - tsm->num_rows());
 150     CompareMatrices(tsm.get(), crsm.get());
 151   }
 152 }
 153
 154 TEST_F(CompressedRowSparseMatrixTest, AppendRows) {
 155   // Clear the row and column blocks as these are purely scalar tests.
 156   crsm->mutable_row_blocks()->clear();
 157   crsm->mutable_col_blocks()->clear();
 158
 159   for (int i = 0; i < num_rows; ++i) {
 160     TripletSparseMatrix tsm_appendage(*tsm);
 161     tsm_appendage.Resize(i, num_cols);
 162
 163     tsm->AppendRows(tsm_appendage);
 164     scoped_ptr<CompressedRowSparseMatrix> crsm_appendage(
 165         CompressedRowSparseMatrix::FromTripletSparseMatrix(tsm_appendage));
 166
 167     crsm->AppendRows(*crsm_appendage);
 168     CompareMatrices(tsm.get(), crsm.get());
 169   }
 170 }
 171
 172 TEST_F(CompressedRowSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) {
 173   int num_diagonal_rows = crsm->num_cols();
 174
 175   scoped_array<double> diagonal(new double[num_diagonal_rows]);
 176   for (int i = 0; i < num_diagonal_rows; ++i) {
 177     diagonal[i] = i;
 178   }
 179
 180   vector<int> row_and_column_blocks;
 181   row_and_column_blocks.push_back(1);
 182   row_and_column_blocks.push_back(2);
 183   row_and_column_blocks.push_back(2);
 184
 185   const vector<int> pre_row_blocks = crsm->row_blocks();
 186   const vector<int> pre_col_blocks = crsm->col_blocks();
 187
 188   scoped_ptr<CompressedRowSparseMatrix> appendage(
 189       CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(
 190           diagonal.get(), row_and_column_blocks));
 191   LOG(INFO) << appendage->row_blocks().size();
 192
 193   crsm->AppendRows(*appendage);
 194
 195   const vector<int> post_row_blocks = crsm->row_blocks();
 196   const vector<int> post_col_blocks = crsm->col_blocks();
 197
 198   vector<int> expected_row_blocks = pre_row_blocks;
 199   expected_row_blocks.insert(expected_row_blocks.end(),
 200                              row_and_column_blocks.begin(),
 201                              row_and_column_blocks.end());
 202
 203   vector<int> expected_col_blocks = pre_col_blocks;
 204
 205   EXPECT_EQ(expected_row_blocks, crsm->row_blocks());
 206   EXPECT_EQ(expected_col_blocks, crsm->col_blocks());
 207
 208   crsm->DeleteRows(num_diagonal_rows);
 209   EXPECT_EQ(crsm->row_blocks(), pre_row_blocks);
 210   EXPECT_EQ(crsm->col_blocks(), pre_col_blocks);
 211 }
 212
 213 TEST_F(CompressedRowSparseMatrixTest, ToDenseMatrix) {
 214   Matrix tsm_dense;
 215   Matrix crsm_dense;
 216
 217   tsm->ToDenseMatrix(&tsm_dense);
 218   crsm->ToDenseMatrix(&crsm_dense);
 219
 220   EXPECT_EQ((tsm_dense - crsm_dense).norm(), 0.0);
 221 }
 222
 223 TEST_F(CompressedRowSparseMatrixTest, ToCRSMatrix) {
 224   CRSMatrix crs_matrix;
 225   crsm->ToCRSMatrix(&crs_matrix);
 226   EXPECT_EQ(crsm->num_rows(), crs_matrix.num_rows);
 227   EXPECT_EQ(crsm->num_cols(), crs_matrix.num_cols);
 228   EXPECT_EQ(crsm->num_rows() + 1, crs_matrix.rows.size());
 229   EXPECT_EQ(crsm->num_nonzeros(), crs_matrix.cols.size());
 230   EXPECT_EQ(crsm->num_nonzeros(), crs_matrix.values.size());
 231
 232   for (int i = 0; i < crsm->num_rows() + 1; ++i) {
 233     EXPECT_EQ(crsm->rows()[i], crs_matrix.rows[i]);
 234   }
 235
 236   for (int i = 0; i < crsm->num_nonzeros(); ++i) {
 237     EXPECT_EQ(crsm->cols()[i], crs_matrix.cols[i]);
 238     EXPECT_EQ(crsm->values()[i], crs_matrix.values[i]);
 239   }
 240 }
 241
 242 TEST(CompressedRowSparseMatrix, CreateBlockDiagonalMatrix) {
 243   vector<int> blocks;
 244   blocks.push_back(1);
 245   blocks.push_back(2);
 246   blocks.push_back(2);
 247
 248   Vector diagonal(5);
 249   for (int i = 0; i < 5; ++i) {
 250     diagonal(i) = i + 1;
 251   }
 252
 253   scoped_ptr<CompressedRowSparseMatrix> matrix(
 254       CompressedRowSparseMatrix::CreateBlockDiagonalMatrix(diagonal.data(),
 255                                                            blocks));
 256
 257   EXPECT_EQ(matrix->num_rows(), 5);
 258   EXPECT_EQ(matrix->num_cols(), 5);
 259   EXPECT_EQ(matrix->num_nonzeros(), 9);
 260   EXPECT_EQ(blocks, matrix->row_blocks());
 261   EXPECT_EQ(blocks, matrix->col_blocks());
 262
 263   Vector x(5);
 264   Vector y(5);
 265
 266   x.setOnes();
 267   y.setZero();
 268   matrix->RightMultiply(x.data(), y.data());
 269   for (int i = 0; i < diagonal.size(); ++i) {
 270     EXPECT_EQ(y[i], diagonal[i]);
 271   }
 272
 273   y.setZero();
 274   matrix->LeftMultiply(x.data(), y.data());
 275   for (int i = 0; i < diagonal.size(); ++i) {
 276     EXPECT_EQ(y[i], diagonal[i]);
 277   }
 278
 279   Matrix dense;
 280   matrix->ToDenseMatrix(&dense);
 281   EXPECT_EQ((dense.diagonal() - diagonal).norm(), 0.0);
 282 }
 283
 284 TEST(CompressedRowSparseMatrix, Transpose) {
 285   //  0  1  0  2  3  0
 286   //  4  6  7  0  0  8
 287   //  9 10  0 11 12  0
 288   // 13  0 14 15  9  0
 289   //  0 16 17  0  0  0
 290
 291   // Block structure:
 292   //  A  A  A  A  B  B
 293   //  A  A  A  A  B  B
 294   //  A  A  A  A  B  B
 295   //  C  C  C  C  D  D
 296   //  C  C  C  C  D  D
 297   //  C  C  C  C  D  D
 298
 299   CompressedRowSparseMatrix matrix(5, 6, 30);
 300   int* rows = matrix.mutable_rows();
 301   int* cols = matrix.mutable_cols();
 302   double* values = matrix.mutable_values();
 303   matrix.mutable_row_blocks()->push_back(3);
 304   matrix.mutable_row_blocks()->push_back(3);
 305   matrix.mutable_col_blocks()->push_back(4);
 306   matrix.mutable_col_blocks()->push_back(2);
 307
 308   rows[0] = 0;
 309   cols[0] = 1;
 310   cols[1] = 3;
 311   cols[2] = 4;
 312
 313   rows[1] = 3;
 314   cols[3] = 0;
 315   cols[4] = 1;
 316   cols[5] = 2;
 317   cols[6] = 5;
 318
 319   rows[2] = 7;
 320   cols[7] = 0;
 321   cols[8] = 1;
 322   cols[9] = 3;
 323   cols[10] = 4;
 324
 325   rows[3] = 11;
 326   cols[11] = 0;
 327   cols[12] = 2;
 328   cols[13] = 3;
 329   cols[14] = 4;
 330
 331   rows[4] = 15;
 332   cols[15] = 1;
 333   cols[16] = 2;
 334   rows[5] = 17;
 335
 336   std::copy(values, values + 17, cols);
 337
 338   scoped_ptr<CompressedRowSparseMatrix> transpose(matrix.Transpose());
 339
 340   ASSERT_EQ(transpose->row_blocks().size(), matrix.col_blocks().size());
 341   for (int i = 0; i < transpose->row_blocks().size(); ++i) {
 342     EXPECT_EQ(transpose->row_blocks()[i], matrix.col_blocks()[i]);
 343   }
 344
 345   ASSERT_EQ(transpose->col_blocks().size(), matrix.row_blocks().size());
 346   for (int i = 0; i < transpose->col_blocks().size(); ++i) {
 347     EXPECT_EQ(transpose->col_blocks()[i], matrix.row_blocks()[i]);
 348   }
 349
 350   Matrix dense_matrix;
 351   matrix.ToDenseMatrix(&dense_matrix);
 352
 353   Matrix dense_transpose;
 354   transpose->ToDenseMatrix(&dense_transpose);
 355   EXPECT_NEAR((dense_matrix - dense_transpose.transpose()).norm(), 0.0, 1e-14);
 356 }
 357
 358 TEST(CompressedRowSparseMatrix, FromTripletSparseMatrix) {
 359   TripletSparseMatrix::RandomMatrixOptions options;
 360   options.num_rows = 5;
 361   options.num_cols = 7;
 362   options.density = 0.5;
 363
 364   const int kNumTrials = 10;
 365   for (int i = 0; i < kNumTrials; ++i) {
 366     scoped_ptr<TripletSparseMatrix> tsm(
 367         TripletSparseMatrix::CreateRandomMatrix(options));
 368     scoped_ptr<CompressedRowSparseMatrix> crsm(
 369         CompressedRowSparseMatrix::FromTripletSparseMatrix(*tsm));
 370
 371     Matrix expected;
 372     tsm->ToDenseMatrix(&expected);
 373     Matrix actual;
 374     crsm->ToDenseMatrix(&actual);
 375     EXPECT_NEAR((expected - actual).norm() / actual.norm(),
 376                 0.0,
 377                 std::numeric_limits<double>::epsilon())
 378         << "\nexpected: \n"
 379         << expected << "\nactual: \n"
 380         << actual;
 381   }
 382 }
 383
 384 TEST(CompressedRowSparseMatrix, FromTripletSparseMatrixTransposed) {
 385   TripletSparseMatrix::RandomMatrixOptions options;
 386   options.num_rows = 5;
 387   options.num_cols = 7;
 388   options.density = 0.5;
 389
 390   const int kNumTrials = 10;
 391   for (int i = 0; i < kNumTrials; ++i) {
 392     scoped_ptr<TripletSparseMatrix> tsm(
 393         TripletSparseMatrix::CreateRandomMatrix(options));
 394     scoped_ptr<CompressedRowSparseMatrix> crsm(
 395         CompressedRowSparseMatrix::FromTripletSparseMatrixTransposed(*tsm));
 396
 397     Matrix tmp;
 398     tsm->ToDenseMatrix(&tmp);
 399     Matrix expected = tmp.transpose();
 400     Matrix actual;
 401     crsm->ToDenseMatrix(&actual);
 402     EXPECT_NEAR((expected - actual).norm() / actual.norm(),
 403                 0.0,
 404                 std::numeric_limits<double>::epsilon())
 405         << "\nexpected: \n"
 406         << expected << "\nactual: \n"
 407         << actual;
 408   }
 409 }
 410
 411 // TODO(sameeragarwal) Add tests for the random matrix creation methods.
 412
 413 }  // namespace internal
 414 }  // namespace ceres