Logo ROOT  
Reference Guide
decomposeQR.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_matrix
3 /// \notebook -nodraw
4 ///
5 /// This tutorial shows how to decompose a matrix A in an orthogonal matrix Q and an uppper triangular matrix R
6 /// uisng QR Householder decomposition with the TDecompQRH class
7 /// The matrix same matrix as in this example is used: https://en.wikipedia.org/wiki/QR_decomposition#Example_2
8 
9 #include <iostream>
10 #include "TMath.h"
11 #include "TDecompQRH.h"
12 
13 void decomposeQR() {
14 
15  const int n = 3;
16 
17 
18  double a[] = {12, -51, 4, 6, 167, -68, -4, 24, -41};
19 
20  TMatrixT<double> A(3, 3, a);
21 
22  std::cout << "initial matrox A " << std::endl;
23 
24  A.Print();
25 
26  TDecompQRH decomp(A);
27 
28  bool ret = decomp.Decompose();
29 
30  std::cout << "Orthogonal Q matrix " << std::endl;
31 
32  // note that decomp.GetQ() returns an intenrnal matrix which is not Q defined as A = QR
33  auto Q = decomp.GetOrthogonalMatrix();
34  Q.Print();
35 
36  std::cout << "Upper Triangular R matrix " << std::endl;
37  auto R = decomp.GetR();
38 
39  R.Print();
40 
41  // check that we have a correct Q-R decomposition
42 
43  TMatrixT<double> compA = Q * R;
44 
45  std::cout << "Computed A matrix from Q * R " << std::endl;
46  compA.Print();
47 
48  for (int i = 0; i < A.GetNrows(); ++i) {
49  for (int j = 0; j < A.GetNcols(); ++j) {
50  if (!TMath::AreEqualAbs( compA(i,j), A(i,j), 1.E-6) )
51  Error("decomposeQR","Reconstrate matrix is not equal to the original : %f different than %f",compA(i,j),A(i,j));
52  }
53  }
54 
55  // chech also that Q is orthogonal (Q^T * Q = I)
56  auto QT = Q;
57  QT.Transpose(Q);
58 
59  auto qtq = QT * Q;
60  for (int i = 0; i < Q.GetNrows(); ++i) {
61  for (int j = 0; j < Q.GetNcols(); ++j) {
62  if ((i == j && !TMath::AreEqualAbs(qtq(i, i), 1., 1.E-6)) ||
63  (i != j && !TMath::AreEqualAbs(qtq(i, j), 0., 1.E-6))) {
64  Error("decomposeQR", "Q matrix is not orthogonal ");
65  qtq.Print();
66  break;
67  }
68  }
69  }
70 }
n
const Int_t n
Definition: legend1.C:16
TDecompQRH.h
ROOT::Math::Cephes::A
static double A[]
Definition: SpecFuncCephes.cxx:170
TMatrixT
TMatrixT.
Definition: TMatrixT.h:39
TDecompQRH::GetR
virtual const TMatrixD & GetR()
Definition: TDecompQRH.h:54
R
#define R(a, b, c, d, e, f, g, h, i)
Definition: RSha256.hxx:110
a
auto * a
Definition: textangle.C:12
TDecompQRH::GetOrthogonalMatrix
TMatrixD GetOrthogonalMatrix() const
For a matrix A(m,n), return the OtrhogonalMatrix Q such as A = Q * R.
Definition: TDecompQRH.cxx:516
decomposeQR
void decomposeQR()
Definition: decomposeQR.C:13
TMath::AreEqualAbs
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:424
TDecompQRH
QR Decomposition class.
Definition: TDecompQRH.h:26
TDecompQRH::Decompose
virtual Bool_t Decompose()
QR decomposition of matrix a by Householder transformations, see Golub & Loan first edition p41 & Sec...
Definition: TDecompQRH.cxx:150
TMath.h
Error
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
ROOT::Math::Cephes::Q
static double Q[]
Definition: SpecFuncCephes.cxx:294