[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Toon-members] TooN Cholesky.h Lapack_Cholesky.h Makefile.in r...
|
From: |
Edward Rosten |
|
Subject: |
[Toon-members] TooN Cholesky.h Lapack_Cholesky.h Makefile.in r... |
|
Date: |
Mon, 19 Oct 2009 14:08:53 +0000 |
CVSROOT: /cvsroot/toon
Module name: TooN
Changes by: Edward Rosten <edrosten> 09/10/19 14:08:53
Modified files:
. : Cholesky.h Lapack_Cholesky.h Makefile.in
Added files:
regressions : chol_lapack.cc chol_lapack.txt chol_toon.cc
chol_toon.txt
Log message:
Patch from Ville Kyrki fixing Cholesky:
- get_L() gets L where L L^T = M.
- Cholesky.h has get_unscaled_L() and get_D() where L D L^T = M.
- Regression tests added.
- Documentation fixed.
CVSWeb URLs:
http://cvs.savannah.gnu.org/viewcvs/TooN/Cholesky.h?cvsroot=toon&r1=1.37&r2=1.38
http://cvs.savannah.gnu.org/viewcvs/TooN/Lapack_Cholesky.h?cvsroot=toon&r1=1.2&r2=1.3
http://cvs.savannah.gnu.org/viewcvs/TooN/Makefile.in?cvsroot=toon&r1=1.16&r2=1.17
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/chol_lapack.cc?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/chol_lapack.txt?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/chol_toon.cc?cvsroot=toon&rev=1.1
http://cvs.savannah.gnu.org/viewcvs/TooN/regressions/chol_toon.txt?cvsroot=toon&rev=1.1
Patches:
Index: Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Cholesky.h,v
retrieving revision 1.37
retrieving revision 1.38
diff -u -b -r1.37 -r1.38
--- Cholesky.h 9 Jun 2009 13:33:22 -0000 1.37
+++ Cholesky.h 19 Oct 2009 14:08:52 -0000 1.38
@@ -38,7 +38,8 @@
/**
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.
+Also can compute the classic A = L*L^T, with L lower triangular. The LDL^T
form is faster to compute than the classical Cholesky decomposition.
+Use get_unscaled_L() and get_D() to access the individual matrices of L*D*L^T
decomposition. Use get_L() to access the lower triangular matrix of the classic
Cholesky decomposition L*L^T.
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:
@@ -122,7 +123,7 @@
/// Compute x = A^-1*v
/// Run time is O(N^2)
template<int Size2, class P2, class B2>
- Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>& v) {
+ Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>& v) const {
int size=my_cholesky.num_rows();
SizeMismatch<Size,Size2>::test(size, v.size());
@@ -157,7 +158,7 @@
/**overload
*/
template<int Size2, int C2, class P2, class B2>
- Matrix<Size, C2, Precision> backsub (const Matrix<Size2, C2, P2, B2>&
m) {
+ Matrix<Size, C2, Precision> backsub (const Matrix<Size2, C2, P2, B2>&
m) const {
int size=my_cholesky.num_rows();
SizeMismatch<Size,Size2>::test(size, m.num_rows());
@@ -206,6 +207,47 @@
return answer;
}
+ template <int Size2, typename P2, typename B2>
+ Precision mahalanobis(const Vector<Size2, P2, B2>& v) const {
+ return v * backsub(v);
+ }
+
+ Matrix<Size,Size,Precision> get_unscaled_L() const {
+ Matrix<Size,Size,Precision> m(my_cholesky.num_rows(),
+ my_cholesky.num_rows());
+ m=Identity;
+ for (int i=1;i<my_cholesky.num_rows();i++) {
+ for (int j=0;j<i;j++) {
+ m(i,j)=my_cholesky(i,j);
+ }
+ }
+ return m;
+ }
+
+ Matrix<Size,Size,Precision> get_D() const {
+ Matrix<Size,Size,Precision> m(my_cholesky.num_rows(),
+ my_cholesky.num_rows());
+ m=Zeros;
+ for (int i=0;i<my_cholesky.num_rows();i++) {
+ m(i,i)=my_cholesky(i,i);
+ }
+ return m;
+ }
+
+ Matrix<Size,Size,Precision> get_L() const {
+ Matrix<Size,Size,Precision> m(my_cholesky.num_rows(),
+ my_cholesky.num_rows());
+ m=Zeros;
+ for (int j=0;j<my_cholesky.num_cols();j++) {
+ Precision sqrtd=sqrt(my_cholesky(j,j));
+ m(j,j)=sqrtd;
+ for (int i=j+1;i<my_cholesky.num_rows();i++) {
+ m(i,j)=my_cholesky(i,j)*sqrtd;
+ }
+ }
+ return m;
+ }
+
private:
Matrix<Size,Size,Precision> my_cholesky;
};
Index: Lapack_Cholesky.h
===================================================================
RCS file: /cvsroot/toon/TooN/Lapack_Cholesky.h,v
retrieving revision 1.2
retrieving revision 1.3
diff -u -b -r1.2 -r1.3
--- Lapack_Cholesky.h 1 Oct 2009 12:36:17 -0000 1.2
+++ Lapack_Cholesky.h 19 Oct 2009 14:08:52 -0000 1.3
@@ -33,6 +33,10 @@
#include <TooN/TooN.h>
+#include <TooN/lapack.h>
+
+#include <assert.h>
+
namespace TooN {
@@ -63,16 +67,20 @@
giving symmetric M = L*D*L.T() where the diagonal of L contains ones
@param Size the size of the matrix
@param Precision the precision of the entries in the matrix and its
decomposition
+
address@hidden Always uses double precision LAPACK routines. On-line checking
of
+precision or template specialization could be used to implement single
+precision if necessary.
**/
-template <int Size, typename Precision>
+template <int Size, typename Precision=DefaultPrecision>
class Lapack_Cholesky {
public:
Lapack_Cholesky(){}
template<class P2, class B2>
- Lapack_Cholesky(const Matrix<Size, Size, P2, B2>& m) {
- : my_cholesky(m) {
+ Lapack_Cholesky(const Matrix<Size, Size, P2, B2>& m)
+ : my_cholesky(m), my_cholesky_lapack(m) {
SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
do_compute();
}
@@ -83,7 +91,7 @@
template<class P2, class B2> void compute(const Matrix<Size, Size, P2,
B2>& m){
SizeMismatch<Size,Size>::test(m.num_rows(), m.num_cols());
SizeMismatch<Size,Size>::test(m.num_rows(),
my_cholesky.num_rows());
- my_cholesky=m;
+ my_cholesky_lapack=m;
do_compute();
}
@@ -92,7 +100,18 @@
void do_compute(){
int N = my_cholesky.num_rows();
int info;
- dpotrf_("L", &N, L.get_data_ptr(), &N, &info);
+ dpotrf_("L", &N, (double*) my_cholesky_lapack.my_data, &N,
&info);
+ for (int i=0;i<N;i++) {
+ int j;
+ for (j=0;j<=i;j++) {
+ my_cholesky[i][j]=my_cholesky_lapack[j][i];
+ }
+ // LAPACK does not set upper triangle to zero,
+ // must be done here
+ for (;j<N;j++) {
+ my_cholesky[i][j]=0;
+ }
+ }
assert(info >= 0);
if (info > 0) {
my_rank = info-1;
@@ -102,64 +121,66 @@
int rank() const { return my_rank; }
template <int Size2, typename P2, typename B2>
- Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>&
v) {
+ Vector<Size, Precision> backsub (const Vector<Size2, P2, B2>&
v) const {
SizeMismatch<Size,Size2>::test(my_cholesky.num_cols(),
v.size());
Vector<Size> result(v);
- int N=L.num_rows();
+ int N=my_cholesky.num_rows();
int NRHS=1;
int info;
- dpotrs_("L", &N, &NRHS, my_cholesky.my_data, &N,
result.my_data, &N, &info);
+ dpotrs_("L", &N, &NRHS, my_cholesky_lapack.my_data, &N,
result.my_data, &N, &info);
assert(info==0);
return result;
}
template <int Size2, int Cols2, typename P2, typename B2>
- Matrix<Size, Cols2, Precision, ColMajor> backsub (const
Matrix<Size2, Cols2, P2, B2>& m) {
+ Matrix<Size, Cols2, Precision, ColMajor> backsub (const
Matrix<Size2, Cols2, P2, B2>& m) const {
SizeMismatch<Size,Size2>::test(my_cholesky.num_cols(),
m.num_rows());
Matrix<Size, Cols2, Precision, ColMajor> result(m);
int N=my_cholesky.num_rows();
int NRHS=m.num_cols();
int info;
- dpotrs_("L", &N, &NRHS, my_cholesky.my_data, &N,
result.my_data, &N, &info);
+ dpotrs_("L", &N, &NRHS, my_cholesky_lapack.my_data, &N,
result.my_data, &N, &info);
assert(info==0);
return result;
}
-
-
-
-
template <int Size2, typename P2, typename B2>
Precision mahalanobis(const Vector<Size2, P2, B2>& v) const {
return v * backsub(v);
}
- const Matrix<>& get_L() const {
+ Matrix<Size,Size,Precision> get_L() const {
return my_cholesky;
}
Precision determinant() const {
- Precision det = L[0][0];
+ Precision det = my_cholesky[0][0];
for (int i=1; i<my_cholesky.num_rows(); i++)
- det *= L[i][i];
+ det *= my_cholesky[i][i];
return det*det;
}
Matrix<> get_inverse() const {
- Matrix<Size> M(my_cholesky);
+ Matrix<Size> M(my_cholesky.num_rows(),my_cholesky.num_rows());
+ M=my_cholesky_lapack;
int N = my_cholesky.num_rows();
int info;
dpotri_("L", &N, M.my_data, &N, &info);
assert(info == 0);
- TooN::Symmetrize(M);
+ for (int i=1;i<N;i++) {
+ for (int j=0;j<i;j++) {
+ M[i][j]=M[j][i];
+ }
+ }
return M;
}
private:
+ Matrix<Size,Size,double> my_cholesky_lapack;
Matrix<Size,Size,Precision> my_cholesky;
- int rank;
+ int my_rank;
};
Index: Makefile.in
===================================================================
RCS file: /cvsroot/toon/TooN/Makefile.in,v
retrieving revision 1.16
retrieving revision 1.17
diff -u -b -r1.16 -r1.17
--- Makefile.in 15 Oct 2009 14:51:50 -0000 1.16
+++ Makefile.in 19 Oct 2009 14:08:52 -0000 1.17
@@ -36,7 +36,7 @@
doxygen
-TESTS=lu slice vector_resize gauss_jordan eigen-sqrt determinant
+TESTS=lu slice vector_resize gauss_jordan eigen-sqrt determinant chol_toon
chol_lapack
TEST_RESULT=$(TESTS:%=regressions/%.result)
Index: regressions/chol_lapack.cc
===================================================================
RCS file: regressions/chol_lapack.cc
diff -N regressions/chol_lapack.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ regressions/chol_lapack.cc 19 Oct 2009 14:08:53 -0000 1.1
@@ -0,0 +1,53 @@
+#include <TooN/Lapack_Cholesky.h>
+
+#include <iostream>
+#include <iomanip>
+
+using namespace std;
+using namespace TooN;
+
+int main(int, char ** ){
+
+ cout << setprecision(10);
+
+ Matrix<3> t = Data(
+ 1, 0.5, 0.5,
+ 0.5, 2, 0.7,
+ 0.5, 0.7, 3);
+
+ Lapack_Cholesky<3> chol(t);
+
+ cout << "Check for determinant\n";
+ cout << chol.determinant() << endl << endl;
+
+ cout << "Static size checks:\n";
+ cout << "Check decomposition, both matrices should be the same.\n";
+ cout << t << endl << chol.get_L()*(chol.get_L().T()) << endl ;
+
+ cout << "Check inverse, third matrix should be close to identity.\n";
+ cout << t << "\n" << chol.get_inverse() << "\n"
+ << t * chol.get_inverse() << endl;
+
+ Matrix<> t2 = t;
+
+ Lapack_Cholesky<Dynamic,float> chol2(t2);
+
+ cout << "Dynamic size, single precision checks:\n";
+ cout << "Check decomposition, both matrices should be the same.\n";
+ cout << t << endl << chol2.get_L()*(chol2.get_L().T()) << endl ;
+
+ cout << "Check inverse, third matrix should be close to identity.\n";
+ cout << t2 << "\n" << chol2.get_inverse() << "\n"
+ << t2 * chol2.get_inverse() << endl;
+
+ Vector<3> bla = makeVector(1,2,3);
+
+ cout << "Check backsub(), the following two vectors should be the same.\n";
+ cout << chol.backsub(bla) << endl;
+ cout << chol.get_inverse() * bla << endl << endl;
+
+ cout << "Check mahalanobis(), result should be zero.\n";
+ cout << chol.mahalanobis(bla) - bla * chol.backsub(bla) << endl;
+
+ return 0;
+}
Index: regressions/chol_lapack.txt
===================================================================
RCS file: regressions/chol_lapack.txt
diff -N regressions/chol_lapack.txt
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ regressions/chol_lapack.txt 19 Oct 2009 14:08:53 -0000 1.1
@@ -0,0 +1,58 @@
+#< t 1e-16
+
+Check for determinant
+4.61
+
+Static size checks:
+Check decomposition, both matrices should be the same.
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+Check inverse, third matrix should be close to identity.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+ 1.1952277657266810e+00 -2.4945770065075915e-01 -1.4099783080260306e-01
+-2.4945770065075915e-01 5.9652928416485884e-01 -9.7613882863340537e-02
+-1.4099783080260306e-01 -9.7613882863340537e-02 3.7960954446854661e-01
+
+1 0 0
+0 1 0
+0 0 1
+
+Dynamic size, single precision checks:
+Check decomposition, both matrices should be the same.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+Check inverse, third matrix should be close to identity.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+ 1.1952277657266810e+00 -2.4945770065075915e-01 -1.4099783080260306e-01
+-2.4945770065075915e-01 5.9652928416485884e-01 -9.7613882863340537e-02
+-1.4099783080260306e-01 -9.7613882863340537e-02 3.7960954446854661e-01
+
+1 0 0
+0 1 0
+0 0 1
+
+Check backsub(), the following two vectors should be the same.
+2.7331887201735361e-01 6.5075921908893697e-01 8.0260303687635570e-01
+2.7331887201735361e-01 6.5075921908893697e-01 8.0260303687635570e-01
+
+Check mahalanobis(), result should be zero.
+0
Index: regressions/chol_toon.cc
===================================================================
RCS file: regressions/chol_toon.cc
diff -N regressions/chol_toon.cc
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ regressions/chol_toon.cc 19 Oct 2009 14:08:53 -0000 1.1
@@ -0,0 +1,57 @@
+#include <TooN/Cholesky.h>
+
+#include <iostream>
+#include <iomanip>
+
+using namespace std;
+using namespace TooN;
+
+int main(int, char ** ){
+
+ cout << setprecision(10);
+
+ Matrix<3> t = Data(
+ 1, 0.5, 0.5,
+ 0.5, 2, 0.7,
+ 0.5, 0.7, 3);
+
+ Cholesky<3> chol(t);
+
+ cout << "Check for determinant\n";
+ cout << chol.determinant() << endl << endl;
+
+ cout << "Static size checks:\n";
+ cout << "Check decomposition, all three matrices should be the same.\n";
+ cout << t << endl << chol.get_L()*(chol.get_L().T()) << endl
+ << chol.get_unscaled_L()*chol.get_D()*(chol.get_unscaled_L().T())
+ << endl;
+
+ cout << "Check inverse, third matrix should be close to identity.\n";
+ cout << t << "\n" << chol.get_inverse() << "\n"
+ << t * chol.get_inverse() << endl;
+
+ Matrix<> t2 = t;
+
+ Cholesky<Dynamic,float> chol2(t2);
+
+ cout << "Dynamic size, single precision checks:\n";
+ cout << "Check decomposition, all three matrices should be the same.\n";
+ cout << t << endl << chol2.get_L()*(chol2.get_L().T()) << endl
+ << chol2.get_unscaled_L()*chol2.get_D()*(chol2.get_unscaled_L().T())
+ << endl;
+
+ cout << "Check inverse, third matrix should be close to identity.\n";
+ cout << t2 << "\n" << chol2.get_inverse() << "\n"
+ << t2 * chol2.get_inverse() << endl;
+
+ Vector<3> bla = makeVector(1,2,3);
+
+ cout << "Check backsub(), the following two vectors should be the same.\n";
+ cout << chol.backsub(bla) << endl;
+ cout << chol.get_inverse() * bla << endl << endl;
+
+ cout << "Check mahalanobis(), result should be zero.\n";
+ cout << chol.mahalanobis(bla) - bla * chol.backsub(bla) << endl;
+
+ return 0;
+}
Index: regressions/chol_toon.txt
===================================================================
RCS file: regressions/chol_toon.txt
diff -N regressions/chol_toon.txt
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ regressions/chol_toon.txt 19 Oct 2009 14:08:53 -0000 1.1
@@ -0,0 +1,65 @@
+#< t 1e-16
+
+Check for determinant
+4.61
+
+Static size checks:
+Check decomposition, all three matrices should be the same.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+Check inverse, third matrix should be close to identity.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+ 1.1952277657266810e+00 -2.4945770065075915e-01 -1.4099783080260306e-01
+-2.4945770065075915e-01 5.9652928416485884e-01 -9.7613882863340537e-02
+-1.4099783080260306e-01 -9.7613882863340537e-02 3.7960954446854661e-01
+
+1 0 0
+0 1 0
+0 0 1
+
+Dynamic size, single precision checks:
+Check decomposition, all three matrices should be the same.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+Check inverse, third matrix should be close to identity.
+1 0.5 0.5
+0.5 2 0.7
+0.5 0.7 3
+
+ 1.1952277657266810e+00 -2.4945770065075915e-01 -1.4099783080260306e-01
+-2.4945770065075915e-01 5.9652928416485884e-01 -9.7613882863340537e-02
+-1.4099783080260306e-01 -9.7613882863340537e-02 3.7960954446854661e-01
+
+1 0 0
+0 1 0
+0 0 1
+
+Check backsub(), the following two vectors should be the same.
+2.7331887201735361e-01 6.5075921908893697e-01 8.0260303687635570e-01
+2.7331887201735361e-01 6.5075921908893697e-01 8.0260303687635570e-01
+
+Check mahalanobis(), result should be zero.
+0
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Toon-members] TooN Cholesky.h Lapack_Cholesky.h Makefile.in r...,
Edward Rosten <=