cholesky.cpp
Go to the documentation of this file.
1 
8 
9 
10 // Copyright (C) 1991,2,3,4: R B Davies
11 
12 #define WANT_MATH
13 //#define WANT_STREAM
14 
15 #include "include.h"
16 
17 #include "newmat.h"
18 #include "newmatrm.h"
19 
20 #ifdef use_namespace
21 namespace NEWMAT {
22 #endif
23 
24 #ifdef DO_REPORT
25 #define REPORT { static ExeCounter ExeCount(__LINE__,14); ++ExeCount; }
26 #else
27 #define REPORT {}
28 #endif
29 
30 /********* Cholesky decomposition of a positive definite matrix *************/
31 
32 // Suppose S is symmetrix and positive definite. Then there exists a unique
33 // lower triangular matrix L such that L L.t() = S;
34 
35 
37 {
38  REPORT
39  Tracer trace("Cholesky");
40  int nr = S.Nrows();
42  Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
43  for (int i=0; i<nr; i++)
44  {
45  Real* tj = t; Real sum; int k;
46  for (int j=0; j<i; j++)
47  {
48  Real* tk = ti; sum = 0.0; k = j;
49  while (k--) { sum += *tj++ * *tk++; }
50  *tk = (*s++ - sum) / *tj++;
51  }
52  sum = 0.0; k = i;
53  while (k--) { sum += square(*ti++); }
54  Real d = *s++ - sum;
55  if (d<=0.0) Throw(NPDException(S));
56  *ti++ = sqrt(d);
57  }
58  T.release(); return T.for_return();
59 }
60 
62 {
63  REPORT
64  Tracer trace("Band-Cholesky");
65  int nr = S.Nrows(); int m = S.lower_val;
66  LowerBandMatrix T(nr,m);
67  Real* s = S.Store(); Real* t = T.Store(); Real* ti = t;
68 
69  for (int i=0; i<nr; i++)
70  {
71  Real* tj = t; Real sum; int l;
72  if (i<m) { REPORT l = m-i; s += l; ti += l; l = i; }
73  else { REPORT t += (m+1); l = m; }
74 
75  for (int j=0; j<l; j++)
76  {
77  Real* tk = ti; sum = 0.0; int k = j; tj += (m-j);
78  while (k--) { sum += *tj++ * *tk++; }
79  *tk = (*s++ - sum) / *tj++;
80  }
81  sum = 0.0;
82  while (l--) { sum += square(*ti++); }
83  Real d = *s++ - sum;
84  if (d<=0.0) Throw(NPDException(S));
85  *ti++ = sqrt(d);
86  }
87 
88  T.release(); return T.for_return();
89 }
90 
91 
92 
93 
94 // Contributed by Nick Bennett of Schlumberger-Doll Research; modified by RBD
95 
96 // The enclosed routines can be used to update the Cholesky decomposition of
97 // a positive definite symmetric matrix. A good reference for this routines
98 // can be found in
99 // LINPACK User's Guide, Chapter 10, Dongarra et. al., SIAM, Philadelphia, 1979
100 
101 // produces the Cholesky decomposition of A + x.t() * x where A = chol.t() * chol
103 {
104  int nc = chol.Nrows();
105  ColumnVector cGivens(nc); cGivens = 0.0;
106  ColumnVector sGivens(nc); sGivens = 0.0;
107 
108  for(int j = 1; j <= nc; ++j) // process the jth column of chol
109  {
110  // apply the previous Givens rotations k = 1,...,j-1 to column j
111  for(int k = 1; k < j; ++k)
112  GivensRotation(cGivens(k), sGivens(k), chol(k,j), x(j));
113 
114  // determine the jth Given's rotation
115  pythag(chol(j,j), x(j), cGivens(j), sGivens(j));
116 
117  // apply the jth Given's rotation
118  {
119  Real tmp0 = cGivens(j) * chol(j,j) + sGivens(j) * x(j);
120  chol(j,j) = tmp0; x(j) = 0.0;
121  }
122 
123  }
124 
125 }
126 
127 
128 // produces the Cholesky decomposition of A - x.t() * x where A = chol.t() * chol
130 {
131  int nRC = chol.Nrows();
132 
133  // solve R^T a = x
134  LowerTriangularMatrix L = chol.t();
135  ColumnVector a(nRC); a = 0.0;
136  int i, j;
137 
138  for (i = 1; i <= nRC; ++i)
139  {
140  // accumulate subtr sum
141  Real subtrsum = 0.0;
142  for(int k = 1; k < i; ++k) subtrsum += a(k) * L(i,k);
143 
144  a(i) = (x(i) - subtrsum) / L(i,i);
145  }
146 
147  // test that l2 norm of a is < 1
148  Real squareNormA = a.SumSquare();
149  if (squareNormA >= 1.0)
150  Throw(ProgramException("downdate_Cholesky() fails", chol));
151 
152  Real alpha = sqrt(1.0 - squareNormA);
153 
154  // compute and apply Givens rotations to the vector a
155  ColumnVector cGivens(nRC); cGivens = 0.0;
156  ColumnVector sGivens(nRC); sGivens = 0.0;
157  for(i = nRC; i >= 1; i--)
158  alpha = pythag(alpha, a(i), cGivens(i), sGivens(i));
159 
160  // apply Givens rotations to the jth column of chol
161  ColumnVector xtilde(nRC); xtilde = 0.0;
162  for(j = nRC; j >= 1; j--)
163  {
164  // only the first j rotations have an affect on chol,0
165  for(int k = j; k >= 1; k--)
166  GivensRotation(cGivens(k), -sGivens(k), chol(k,j), xtilde(j));
167  }
168 }
169 
170 
171 
172 // produces the Cholesky decomposition of EAE where A = chol.t() * chol
173 // and E produces a RIGHT circular shift of the rows and columns from
174 // 1,...,k-1,k,k+1,...l,l+1,...,p to
175 // 1,...,k-1,l,k,k+1,...l-1,l+1,...p
177 {
178  int nRC = chol.Nrows();
179  int i, j;
180 
181  // I. compute shift of column l to the kth position
182  Matrix cholCopy = chol;
183  // a. grab column l
184  ColumnVector columnL = cholCopy.Column(l);
185  // b. shift columns k,...l-1 to the RIGHT
186  for(j = l-1; j >= k; --j)
187  cholCopy.Column(j+1) = cholCopy.Column(j);
188  // c. copy the top k-1 elements of columnL into the kth column of cholCopy
189  cholCopy.Column(k) = 0.0;
190  for(i = 1; i < k; ++i) cholCopy(i,k) = columnL(i);
191 
192  // II. determine the l-k Given's rotations
193  int nGivens = l-k;
194  ColumnVector cGivens(nGivens); cGivens = 0.0;
195  ColumnVector sGivens(nGivens); sGivens = 0.0;
196  for(i = l; i > k; i--)
197  {
198  int givensIndex = l-i+1;
199  columnL(i-1) = pythag(columnL(i-1), columnL(i),
200  cGivens(givensIndex), sGivens(givensIndex));
201  columnL(i) = 0.0;
202  }
203  // the kth entry of columnL is the new diagonal element in column k of cholCopy
204  cholCopy(k,k) = columnL(k);
205 
206  // III. apply these Given's rotations to subsequent columns
207  // for columns k+1,...,l-1 we only need to apply the last nGivens-(j-k) rotations
208  for(j = k+1; j <= nRC; ++j)
209  {
210  ColumnVector columnJ = cholCopy.Column(j);
211  int imin = nGivens - (j-k) + 1; if (imin < 1) imin = 1;
212  for(int gIndex = imin; gIndex <= nGivens; ++gIndex)
213  {
214  // apply gIndex Given's rotation
215  int topRowIndex = k + nGivens - gIndex;
216  GivensRotationR(cGivens(gIndex), sGivens(gIndex),
217  columnJ(topRowIndex), columnJ(topRowIndex+1));
218  }
219  cholCopy.Column(j) = columnJ;
220  }
221 
222  chol << cholCopy;
223 }
224 
225 
226 
227 // produces the Cholesky decomposition of EAE where A = chol.t() * chol
228 // and E produces a LEFT circular shift of the rows and columns from
229 // 1,...,k-1,k,k+1,...l,l+1,...,p to
230 // 1,...,k-1,k+1,...l,k,l+1,...,p to
232 {
233  int nRC = chol.Nrows();
234  int i, j;
235 
236  // I. compute shift of column k to the lth position
237  Matrix cholCopy = chol;
238  // a. grab column k
239  ColumnVector columnK = cholCopy.Column(k);
240  // b. shift columns k+1,...l to the LEFT
241  for(j = k+1; j <= l; ++j)
242  cholCopy.Column(j-1) = cholCopy.Column(j);
243  // c. copy the elements of columnK into the lth column of cholCopy
244  cholCopy.Column(l) = 0.0;
245  for(i = 1; i <= k; ++i)
246  cholCopy(i,l) = columnK(i);
247 
248  // II. apply and compute Given's rotations
249  int nGivens = l-k;
250  ColumnVector cGivens(nGivens); cGivens = 0.0;
251  ColumnVector sGivens(nGivens); sGivens = 0.0;
252  for(j = k; j <= nRC; ++j)
253  {
254  ColumnVector columnJ = cholCopy.Column(j);
255 
256  // apply the previous Givens rotations to columnJ
257  int imax = j - k; if (imax > nGivens) imax = nGivens;
258  for(int i = 1; i <= imax; ++i)
259  {
260  int gIndex = i;
261  int topRowIndex = k + i - 1;
262  GivensRotationR(cGivens(gIndex), sGivens(gIndex),
263  columnJ(topRowIndex), columnJ(topRowIndex+1));
264  }
265 
266  // compute a new Given's rotation when j < l
267  if(j < l)
268  {
269  int gIndex = j-k+1;
270  columnJ(j) = pythag(columnJ(j), columnJ(j+1), cGivens(gIndex),
271  sGivens(gIndex));
272  columnJ(j+1) = 0.0;
273  }
274 
275  cholCopy.Column(j) = columnJ;
276  }
277 
278  chol << cholCopy;
279 
280 }
281 
282 
283 
284 
285 #ifdef use_namespace
286 }
287 #endif
288 
void downdate_Cholesky(UpperTriangularMatrix &chol, RowVector x)
Definition: cholesky.cpp:129
void GivensRotation(Real cGivens, Real sGivens, Real &x, Real &y)
Definition: newmatrm.h:126
void GivensRotationR(Real cGivens, Real sGivens, Real &x, Real &y)
Definition: newmatrm.h:134
Miscellaneous exception (details in character string).
Definition: newmat.h:1947
ReturnMatrix for_return() const
Definition: newmat4.cpp:249
double Real
Definition: include.h:307
Real pythag(Real f, Real g, Real &c, Real &s)
Definition: newmatrm.cpp:189
Real square(Real x)
Definition: bandmat.cpp:672
int Nrows() const
Definition: newmat.h:494
Upper triangular matrix.
Definition: newmat.h:799
void update_Cholesky(UpperTriangularMatrix &chol, RowVector x)
Definition: cholesky.cpp:102
GetSubMatrix Column(int f) const
Definition: newmat.h:2152
FloatVector * d
#define REPORT
Definition: cholesky.cpp:27
TransposedMatrix t() const
Definition: newmat6.cpp:320
Real * Store() const
Definition: newmat.h:497
#define Throw(E)
Definition: myexcept.h:191
The usual rectangular matrix.
Definition: newmat.h:625
void left_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
Definition: cholesky.cpp:231
Real SumSquare() const
Definition: newmat.h:346
ReturnMatrix Cholesky(const SymmetricMatrix &S)
Definition: cholesky.cpp:36
FloatVector FloatVector * a
Real trace(const BaseMatrix &B)
Definition: newmat.h:2099
Lower triangular matrix.
Definition: newmat.h:848
void right_circular_update_Cholesky(UpperTriangularMatrix &chol, int k, int l)
Definition: cholesky.cpp:176
Symmetric band matrix.
Definition: newmat.h:1245
Lower triangular band matrix.
Definition: newmat.h:1206
FloatVector FloatVector FloatVector * alpha
Real sum(const BaseMatrix &B)
Definition: newmat.h:2105
Not positive definite exception.
Definition: newmat.h:1914
Row vector.
Definition: newmat.h:953
Column vector.
Definition: newmat.h:1008
void release()
Definition: newmat.h:517
Symmetric matrix.
Definition: newmat.h:753


kni
Author(s): Martin Günther
autogenerated on Fri Jan 3 2020 04:01:16