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