mirror of
https://github.com/ceres-solver/ceres-solver.git
synced 2026-08-29 16:40:38 +08:00
Sameer Agarwal
08e60379ba
Despite its relative size, this is very significant change to Ceres. Why === Up till now, when the user chose SPARSE_NORMAL_CHOLESKY, the Jacobian was evaluated in a CompressedRowSparseMatrix, which was then use to compute the normal equations which were passed to a sparse linear algebra library for factorization. The reason to do this was because in the case of SuiteSparse, we were able to pass the Jacobian matrix directly without computing the normal equations and SuiteSparse/CHOLMOD did the normal equation computation. This turned out to be slow, so Cheng Wang implemented a high performance version of the matrix-matrix multiply to compute the normal equations, and all the sparse linear algebra libraries now are passed the normal equations. So that raises the question, as to what the best representation of the Jacobian which is suitable for the normal equation computation. Turns out BlockSparseMatrix is ideal. It brings two advantages. 1. Jacobian evaluation into a BlockSparseMatrix is considerably faster when using a BlockSparseMatrix than CompressedRowSparseMatrix. This is because we save on a bunch of memory copies. 2. To make the matrix multiplication fast and use the block structure Cheng Wang had to essentially make the CompressedRowSparseMatrix carry a bunch of sidecar information about the block sparsity, essentially making it behave like a BlockSparseMatrix. The resulting code had fairly complicated indexing and complicated the semantics of CompressedRowSparseMatrix. The new InnerProductComputer class does away with all that and once this CL goes in, I will be able to remove all that code and simplify the semantics of CompressedRowSparseMatrix. Changes ======= 1. Use InnerProductComputer in SparseNormalCholeskySolver. 2. Change the evaluator instantiated for SPARSE_NORMAL_CHOLESKY with static sparsity inside evaluator.cc 3. The former change necessitates that we change ProblemImpl::Evaluate to create the evaluate it needs on its own, because it was depending on passing "SPARSE_NORMAL_CHOLESKY" as linear solver type to the evaluator factor to get an Evaluator which can use CompressedRowSparseMatrix objects for storing the Jacobian. 4. Update the tests for SparseNormalCholeskySolver. 5. Separate out the tests for DynamicSparseNormalCholeskySolver into its own file. Change-Id: I2ef7ef8fbfbb4967d0c1ec2068c1c778248fdf5b
111 lines
4.1 KiB
C++
111 lines
4.1 KiB
C++
// 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/sparse_normal_cholesky_solver.h"
|
|
|
|
#include <algorithm>
|
|
#include <cstring>
|
|
#include <ctime>
|
|
|
|
#include "ceres/block_sparse_matrix.h"
|
|
#include "ceres/inner_product_computer.h"
|
|
#include "ceres/internal/eigen.h"
|
|
#include "ceres/internal/scoped_ptr.h"
|
|
#include "ceres/linear_solver.h"
|
|
#include "ceres/sparse_cholesky.h"
|
|
#include "ceres/triplet_sparse_matrix.h"
|
|
#include "ceres/types.h"
|
|
#include "ceres/wall_time.h"
|
|
|
|
namespace ceres {
|
|
namespace internal {
|
|
|
|
SparseNormalCholeskySolver::SparseNormalCholeskySolver(
|
|
const LinearSolver::Options& options)
|
|
: options_(options) {
|
|
sparse_cholesky_.reset(
|
|
SparseCholesky::Create(options_.sparse_linear_algebra_library_type,
|
|
options_.use_postordering ? AMD : NATURAL));
|
|
}
|
|
|
|
SparseNormalCholeskySolver::~SparseNormalCholeskySolver() {}
|
|
|
|
LinearSolver::Summary SparseNormalCholeskySolver::SolveImpl(
|
|
BlockSparseMatrix* A,
|
|
const double* b,
|
|
const LinearSolver::PerSolveOptions& per_solve_options,
|
|
double* x) {
|
|
EventLogger event_logger("SparseNormalCholeskySolver::Solve");
|
|
LinearSolver::Summary summary;
|
|
summary.num_iterations = 1;
|
|
summary.termination_type = LINEAR_SOLVER_SUCCESS;
|
|
summary.message = "Success.";
|
|
|
|
const int num_cols = A->num_cols();
|
|
VectorRef(x, num_cols).setZero();
|
|
A->LeftMultiply(b, x);
|
|
event_logger.AddEvent("Compute RHS");
|
|
|
|
if (per_solve_options.D != NULL) {
|
|
// Temporarily append a diagonal block to the A matrix, but undo
|
|
// it before returning the matrix to the user.
|
|
scoped_ptr<BlockSparseMatrix> regularizer;
|
|
regularizer.reset(BlockSparseMatrix::CreateDiagonalMatrix(
|
|
per_solve_options.D, A->block_structure()->cols));
|
|
event_logger.AddEvent("Diagonal");
|
|
A->AppendRows(*regularizer);
|
|
event_logger.AddEvent("Append");
|
|
}
|
|
event_logger.AddEvent("Append Rows");
|
|
|
|
if (inner_product_computer_.get() == NULL) {
|
|
inner_product_computer_.reset(
|
|
InnerProductComputer::Create(*A, sparse_cholesky_->StorageType()));
|
|
|
|
event_logger.AddEvent("InnerProductComputer::Create");
|
|
}
|
|
|
|
inner_product_computer_->Compute();
|
|
event_logger.AddEvent("InnerProductComputer::Compute");
|
|
|
|
// TODO(sameeragarwal):
|
|
|
|
if (per_solve_options.D != NULL) {
|
|
A->DeleteRowBlocks(A->block_structure()->cols.size());
|
|
}
|
|
summary.termination_type = sparse_cholesky_->FactorAndSolve(
|
|
inner_product_computer_->mutable_result(), x, x, &summary.message);
|
|
event_logger.AddEvent("Factor & Solve");
|
|
return summary;
|
|
}
|
|
|
|
} // namespace internal
|
|
} // namespace ceres
|