src/math/jama_qr.hpp

#ifndef JAMA_QR_H
#define JAMA_QR_H

#include "tnt_array1d.hpp"
#include "tnt_array2d.hpp"
#include "tnt_math_utils.hpp"

namespace JAMA
{

/** 
<p>
        Classical QR Decompisition:
   for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
   orthogonal matrix Q and an n-by-n upper triangular matrix R so that
   A = Q*R.
<P>
   The QR decompostion always exists, even if the matrix does not have
   full rank, so the constructor will never fail.  The primary use of the
   QR decomposition is in the least squares solution of nonsquare systems
   of simultaneous linear equations.  This will fail if isFullRank()
   returns 0 (false).

<p>
        The Q and R factors can be retrived via the getQ() and getR()
        methods. Furthermore, a solve() method is provided to find the
        least squares solution of Ax=b using the QR factors.  

   <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 QR {


   /** Array for internal storage of decomposition.
   @serial internal array storage.
   */
   
   TNT::Array2D<Real> QR_;

   /** Row and column dimensions.
   @serial column dimension.
   @serial row dimension.
   */
   int m, n;

   /** Array for internal storage of diagonal of R.
   @serial diagonal of R.
   */
   TNT::Array1D<Real> Rdiag;


public:
        
/**
        Create a QR factorization object for A.

        @param A rectangular (m>=n) matrix.
*/
        QR(const TNT::Array2D<Real> &A)         /* constructor */
        {
      QR_ = A.copy();
      m = A.dim1();
      n = A.dim2();
      Rdiag = TNT::Array1D<Real>(n);
          int i=0, j=0, k=0;

      // Main loop.
      for (k = 0; k < n; k++) {
         // Compute 2-norm of k-th column without under/overflow.
         Real nrm = 0;
         for (i = k; i < m; i++) {
            nrm = TNT::hypot(nrm,QR_[i][k]);
         }

         if (nrm != 0.0) {
            // Form k-th Householder vector.
            if (QR_[k][k] < 0) {
               nrm = -nrm;
            }
            for (i = k; i < m; i++) {
               QR_[i][k] /= nrm;
            }
            QR_[k][k] += 1.0;

            // Apply transformation to remaining columns.
            for (j = k+1; j < n; j++) {
               Real s = 0.0; 
               for (i = k; i < m; i++) {
                  s += QR_[i][k]*QR_[i][j];
               }
               s = -s/QR_[k][k];
               for (i = k; i < m; i++) {
                  QR_[i][j] += s*QR_[i][k];
               }
            }
         }
         Rdiag[k] = -nrm;
      }
   }


/**
        Flag to denote the matrix is of full rank.

        @return 1 if matrix is full rank, 0 otherwise.
*/
        int isFullRank() const          
        {
      for (int j = 0; j < n; j++) 
          {
         if (Rdiag[j] == 0)
            return 0;
      }
      return 1;
        }
        
        
   /** 
   
   Retreive the Householder vectors from QR factorization
   @returns lower trapezoidal matrix whose columns define the reflections
   */

   TNT::Array2D<Real> getHouseholder (void)  const
   {
          TNT::Array2D<Real> H(m,n);

          /* note: H is completely filled in by algorithm, so
             initializaiton of H is not necessary.
          */
      for (int i = 0; i < m; i++) 
          {
         for (int j = 0; j < n; j++) 
                 {
            if (i >= j) {
               H[i][j] = QR_[i][j];
            } else {
               H[i][j] = 0.0;
            }
         }
      }
          return H;
   }


   /** Return the upper triangular factor, R, of the QR factorization
   @return     R
   */

        TNT::Array2D<Real> getR() const
        {
      TNT::Array2D<Real> R(n,n);
      for (int i = 0; i < n; i++) {
         for (int j = 0; j < n; j++) {
            if (i < j) {
               R[i][j] = QR_[i][j];
            } else if (i == j) {
               R[i][j] = Rdiag[i];
            } else {
               R[i][j] = 0.0;
            }
         }
      }
          return R;
   }
        
        
   /** 
        Generate and return the (economy-sized) orthogonal factor
   @param     Q the (ecnomy-sized) orthogonal factor (Q*R=A).
   */

