src/math/jama_cholesky.hpp

#ifndef JAMA_CHOLESKY_H
#define JAMA_CHOLESKY_H

#include <cmath>
        /* needed for sqrt() below. */


namespace JAMA
{

using namespace TNT;

/** 
   <P>
   For a symmetric, positive definite matrix A, this function
   computes the Cholesky factorization, i.e. it computes a lower 
   triangular matrix L such that A = L*L'.
   If the matrix is not symmetric or positive definite, the function
   computes only a partial decomposition.  This can be tested with
   the is_spd() flag.

   <p>Typical usage looks like:
   <pre>
        Array2D<double> A(n,n);
        Array2D<double> L;

         ... 

        Cholesky<double> chol(A);

        if (chol.is_spd())
                L = chol.getL();
                
        else
                cout << "factorization was not complete.\n";

        </pre>


   <p>
        (Adapted from JAMA, a Java Matrix Library, developed by jointly 
        by the Mathworks and NIST; see  http://math.nist.gov/javanumerics/jama).

   */

template <class Real>
class Cholesky
{
        Array2D<Real> L_;               // lower triangular factor
        int isspd;                              // 1 if matrix to be factored was SPD

public:

        Cholesky();
        Cholesky(const Array2D<Real> &A);
        Array2D<Real> getL() const;
        Array1D<Real> solve(const Array1D<Real> &B);
        Array2D<Real> solve(const Array2D<Real> &B);
        int is_spd() const;

};

template <class Real>
Cholesky<Real>::Cholesky() : L_(0,0), isspd(0) {}

/**
        @return 1, if original matrix to be factored was symmetric 
                positive-definite (SPD).
*/
template <class Real>
int Cholesky<Real>::is_spd() const
{
        return isspd;
}

/**
        @return the lower triangular factor, L, such that L*L'=A.
*/
template <class Real>
Array2D<Real> Cholesky<Real>::getL() const
{
        return L_;
}

/**
        Constructs a lower triangular matrix L, such that L*L'= A.
        If A is not symmetric positive-definite (SPD), only a
        partial factorization is performed.  If is_spd()
        evalutate true (1) then the factorizaiton was successful.
*/
template <class Real>
Cholesky<Real>::Cholesky(const Array2D<Real> &A)
{


        int m = A.dim1();
        int n = A.dim2();
        
        isspd = (m == n);

        if (m != n)
        {
                L_ = Array2D<Real>(0,0);
                return;
        }

        L_ = Array2D<Real>(n,n);


      // Main loop.
     for (int j = 0; j < n; j++) 
         {
        Real d(0.0);
        for (int k = 0; k < j; k++) 
                {
            Real s(0.0);
            for (int i = 0; i < k; i++) 
                        {
               s += L_[k][i]*L_[j][i];
            }
            L_[j][k] = s = (A[j][k] - s)/L_[k][k];
            d = d + s*s;
            isspd = isspd && (A[k][j] == A[j][k]); 
         }
         d = A[j][j] - d;
         isspd = isspd && (d > 0.0);
         L_[j][j] = sqrt(d > 0.0 ? d : 0.0);
         for (int k = j+1; k < n; k++) 
                 {
            L_[j][k] = 0.0;
         }
        }
}

/**

        Solve a linear system A*x = b, using the previously computed
        cholesky factorization of A: L*L'.

   @param  B   A Matrix with as many rows as A and any number of columns.
   @return     x so that L*L'*x = b.  If b is nonconformat, or if A
                                was not symmetric posidtive definite, a null (0x0)
                                                array is returned.
*/
template <class Real>
Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
{
        int n = L_.dim1();
        if (b.dim1() != n)
                return Array1D<Real>();


        Array1D<Real> x = b.copy();


      // Solve L*y = b;
      for (int k = 0; k < n; k++) 
          {
         for (int i = 0; i < k; i++) 
               x[k] -= x[i]*L_[k][i];
                 x[k] /= L_[k][k];
                
      }

      // Solve L'*X = Y;
      for (int k = n-1; k >= 0; k--) 
          {
         for (int i = k+1; i < n; i++) 
               x[k] -= x[i]*L_[i][k];
         x[k] /= L_[k][k];
      }

        return x;
}


/**

        Solve a linear system A*X = B, using the previously computed
        cholesky factorization of A: L*L'.

   @param  B   A Matrix with as many rows as A and any number of columns.
   @return     X so that L*L'*X = B.  If B is nonconformat, or if A
                                was not symmetric posidtive definite, a null (0x0)
                                                array is returned.
*/
template <class Real>
Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
{
        int n = L_.dim1();
        if (B.dim1() != n)
                return Array2D<Real>();


        Array2D<Real> X = B.copy();
        int nx = B.dim2();

// Cleve's original code
#if 0
      // Solve L*Y = B;
      for (int k = 0; k < n; k++) {
         for (int i = k+1; i < n; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*L_[k][i];
            }
         }
         for (int j = 0; j < nx; j++) {
            X[k][j] /= L_[k][k];
         }
      }

