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, 1 month ago) by gezelter
File size:	7361 byte(s)
Log Message:	adding jama and tnt libraries for various linear algebra routines
#	User	Rev	Content
1	gezelter	1336	#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