diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt index cb284b0c6..e4dca2467 100644 --- a/internal/ceres/CMakeLists.txt +++ b/internal/ceres/CMakeLists.txt @@ -299,6 +299,7 @@ if (BUILD_TESTING AND GFLAGS) ceres_test(covariance) ceres_test(cubic_interpolation) ceres_test(detect_structure) + ceres_test(dense_linear_solver) ceres_test(dense_sparse_matrix) ceres_test(dynamic_autodiff_cost_function) ceres_test(dynamic_compressed_row_sparse_matrix) @@ -343,6 +344,7 @@ if (BUILD_TESTING AND GFLAGS) ceres_test(small_blas) ceres_test(solver) ceres_test(sparse_cholesky) + ceres_test(sparse_normal_cholesky_solver) # TODO(sameeragarwal): This test should ultimately be made # independent of SuiteSparse. @@ -356,7 +358,6 @@ if (BUILD_TESTING AND GFLAGS) ceres_test(triplet_sparse_matrix) ceres_test(trust_region_minimizer) ceres_test(trust_region_preprocessor) - ceres_test(unsymmetric_linear_solver) ceres_test(visibility) ceres_test(visibility_based_preconditioner) diff --git a/internal/ceres/dense_linear_solver_test.cc b/internal/ceres/dense_linear_solver_test.cc new file mode 100644 index 000000000..770a2d3bf --- /dev/null +++ b/internal/ceres/dense_linear_solver_test.cc @@ -0,0 +1,123 @@ +// 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/casts.h" +#include "ceres/internal/scoped_ptr.h" +#include "ceres/linear_least_squares_problems.h" +#include "ceres/linear_solver.h" +#include "ceres/triplet_sparse_matrix.h" +#include "ceres/types.h" +#include "glog/logging.h" +#include "gtest/gtest.h" + +namespace ceres { +namespace internal { + +typedef ::testing:: + tuple + Param; + +std::string ParamInfoToString(testing::TestParamInfo info) { + Param param = info.param; + std::stringstream ss; + ss << LinearSolverTypeToString(::testing::get<0>(param)) << "_" + << DenseLinearAlgebraLibraryTypeToString(::testing::get<1>(param)) << "_" + << (::testing::get<2>(param) ? "Regularized" : "Unregularized") << "_" + << ::testing::get<3>(param); + return ss.str(); +} + +class DenseLinearSolverTest : public ::testing::TestWithParam {}; + +TEST_P(DenseLinearSolverTest, _) { + Param param = GetParam(); + const bool regularized = testing::get<2>(param); + + scoped_ptr problem( + CreateLinearLeastSquaresProblemFromId(testing::get<3>(param))); + DenseSparseMatrix lhs(*down_cast(problem->A.get())); + + const int num_cols = lhs.num_cols(); + const int num_rows = lhs.num_rows(); + + Vector rhs = Vector::Zero(num_rows + num_cols); + rhs.head(num_rows) = ConstVectorRef(problem->b.get(), num_rows); + + LinearSolver::Options options; + options.type = ::testing::get<0>(param); + options.dense_linear_algebra_library_type = ::testing::get<1>(param); + scoped_ptr solver(LinearSolver::Create(options)); + + LinearSolver::PerSolveOptions per_solve_options; + if (regularized) { + per_solve_options.D = problem->D.get(); + } + + Vector solution(num_cols); + LinearSolver::Summary summary = + solver->Solve(&lhs, rhs.data(), per_solve_options, solution.data()); + EXPECT_EQ(summary.termination_type, LINEAR_SOLVER_SUCCESS); + + // If solving for the regularized solution, add the diagonal to the + // matrix. This makes subsequent computations simpler. + if (testing::get<2>(param)) { + lhs.AppendDiagonal(problem->D.get()); + }; + + Vector tmp = Vector::Zero(num_rows + num_cols); + lhs.RightMultiply(solution.data(), tmp.data()); + Vector actual_normal_rhs = Vector::Zero(num_cols); + lhs.LeftMultiply(tmp.data(), actual_normal_rhs.data()); + + Vector expected_normal_rhs = Vector::Zero(num_cols); + lhs.LeftMultiply(rhs.data(), expected_normal_rhs.data()); + const double residual = (expected_normal_rhs - actual_normal_rhs).norm() / + expected_normal_rhs.norm(); + + EXPECT_NEAR(residual, 0.0, 10 * std::numeric_limits::epsilon()); +} + +// TODO(sameeragarwal): Should we move away from hard coded linear +// least squares problem to randomly generated ones? +INSTANTIATE_TEST_CASE_P( + DenseLinearSolver, + DenseLinearSolverTest, + ::testing::Combine(::testing::Values(DENSE_QR, DENSE_NORMAL_CHOLESKY), +#ifdef CERES_NO_LAPACK + ::testing::Values(EIGEN), +#else + ::testing::Values(EIGEN, LAPACK), +#endif + ::testing::Values(true, false), + ::testing::Values(0, 1)), + ParamInfoToString); + +} // namespace internal +} // namespace ceres diff --git a/internal/ceres/unsymmetric_linear_solver_test.cc b/internal/ceres/sparse_normal_cholesky_solver_test.cc similarity index 68% rename from internal/ceres/unsymmetric_linear_solver_test.cc rename to internal/ceres/sparse_normal_cholesky_solver_test.cc index 401db1e50..c4a1c6972 100644 --- a/internal/ceres/unsymmetric_linear_solver_test.cc +++ b/internal/ceres/sparse_normal_cholesky_solver_test.cc @@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2015 Google Inc. All rights reserved. +// Copyright 2017 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -41,9 +41,15 @@ namespace ceres { namespace internal { -// TODO(sameeragarwal): Refactor and expand these tests. -class UnsymmetricLinearSolverTest : public ::testing::Test { - protected : +// TODO(sameeragarwal): This tests needs to be re-written, since +// SparseNormalCholeskySolver is a composition of two classes now, +// OuterProduct and SparseCholesky. +// +// So the test should exercise the composition, rather than the +// numerics of the solver, which are well covered by tests for those +// classes. +class SparseNormalCholeskyLinearSolverTest : public ::testing::Test { + protected: virtual void SetUp() { scoped_ptr problem( CreateLinearLeastSquaresProblemFromId(0)); @@ -65,51 +71,42 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { scoped_ptr transformed_A; - if (options.type == DENSE_QR || - options.type == DENSE_NORMAL_CHOLESKY) { - transformed_A.reset(new DenseSparseMatrix(*A_)); - } else if (options.type == SPARSE_NORMAL_CHOLESKY) { - CompressedRowSparseMatrix* crsm = - CompressedRowSparseMatrix::FromTripletSparseMatrix(*A_); - // Add row/column blocks structure. - for (int i = 0; i < A_->num_rows(); ++i) { - crsm->mutable_row_blocks()->push_back(1); - } - - for (int i = 0; i < A_->num_cols(); ++i) { - crsm->mutable_col_blocks()->push_back(1); - } - - // With all blocks of size 1, crsb_rows and crsb_cols are equivalent to - // rows and cols. - std::copy(crsm->rows(), crsm->rows() + crsm->num_rows() + 1, - std::back_inserter(*crsm->mutable_crsb_rows())); - - std::copy(crsm->cols(), crsm->cols() + crsm->num_nonzeros(), - std::back_inserter(*crsm->mutable_crsb_cols())); - - transformed_A.reset(crsm); - } else { - LOG(FATAL) << "Unknown linear solver : " << options.type; + CompressedRowSparseMatrix* crsm = + CompressedRowSparseMatrix::FromTripletSparseMatrix(*A_); + // Add row/column blocks structure. + for (int i = 0; i < A_->num_rows(); ++i) { + crsm->mutable_row_blocks()->push_back(1); } + for (int i = 0; i < A_->num_cols(); ++i) { + crsm->mutable_col_blocks()->push_back(1); + } + + // With all blocks of size 1, crsb_rows and crsb_cols are equivalent to + // rows and cols. + std::copy(crsm->rows(), + crsm->rows() + crsm->num_rows() + 1, + std::back_inserter(*crsm->mutable_crsb_rows())); + + std::copy(crsm->cols(), + crsm->cols() + crsm->num_nonzeros(), + std::back_inserter(*crsm->mutable_crsb_cols())); + + transformed_A.reset(crsm); + // Unregularized scoped_ptr solver(LinearSolver::Create(options)); - unregularized_solve_summary = - solver->Solve(transformed_A.get(), - b_.get(), - per_solve_options, - x_unregularized.data()); + unregularized_solve_summary = solver->Solve(transformed_A.get(), + b_.get(), + per_solve_options, + x_unregularized.data()); // Sparsity structure is changing, reset the solver. solver.reset(LinearSolver::Create(options)); // Regularized solution per_solve_options.D = D_.get(); - regularized_solve_summary = - solver->Solve(transformed_A.get(), - b_.get(), - per_solve_options, - x_regularized.data()); + regularized_solve_summary = solver->Solve( + transformed_A.get(), b_.get(), per_solve_options, x_regularized.data()); EXPECT_EQ(unregularized_solve_summary.termination_type, LINEAR_SOLVER_SUCCESS); @@ -117,8 +114,8 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { for (int i = 0; i < A_->num_cols(); ++i) { EXPECT_NEAR(sol_unregularized_[i], x_unregularized[i], 1e-8) << "\nExpected: " - << ConstVectorRef(sol_unregularized_.get(), - A_->num_cols()).transpose() + << ConstVectorRef(sol_unregularized_.get(), A_->num_cols()) + .transpose() << "\nActual: " << x_unregularized.transpose(); } @@ -139,38 +136,8 @@ class UnsymmetricLinearSolverTest : public ::testing::Test { scoped_array sol_regularized_; }; -TEST_F(UnsymmetricLinearSolverTest, EigenDenseQR) { - LinearSolver::Options options; - options.type = DENSE_QR; - options.dense_linear_algebra_library_type = EIGEN; - TestSolver(options); -} - -TEST_F(UnsymmetricLinearSolverTest, EigenDenseNormalCholesky) { - LinearSolver::Options options; - options.dense_linear_algebra_library_type = EIGEN; - options.type = DENSE_NORMAL_CHOLESKY; - TestSolver(options); -} - -#ifndef CERES_NO_LAPACK -TEST_F(UnsymmetricLinearSolverTest, LAPACKDenseQR) { - LinearSolver::Options options; - options.type = DENSE_QR; - options.dense_linear_algebra_library_type = LAPACK; - TestSolver(options); -} - -TEST_F(UnsymmetricLinearSolverTest, LAPACKDenseNormalCholesky) { - LinearSolver::Options options; - options.dense_linear_algebra_library_type = LAPACK; - options.type = DENSE_NORMAL_CHOLESKY; - TestSolver(options); -} -#endif - #ifndef CERES_NO_SUITESPARSE -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingSuiteSparsePreOrdering) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = SUITE_SPARSE; @@ -179,7 +146,7 @@ TEST_F(UnsymmetricLinearSolverTest, TestSolver(options); } -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingSuiteSparsePostOrdering) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = SUITE_SPARSE; @@ -188,7 +155,7 @@ TEST_F(UnsymmetricLinearSolverTest, TestSolver(options); } -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingSuiteSparseDynamicSparsity) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = SUITE_SPARSE; @@ -199,7 +166,7 @@ TEST_F(UnsymmetricLinearSolverTest, #endif #ifndef CERES_NO_CXSPARSE -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingCXSparsePreOrdering) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = CX_SPARSE; @@ -208,7 +175,7 @@ TEST_F(UnsymmetricLinearSolverTest, TestSolver(options); } -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingCXSparsePostOrdering) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = CX_SPARSE; @@ -217,7 +184,7 @@ TEST_F(UnsymmetricLinearSolverTest, TestSolver(options); } -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingCXSparseDynamicSparsity) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = CX_SPARSE; @@ -228,7 +195,7 @@ TEST_F(UnsymmetricLinearSolverTest, #endif #ifdef CERES_USE_EIGEN_SPARSE -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingEigenPreOrdering) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = EIGEN_SPARSE; @@ -237,7 +204,7 @@ TEST_F(UnsymmetricLinearSolverTest, TestSolver(options); } -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingEigenPostOrdering) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = EIGEN_SPARSE; @@ -246,7 +213,7 @@ TEST_F(UnsymmetricLinearSolverTest, TestSolver(options); } -TEST_F(UnsymmetricLinearSolverTest, +TEST_F(SparseNormalCholeskyLinearSolverTest, SparseNormalCholeskyUsingEigenDynamicSparsity) { LinearSolver::Options options; options.sparse_linear_algebra_library_type = EIGEN_SPARSE;