{ "cells": [ { "cell_type": "markdown", "id": "7ce92596", "metadata": {}, "source": [ "# solveLinear\n", "This macro shows several ways to perform a linear least-squares\n", "analysis . To keep things simple we fit a straight line to 4\n", "data points\n", "The first 4 methods use the linear algebra package to find\n", " x such that min $ (A x - b)^T (A x - b) $ where A and b\n", " are calculated with the data points and the functional expression :\n", "\n", " 1. Normal equations:\n", " Expanding the expression $ (A x - b)^T (A x - b) $ and taking the\n", " derivative wrt x leads to the \"Normal Equations\":\n", " $ A^T A x = A^T b $ where $ A^T A $ is a positive definite matrix. Therefore,\n", " a Cholesky decomposition scheme can be used to calculate its inverse .\n", " This leads to the solution $ x = (A^T A)^-1 A^T b $ . All this is done in\n", " routine NormalEqn . We made it a bit more complicated by giving the\n", " data weights .\n", " Numerically this is not the best way to proceed because effectively the\n", " condition number of $ A^T A $ is twice as large as that of A, making inversion\n", " more difficult\n", "\n", " 2. SVD\n", " One can show that solving $ A x = b $ for x with A of size $ (m x n) $\n", " and $ m > n $ through a Singular Value Decomposition is equivalent to minimizing\n", " $ (A x - b)^T (A x - b) $ Numerically , this is the most stable method of all 5\n", "\n", " 3. Pseudo Inverse\n", " Here we calculate the generalized matrix inverse (\"pseudo inverse\") by\n", " solving $ A X = Unit $ for matrix $ X $ through an SVD . The formal expression for\n", " is $ X = (A^T A)^-1 A^T $ . Then we multiply it by $ b $ .\n", " Numerically, not as good as 2 and not as fast . In general it is not a\n", " good idea to solve a set of linear equations with a matrix inversion .\n", "\n", " 4. Pseudo Inverse , brute force\n", " The pseudo inverse is calculated brute force through a series of matrix\n", " manipulations . It shows nicely some operations in the matrix package,\n", " but is otherwise a big \"no no\" .\n", "\n", " 5. Least-squares analysis with Minuit\n", " An objective function L is minimized by Minuit, where\n", " $ L = sum_i { (y - c_0 -c_1 * x / e)^2 } $\n", " Minuit will calculate numerically the derivative of L wrt c_0 and c_1 .\n", " It has not been told that these derivatives are linear in the parameters\n", " c_0 and c_1 .\n", " For ill-conditioned linear problems it is better to use the fact it is\n", " a linear fit as in 2 .\n", "\n", "Another interesting thing is the way we assign data to the vectors and\n", "matrices through adoption .\n", "This allows data assignment without physically moving bytes around .\n", "\n", " #### USAGE\n", "\n", "This macro can be executed via CLING or via ACLIC\n", "- via the interpreter, do\n", "```cpp\n", " root > .x solveLinear.C\n", "```\n", "- via ACLIC\n", "```cpp\n", " root > gSystem->Load(\"libMatrix\");\n", " root > gSystem->Load(\"libGpad\");\n", " root > .x solveLinear.C+\n", "```\n", "\n", "\n", "\n", "\n", "**Author:** Eddy Offermann \n", "This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Tuesday, May 19, 2026 at 08:26 PM." ] }, { "cell_type": "markdown", "id": "a0e3acce", "metadata": {}, "source": [ " Arguments are defined. " ] }, { "cell_type": "code", "execution_count": null, "id": "0a4a97cc", "metadata": { "collapsed": false }, "outputs": [], "source": [ "Double_t eps = 1.e-12;" ] }, { "cell_type": "code", "execution_count": null, "id": "fb6e242b", "metadata": { "collapsed": false }, "outputs": [], "source": [ "cout << \"Perform the fit y = c0 + c1 * x in four different ways\" << endl;\n", "\n", "const Int_t nrVar = 2;\n", "const Int_t nrPnts = 4;\n", "\n", "Double_t ax[] = {0.0,1.0,2.0,3.0};\n", "Double_t ay[] = {1.4,1.5,3.7,4.1};\n", "Double_t ae[] = {0.5,0.2,1.0,0.5};" ] }, { "cell_type": "markdown", "id": "faa14261", "metadata": {}, "source": [ "Make the vectors 'Use\" the data : they are not copied, the vector data\n", "pointer is just set appropriately" ] }, { "cell_type": "code", "execution_count": null, "id": "882555ff", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TVectorD x; x.Use(nrPnts,ax);\n", "TVectorD y; y.Use(nrPnts,ay);\n", "TVectorD e; e.Use(nrPnts,ae);\n", "\n", "TMatrixD A(nrPnts,nrVar);\n", "TMatrixDColumn(A,0) = 1.0;\n", "TMatrixDColumn(A,1) = x;\n", "\n", "cout << \" - 1. solve through Normal Equations\" << endl;\n", "\n", "const TVectorD c_norm = NormalEqn(A,y,e);\n", "\n", "cout << \" - 2. solve through SVD\" << endl;" ] }, { "cell_type": "markdown", "id": "3c87069b", "metadata": {}, "source": [ "numerically preferred method" ] }, { "cell_type": "markdown", "id": "2ea58130", "metadata": {}, "source": [ "first bring the weights in place" ] }, { "cell_type": "code", "execution_count": null, "id": "52f12aad", "metadata": { "collapsed": false }, "outputs": [], "source": [ "TMatrixD Aw = A;\n", "TVectorD yw = y;\n", "for (Int_t irow = 0; irow < A.GetNrows(); irow++) {\n", " TMatrixDRow(Aw,irow) *= 1/e(irow);\n", " yw(irow) /= e(irow);\n", "}\n", "\n", "TDecompSVD svd(Aw);\n", "Bool_t ok;\n", "const TVectorD c_svd = svd.Solve(yw,ok);\n", "\n", "cout << \" - 3. solve with pseudo inverse\" << endl;\n", "\n", "const TMatrixD pseudo1 = svd.Invert();\n", "TVectorD c_pseudo1 = yw;\n", "c_pseudo1 *= pseudo1;\n", "\n", "cout << \" - 4. solve with pseudo inverse, calculated brute force\" << endl;\n", "\n", "TMatrixDSym AtA(TMatrixDSym::kAtA,Aw);\n", "const TMatrixD pseudo2 = AtA.Invert() * Aw.T();\n", "TVectorD c_pseudo2 = yw;\n", "c_pseudo2 *= pseudo2;\n", "\n", "cout << \" - 5. Minuit through TGraph\" << endl;\n", "\n", "TGraphErrors *gr = new TGraphErrors(nrPnts,ax,ay,0,ae);\n", "TF1 *f1 = new TF1(\"f1\",\"pol1\",0,5);\n", "gr->Fit(\"f1\",\"Q\");\n", "TVectorD c_graph(nrVar);\n", "c_graph(0) = f1->GetParameter(0);\n", "c_graph(1) = f1->GetParameter(1);" ] }, { "cell_type": "markdown", "id": "1ff27a23", "metadata": {}, "source": [ "Check that all 4 answers are identical within a certain\n", "tolerance . The 1e-12 is somewhat arbitrary . It turns out that\n", "the TGraph fit is different by a few times 1e-13." ] }, { "cell_type": "code", "execution_count": null, "id": "c684dc1e", "metadata": { "collapsed": false }, "outputs": [], "source": [ "Bool_t same = kTRUE;\n", "same &= VerifyVectorIdentity(c_norm,c_svd,0,eps);\n", "same &= VerifyVectorIdentity(c_norm,c_pseudo1,0,eps);\n", "same &= VerifyVectorIdentity(c_norm,c_pseudo2,0,eps);\n", "same &= VerifyVectorIdentity(c_norm,c_graph,0,eps);\n", "if (same)\n", " cout << \" All solutions are the same within tolerance of \" << eps << endl;\n", "else\n", " cout << \" Some solutions differ more than the allowed tolerance of \" << eps << endl;" ] } ], "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 }