{
"cells": [
{
"cell_type": "markdown",
"id": "4b6e0aeb",
"metadata": {},
"source": [
"# decomposeQR\n",
"This tutorial shows how to decompose a matrix A in an orthogonal matrix Q and an upper\n",
"triangular matrix R using QR Householder decomposition with the TDecompQRH class.\n",
"The matrix same matrix as in this example is used: https://en.wikipedia.org/wiki/QR_decomposition#Example_2\n",
"\n",
"\n",
"**Author:** \n",
"This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, March 19, 2024 at 07:13 PM."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "99bbc34c",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:13:08.468597Z",
"iopub.status.busy": "2024-03-19T19:13:08.468103Z",
"iopub.status.idle": "2024-03-19T19:13:09.423855Z",
"shell.execute_reply": "2024-03-19T19:13:09.391250Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"initial matrox A \n",
"\n",
"3x3 matrix is as follows\n",
"\n",
" | 0 | 1 | 2 |\n",
"--------------------------------------------\n",
" 0 | 12 -51 4 \n",
" 1 | 6 167 -68 \n",
" 2 | -4 24 -41 \n",
"\n",
"Orthogonal Q matrix \n"
]
}
],
"source": [
"const int n = 3;\n",
"\n",
"\n",
"double a[] = {12, -51, 4, 6, 167, -68, -4, 24, -41};\n",
"\n",
"TMatrixT A(3, 3, a);\n",
"\n",
"std::cout << \"initial matrox A \" << std::endl;\n",
"\n",
"A.Print();\n",
"\n",
"TDecompQRH decomp(A);\n",
"\n",
"bool ret = decomp.Decompose();\n",
"\n",
"std::cout << \"Orthogonal Q matrix \" << std::endl;"
]
},
{
"cell_type": "markdown",
"id": "bd20e06b",
"metadata": {},
"source": [
"note that decomp.GetQ() returns an intenrnal matrix which is not Q defined as A = QR"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "26fd5756",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:13:09.435066Z",
"iopub.status.busy": "2024-03-19T19:13:09.434644Z",
"iopub.status.idle": "2024-03-19T19:13:09.666741Z",
"shell.execute_reply": "2024-03-19T19:13:09.665368Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"3x3 matrix is as follows\n",
"\n",
" | 0 | 1 | 2 |\n",
"--------------------------------------------\n",
" 0 | -0.8571 0.3943 0.3314 \n",
" 1 | -0.4286 -0.9029 -0.03429 \n",
" 2 | 0.2857 -0.1714 0.9429 \n",
"\n",
"Upper Triangular R matrix \n",
"\n",
"3x3 matrix is as follows\n",
"\n",
" | 0 | 1 | 2 |\n",
"--------------------------------------------\n",
" 0 | -14 -21 14 \n",
" 1 | 0 -175 70 \n",
" 2 | 0 0 -35 \n",
"\n"
]
}
],
"source": [
"auto Q = decomp.GetOrthogonalMatrix();\n",
"Q.Print();\n",
"\n",
"std::cout << \"Upper Triangular R matrix \" << std::endl;\n",
"auto R = decomp.GetR();\n",
"\n",
"R.Print();"
]
},
{
"cell_type": "markdown",
"id": "974c71a5",
"metadata": {},
"source": [
"check that we have a correct Q-R decomposition"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "c115fb0b",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:13:09.677078Z",
"iopub.status.busy": "2024-03-19T19:13:09.676672Z",
"iopub.status.idle": "2024-03-19T19:13:09.888499Z",
"shell.execute_reply": "2024-03-19T19:13:09.887354Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Computed A matrix from Q * R \n",
"\n",
"3x3 matrix is as follows\n",
"\n",
" | 0 | 1 | 2 |\n",
"--------------------------------------------\n",
" 0 | 12 -51 4 \n",
" 1 | 6 167 -68 \n",
" 2 | -4 24 -41 \n",
"\n"
]
}
],
"source": [
"TMatrixT compA = Q * R;\n",
"\n",
"std::cout << \"Computed A matrix from Q * R \" << std::endl;\n",
"compA.Print();\n",
"\n",
"for (int i = 0; i < A.GetNrows(); ++i) {\n",
" for (int j = 0; j < A.GetNcols(); ++j) {\n",
" if (!TMath::AreEqualAbs( compA(i,j), A(i,j), 1.E-6) )\n",
" Error(\"decomposeQR\",\"Reconstrate matrix is not equal to the original : %f different than %f\",compA(i,j),A(i,j));\n",
" }\n",
"}"
]
},
{
"cell_type": "markdown",
"id": "b4bfae56",
"metadata": {},
"source": [
"chech also that Q is orthogonal (Q^T * Q = I)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "3ed43f97",
"metadata": {
"collapsed": false,
"execution": {
"iopub.execute_input": "2024-03-19T19:13:09.894612Z",
"iopub.status.busy": "2024-03-19T19:13:09.894169Z",
"iopub.status.idle": "2024-03-19T19:13:10.155501Z",
"shell.execute_reply": "2024-03-19T19:13:10.139320Z"
}
},
"outputs": [],
"source": [
"auto QT = Q;\n",
"QT.Transpose(Q);\n",
"\n",
"auto qtq = QT * Q;\n",
"for (int i = 0; i < Q.GetNrows(); ++i) {\n",
" for (int j = 0; j < Q.GetNcols(); ++j) {\n",
" if ((i == j && !TMath::AreEqualAbs(qtq(i, i), 1., 1.E-6)) ||\n",
" (i != j && !TMath::AreEqualAbs(qtq(i, j), 0., 1.E-6))) {\n",
" Error(\"decomposeQR\", \"Q matrix is not orthogonal \");\n",
" qtq.Print();\n",
" break;\n",
" }\n",
" }\n",
"}"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "ROOT C++",
"language": "c++",
"name": "root"
},
"language_info": {
"codemirror_mode": "text/x-c++src",
"file_extension": ".C",
"mimetype": " text/x-c++src",
"name": "c++"
}
},
"nbformat": 4,
"nbformat_minor": 5
}