      // Solve L'*X = Y;
      for (int k = n-1; k >= 0; k--) {
         for (int j = 0; j < nx; j++) {
            X[k][j] /= L_[k][k];
         }
         for (int i = 0; i < k; i++) {
            for (int j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*L_[k][i];
            }
         }
      }
#endif


      // Solve L*y = b;
          for (int j=0; j< nx; j++)
          {
        for (int k = 0; k < n; k++) 
                {
                        for (int i = 0; i < k; i++) 
               X[k][j] -= X[i][j]*L_[k][i];
                    X[k][j] /= L_[k][k];
                 }
      }

      // Solve L'*X = Y;
     for (int j=0; j<nx; j++)
         {
        for (int k = n-1; k >= 0; k--) 
                {
                for (int i = k+1; i < n; i++) 
               X[k][j] -= X[i][j]*L_[i][k];
                X[k][j] /= L_[k][k];
                }
      }


        return X;
}


}
// namespace JAMA

#endif
// JAMA_CHOLESKY_H
Revision:	1336
Committed:	Tue Apr 7 20:16:26 2009 UTC (16 years, 1 month ago) by gezelter
File size:	5232 byte(s)
Log Message:	adding jama and tnt libraries for various linear algebra routines
#	User	Rev	Content
1	gezelter	1336	#ifndef JAMA_CHOLESKY_H
2			#define JAMA_CHOLESKY_H
3
4			#include <cmath>
5			/* needed for sqrt() below. */
6
7
8			namespace JAMA
9			{
10
11			using namespace TNT;
12
13			/**
14			<P>
15			For a symmetric, positive definite matrix A, this function
16			computes the Cholesky factorization, i.e. it computes a lower
17			triangular matrix L such that A = L*L'.
18			If the matrix is not symmetric or positive definite, the function
19			computes only a partial decomposition. This can be tested with
20			the is_spd() flag.
21
22			<p>Typical usage looks like:
23			<pre>
24			Array2D<double> A(n,n);
25			Array2D<double> L;
26
27			...
28
29			Cholesky<double> chol(A);
30
31			if (chol.is_spd())
32			L = chol.getL();
33
34			else
35			cout << "factorization was not complete.\n";
36
37			</pre>
38
39
40			<p>
41			(Adapted from JAMA, a Java Matrix Library, developed by jointly
42			by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
43
44			*/
45
46			template <class Real>
47			class Cholesky
48			{
49			Array2D<Real> L_; // lower triangular factor
50			int isspd; // 1 if matrix to be factored was SPD
51
52			public:
53
54			Cholesky();
55			Cholesky(const Array2D<Real> &A);
56			Array2D<Real> getL() const;
57			Array1D<Real> solve(const Array1D<Real> &B);
58			Array2D<Real> solve(const Array2D<Real> &B);
59			int is_spd() const;
60
61			};
62
63			template <class Real>
64			Cholesky<Real>::Cholesky() : L_(0,0), isspd(0) {}
65
66			/**
67			@return 1, if original matrix to be factored was symmetric
68			positive-definite (SPD).
69			*/
70			template <class Real>
71			int Cholesky<Real>::is_spd() const
72			{
73			return isspd;
74			}
75
76			/**
77			@return the lower triangular factor, L, such that L*L'=A.
78			*/
79			template <class Real>
80			Array2D<Real> Cholesky<Real>::getL() const
81			{
82			return L_;
83			}
84
85			/**
86			Constructs a lower triangular matrix L, such that L*L'= A.
87			If A is not symmetric positive-definite (SPD), only a
88			partial factorization is performed. If is_spd()
89			evalutate true (1) then the factorizaiton was successful.
90			*/
91			template <class Real>
92			Cholesky<Real>::Cholesky(const Array2D<Real> &A)
93			{
94
95
96			int m = A.dim1();
97			int n = A.dim2();
98
99			isspd = (m == n);
100
101			if (m != n)
102			{
103			L_ = Array2D<Real>(0,0);
104			return;
105			}
106
107			L_ = Array2D<Real>(n,n);
108
109
110			// Main loop.
111			for (int j = 0; j < n; j++)
112			{
113			Real d(0.0);
114			for (int k = 0; k < j; k++)
115			{
116			Real s(0.0);
117			for (int i = 0; i < k; i++)
118			{
119			s += L_[k][i]*L_[j][i];
120			}
121			L_[j][k] = s = (A[j][k] - s)/L_[k][k];
122			d = d + s*s;
123			isspd = isspd && (A[k][j] == A[j][k]);
124			}
125			d = A[j][j] - d;
126			isspd = isspd && (d > 0.0);
127			L_[j][j] = sqrt(d > 0.0 ? d : 0.0);
128			for (int k = j+1; k < n; k++)
129			{
130			L_[j][k] = 0.0;
131			}
132			}
133			}
134
135			/**
136
137			Solve a linear system A*x = b, using the previously computed
138			cholesky factorization of A: L*L'.
139
140			@param B A Matrix with as many rows as A and any number of columns.
141			@return x so that LL'x = b. If b is nonconformat, or if A
142			was not symmetric posidtive definite, a null (0x0)
143			array is returned.
144			*/
145			template <class Real>
146			Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
147			{
148			int n = L_.dim1();
149			if (b.dim1() != n)
150			return Array1D<Real>();
151
152
153			Array1D<Real> x = b.copy();
154
155
156			// Solve L*y = b;
157			for (int k = 0; k < n; k++)
158			{
159			for (int i = 0; i < k; i++)
160			x[k] -= x[i]*L_[k][i];
161			x[k] /= L_[k][k];
162
163			}
164
165			// Solve L'*X = Y;
166			for (int k = n-1; k >= 0; k--)
167			{
168			for (int i = k+1; i < n; i++)
169			x[k] -= x[i]*L_[i][k];
170			x[k] /= L_[k][k];
171			}
172
173			return x;
174			}
175
176
177			/**
178
179			Solve a linear system A*X = B, using the previously computed
180			cholesky factorization of A: L*L'.
181
182			@param B A Matrix with as many rows as A and any number of columns.
183			@return X so that LL'X = B. If B is nonconformat, or if A
184			was not symmetric posidtive definite, a null (0x0)
185			array is returned.
186			*/
187			template <class Real>
188			Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
189			{
190			int n = L_.dim1();
191			if (B.dim1() != n)
192			return Array2D<Real>();
193
194
195			Array2D<Real> X = B.copy();
196			int nx = B.dim2();
197
198			// Cleve's original code
199			#if 0
200			// Solve L*Y = B;
201			for (int k = 0; k < n; k++) {
202			for (int i = k+1; i < n; i++) {
203			for (int j = 0; j < nx; j++) {
204			X[i][j] -= X[k][j]*L_[k][i];
205			}
206			}
207			for (int j = 0; j < nx; j++) {
208			X[k][j] /= L_[k][k];
209			}
210			}
211
212			// Solve L'*X = Y;
213			for (int k = n-1; k >= 0; k--) {
214			for (int j = 0; j < nx; j++) {
215			X[k][j] /= L_[k][k];
216			}
217			for (int i = 0; i < k; i++) {
218			for (int j = 0; j < nx; j++) {
219			X[i][j] -= X[k][j]*L_[k][i];
220			}
221			}
222			}
223			#endif
224
225
226			// Solve L*y = b;
227			for (int j=0; j< nx; j++)
228			{
229			for (int k = 0; k < n; k++)
230			{
231			for (int i = 0; i < k; i++)
232			X[k][j] -= X[i][j]*L_[k][i];
233			X[k][j] /= L_[k][k];
234			}
235			}
236
237			// Solve L'*X = Y;
238			for (int j=0; j<nx; j++)
239			{
240			for (int k = n-1; k >= 0; k--)
241			{
242			for (int i = k+1; i < n; i++)
243			X[k][j] -= X[i][j]*L_[i][k];
244			X[k][j] /= L_[k][k];
245			}
246			}
247
248
249
250			return X;
251			}
252
253
254			}
255			// namespace JAMA
256
257			#endif
258			// JAMA_CHOLESKY_H