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 (at) 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 39 namespace ceres { 40 namespace internal { 41 42 class 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. 80 class 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
