1 // This is core/vnl/algo/vnl_cholesky.h 2 #ifndef vnl_cholesky_h_ 3 #define vnl_cholesky_h_ 4 //: 5 // \file 6 // \brief Decomposition of symmetric matrix 7 // \author Andrew W. Fitzgibbon, Oxford RRG 8 // \date 08 Dec 96 9 // 10 // \verbatim 11 // Modifications 12 // Peter Vanroose, Leuven, Apr 1998: added L() (return decomposition matrix) 13 // dac (Manchester) 26/03/2001: tidied up documentation 14 // Feb.2002 - Peter Vanroose - brief doxygen comment placed on single line 15 // \endverbatim 16 17 #include <vnl/vnl_vector.h> 18 #include <vnl/vnl_matrix.h> 19 #include <vnl/algo/vnl_algo_export.h> 20 21 //: Decomposition of symmetric matrix. 22 // A class to hold the Cholesky decomposition of a symmetric matrix and 23 // use that to solve linear systems, compute determinants and inverses. 24 // The cholesky decomposition decomposes symmetric A = L*L.transpose() 25 // where L is lower triangular 26 // 27 // To check that the decomposition can be used safely for solving a linear 28 // equation it is wise to construct with mode==estimate_condition and 29 // check that rcond()>sqrt(machine precision). If this is not the case 30 // it might be a good idea to use vnl_svd instead. 31 class VNL_ALGO_EXPORT vnl_cholesky 32 { 33 public: 34 //: Modes of computation. See constructor for details. 35 enum Operation { 36 quiet, 37 verbose, 38 estimate_condition 39 }; 40 41 //: Make cholesky decomposition of M optionally computing the reciprocal condition number. 42 vnl_cholesky(vnl_matrix<double> const& M, Operation mode = verbose); 43 ~vnl_cholesky() = default; 44 45 //: Solve LS problem M x = b 46 vnl_vector<double> solve(vnl_vector<double> const& b) const; 47 48 //: Solve LS problem M x = b 49 void solve(vnl_vector<double> const& b, vnl_vector<double>* x) const; 50 51 //: Compute determinant 52 double determinant() const; 53 54 // Compute inverse. Not efficient. 55 // It's broken, I don't have time to fix it. 56 // Mail awf@robots if you need it and I'll tell you as much as I can 57 // to fix it. 58 vnl_matrix<double> inverse() const; 59 60 //: Return lower-triangular factor. 61 vnl_matrix<double> lower_triangle() const; 62 63 //: Return upper-triangular factor. 64 vnl_matrix<double> upper_triangle() const; 65 66 //: Return the decomposition matrix L_badly_named_method()67 vnl_matrix<double> const& L_badly_named_method() const { return A_; } 68 69 //: A Success/failure flag rank_deficiency()70 int rank_deficiency() const { return num_dims_rank_def_; } 71 72 //: Return reciprocal condition number (smallest/largest singular values). 73 // As long as rcond()>sqrt(precision) the decomposition can be used for 74 // solving equations safely. 75 // Not calculated unless Operation mode at construction was estimate_condition. rcond()76 double rcond() const { return rcond_; } 77 78 //: Return computed nullvector. 79 // Not calculated unless Operation mode at construction was estimate_condition. nullvector()80 vnl_vector<double> & nullvector() { return nullvector_; } nullvector()81 vnl_vector<double> const& nullvector() const { return nullvector_; } 82 83 protected: 84 // Data Members-------------------------------------------------------------- 85 vnl_matrix<double> A_; 86 double rcond_; 87 long num_dims_rank_def_; 88 vnl_vector<double> nullvector_; 89 90 private: 91 //: Copy constructor - privatised to avoid it being used 92 vnl_cholesky(vnl_cholesky const & that) = delete; 93 //: Assignment operator - privatised to avoid it being used 94 vnl_cholesky& operator=(vnl_cholesky const & that) = delete; 95 }; 96 97 #endif // vnl_cholesky_h_ 98