qr.h

/*
This file is distributed as part of the BIAS library (Basic ImageAlgorithmS)
but it has not been developed by the authors of BIAS.

For copyright, author and license information see below.
*/


/*
*
* Template Numerical Toolkit (TNT): Linear Algebra Module
*
* Mathematical and Computational Sciences Division
* National Institute of Technology,
* Gaithersburg, MD USA
*
*
* This software was developed at the National Institute of Standards and
* Technology (NIST) by employees of the Federal Government in the course
* of their official duties. Pursuant to title 17 Section 105 of the
* United States Code, this software is not subject to copyright protection
* and is in the public domain.  The Template Numerical Toolkit (TNT) is
* an experimental system.  NIST assumes no responsibility whatsoever for
* its use by other parties, and makes no guarantees, expressed or implied,
* about its quality, reliability, or any other characteristic.
*
* BETA VERSION INCOMPLETE AND SUBJECT TO CHANGE
* see http://math.nist.gov/tnt for latest updates.
*
*/


#ifndef QR_H
#define QR_H

// Classical QR factorization example, based on Stewart[1973].
//
//
// This algorithm computes the factorization of a matrix A
// into a product of an orthognal matrix (Q) and an upper triangular 
// matrix (R), such that QR = A.
//
// Parameters:
//
//  A   (in):       Matrix(1:N, 1:N)
//
//  Q   (output):   Matrix(1:N, 1:N), collection of Householder
//                      column vectors Q1, Q2, ... QN
//
//  R   (output):   upper triangular Matrix(1:N, 1:N)
//
// Returns:  
//
//  0 if successful, 1 if A is detected to be singular
//


#include <cmath>      //for sqrt() & fabs()
#include "tntmath.h"  // for sign()

// Classical QR factorization, based on Stewart[1973].
//
//
// This algorithm computes the factorization of a matrix A
// into a product of an orthognal matrix (Q) and an upper triangular 
// matrix (R), such that QR = A.
//
// Parameters:
//
//  A   (in/out):  On input, A is square, Matrix(1:N, 1:N), that represents
//                  the matrix to be factored.
//
//                 On output, Q and R is encoded in the same Matrix(1:N,1:N)
//                 in the following manner:
//
//                  R is contained in the upper triangular section of A,
//                  except that R's main diagonal is in D.  The lower
//                  triangular section of A represents Q, where each
//                  column j is the vector  Qj = I - uj*uj'/pi_j.
//
//  C  (output):    vector of Pi[j]
//  D  (output):    main diagonal of R, i.e. D(i) is R(i,i)
//
// Returns:  
//
//  0 if successful, 1 if A is detected to be singular
//

namespace TNT
{

  template <class MaTRiX, class Vector>
    int QR_factor(MaTRiX &A, Vector& C, Vector &D)
  {
    BIASASSERT(A.lbound() == 1);        // ensure these are all 
    BIASASSERT(C.lbound() == 1);        // 1-based arrays and vectors
    BIASASSERT(D.lbound() == 1);

    Subscript M = A.num_rows();
    Subscript N = A.num_cols(); 

    BIASASSERT(M == N);                 // make sure A is square

    Subscript i,j,k;
    typename MaTRiX::element_type eta, sigma, sum;

    // adjust the shape of C and D, if needed...

    if (N != C.size())  C.newsize(N);
    if (N != D.size())  D.newsize(N);

    for (k=1; k<N; k++)
    {
      // eta = max |M(i,k)|,  for k <= i <= n
      //
      eta = 0;
      for (i=k; i<=N; i++)
      {
        double absA = fabs(A(i,k));
        eta = ( absA >  eta ? absA : eta ); 
      }

      if (eta == 0)           // matrix is singular
      {
        cerr << "QR: k=" << k << "\n";
        return 1;
      }

      // form Qk and premiltiply M by it
      //
      for(i=k; i<=N; i++)
        A(i,k)  = A(i,k) / eta;

      sum = 0;
      for (i=k; i<=N; i++)
        sum = sum + A(i,k)*A(i,k);
      sigma = sign(A(k,k)) *  sqrt(sum);


      A(k,k) = A(k,k) + sigma;
      C(k) = sigma * A(k,k);
      D(k) = -eta * sigma;

      for (j=k+1; j<=N; j++)
      {
        sum = 0;
        for (i=k; i<=N; i++)
          sum = sum + A(i,k)*A(i,j);
        sum = sum / C(k);

        for (i=k; i<=N; i++)
          A(i,j) = A(i,j) - sum * A(i,k);
      }

      D(N) = A(N,N);
    }

    return 0;
  }

  // modified form of upper triangular solve, except that the main diagonal
  // of R (upper portion of A) is in D.
  //
  template <class MaTRiX, class Vector>
    int R_solve(const MaTRiX &A, /*const*/ Vector &D, Vector &b)
  {
    BIASASSERT(A.lbound() == 1);        // ensure these are all 
    BIASASSERT(D.lbound() == 1);        // 1-based arrays and vectors
    BIASASSERT(b.lbound() == 1);

    Subscript i,j;
    Subscript N = A.num_rows();

    BIASASSERT(N == A.num_cols());
    BIASASSERT(N == D.dim());
    BIASASSERT(N == b.dim());

    typename MaTRiX::element_type sum;

    if (D(N) == 0)
      return 1;

    b(N) = b(N) / 
      D(N);

    for (i=N-1; i>=1; i--)
    {
      if (D(i) == 0)
        return 1;
      sum = 0;
      for (j=i+1; j<=N; j++)
        sum = sum + A(i,j)*b(j);
      b(i) = ( b(i) - sum ) / 
        D(i);
    }

    return 0;
  }


  template <class MaTRiX, class Vector>
    int QR_solve(const MaTRiX &A, const Vector &c, /*const*/ Vector &d, 
    Vector &b)
  {
    BIASASSERT(A.lbound() == 1);        // ensure these are all 
    BIASASSERT(c.lbound() == 1);        // 1-based arrays and vectors
    BIASASSERT(d.lbound() == 1);

    Subscript N=A.num_rows();

    BIASASSERT(N == A.num_cols());
    BIASASSERT(N == c.dim());
    BIASASSERT(N == d.dim());
    BIASASSERT(N == b.dim());

    Subscript i,j;
    typename MaTRiX::element_type sum, tau;

    for (j=1; j<N; j++)
    {
      // form Q'*b
      sum = 0;
      for (i=j; i<=N; i++)
        sum = sum + A(i,j)*b(i);
      if (c(j) == 0)
        return 1;
      tau = sum / c(j);
      for (i=j; i<=N; i++)
        b(i) = b(i) - tau * A(i,j);
    }
    return R_solve(A, d, b);        // solve Rx = Q'b
  }

} // namespace TNT

#endif
// QR_H