        TNT::Array2D<Real> getQ() const
        {
          int i=0, j=0, k=0;

          TNT::Array2D<Real> Q(m,n);
      for (k = n-1; k >= 0; k--) {
         for (i = 0; i < m; i++) {
            Q[i][k] = 0.0;
         }
         Q[k][k] = 1.0;
         for (j = k; j < n; j++) {
            if (QR_[k][k] != 0) {
               Real s = 0.0;
               for (i = k; i < m; i++) {
                  s += QR_[i][k]*Q[i][j];
               }
               s = -s/QR_[k][k];
               for (i = k; i < m; i++) {
                  Q[i][j] += s*QR_[i][k];
               }
            }
         }
      }
          return Q;
   }


   /** Least squares solution of A*x = b
   @param B     m-length array (vector).
   @return x    n-length array (vector) that minimizes the two norm of Q*R*X-B.
                If B is non-conformant, or if QR.isFullRank() is false,
                                                the routine returns a null (0-length) vector.
   */

   TNT::Array1D<Real> solve(const TNT::Array1D<Real> &b) const
   {
          if (b.dim1() != m)            /* arrays must be conformant */
                return TNT::Array1D<Real>();

          if ( !isFullRank() )          /* matrix is rank deficient */
          {
                return TNT::Array1D<Real>();
          }

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

      // Compute Y = transpose(Q)*b
      for (int k = 0; k < n; k++) 
          {
            Real s = 0.0; 
            for (int i = k; i < m; i++) 
                        {
               s += QR_[i][k]*x[i];
            }
            s = -s/QR_[k][k];
            for (int i = k; i < m; i++) 
                        {
               x[i] += s*QR_[i][k];
            }
      }
      // Solve R*X = Y;
      for (int k = n-1; k >= 0; k--) 
          {
         x[k] /= Rdiag[k];
         for (int i = 0; i < k; i++) {
               x[i] -= x[k]*QR_[i][k];
         }
      }


          /* return n x nx portion of X */
          TNT::Array1D<Real> x_(n);
          for (int i=0; i<n; i++)
                        x_[i] = x[i];

          return x_;
   }

   /** Least squares solution of A*X = B
   @param B     m x k Array (must conform).
   @return X     n x k Array that minimizes the two norm of Q*R*X-B. If
                                                B is non-conformant, or if QR.isFullRank() is false,
                                                the routine returns a null (0x0) array.
   */

   TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
   {
          if (B.dim1() != m)            /* arrays must be conformant */
                return TNT::Array2D<Real>(0,0);

          if ( !isFullRank() )          /* matrix is rank deficient */
          {
                return TNT::Array2D<Real>(0,0);
          }

      int nx = B.dim2(); 
          TNT::Array2D<Real> X = B.copy();
          int i=0, j=0, k=0;

      // Compute Y = transpose(Q)*B
      for (k = 0; k < n; k++) {
         for (j = 0; j < nx; j++) {
            Real s = 0.0; 
            for (i = k; i < m; i++) {
               s += QR_[i][k]*X[i][j];
            }
            s = -s/QR_[k][k];
            for (i = k; i < m; i++) {
               X[i][j] += s*QR_[i][k];
            }
         }
      }
      // Solve R*X = Y;
      for (k = n-1; k >= 0; k--) {
         for (j = 0; j < nx; j++) {
            X[k][j] /= Rdiag[k];
         }
         for (i = 0; i < k; i++) {
            for (j = 0; j < nx; j++) {
               X[i][j] -= X[k][j]*QR_[i][k];
            }
         }
      }


          /* return n x nx portion of X */
          TNT::Array2D<Real> X_(n,nx);
          for (i=0; i<n; i++)
                for (j=0; j<nx; j++)
                        X_[i][j] = X[i][j];

          return X_;
   }


};


}
// namespace JAMA

#endif
// JAMA_QR__H

