| // Ceres Solver - A fast non-linear least squares minimizer |
| // Copyright 2017 Google Inc. All rights reserved. |
| // http://ceres-solver.org/ |
| // |
| // Redistribution and use in source and binary forms, with or without |
| // modification, are permitted provided that the following conditions are met: |
| // |
| // * Redistributions of source code must retain the above copyright notice, |
| // this list of conditions and the following disclaimer. |
| // * Redistributions in binary form must reproduce the above copyright notice, |
| // this list of conditions and the following disclaimer in the documentation |
| // and/or other materials provided with the distribution. |
| // * Neither the name of Google Inc. nor the names of its contributors may be |
| // used to endorse or promote products derived from this software without |
| // specific prior written permission. |
| // |
| // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" |
| // AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE |
| // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE |
| // ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE |
| // LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR |
| // CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF |
| // SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS |
| // INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN |
| // CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) |
| // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
| // POSSIBILITY OF SUCH DAMAGE. |
| // |
| // Author: sameeragarwal@google.com (Sameer Agarwal) |
| |
| #include "ceres/subset_preconditioner.h" |
| #include "Eigen/Dense" |
| #include "Eigen/SparseCore" |
| #include "ceres/block_sparse_matrix.h" |
| #include "ceres/compressed_row_sparse_matrix.h" |
| #include "ceres/inner_product_computer.h" |
| #include "ceres/internal/eigen.h" |
| #include "ceres/internal/scoped_ptr.h" |
| #include "glog/logging.h" |
| #include "gtest/gtest.h" |
| |
| namespace ceres { |
| namespace internal { |
| |
| // TODO(sameeragarwal): Refactor the following two functions out of |
| // here and sparse_cholesky_test.cc into a more suitable place. |
| template <int UpLoType> |
| bool SolveLinearSystemUsingEigen(const Matrix& lhs, |
| const Vector rhs, |
| Vector* solution) { |
| Eigen::LLT<Matrix, UpLoType> llt = lhs.selfadjointView<UpLoType>().llt(); |
| if (llt.info() != Eigen::Success) { |
| return false; |
| } |
| *solution = llt.solve(rhs); |
| return (llt.info() == Eigen::Success); |
| } |
| |
| // Use Eigen's Dense Cholesky solver to compute the solution to a |
| // sparse linear system. |
| bool ComputeExpectedSolution(const CompressedRowSparseMatrix& lhs, |
| const Vector& rhs, |
| Vector* solution) { |
| Matrix dense_triangular_lhs; |
| lhs.ToDenseMatrix(&dense_triangular_lhs); |
| if (lhs.storage_type() == CompressedRowSparseMatrix::UPPER_TRIANGULAR) { |
| Matrix full_lhs = dense_triangular_lhs.selfadjointView<Eigen::Upper>(); |
| return SolveLinearSystemUsingEigen<Eigen::Upper>(full_lhs, rhs, solution); |
| } |
| return SolveLinearSystemUsingEigen<Eigen::Lower>( |
| dense_triangular_lhs, rhs, solution); |
| } |
| |
| typedef ::testing::tuple<SparseLinearAlgebraLibraryType, bool> Param; |
| |
| std::string ParamInfoToString(testing::TestParamInfo<Param> info) { |
| Param param = info.param; |
| std::stringstream ss; |
| ss << SparseLinearAlgebraLibraryTypeToString(::testing::get<0>(param)) << "_" |
| << (::testing::get<1>(param) ? "Diagonal" : "NoDiagonal"); |
| return ss.str(); |
| } |
| |
| class SubsetPreconditionerTest : public ::testing::TestWithParam<Param> { |
| protected: |
| virtual void SetUp() { |
| BlockSparseMatrix::RandomMatrixOptions options; |
| options.num_col_blocks = 4; |
| options.min_col_block_size = 1; |
| options.max_col_block_size = 4; |
| options.num_row_blocks = 8; |
| options.min_row_block_size = 1; |
| options.max_row_block_size = 4; |
| options.block_density = 0.9; |
| |
| m_.reset(BlockSparseMatrix::CreateRandomMatrix(options)); |
| start_row_block_ = m_->block_structure()->rows.size(); |
| |
| // Ensure that the bottom part of the matrix has the same column |
| // block structure. |
| options.col_blocks = m_->block_structure()->cols; |
| b_.reset(BlockSparseMatrix::CreateRandomMatrix(options)); |
| m_->AppendRows(*b_); |
| |
| // Create a Identity block diagonal matrix with the same column |
| // block structure. |
| diagonal_ = Vector::Ones(m_->num_cols()); |
| block_diagonal_.reset(BlockSparseMatrix::CreateDiagonalMatrix( |
| diagonal_.data(), b_->block_structure()->cols)); |
| |
| // Unconditionally add the block diagonal to the matrix b_, |
| // because either it is either part of b_ to make it full rank, or |
| // we pass the same diagonal matrix later as the parameter D. In |
| // either case the preconditioner matrix is b_' b + D'D. |
| b_->AppendRows(*block_diagonal_); |
| inner_product_computer_.reset(InnerProductComputer::Create( |
| *b_, CompressedRowSparseMatrix::UPPER_TRIANGULAR)); |
| inner_product_computer_->Compute(); |
| } |
| |
| scoped_ptr<BlockSparseMatrix> m_; |
| scoped_ptr<BlockSparseMatrix> b_; |
| scoped_ptr<BlockSparseMatrix> block_diagonal_; |
| scoped_ptr<InnerProductComputer> inner_product_computer_; |
| scoped_ptr<Preconditioner> preconditioner_; |
| Vector diagonal_; |
| int start_row_block_; |
| }; |
| |
| TEST_P(SubsetPreconditionerTest, foo) { |
| Param param = GetParam(); |
| Preconditioner::Options options; |
| options.subset_preconditioner_start_row_block = start_row_block_; |
| options.sparse_linear_algebra_library_type = ::testing::get<0>(param); |
| preconditioner_.reset(new SubsetPreconditioner(options, *m_)); |
| |
| const bool with_diagonal = ::testing::get<1>(param); |
| if (!with_diagonal) { |
| m_->AppendRows(*block_diagonal_); |
| } |
| |
| EXPECT_TRUE( |
| preconditioner_->Update(*m_, with_diagonal ? diagonal_.data() : NULL)); |
| |
| // Repeatedly apply the preconditioner to random vectors and check |
| // that the preconditioned value is the same as one obtained by |
| // solving the linear system directly. |
| for (int i = 0; i < 5; ++i) { |
| CompressedRowSparseMatrix* lhs = inner_product_computer_->mutable_result(); |
| Vector rhs = Vector::Random(lhs->num_rows()); |
| Vector expected(lhs->num_rows()); |
| EXPECT_TRUE(ComputeExpectedSolution(*lhs, rhs, &expected)); |
| |
| Vector actual(lhs->num_rows()); |
| preconditioner_->RightMultiply(rhs.data(), actual.data()); |
| |
| Matrix eigen_lhs; |
| lhs->ToDenseMatrix(&eigen_lhs); |
| EXPECT_NEAR((actual - expected).norm() / actual.norm(), |
| 0.0, |
| std::numeric_limits<double>::epsilon() * 10) |
| << "\n" |
| << eigen_lhs << "\n" |
| << expected.transpose() << "\n" |
| << actual.transpose(); |
| } |
| } |
| |
| #ifndef CERES_NO_SUITESPARSE |
| INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithSuiteSparse, |
| SubsetPreconditionerTest, |
| ::testing::Combine(::testing::Values(SUITE_SPARSE), |
| ::testing::Values(true, false)), |
| ParamInfoToString); |
| #endif |
| |
| #ifndef CERES_NO_CXSPARSE |
| INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithCXSparse, |
| SubsetPreconditionerTest, |
| ::testing::Combine(::testing::Values(CX_SPARSE), |
| ::testing::Values(true, false)), |
| ParamInfoToString); |
| #endif |
| |
| #ifdef CERES_USE_EIGEN_SPARSE |
| INSTANTIATE_TEST_CASE_P(SubsetPreconditionerWithEigenSparse, |
| SubsetPreconditionerTest, |
| ::testing::Combine(::testing::Values(EIGEN_SPARSE), |
| ::testing::Values(true, false)), |
| ParamInfoToString); |
| #endif |
| |
| } // namespace internal |
| } // namespace ceres |