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