# Decompose Q R¶

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

Author:
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Sunday, September 20, 2020 at 08:20 AM.

In [1]:
%%cpp -d
#include <iostream>
#include "TMath.h"
#include "TDecompQRH.h"

In [2]:
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;

initial matrox A

3x3 matrix is as follows

|      0    |      1    |      2    |
--------------------------------------------
0 |         12         -51           4
1 |          6         167         -68
2 |         -4          24         -41

Orthogonal Q matrix


Note that decomp.getq() returns an intenrnal matrix which is not q defined as a = qr

In [3]:
auto Q = decomp.GetOrthogonalMatrix();
Q.Print();

std::cout << "Upper Triangular R matrix " << std::endl;
auto R = decomp.GetR();

R.Print();

3x3 matrix is as follows

|      0    |      1    |      2    |
--------------------------------------------
0 |    -0.8571      0.3943      0.3314
1 |    -0.4286     -0.9029    -0.03429
2 |     0.2857     -0.1714      0.9429

Upper Triangular R matrix

3x3 matrix is as follows

|      0    |      1    |      2    |
--------------------------------------------
0 |        -14         -21          14
1 |          0        -175          70
2 |          0           0         -35



Check that we have a correct q-r decomposition

In [4]:
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));
}
}

Computed A matrix from Q * R

3x3 matrix is as follows

|      0    |      1    |      2    |
--------------------------------------------
0 |         12         -51           4
1 |          6         167         -68
2 |         -4          24         -41



Chech also that q is orthogonal (q^t * q = i)

In [5]:
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;
}
}
}