qrdec.h

Go to the documentation of this file.

#ifndef __QRDEC_H__
#define __QRDEC_H__

#include <vector>
#include <iostream>

#ifdef __MPF_BUILD__
namespace mpf_sdp {
#else
namespace sdp {
#endif

template <class T>
class QRDecomposition {
  int _m, _n;             // Number of rows and columns
  T sfactor;              // Scaling factor
  std::vector<T> QR;      // Matrix containing Q and R
  std::vector<T> sigma;   // Reflector parameters & diagonal of R
  std::vector<T> msigma;  // The opposites of the sigma numbers
  std::vector<int> p;     // Pivot vector

public:
  QRDecomposition(const std::vector<T>& A, int m, int n, double sing_eps = 1e-10)
    : _m(m), _n(n), QR(A), p(n)
  {
    if (m < n)
      throw ValueError(ERROR_MSG("Number of rows smaller than of columns"));
    
    /* Finds the scaling factor and scales the matrix */
    sfactor = 0;
    for (int i = 0; i < QR.size(); i++)
      if (NumberOperations::abs(QR[i]) > sfactor)
        sfactor = NumberOperations::abs(QR[i]);

    for (int i = 0; i < QR.size(); i++) QR[i] /= sfactor;

    /* Computes column squared norms */
    std::vector<T> cnorms(n, 0);

    for (int j = 0; j < n; j++)
      for (int i = 0; i < m; i++)
        cnorms[j] += QR[i + j * m] * QR[i + j * m];

    /* Does the QR decomposition with column pivoting */
    sigma.reserve(n);
    for (int j = 0; j < n; j++) p[j] = j;

    for (int k = 0; k < n; k++) {
      /* Pivoting step */
      int pivot = k;
      for (int j = k + 1; j < n; j++)
        if (cnorms[p[j]] > cnorms[p[pivot]])
          pivot = j;

      if (cnorms[p[pivot]] <= sing_eps) break;  // We are finished

      { int tmp = p[k];
        p[k] = p[pivot];
        p[pivot] = tmp;
      }

      /* Computes the reflector */
      T s = NumberOperations::sqrt(cnorms[p[k]]);

      if (QR[k + p[k] * m] < 0) s *= -1;
      QR[k + p[k] * m] += s;

      T gamma_inv = s * QR[k + p[k] * m];
      sigma.push_back(s);

      /* Applies reflector to other columns */
      for (int j = k + 1; j < n; j++) {
        T tau = 0;
        for (int i = k; i < m; i++)
          tau += QR[i + p[k] * m] * QR[i + p[j] * m];

        tau /= gamma_inv;
        for (int i = k; i < m; i++)
          QR[i + p[j] * m] -= tau * QR[i + p[k] * m];
      }

      /* Updates column norms */
      for (int j = k + 1; j < n; j++)
        cnorms[p[j]] -= QR[k + p[j] * m] * QR[k + p[j] * m];
    }

    /* And we fill msigma */
    msigma = sigma;
    for (int i = 0; i < msigma.size(); i++)
      msigma[i] *= -1;
  }
  
  int rank() const { return sigma.size(); }

  int perm(int j) const { return p[j]; }

  const T& scaling_factor() const { return sfactor; }

  const T& R(int i, int j) const 
  { 
    if (i == j) return msigma[i];

    return QR[i + p[j] * _m]; 
  }

  void Q_mul(std::vector<T>& x) const
  {
    for (int k = sigma.size() - 1; k >= 0; k--) {
      T tau = 0;

      for (int i = k; i < _m; i++)
        tau += QR[i + p[k] * _m] * x[i];

      tau /= sigma[k] * QR[k + p[k] * _m];

      for (int i = k; i < _m; i++)
        x[i] -= tau * QR[i + p[k] * _m];
    }
  }

  void Q_transpose_mul(std::vector<T>& x) const
  {
    for (int k = 0; k < sigma.size(); k++) {
      T tau = 0;

      for (int i = k; i < _m; i++)
        tau += QR[i + p[k] * _m] * x[i];

      tau /= sigma[k] * QR[k + p[k] * _m];

      for (int i = k; i < _m; i++)
        x[i] -= tau * QR[i + p[k] * _m];
    }
  }
};

}  // namespace

#endif

// Local variables:
// mode: c++
// End: