Tools/Math/LA_Qr_decomp.h

Go to the documentation of this file.

//------------------------------------------------------------------------------
/**
 * @file LA_Qr_decomp.h
 * Contains class LinAlg::QR
 *
 * @author <a href="mailto:stefanuhrig@gmx.net">Stefan Uhrig</a>
 */
//------------------------------------------------------------------------------
#ifndef LINALG_QR_H_INCLUDED
#define LINALG_QR_H_INCLUDED
//------------------------------------------------------------------------------
#include "LA_Vector.h"
#include "LA_Matrix.h"
//------------------------------------------------------------------------------
namespace LA
{
//------------------------------------------------------------------------------
/** 
  <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).
 */
class QR
{
public:
  /**
   * Create a QR factorization object for A.
   * @param A rectangular (rows>=columns) matrix.
   */
  QR(const Matrix& A);

public:
  /**
   * Flag to denote the matrix is of full rank.
   * @return true if matrix is full rank, false otherwise.
   */
  bool isFullRank() const;
  
  /** 
   * Retreive the Householder vectors from QR factorization
   * @returns lower trapezoidal matrix whose columns define the reflections
   */
  Matrix getHouseholder() const;
  
  /**
   * Return the upper triangular factor, R, of the QR factorization
   * @return R
   */
  Matrix getR() const;
  
  /** 
   * Generate and return the (economy-sized) orthogonal factor
   * @return Q the (ecnomy-sized) orthogonal factor (Q*R=A).
   */
  Matrix getQ() const;

  /**
   * 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.
   */
  Vector solve(const Vector& b) const;


private:
  Matrix QR_;
  int m;
  int n;
  Vector Rdiag;
};
//------------------------------------------------------------------------------
}
//------------------------------------------------------------------------------
#endif
//------------------------------------------------------------------------------

---