1// Ceres Solver - A fast non-linear least squares minimizer 2// Copyright 2012 Google Inc. All rights reserved. 3// http://code.google.com/p/ceres-solver/ 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: keir@google.com (Keir Mierle) 30 31#ifndef CERES_INTERNAL_CGNR_LINEAR_OPERATOR_H_ 32#define CERES_INTERNAL_CGNR_LINEAR_OPERATOR_H_ 33 34#include <algorithm> 35#include "ceres/linear_operator.h" 36#include "ceres/internal/scoped_ptr.h" 37#include "ceres/internal/eigen.h" 38 39namespace ceres { 40namespace internal { 41 42class SparseMatrix; 43 44// A linear operator which takes a matrix A and a diagonal vector D and 45// performs products of the form 46// 47// (A^T A + D^T D)x 48// 49// This is used to implement iterative general sparse linear solving with 50// conjugate gradients, where A is the Jacobian and D is a regularizing 51// parameter. A brief proof that D^T D is the correct regularizer: 52// 53// Given a regularized least squares problem: 54// 55// min ||Ax - b||^2 + ||Dx||^2 56// x 57// 58// First expand into matrix notation: 59// 60// (Ax - b)^T (Ax - b) + xD^TDx 61// 62// Then multiply out to get: 63// 64// = xA^TAx - 2b^T Ax + b^Tb + xD^TDx 65// 66// Take the derivative: 67// 68// 0 = 2A^TAx - 2A^T b + 2 D^TDx 69// 0 = A^TAx - A^T b + D^TDx 70// 0 = (A^TA + D^TD)x - A^T b 71// 72// Thus, the symmetric system we need to solve for CGNR is 73// 74// Sx = z 75// 76// with S = A^TA + D^TD 77// and z = A^T b 78// 79// Note: This class is not thread safe, since it uses some temporary storage. 80class CgnrLinearOperator : public LinearOperator { 81 public: 82 CgnrLinearOperator(const LinearOperator& A, const double *D) 83 : A_(A), D_(D), z_(new double[A.num_rows()]) { 84 } 85 virtual ~CgnrLinearOperator() {} 86 87 virtual void RightMultiply(const double* x, double* y) const { 88 std::fill(z_.get(), z_.get() + A_.num_rows(), 0.0); 89 90 // z = Ax 91 A_.RightMultiply(x, z_.get()); 92 93 // y = y + Atz 94 A_.LeftMultiply(z_.get(), y); 95 96 // y = y + DtDx 97 if (D_ != NULL) { 98 int n = A_.num_cols(); 99 VectorRef(y, n).array() += ConstVectorRef(D_, n).array().square() * 100 ConstVectorRef(x, n).array(); 101 } 102 } 103 104 virtual void LeftMultiply(const double* x, double* y) const { 105 RightMultiply(x, y); 106 } 107 108 virtual int num_rows() const { return A_.num_cols(); } 109 virtual int num_cols() const { return A_.num_cols(); } 110 111 private: 112 const LinearOperator& A_; 113 const double* D_; 114 scoped_array<double> z_; 115}; 116 117} // namespace internal 118} // namespace ceres 119 120#endif // CERES_INTERNAL_CGNR_LINEAR_OPERATOR_H_ 121