[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN Cholesky.h doc/Choleskydoc.h
|
From: |
Ethan Eade |
|
Subject: |
[Toon-members] TooN Cholesky.h doc/Choleskydoc.h |
|
Date: |
Mon, 15 Jan 2007 18:00:54 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Ethan Eade <ethaneade> 07/01/15 18:00:54
Modified files:
. : Cholesky.h
Added files:
doc : Choleskydoc.h
Log message:
Cholesky now computes LDL^T instead of LL^T. It is much faster, and
preserves the old interface. Any code that used it before should still
work, however the get_L() method is deprecated in favor of others.
Also added new methods for transforming the inverse.
Algorithms for fixed-size matrices are now implemented in compile-time
generated code instead of externally generated code. The quality of
generated code is identical on gcc-4.1 (the only one I checked).
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.12&r2=1.13
http://cvs.savannah.gnu.org/viewcvs/TooN/doc/Choleskydoc.h?cvsroot=toon&rev=1.1
Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.12
retrieving revision 1.13
diff -u -b -r1.12 -r1.13
--- Cholesky.h 23 Jul 2006 10:17:11 -0000 1.12
+++ Cholesky.h 15 Jan 2007 18:00:53 -0000 1.13
@@ -17,233 +17,381 @@
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
USA
-#ifndef __CHOLESKY_H
-#define __CHOLESKY_H
+#ifndef CHOLESKY_H
+#define CHOLESKY_H
#include <iostream>
#include <TooN/lapack.h>
#include <TooN/TooN.h>
-#include <TooN/util.h>
+#include <TooN/helpers.h>
#include <limits>
#ifndef TOON_NO_NAMESPACE
namespace TooN {
#endif
namespace util {
- template <int Size, class A1, class A2, class A3, class A4> inline
- void cholesky_backsub(const FixedMatrix<Size,Size,A1>& L, const
FixedVector<Size,A2>& invdiag, const FixedVector<Size,A3>& v,
FixedVector<Size,A4>& x)
- {
- Vector<Size> t;
- // forward substitution
- for (int i=0; i<Size; i++) {
- double b = v[i];
- for (int j=0; j<i;j++)
- b -= L[i][j]*t[j];
- t[i] = b*invdiag[i];
- }
- // back substitution
- for (int i=Size-1; i >=0; i--) {
- double b = t[i];
- for (int j=i+1; j<Size;j++)
- b -= L[j][i]*x[j];
- x[i] = b*invdiag[i];
+
+ template <int B, int E> struct Dot3 {
+ template <class V1, class V2, class V3> static inline double
eval(const V1& a, const V2& b, const V3& c) { return a[B]*b[B]*c[B] +
Dot3<B+1,E>::eval(a,b,c); }
+ };
+ template <int B> struct Dot3<B,B> {
+ template <class V1, class V2, class V3> static inline double
eval(const V1& a, const V2& b, const V3& c) { return a[B]*b[B]*c[B]; }
+ };
+
+ //
+ // Forward Sub
+ //
+ template <int N, int I=0> struct Forwardsub_L {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v,
FixedVector<N,A3>& x) {
+ x[I] = v[I] - Dot<0,I-1>::eval(L[I], x);
+ Forwardsub_L<N,I+1>::eval(L, v, x);
}
+ };
+ template <int N> struct Forwardsub_L<N,0> {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v,
FixedVector<N,A3>& x) {
+ x[0] = v[0];
+ Forwardsub_L<N,1>::eval(L, v, x);
}
+ };
+ template <int N> struct Forwardsub_L<N,N> {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v,
FixedVector<N,A3>& x) {}
+ };
- inline double dynamic_dot(const double* a, const double* b, int size)
- {
- double sum = 0;
- const double* const end = a+size;
- for (;a != end; ++a, ++b)
- sum += *a * *b;
- return sum;
+ //
+ // Back Sub
+ //
+ template <int N, int I=N> struct Backsub_LT {
+ template <class A1, class A2, class A3, class A4> static inline
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v,
+
const FixedVector<N,A3>& invdiag, FixedVector<N,A4>& x) {
+ x[I] = v[I]*invdiag[I] - Dot<I+1,N-1>::eval(L.T()[I], x);
+ Backsub_LT<N,I-1>::eval(L, v, invdiag, x);
}
-
-
- template <class A1, class A2, class V, class A4>
- void cholesky_backsub(const DynamicMatrix<A1>& L, const
DynamicVector<A2>& invdiag, const V& v, DynamicVector<A4>& x)
- {
- const int Size = L.num_rows();
- Vector<> t(Size);
- // forward substitution
- for (int i=0; i<Size; i++) {
- const double* Li = &L(i,0);
- t[i] = invdiag[i] * (v[i] - dynamic_dot(Li, t.get_data_ptr(),
i));
- }
- // back substitution
- for (int i=Size-1; i >=0; i--) {
- double b = t[i];
- const double* Lji = &L(i+1,i);
- for (int j=i+1; j<Size;j++, Lji += Size)
- b -= *Lji * x[j];
- x[i] = b*invdiag[i];
+ };
+ template <int N> struct Backsub_LT<N,N> {
+ template <class A1, class A2, class A3, class A4> static inline
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v,
+
const FixedVector<N,A3>& invdiag, FixedVector<N,A4>& x) {
+ x[N-1] = v[N-1]*invdiag[N-1];
+ Backsub_LT<N,N-2>::eval(L, v, invdiag, x);
}
+ };
+ template <int N> struct Backsub_LT<N,0> {
+ template <class A1, class A2, class A3, class A4> static inline
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& v,
+
const FixedVector<N,A3>& invdiag, FixedVector<N,A4>& x) {
+ x[0] = v[0]*invdiag[0] - Dot<1,N-1>::eval(L.T()[0], x);
}
+ };
- template <int Size, class A1, class A2, class A3> inline int
cholesky_compute(const FixedMatrix<Size,Size,A1>& M, FixedMatrix<Size,Size,A2>&
L, FixedVector<Size,A3>& invdiag)
+ template <int N, class A1, class A2, class A3, class A4>
+ void cholesky_solve(const FixedMatrix<N,N,A1>& L, const
FixedVector<N,A2>& invdiag, const FixedVector<N,A3>& v, FixedVector<N,A4>& x)
{
- int rank = Size;
- const double eps = std::numeric_limits<double>::min();
- for(int i=0;i<Size;i++) {
- double a=M[i][i];
- for(int k=0;k<i;k++)
- a-=L[i][k]*L[i][k];
- if (a < eps) {
- --rank;
- L[i][i] = invdiag[i] = 0;
- } else {
- L[i][i]=sqrt(a);
- invdiag[i] = 1.0/L[i][i];
- }
- for(int j=i+1;j<Size;j++) {
- L[i][j] = 0;
- a=M[i][j];
- for(int k=0;k<i;k++)
- a-=L[i][k]*L[j][k];
- L[j][i]=a*invdiag[i];
- }
+ Vector<N> t;
+ Forwardsub_L<N>::eval(L, v, t);
+ Backsub_LT<N>::eval(L, t, invdiag, x);
}
- return rank;
+
+ //
+ // Compute cholesky
+ //
+ template <int N, int I, int J=I+1> struct CholeskyInner {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, const
FixedVector<N,A3>& invdiag) {
+ const double a = M[I][J] - Dot<0,I-1>::eval(L[I], L.T()[J]);
+ L[I][J] = a;
+ L[J][I] = a * invdiag[I];
+ CholeskyInner<N, I, J+1>::eval(M, L, invdiag);
}
+ };
- template <class A1, class A2, class A3> inline int
cholesky_compute(const DynamicMatrix<A1>& M, DynamicMatrix<A2>& L,
DynamicVector<A3>& invdiag)
- {
- const int Size = M.num_rows();
- int rank = Size;
- const double eps = std::numeric_limits<double>::min();
- for(int i=0;i<Size;i++) {
- double* const Li = &L(i,0);
- double a = M(i,i);
- for (int j=0; j<i; ++j) {
- double sum = invdiag[j] * (M(j,i) - dynamic_dot(Li,
&L(j,0), j));
- Li[j] = sum;
- a -= sum*sum;
- }
- if (eps < a) {
- const double s = sqrt(a);
- Li[i]= s;
- invdiag[i] = 1.0/s;
+ template <int N, int I> struct CholeskyInner<N,I,N> {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, const
FixedVector<N,A3>& invdiag) {}
+ };
+ template <int N, int I=0> struct CholeskyOuter {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, FixedVector<N,A3>&
invdiag, int& rank) {
+ const double a = M[I][I] - Dot<0,I-1>::eval(L[I], L.T()[I]);
+ if (0 < a) {
+ L[I][I] = a;
+ invdiag[I] = 1.0/a;
+ ++rank;
} else {
- --rank;
- Li[i] = invdiag[i] = 0;
+ L[I][I] = invdiag[I] = 0;
+ return;
}
+ CholeskyInner<N,I>::eval(M,L,invdiag);
+ CholeskyOuter<N,I+1>::eval(M,L,invdiag,rank);
}
+ };
+ template <int N> struct CholeskyOuter<N,N> {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L, FixedVector<N,A3>&
invdiag, int& rank) {}
+ };
+
+
+ template <int N, class A1, class A2, class A3> int
cholesky_compute(const FixedMatrix<N,N,A1>& M, FixedMatrix<N,N,A2>& L,
FixedVector<N,A3>& invdiag)
+ {
+ int rank=0;
+ CholeskyOuter<N>::eval(M, L, invdiag, rank);
return rank;
}
+ //
+ // Compute inverse
+ //
+ template <int N, int Col, int I=Col+1> struct InverseForward {
+ template <class A1, class A2> static inline void eval(const
FixedMatrix<N,N,A1>& L, FixedVector<N,A2>& t) {
+ t[I] = -(L[I][Col] + Dot<Col+1,I-1>::eval(L[I], t));
+ InverseForward<N,Col,I+1>::eval(L,t);
+ }
+ };
+ template <int N, int Col> struct InverseForward<N,Col,N> {
+ template <class A1, class A2> static inline void eval(const
FixedMatrix<N,N,A1>& L, FixedVector<N,A2>& t) {}
+ };
- template <int S, class A1, class A2, class A3> inline void
cholesky_inverse(const FixedMatrix<S,S,A1>& L, const FixedVector<S,A2>&
invdiag, FixedMatrix<S,S,A3>& I)
- {
- for (int col = 0; col<S; col++) {
- Vector<S> t,x;
- t[col] = invdiag[col];
- for (int i=col+1; i<S; i++) {
- double psum = 0;
- for (int j=col; j<i; j++)
- psum -= L[i][j]*t[j];
- t[i] = invdiag[i] * psum;
- }
- for (int i=S-1; i>col; i--) {
- double psum = t[i];
- for (int j=i+1; j<S; j++)
- psum -= L[j][i]*x[j];
- I[i][col] = I[col][i] = x[i] = invdiag[i] * psum;
- }
- double psum = t[col];
- for (int j=col+1; j<S; j++)
- psum -= L[j][col]*x[j];
- I[col][col] = invdiag[col]*psum;
+ template <int N, int Col, int I=N-1> struct InverseBack {
+ template <class A1, class A2, class A3, class A4> static inline
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& t,
+
const FixedVector<N,A3>& invdiag, FixedMatrix<N,N,A4>& Inv) {
+ Inv[I][Col] = Inv[Col][I] = invdiag[I]*t[I] -
Dot<I+1,N-1>::eval(L.T()[I], Inv[Col]);
+ InverseBack<N,Col,I-1>::eval(L, t, invdiag, Inv);
+ }
+ };
+ template <int N, int Col> struct InverseBack<N,Col,Col> {
+ template <class A1, class A2, class A3, class A4> static inline
void eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& t,
+
const FixedVector<N,A3>& invdiag, FixedMatrix<N,N,A4>& Inv) {
+ Inv[Col][Col] = invdiag[Col] - Dot<Col+1,N-1>::eval(L.T()[Col],
Inv[Col]);
}
+ };
+
+ template <int N, int Col=0> struct InverseOuter {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& invdiag,
FixedMatrix<N,N,A3>& Inv) {
+ Vector<N> t;
+ InverseForward<N,Col>::eval(L, t);
+ InverseBack<N,Col>::eval(L, t, invdiag, Inv);
+ InverseOuter<N,Col+1>::eval(L,invdiag,Inv);
}
+ };
+ template <int N> struct InverseOuter<N,N> {
+ template <class A1, class A2, class A3> static inline void
eval(const FixedMatrix<N,N,A1>& L, const FixedVector<N,A2>& invdiag,
FixedMatrix<N,N,A3>& Inv) {}
+ };
- template <class A1, class A2, class A3> inline void
cholesky_inverse(const DynamicMatrix<A1>& L, const DynamicVector<A2>& invdiag,
DynamicMatrix<A3>& I)
+ template <int S, class A1, class A2, class A3> void
cholesky_inverse(const FixedMatrix<S,S,A1>& L, const FixedVector<S,A2>&
invdiag, FixedMatrix<S,S,A3>& I)
{
- int S = L.num_rows();
- for (int col = 0; col<S; col++) {
- Vector<> t(S),x(S);
- t[col] = invdiag[col];
- for (int i=col+1; i<S; i++) {
- const double* Li = &L(i,0);
- double psum = 0;
- for (int j=col; j<i; j++)
- psum -= Li[j]*t[j];
- t[i] = invdiag[i] * psum;
- }
- for (int i=S-1; i>col; i--) {
- const double* Lji = &L(i+1,i);
- double psum = t[i];
- for (int j=i+1; j<S; j++, Lji += S)
- psum -= *Lji * x[j];
- I(i,col) = I(col,i) = x[i] = invdiag[i] * psum;
+ InverseOuter<S>::eval(L, invdiag, I);
}
- double psum = t[col];
- const double* Ljc = &L(col+1,col);
- for (int j=col+1; j<S; j++, Ljc += S)
- psum -= *Ljc * x[j];
- I(col,col) = invdiag[col]*psum;
+
+ //
+ // Cholesky update
+ //
+ template <int N, int I, int J=I-1> struct UpdateInner {
+ template <class LT, class PT, class FT, class ST> static inline
void eval(LT& L, const PT& p, const FT& F, const ST sigma) {
+ const double old = L[I][J];
+ L[I][J] = old + (sigma * p[J]) * F[J];
+ UpdateInner<N,I,J-1>::eval(L, p, F, sigma + old * p[J]);
}
+ };
+
+ template <int N, int I> struct UpdateInner<N,I,-1> {
+ template <class LT, class PT, class FT, class ST> static inline
void eval(LT& L, const PT& p, const FT& F, const ST sigma) {}
+ };
+
+ template <int N, int I=1> struct UpdateOuter {
+ template <class LT, class PT, class FT> static inline void eval(LT&
L, const PT& p, const FT& F) {
+ UpdateInner<N,I>::eval(L, p, F, p[I]);
+ UpdateOuter<N,I+1>::eval(L, p ,F);
}
+ };
+ template <int N> struct UpdateOuter<N,N> {
+ template <class LT, class PT, class FT> static inline void eval(LT&
L, const PT& p, const FT& F) {}
+ };
+ template <class P, int N, class A1, class A2> void
cholesky_update(TooN::FixedMatrix<N,N,A1>& L, TooN::FixedVector<N,A2>& invdiag,
const P& p)
+ {
+ double F[N-1];
+ double alpha = 1;
+ for (int i=0; i<N-1; ++i) {
+ const double p2 = p[i]*p[i];
+ const double di = L[i][i];
+ L[i][i] += alpha * p2;
+ invdiag[i] = 1.0/L[i][i];
+ const double a_over_d = alpha*invdiag[i];
+ F[i] = a_over_d;
+ alpha = di * a_over_d;
+ }
+ L[N-1][N-1] += alpha * (p[N-1]*p[N-1]);
+ invdiag[N-1] = 1.0/L[N-1][N-1];
+
+ UpdateOuter<N>::eval(L, p, F);
+ }
}
- template <int Size=-1>
+ template <int N=-1>
class Cholesky {
public:
template<class Accessor>
- Cholesky(const FixedMatrix<Size,Size,Accessor>& m){
+ Cholesky(const FixedMatrix<N,N,Accessor>& m){
compute(m);
}
template<class Accessor>
- void compute(const FixedMatrix<Size,Size,Accessor>& m){
+ void compute(const FixedMatrix<N,N,Accessor>& m){
rank = util::cholesky_compute(m,L,invdiag);
}
int get_rank() const { return rank; }
double get_determinant() const {
double det = L[0][0];
- for (int i=1; i<Size; i++)
+ for (int i=1; i<N; i++)
det *= L[i][i];
- return det*det;
+ return det;
+ }
+
+ const Matrix<N>& get_L_D() const { return L; }
+
+ template <class A1, class A2> static void sqrt(const
FixedMatrix<N,N,A1>& A, FixedMatrix<N,N,A2>& L) {
+ for (int i=0; i<N; ++i) {
+ double a = A[i][i];
+ for (int k=0; k<i; ++k)
+ a -= L[i][k]*L[i][k];
+ if (0<a) {
+ L[i][i] = ::sqrt(a);
+ } else {
+ Zero(L.slice(i,i,N-i,N-i));
+ return;
+ }
+ const double id = 1.0/L[i][i];
+ for (int j=i+1; j<N; ++j) {
+ L[i][j] = 0;
+ double a = A[i][j];
+ for (int k=0; k<i; ++k)
+ a -= L[i][k]*L[j][k];
+ L[j][i] = a * id;
+ }
+ }
}
- const Matrix<Size>& get_L() const {
+ template <class A1> static Matrix<N> sqrt(const FixedMatrix<N,N,A1>& A)
{
+ Matrix<N> L;
+ Cholesky<N>::sqrt(A,L);
return L;
}
+ template <class A> void get_sqrt(FixedMatrix<N,N,A>& M) const {
+ for (int i=0; i<N; ++i) {
+ const double root_d = ::sqrt(L[i][i]);
+ M[i][i] = root_d;
+ for (int j=i+1; j<N; ++j) {
+ M[j][i] = L[j][i]*root_d;
+ M[i][j] = 0;
+ }
+ }
+ }
+
+ Matrix<N> get_sqrt() const {
+ Matrix<N> S;
+ get_sqrt(S);
+ return S;
+ }
+
+ Matrix<N> get_L() const __attribute__ ((deprecated)) { return
get_sqrt(); }
+
+ template <class A> void get_inv_sqrt(FixedMatrix<N,N,A>& M) const {
+ Vector<N> root;
+ for (int i=0; i<N; ++i)
+ root[i] = ::sqrt(invdiag[i]);
+ for (int j=0; j<N; ++j) {
+ for (int i=j+1; i<N; ++i) {
+ double sum = L[i][j];
+ for (int k=j+1; k<i; ++k)
+ sum += L[i][k]*M[j][k];
+ M[j][i] = -sum;
+ M[i][j] = 0;
+ }
+ M[j][j] = root[j];
+ for (int i=j+1; i<N; ++i)
+ M[j][i] *= root[i];
+ }
+ }
+
+ template <class A> double mahalanobis(const FixedVector<N,A>& v) const {
+ Vector<N> L_inv_v;
+ util::Forwardsub_L<N>::eval(L, v, L_inv_v);
+ return util::Dot3<0,N-1>::eval(L_inv_v, L_inv_v, invdiag);
+ }
+
+ template <class F, int M, class A1, class A2> void
transform_inverse(const FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const {
+ Matrix<M,N> L_inv_JT;
+ for (int i=0; i<M; ++i)
+ util::Forwardsub_L<N>::eval(L, J[i], L_inv_JT[i]);
+ for (int i=0; i<M; ++i) {
+ F::eval(T[i][i],util::Dot3<0,N-1>::eval(L_inv_JT[i],
L_inv_JT[i], invdiag));
+ for (int j=i+1; j<M; ++j) {
+ const double x = util::Dot3<0,N-1>::eval(L_inv_JT[i],
L_inv_JT[j], invdiag);
+ F::eval(T[i][j],x);
+ F::eval(T[j][i],x);
+ }
+ }
+ }
+ template <int M, class A1, class A2> void transform_inverse(const
FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const {
+ transform_inverse<util::Assign>(J,T);
+ }
+ template <int M, class A> Matrix<M> transform_inverse(const
FixedMatrix<M,N,A>& J) const {
+ Matrix<M> T;
+ transform_inverse(J,T);
+ return T;
+ }
+
+ template <class A1, class A2> inline
+ void inverse_times(const FixedVector<N,A1>& v, FixedVector<N,A2>& x)
const
+ {
+ util::cholesky_solve(L, invdiag, v, x);
+ }
+
template <class Accessor> inline
- Vector<Size> backsub(const FixedVector<Size,Accessor>& v) const
+ Vector<N> inverse_times(const FixedVector<N,Accessor>& v) const
{
- Vector<Size> x;
- util::cholesky_backsub(L, invdiag, v, x);
+ Vector<N> x;
+ inverse_times(v, x);
return x;
}
- template <class Accessor, int Cols> inline Matrix<Size,Cols>
inverse_times(const FixedMatrix<Size,Cols,Accessor>& m)
+ template <class Accessor> inline
+ Vector<N> backsub(const FixedVector<N,Accessor>& v) const { return
inverse_times(v); }
+
+ template <class A, int M> inline Matrix<N,M> inverse_times(const
FixedMatrix<N,M,A>& m)
{
- Matrix<Cols, Size> result;
- for (int i=0; i<Cols; i++)
- result[i] = backsub(m.T()[i]);
- return result.T();
+ Matrix<N,M> result;
+ for (int i=0; i<M; i++)
+ inverse_times(m.T()[i], result.T()[i]);
+ return result;
}
- template <class A> void get_inverse(FixedMatrix<Size,Size,A>& M) const {
+ template <class A> void get_inverse(FixedMatrix<N,N,A>& M) const {
util::cholesky_inverse(L, invdiag, M);
}
- Matrix<Size,Size> get_inverse() const {
- Matrix<Size,Size> M;
+ Matrix<N> get_inverse() const {
+ Matrix<N> M;
get_inverse(M);
return M;
}
+ template <int M, class A> void update(const FixedMatrix<N,M,A>& V) {
+ for (int i=0; i<M; ++i) {
+ Vector<N> p;
+ util::Forwardsub_L<N>::eval(L, V.T()[i], p);
+ util::cholesky_update(L, invdiag, p);
+ }
+ }
+
+ template <class A> void update(const FixedVector<N,A>& v) {
+ Vector<N> p;
+ util::Forwardsub_L<N>::eval(L, v, p);
+ util::cholesky_update(L, invdiag, p);
+ }
+
private:
- Matrix<Size,Size,RowMajor> L;
- Vector<Size> invdiag;
+ Matrix<N> L;
+ Vector<N> invdiag;
int rank;
};
@@ -252,14 +400,20 @@
public:
template<class Accessor>
- Cholesky(const DynamicMatrix<Accessor>& m) : L(m.num_rows(),
m.num_cols()), invdiag(m.num_rows()) {
- assert(m.num_rows() == m.num_cols());
+ Cholesky(const DynamicMatrix<Accessor>& m) {
compute(m);
}
template<class Accessor>
void compute(const DynamicMatrix<Accessor>& m){
- rank = util::cholesky_compute(m,L,invdiag);
+ assert(m.num_rows() == m.num_cols());
+ L.assign(m);
+ int N = L.num_rows();
+ int info;
+ dpotrf_("L", &N, L.get_data_ptr(), &N, &info);
+ assert(info >= 0);
+ if (info > 0)
+ rank = info-1;
}
int get_rank() const { return rank; }
@@ -267,8 +421,12 @@
Vector<> backsub(const V& v) const
{
assert(v.size() == L.num_rows());
- Vector<> x(L.num_rows());
- util::cholesky_backsub(L, invdiag, v, x);
+ Vector<> x = v;
+ int N=L.num_rows();
+ int NRHS=1;
+ int info;
+ dpotrs_("L", &N, &NRHS, L.get_data_ptr(), &N, x.get_data_ptr(), &N,
&info);
+ assert(info==0);
return x;
}
@@ -284,8 +442,12 @@
}
template <class Mat> void get_inverse(Mat& M) const {
- assert(M.num_rows() == M.num_cols() && M.num_rows() ==
L.num_rows());
- util::cholesky_inverse(L, invdiag, M);
+ M = L;
+ int N = L.num_rows();
+ int info;
+ dpotri_("L", &N, M.get_data_ptr(), &N, &info);
+ assert(info == 0);
+ TooN::Symmetrize(M);
}
Matrix<> get_inverse() const {
@@ -296,7 +458,6 @@
private:
Matrix<> L;
- Vector<> invdiag;
int rank;
};
Index: doc/Choleskydoc.h
===================================================================
RCS file: doc/Choleskydoc.h
diff -N doc/Choleskydoc.h
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ doc/Choleskydoc.h 15 Jan 2007 18:00:54 -0000 1.1
@@ -0,0 +1,168 @@
+/*
+ Copyright (c) 2005 Paul Smith
+
+ Permission is granted to copy, distribute and/or modify this document under
+ the terms of the GNU Free Documentation License, Version 1.2 or any later
+ version published by the Free Software Foundation; with no Invariant
+ Sections, no Front-Cover Texts, and no Back-Cover Texts.
+
+ You should have received a copy of the GNU Free Documentation License
+ License along with this library; if not, write to the Free Software
+ Foundation, Inc.
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+
+*/
+// A proxy version of the LU class,
+// cleaned up to present a comprehensible
+// version of the LU interface
+
+#ifdef DOXYGEN_INCLUDE_ONLY_FOR_DOCS
+
+
+#include <iostream>
+
+#include <TooN/lapack.h>
+
+#include <TooN/TooN.h>
+#include <TooN/helpers.h>
+#include <limits>
+
+namespace TooN {
+/// @class Cholesky Choleskydoc.h TooN/Cholesky.h
+/// Decomposes a positive-semidefinite symmetric matrix A (such as a
covariance) into L*D*L^T, where L is lower-triangular and D is diagonal.
+/// Also can compute A = S*S^T, with S lower triangular. The LDL^T form is
faster to compute than the class Cholesky decomposition.
+/// The decomposition can be used to compute A^-1*x, A^-1*M, M*A^-1*M^T, and
A^-1 itself, though the latter rarely needs to be explicitly represented.
+/// Also efficiently computes det(A) and rank(A).
+/// It can be used as follows:
+/// @code
+/// // Declare some matrices.
+/// Matrix<3> A = ...; // we'll pretend it is pos-def
+/// Matrix<2,3> M;
+/// Matrix<2> B;
+/// Vector<3> y = (make_Vector, 2,3,4);
+/// // create the Cholesky decomposition of A
+/// Cholesky<3> chol(A);
+/// // compute x = A^-1 * y
+/// Vector<3> x = cholA.inverse_times(y);
+/// Identical to above
+/// x = cholA.backsub(y);
+/// // compute B = M*A^-1*M^T
+/// B = cholA.transform_inverse(M);
+/// //compute A^-1
+/// Matrix<3> Ainv = cholA.get_inverse();
+/// Matrix<3> C = ... // again, C is pos-def
+/// //compute the 'square-root' of C
+/// Matrix<3> L = Cholesky<3>::sqrt(C);
+/// @endcode
+/// @ingroup gDecomps
+template <int N>
+class Cholesky {
+public:
+
+ /// Construct the Cholesky-ish decomposition of a matrix. This initialises
the class, and
+ /// performs the decomposition immediately.
+ /// Run time is O(N^3)
+ template<class Accessor> Cholesky(const FixedMatrix<N,N,Accessor>& m);
+
+ /// Compute the LDL^T decomposition of another matrix.
+ /// Run time is O(N^3)
+ template<class Accessor> void compute(const FixedMatrix<N,N,Accessor>& m);
+
+ /// Get the rank of the matrix
+ int get_rank() const;
+
+ /// Get the determinant
+ double get_determinant() const;
+
+ /// Get the internal representation of L and D.
+ /// L has 1's on the diagonal, so here D is stored in its diagonal.
+ /// L is lower triangular; the upper-triangular contents are not
initialized.
+ /// Not yet present in Cholesky<-1>
+ const Matrix<N>& get_L_D() const { return L; }
+
+ /// Compute the lower-triangular 'square-root' L of A, s.t. A = L*L^T
+ /// Run time is O(N^3)
+ /// Not yet present in Cholesky<-1>
+ template <class A1, class A2> static void sqrt(const FixedMatrix<N,N,A1>&
A, FixedMatrix<N,N,A2>& L);
+
+ /// Compute the lower-triangular 'square-root' L, of A, s.t. A = L*L^T
+ /// Run time is O(N^3)
+ /// Not yet present in Cholesky<-1>
+ template <class A1> static Matrix<N> sqrt(const FixedMatrix<N,N,A1>& A);
+
+ /// Compute the lower-triangular 'square-root' L of the decomposed matrix,
s.t. A = L*L^T
+ /// L is stored in M
+ /// Not yet present in Cholesky<-1>
+ template <class A> void get_sqrt(FixedMatrix<N,N,A>& M) const;
+
+ /// Compute the lower-triangular 'square-root' L of the decomposed matrix,
s.t. A = L*L^T
+ /// Not yet present in Cholesky<-1>
+ Matrix<N> get_sqrt() const;
+
+ /// Compute the lower-triangular 'square-root' L of the decomposed matrix,
s.t. A = L*L^T
+ /// @deprecated Use get_sqrt, or the static sqrt function
+ Matrix<N> get_L() const __attribute__ ((deprecated));
+
+ /// Compute the lower-triangular 'square-root' L of the inverse of the
decomposed matrix, s.t. A^-1 = L*L^T
+ /// Not yet present in Cholesky<-1>
+ template <class A> void get_inv_sqrt(FixedMatrix<N,N,A>& M) const;
+
+ /// Compute the Mahalanobis squared magnitude of x, d = x*A^-1*x
+ /// Run time is O(N^2)
+ /// Not yet present in Cholesky<-1>
+ template <class A> double mahalanobis(const FixedVector<N,A>& v) const;
+
+ /// Compute J*A^-1*J^T, and store in T using F::eval (for instance, F
might be util::PlusAssign to add the result to T)
+ /// Run time is O(MN^2 + NM^2)
+ /// Not yet present in Cholesky<-1>
+ template <class F, int M, class A1, class A2> void transform_inverse(const
FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const;
+
+ /// Compute J*A^-1*J^T, and store in T
+ /// Run time is O(MN^2 + NM^2)
+ /// Not yet present in Cholesky<-1>
+ template <int M, class A1, class A2> void transform_inverse(const
FixedMatrix<M,N,A1>& J, FixedMatrix<M,M,A2>& T) const;
+
+ /// Compute J*A^-1*J^T
+ /// Run time is O(MN^2 + NM^2)
+ /// Not yet present in Cholesky<-1>
+ template <int M, class A> Matrix<M> transform_inverse(const
FixedMatrix<M,N,A>& J) const;
+
+ /// Compute x = A^-1*v
+ /// Run time is O(N^2)
+ /// Not yet present in Cholesky<-1>
+ template <class A1, class A2> inline void inverse_times(const
FixedVector<N,A1>& v, FixedVector<N,A2>& x) const;
+
+ /// Compute A^-1*v
+ /// Run time is O(N^2)
+ /// Not yet present in Cholesky<-1>
+ template <class Accessor> inline Vector<N> inverse_times(const
FixedVector<N,Accessor>& v) const;
+
+ /// Compute x = A^-1*v
+ /// Run time is O(N^2)
+ template <class Accessor> inline Vector<N> backsub(const
FixedVector<N,Accessor>& v) const;
+
+ /// Compute A^-1*B
+ /// Run time is O(MN^2)
+ /// Not yet present in Cholesky<-1>
+ template <class A, int M> inline Matrix<N,M> inverse_times(const
FixedMatrix<N,M,A>& B);
+
+ /// Compute A^-1 and store in M
+ /// Run time is O(N^3)
+ template <class A> void get_inverse(FixedMatrix<N,N,A>& M) const;
+
+ /// Compute A^-1
+ /// Run time is O(N^3)
+ Matrix<N> get_inverse() const;
+
+ /// Update the decomposition for A' = A + V*V^T
+ /// Run time is O(MN^2), but for small N, it's faster to compute A' and
recompute the decomposition
+ /// Not yet present in Cholesky<-1>
+ template <int M, class A> void update(const FixedMatrix<N,M,A>& V);
+
+ /// Update the decomposition for A' = A + v*v^T
+ /// Run time is O(N^2), but for small N, it's faster to compute A' and
recompute the decomposition
+ /// Not yet present in Cholesky<-1>
+ template <class A> void update(const FixedVector<N,A>& v);
+};
+
+#endif
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Cholesky.h doc/Choleskydoc.h,
Ethan Eade <=