decomposeQR.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_matrix
/// \notebook -nodraw
///
/// This tutorial shows how to decompose a matrix A in an orthogonal matrix Q and an uppper triangular matrix R
/// uisng QR Householder decomposition with the TDecompQRH class
/// The matrix same matrix as in this example is used: https://en.wikipedia.org/wiki/QR_decomposition#Example_2
 
#include <iostream>
#include "TMath.h"
#include "TDecompQRH.h"
 
void decomposeQR() {
 
   const int n = 3;
 
 
   double a[] = {12, -51, 4, 6, 167, -68, -4, 24, -41};
 
   TMatrixT<double> A(3, 3, a);
 
   std::cout << "initial matrox A " << std::endl;
 
   A.Print();
 
   TDecompQRH decomp(A);
 
   bool ret = decomp.Decompose();
 
   std::cout << "Orthogonal Q matrix " << std::endl;
 
   // note that decomp.GetQ()  returns an intenrnal matrix which is not Q defined as A = QR
   auto Q = decomp.GetOrthogonalMatrix();
   Q.Print();
 
   std::cout << "Upper Triangular R matrix " << std::endl;
   auto R = decomp.GetR();
 
   R.Print();
 
   // check that we have a correct Q-R decomposition
 
   TMatrixT<double> compA = Q * R;
 
   std::cout << "Computed A matrix from Q * R " << std::endl;
   compA.Print();
 
   for (int i = 0; i < A.GetNrows(); ++i) {
      for (int j = 0; j < A.GetNcols(); ++j) {
         if (!TMath::AreEqualAbs( compA(i,j), A(i,j), 1.E-6) )
            Error("decomposeQR","Reconstrate matrix is not equal to the original : %f different than %f",compA(i,j),A(i,j));
      }
   }
 
   // chech also that Q is orthogonal (Q^T * Q = I)
   auto QT = Q;
   QT.Transpose(Q);
 
   auto qtq = QT * Q;
   for (int i = 0; i < Q.GetNrows(); ++i) {
      for (int j = 0; j < Q.GetNcols(); ++j) {
         if ((i == j && !TMath::AreEqualAbs(qtq(i, i), 1., 1.E-6)) ||
             (i != j && !TMath::AreEqualAbs(qtq(i, j), 0., 1.E-6))) {
            Error("decomposeQR", "Q matrix is not orthogonal ");
            qtq.Print();
            break;
         }
      }
   }
}