1 // This is core/vnl/algo/vnl_ldl_cholesky.h 2 #ifndef vnl_ldl_cholesky_h_ 3 #define vnl_ldl_cholesky_h_ 4 //: 5 // \file 6 // \brief Updateable Cholesky decomposition: A=LDL' 7 // \author Tim Cootes 8 // \date 29 Mar 2006 9 10 #include <vnl/vnl_vector.h> 11 #include <vnl/vnl_matrix.h> 12 #include <vnl/algo/vnl_algo_export.h> 13 14 //: Updateable Cholesky decomposition: A=LDL' 15 // A class to hold the Cholesky decomposition of a positive definite 16 // symmetric matrix, using the form A=LDL', where L is lower triangular 17 // with ones on the leading diagonal, and D is a diagonal matrix. 18 // This differs from vnl_cholesky, which decomposes as A=LL', without 19 // a constraint on the diagonal. The advantage of the LDL' form is that 20 // it can be efficiently updated. Thus given the decomposition of A, 21 // we can compute the decomposition of (A+v.v') in O(n^2) time. 22 // 23 // This can be used to solve linear systems, compute determinants and inverses. 24 // 25 // To check that the decomposition can be used safely for solving a linear 26 // equation it is wise to construct with mode==estimate_condition and 27 // check that rcond()>sqrt(machine precision). If this is not the case 28 // it might be a good idea to use vnl_svd instead. 29 class VNL_ALGO_EXPORT vnl_ldl_cholesky 30 { 31 public: 32 //: Modes of computation. See constructor for details. 33 enum Operation { 34 quiet, 35 verbose, 36 estimate_condition 37 }; 38 39 //: Make cholesky decomposition of M optionally computing the reciprocal condition number. 40 vnl_ldl_cholesky(vnl_matrix<double> const& M, Operation mode = verbose); 41 ~vnl_ldl_cholesky() = default; 42 43 //: Solve LS problem M x = b 44 vnl_vector<double> solve(vnl_vector<double> const& b) const; 45 46 //: Solve LS problem M x = b 47 void solve(vnl_vector<double> const& b, vnl_vector<double>* x) const; 48 49 //: Solve equation of form Lx=y (in-place) 50 // x is overwritten with solution 51 void solve_lx(vnl_vector<double>& y); 52 53 //: Compute determinant 54 double determinant() const; 55 56 //: Compute rank-1 update, ie the decomposition of (M+v.v') 57 // If the initial state is the decomposition of M, then 58 // L and D are updated so that on exit LDL'=M+v.v' 59 void rank1_update(const vnl_vector<double>& v); 60 61 //: Multi-rank update, ie the decomposition of (M+W.W') 62 // If the initial state is the decomposition of M, then 63 // L and D are updated so that on exit LDL'=M+W.W' 64 void update(const vnl_matrix<double>& W); 65 66 //: Compute inverse. Not efficient. 67 // Note that you rarely need the inverse - backsubstitution 68 // is faster and less prone to rounding errors. 69 vnl_matrix<double> inverse() const; 70 71 //: Return lower-triangular factor. lower_triangle()72 const vnl_matrix<double>& lower_triangle() const { return L_; } 73 74 //: Return upper-triangular factor. upper_triangle()75 vnl_matrix<double> upper_triangle() const { return L_.transpose(); } 76 77 //: Return elements of diagonal matrix D in LDL' diagonal()78 const vnl_vector<double>& diagonal() const { return d_; } 79 80 //: Efficient computation of x' * inv(M) * x 81 // Useful when M is a covariance matrix! 82 double xt_m_inv_x(const vnl_vector<double>& x) const; 83 84 //: Efficient computation of x' * M * x 85 // Twice as fast as explicitly computing x' * M * x 86 double xt_m_x(const vnl_vector<double>& x) const; 87 88 //: A Success/failure flag rank_deficiency()89 int rank_deficiency() const { return num_dims_rank_def_; } 90 91 //: Return reciprocal condition number (smallest/largest singular values). 92 // As long as rcond()>sqrt(precision) the decomposition can be used for 93 // solving equations safely. 94 // Not calculated unless Operation mode at construction was estimate_condition. rcond()95 double rcond() const { return rcond_; } 96 97 //: Return computed nullvector. 98 // Not calculated unless Operation mode at construction was estimate_condition. nullvector()99 vnl_vector<double> & nullvector() { return nullvector_; } nullvector()100 vnl_vector<double> const& nullvector() const { return nullvector_; } 101 102 protected: 103 // Data Members-------------------------------------------------------------- 104 105 //: Lower triangular matrix 106 vnl_matrix<double> L_; 107 108 //: Elements of diagonal matrix 109 vnl_vector<double> d_; 110 111 //: 1/(condition number) 112 double rcond_; 113 long num_dims_rank_def_; 114 vnl_vector<double> nullvector_; 115 116 private: 117 //: Copy constructor - privatised to avoid it being used 118 vnl_ldl_cholesky(vnl_ldl_cholesky const & that) = delete; 119 //: Assignment operator - privatised to avoid it being used 120 vnl_ldl_cholesky& operator=(vnl_ldl_cholesky const & that) = delete; 121 122 //: Solve Mx=b, overwriting input vector with the solution. 123 // x points to beginning of an n-element vector containing b 124 // On exit, x[i] filled with solution vector. 125 void inplace_solve(double* x) const; 126 }; 127 128 #endif // vnl_ldl_cholesky_h_ 129