internal/ceres/block_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/block_sparse_matrix.h"
  32
  33 #include <string>
  34 #include "ceres/casts.h"
  35 #include "ceres/internal/eigen.h"
  36 #include "ceres/internal/scoped_ptr.h"
  37 #include "ceres/linear_least_squares_problems.h"
  38 #include "ceres/triplet_sparse_matrix.h"
  39 #include "glog/logging.h"
  40 #include "gtest/gtest.h"
  41
  42 namespace ceres {
  43 namespace internal {
  44
  45 class BlockSparseMatrixTest : public ::testing::Test {
  46  protected :
  47   virtual void SetUp() {
  48     scoped_ptr<LinearLeastSquaresProblem> problem(
  49         CreateLinearLeastSquaresProblemFromId(2));
  50     CHECK_NOTNULL(problem.get());
  51     A_.reset(down_cast<BlockSparseMatrix*>(problem->A.release()));
  52
  53     problem.reset(CreateLinearLeastSquaresProblemFromId(1));
  54     CHECK_NOTNULL(problem.get());
  55     B_.reset(down_cast<TripletSparseMatrix*>(problem->A.release()));
  56
  57     CHECK_EQ(A_->num_rows(), B_->num_rows());
  58     CHECK_EQ(A_->num_cols(), B_->num_cols());
  59     CHECK_EQ(A_->num_nonzeros(), B_->num_nonzeros());
  60   }
  61
  62   scoped_ptr<BlockSparseMatrix> A_;
  63   scoped_ptr<TripletSparseMatrix> B_;
  64 };
  65
  66 TEST_F(BlockSparseMatrixTest, SetZeroTest) {
  67   A_->SetZero();
  68   EXPECT_EQ(13, A_->num_nonzeros());
  69 }
  70
  71 TEST_F(BlockSparseMatrixTest, RightMultiplyTest) {
  72   Vector y_a = Vector::Zero(A_->num_rows());
  73   Vector y_b = Vector::Zero(A_->num_rows());
  74   for (int i = 0; i < A_->num_cols(); ++i) {
  75     Vector x = Vector::Zero(A_->num_cols());
  76     x[i] = 1.0;
  77     A_->RightMultiply(x.data(), y_a.data());
  78     B_->RightMultiply(x.data(), y_b.data());
  79     EXPECT_LT((y_a - y_b).norm(), 1e-12);
  80   }
  81 }
  82
  83 TEST_F(BlockSparseMatrixTest, LeftMultiplyTest) {
  84   Vector y_a = Vector::Zero(A_->num_cols());
  85   Vector y_b = Vector::Zero(A_->num_cols());
  86   for (int i = 0; i < A_->num_rows(); ++i) {
  87     Vector x = Vector::Zero(A_->num_rows());
  88     x[i] = 1.0;
  89     A_->LeftMultiply(x.data(), y_a.data());
  90     B_->LeftMultiply(x.data(), y_b.data());
  91     EXPECT_LT((y_a - y_b).norm(), 1e-12);
  92   }
  93 }
  94
  95 TEST_F(BlockSparseMatrixTest, SquaredColumnNormTest) {
  96   Vector y_a = Vector::Zero(A_->num_cols());
  97   Vector y_b = Vector::Zero(A_->num_cols());
  98   A_->SquaredColumnNorm(y_a.data());
  99   B_->SquaredColumnNorm(y_b.data());
 100   EXPECT_LT((y_a - y_b).norm(), 1e-12);
 101 }
 102
 103 TEST_F(BlockSparseMatrixTest, ToDenseMatrixTest) {
 104   Matrix m_a;
 105   Matrix m_b;
 106   A_->ToDenseMatrix(&m_a);
 107   B_->ToDenseMatrix(&m_b);
 108   EXPECT_LT((m_a - m_b).norm(), 1e-12);
 109 }
 110
 111 TEST_F(BlockSparseMatrixTest, AppendRows) {
 112   scoped_ptr<LinearLeastSquaresProblem> problem(
 113       CreateLinearLeastSquaresProblemFromId(2));
 114   scoped_ptr<BlockSparseMatrix> m(
 115       down_cast<BlockSparseMatrix*>(problem->A.release()));
 116   A_->AppendRows(*m);
 117   EXPECT_EQ(A_->num_rows(), 2 * m->num_rows());
 118   EXPECT_EQ(A_->num_cols(), m->num_cols());
 119
 120   problem.reset(CreateLinearLeastSquaresProblemFromId(1));
 121   scoped_ptr<TripletSparseMatrix> m2(
 122       down_cast<TripletSparseMatrix*>(problem->A.release()));
 123   B_->AppendRows(*m2);
 124
 125   Vector y_a = Vector::Zero(A_->num_rows());
 126   Vector y_b = Vector::Zero(A_->num_rows());
 127   for (int i = 0; i < A_->num_cols(); ++i) {
 128     Vector x = Vector::Zero(A_->num_cols());
 129     x[i] = 1.0;
 130     y_a.setZero();
 131     y_b.setZero();
 132
 133     A_->RightMultiply(x.data(), y_a.data());
 134     B_->RightMultiply(x.data(), y_b.data());
 135     EXPECT_LT((y_a - y_b).norm(), 1e-12);
 136   }
 137 }
 138
 139 TEST_F(BlockSparseMatrixTest, AppendAndDeleteBlockDiagonalMatrix) {
 140   const std::vector<Block>& column_blocks = A_->block_structure()->cols;
 141   const int num_cols =
 142       column_blocks.back().size + column_blocks.back().position;
 143   Vector diagonal(num_cols);
 144   for (int i = 0; i < num_cols; ++i) {
 145     diagonal(i) = 2 * i * i + 1;
 146   }
 147   scoped_ptr<BlockSparseMatrix> appendage(
 148       BlockSparseMatrix::CreateDiagonalMatrix(diagonal.data(), column_blocks));
 149
 150   A_->AppendRows(*appendage);
 151   Vector y_a, y_b;
 152   y_a.resize(A_->num_rows());
 153   y_b.resize(A_->num_rows());
 154   for (int i = 0; i < A_->num_cols(); ++i) {
 155     Vector x = Vector::Zero(A_->num_cols());
 156     x[i] = 1.0;
 157     y_a.setZero();
 158     y_b.setZero();
 159
 160     A_->RightMultiply(x.data(), y_a.data());
 161     B_->RightMultiply(x.data(), y_b.data());
 162     EXPECT_LT((y_a.head(B_->num_rows()) - y_b.head(B_->num_rows())).norm(), 1e-12);
 163     Vector expected_tail = Vector::Zero(A_->num_cols());
 164     expected_tail(i) = diagonal(i);
 165     EXPECT_LT((y_a.tail(A_->num_cols()) - expected_tail).norm(), 1e-12);
 166   }
 167
 168
 169   A_->DeleteRowBlocks(column_blocks.size());
 170   EXPECT_EQ(A_->num_rows(), B_->num_rows());
 171   EXPECT_EQ(A_->num_cols(), B_->num_cols());
 172
 173   y_a.resize(A_->num_rows());
 174   y_b.resize(A_->num_rows());
 175   for (int i = 0; i < A_->num_cols(); ++i) {
 176     Vector x = Vector::Zero(A_->num_cols());
 177     x[i] = 1.0;
 178     y_a.setZero();
 179     y_b.setZero();
 180
 181     A_->RightMultiply(x.data(), y_a.data());
 182     B_->RightMultiply(x.data(), y_b.data());
 183     EXPECT_LT((y_a - y_b).norm(), 1e-12);
 184   }
 185 }
 186
 187 TEST(BlockSparseMatrix, CreateDiagonalMatrix) {
 188   std::vector<Block> column_blocks;
 189   column_blocks.push_back(Block(2, 0));
 190   column_blocks.push_back(Block(1, 2));
 191   column_blocks.push_back(Block(3, 3));
 192   const int num_cols =
 193       column_blocks.back().size + column_blocks.back().position;
 194   Vector diagonal(num_cols);
 195   for (int i = 0; i < num_cols; ++i) {
 196     diagonal(i) = 2 * i * i + 1;
 197   }
 198
 199   scoped_ptr<BlockSparseMatrix> m(
 200       BlockSparseMatrix::CreateDiagonalMatrix(diagonal.data(), column_blocks));
 201   const CompressedRowBlockStructure* bs = m->block_structure();
 202   EXPECT_EQ(bs->cols.size(), column_blocks.size());
 203   for (int i = 0; i < column_blocks.size(); ++i) {
 204     EXPECT_EQ(bs->cols[i].size, column_blocks[i].size);
 205     EXPECT_EQ(bs->cols[i].position, column_blocks[i].position);
 206   }
 207   EXPECT_EQ(m->num_rows(), m->num_cols());
 208   Vector x = Vector::Ones(num_cols);
 209   Vector y = Vector::Zero(num_cols);
 210   m->RightMultiply(x.data(), y.data());
 211   for (int i = 0; i < num_cols; ++i) {
 212     EXPECT_NEAR(y[i], diagonal[i], std::numeric_limits<double>::epsilon());
 213   }
 214 }
 215
 216
 217 }  // namespace internal
 218 }  // namespace ceres