Revision:	1336
Committed:	Tue Apr 7 20:16:26 2009 UTC (16 years, 6 months ago) by gezelter
File size:	7361 byte(s)
Log Message:	adding jama and tnt libraries for various linear algebra routines
#	Content
1	#ifndef JAMA_QR_H
2	#define JAMA_QR_H
3
4	#include "tnt_array1d.hpp"
5	#include "tnt_array2d.hpp"
6	#include "tnt_math_utils.hpp"
7
8	namespace JAMA
9	{
10
11	/**
12	<p>
13	Classical QR Decompisition:
14	for an m-by-n matrix A with m >= n, the QR decomposition is an m-by-n
15	orthogonal matrix Q and an n-by-n upper triangular matrix R so that
16	A = Q*R.
17	<P>
18	The QR decompostion always exists, even if the matrix does not have
19	full rank, so the constructor will never fail. The primary use of the
20	QR decomposition is in the least squares solution of nonsquare systems
21	of simultaneous linear equations. This will fail if isFullRank()
22	returns 0 (false).
23
24	<p>
25	The Q and R factors can be retrived via the getQ() and getR()
26	methods. Furthermore, a solve() method is provided to find the
27	least squares solution of Ax=b using the QR factors.
28
29	<p>
30	(Adapted from JAMA, a Java Matrix Library, developed by jointly
31	by the Mathworks and NIST; see http://math.nist.gov/javanumerics/jama).
32	*/
33
34	template <class Real>
35	class QR {
36
37
38	/** Array for internal storage of decomposition.
39	@serial internal array storage.
40	*/
41
42	TNT::Array2D<Real> QR_;
43
44	/** Row and column dimensions.
45	@serial column dimension.
46	@serial row dimension.
47	*/
48	int m, n;
49
50	/** Array for internal storage of diagonal of R.
51	@serial diagonal of R.
52	*/
53	TNT::Array1D<Real> Rdiag;
54
55
56	public:
57
58	/**
59	Create a QR factorization object for A.
60
61	@param A rectangular (m>=n) matrix.
62	*/
63	QR(const TNT::Array2D<Real> &A) /* constructor */
64	{
65	QR_ = A.copy();
66	m = A.dim1();
67	n = A.dim2();
68	Rdiag = TNT::Array1D<Real>(n);
69	int i=0, j=0, k=0;
70
71	// Main loop.
72	for (k = 0; k < n; k++) {
73	// Compute 2-norm of k-th column without under/overflow.
74	Real nrm = 0;
75	for (i = k; i < m; i++) {
76	nrm = TNT::hypot(nrm,QR_[i][k]);
77	}
78
79	if (nrm != 0.0) {
80	// Form k-th Householder vector.
81	if (QR_[k][k] < 0) {
82	nrm = -nrm;
83	}
84	for (i = k; i < m; i++) {
85	QR_[i][k] /= nrm;
86	}
87	QR_[k][k] += 1.0;
88
89	// Apply transformation to remaining columns.
90	for (j = k+1; j < n; j++) {
91	Real s = 0.0;
92	for (i = k; i < m; i++) {
93	s += QR_[i][k]*QR_[i][j];
94	}
95	s = -s/QR_[k][k];
96	for (i = k; i < m; i++) {
97	QR_[i][j] += s*QR_[i][k];
98	}
99	}
100	}
101	Rdiag[k] = -nrm;
102	}
103	}
104
105
106	/**
107	Flag to denote the matrix is of full rank.
108
109	@return 1 if matrix is full rank, 0 otherwise.
110	*/
111	int isFullRank() const
112	{
113	for (int j = 0; j < n; j++)
114	{
115	if (Rdiag[j] == 0)
116	return 0;
117	}
118	return 1;
119	}
120
121
122
123
124	/**
125
126	Retreive the Householder vectors from QR factorization
127	@returns lower trapezoidal matrix whose columns define the reflections
128	*/
129
130	TNT::Array2D<Real> getHouseholder (void) const
131	{
132	TNT::Array2D<Real> H(m,n);
133
134	/* note: H is completely filled in by algorithm, so
135	initializaiton of H is not necessary.
136	*/
137	for (int i = 0; i < m; i++)
138	{
139	for (int j = 0; j < n; j++)
140	{
141	if (i >= j) {
142	H[i][j] = QR_[i][j];
143	} else {
144	H[i][j] = 0.0;
145	}
146	}
147	}
148	return H;
149	}
150
151
152
153	/** Return the upper triangular factor, R, of the QR factorization
154	@return R
155	*/
156
157	TNT::Array2D<Real> getR() const
158	{
159	TNT::Array2D<Real> R(n,n);
160	for (int i = 0; i < n; i++) {
161	for (int j = 0; j < n; j++) {
162	if (i < j) {
163	R[i][j] = QR_[i][j];
164	} else if (i == j) {
165	R[i][j] = Rdiag[i];
166	} else {
167	R[i][j] = 0.0;
168	}
169	}
170	}
171	return R;
172	}
173
174
175
176
177
178	/**
179	Generate and return the (economy-sized) orthogonal factor
180	@param Q the (ecnomy-sized) orthogonal factor (Q*R=A).
181	*/
182
183	TNT::Array2D<Real> getQ() const
184	{
185	int i=0, j=0, k=0;
186
187	TNT::Array2D<Real> Q(m,n);
188	for (k = n-1; k >= 0; k--) {
189	for (i = 0; i < m; i++) {
190	Q[i][k] = 0.0;
191	}
192	Q[k][k] = 1.0;
193	for (j = k; j < n; j++) {
194	if (QR_[k][k] != 0) {
195	Real s = 0.0;
196	for (i = k; i < m; i++) {
197	s += QR_[i][k]*Q[i][j];
198	}
199	s = -s/QR_[k][k];
200	for (i = k; i < m; i++) {
201	Q[i][j] += s*QR_[i][k];
202	}
203	}
204	}
205	}
206	return Q;
207	}
208
209
210	/** Least squares solution of A*x = b
211	@param B m-length array (vector).
212	@return x n-length array (vector) that minimizes the two norm of QRX-B.
213	If B is non-conformant, or if QR.isFullRank() is false,
214	the routine returns a null (0-length) vector.
215	*/
216
217	TNT::Array1D<Real> solve(const TNT::Array1D<Real> &b) const
218	{
219	if (b.dim1() != m) /* arrays must be conformant */
220	return TNT::Array1D<Real>();
221
222	if ( !isFullRank() ) /* matrix is rank deficient */
223	{
224	return TNT::Array1D<Real>();
225	}
226
227	TNT::Array1D<Real> x = b.copy();
228
229	// Compute Y = transpose(Q)*b
230	for (int k = 0; k < n; k++)
231	{
232	Real s = 0.0;
233	for (int i = k; i < m; i++)
234	{
235	s += QR_[i][k]*x[i];
236	}
237	s = -s/QR_[k][k];
238	for (int i = k; i < m; i++)
239	{
240	x[i] += s*QR_[i][k];
241	}
242	}
243	// Solve R*X = Y;
244	for (int k = n-1; k >= 0; k--)
245	{
246	x[k] /= Rdiag[k];
247	for (int i = 0; i < k; i++) {
248	x[i] -= x[k]*QR_[i][k];
249	}
250	}
251
252
253	/* return n x nx portion of X */
254	TNT::Array1D<Real> x_(n);
255	for (int i=0; i<n; i++)
256	x_[i] = x[i];
257
258	return x_;
259	}
260
261	/** Least squares solution of A*X = B
262	@param B m x k Array (must conform).
263	@return X n x k Array that minimizes the two norm of QRX-B. If
264	B is non-conformant, or if QR.isFullRank() is false,
265	the routine returns a null (0x0) array.
266	*/
267
268	TNT::Array2D<Real> solve(const TNT::Array2D<Real> &B) const
269	{
270	if (B.dim1() != m) /* arrays must be conformant */
271	return TNT::Array2D<Real>(0,0);
272
273	if ( !isFullRank() ) /* matrix is rank deficient */
274	{
275	return TNT::Array2D<Real>(0,0);
276	}
277
278	int nx = B.dim2();
279	TNT::Array2D<Real> X = B.copy();
280	int i=0, j=0, k=0;
281
282	// Compute Y = transpose(Q)*B
283	for (k = 0; k < n; k++) {
284	for (j = 0; j < nx; j++) {
285	Real s = 0.0;
286	for (i = k; i < m; i++) {
287	s += QR_[i][k]*X[i][j];
288	}
289	s = -s/QR_[k][k];
290	for (i = k; i < m; i++) {
291	X[i][j] += s*QR_[i][k];
292	}
293	}
294	}
295	// Solve R*X = Y;
296	for (k = n-1; k >= 0; k--) {
297	for (j = 0; j < nx; j++) {
298	X[k][j] /= Rdiag[k];
299	}
300	for (i = 0; i < k; i++) {
301	for (j = 0; j < nx; j++) {
302	X[i][j] -= X[k][j]*QR_[i][k];
303	}
304	}
305	}
306
307
308	/* return n x nx portion of X */
309	TNT::Array2D<Real> X_(n,nx);
310	for (i=0; i<n; i++)
311	for (j=0; j<nx; j++)
312	X_[i][j] = X[i][j];
313
314	return X_;
315	}
316
317
318	};
319
320
321	}
322	// namespace JAMA
323
324	#endif
325	// JAMA_QR__H
326