ROOT   Reference Guide
solveLinear.C File Reference

## Detailed Description

This macro shows several ways to perform a linear least-squares analysis .

To keep things simple we fit a straight line to 4 data points The first 4 methods use the linear algebra package to find x such that min $$(A x - b)^T (A x - b)$$ where A and b are calculated with the data points and the functional expression :

1. Normal equations: Expanding the expression $$(A x - b)^T (A x - b)$$ and taking the derivative wrt x leads to the "Normal Equations": $$A^T A x = A^T b$$ where $$A^T A$$ is a positive definite matrix. Therefore, a Cholesky decomposition scheme can be used to calculate its inverse . This leads to the solution $$x = (A^T A)^-1 A^T b$$ . All this is done in routine NormalEqn . We made it a bit more complicated by giving the data weights . Numerically this is not the best way to proceed because effectively the condition number of $$A^T A$$ is twice as large as that of A, making inversion more difficult
2. SVD One can show that solving $$A x = b$$ for x with A of size $$(m x n)$$ and $$m > n$$ through a Singular Value Decomposition is equivalent to minimizing $$(A x - b)^T (A x - b)$$ Numerically , this is the most stable method of all 5
3. Pseudo Inverse Here we calculate the generalized matrix inverse ("pseudo inverse") by solving $$A X = Unit$$ for matrix $$X$$ through an SVD . The formal expression for is $$X = (A^T A)^-1 A^T$$ . Then we multiply it by $$b$$ . Numerically, not as good as 2 and not as fast . In general it is not a good idea to solve a set of linear equations with a matrix inversion .
4. Pseudo Inverse , brute force The pseudo inverse is calculated brute force through a series of matrix manipulations . It shows nicely some operations in the matrix package, but is otherwise a big "no no" .
5. Least-squares analysis with Minuit An objective function L is minimized by Minuit, where $$L = sum_i { (y - c_0 -c_1 * x / e)^2 }$$ Minuit will calculate numerically the derivative of L wrt c_0 and c_1 . It has not been told that these derivatives are linear in the parameters c_0 and c_1 . For ill-conditioned linear problems it is better to use the fact it is a linear fit as in 2 .

Another interesting thing is the way we assign data to the vectors and matrices through adoption . This allows data assignment without physically moving bytes around .

#### USAGE

This macro can be executed via CINT or via ACLIC

• via the interpretor, do
root > .x solveLinear.C
• via ACLIC
root > .x solveLinear.C+
Perform the fit y = c0 + c1 * x in four different ways
- 1. solve through Normal Equations
- 2. solve through SVD
- 3. solve with pseudo inverse
- 4. solve with pseudo inverse, calculated brute force
- 5. Minuit through TGraph
All solutions are the same within tolerance of 1e-12
#include "Riostream.h"
#include "TMatrixD.h"
#include "TVectorD.h"
#include "TGraphErrors.h"
#include "TDecompChol.h"
#include "TDecompSVD.h"
#include "TF1.h"
void solveLinear(Double_t eps = 1.e-12)
{
cout << "Perform the fit y = c0 + c1 * x in four different ways" << endl;
const Int_t nrVar = 2;
const Int_t nrPnts = 4;
Double_t ax[] = {0.0,1.0,2.0,3.0};
Double_t ay[] = {1.4,1.5,3.7,4.1};
Double_t ae[] = {0.5,0.2,1.0,0.5};
// Make the vectors 'Use" the data : they are not copied, the vector data
// pointer is just set appropriately
TVectorD x; x.Use(nrPnts,ax);
TVectorD y; y.Use(nrPnts,ay);
TVectorD e; e.Use(nrPnts,ae);
TMatrixD A(nrPnts,nrVar);
TMatrixDColumn(A,0) = 1.0;
cout << " - 1. solve through Normal Equations" << endl;
const TVectorD c_norm = NormalEqn(A,y,e);
cout << " - 2. solve through SVD" << endl;
// numerically preferred method
// first bring the weights in place
TMatrixD Aw = A;
TVectorD yw = y;
for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
TMatrixDRow(Aw,irow) *= 1/e(irow);
yw(irow) /= e(irow);
}
TDecompSVD svd(Aw);
Bool_t ok;
const TVectorD c_svd = svd.Solve(yw,ok);
cout << " - 3. solve with pseudo inverse" << endl;
const TMatrixD pseudo1 = svd.Invert();
TVectorD c_pseudo1 = yw;
c_pseudo1 *= pseudo1;
cout << " - 4. solve with pseudo inverse, calculated brute force" << endl;
const TMatrixD pseudo2 = AtA.Invert() * Aw.T();
TVectorD c_pseudo2 = yw;
c_pseudo2 *= pseudo2;
cout << " - 5. Minuit through TGraph" << endl;
TGraphErrors *gr = new TGraphErrors(nrPnts,ax,ay,0,ae);
TF1 *f1 = new TF1("f1","pol1",0,5);
gr->Fit("f1","Q");
TVectorD c_graph(nrVar);
c_graph(0) = f1->GetParameter(0);
c_graph(1) = f1->GetParameter(1);
// Check that all 4 answers are identical within a certain
// tolerance . The 1e-12 is somewhat arbitrary . It turns out that
// the TGraph fit is different by a few times 1e-13.
Bool_t same = kTRUE;
same &= VerifyVectorIdentity(c_norm,c_svd,0,eps);
same &= VerifyVectorIdentity(c_norm,c_pseudo1,0,eps);
same &= VerifyVectorIdentity(c_norm,c_pseudo2,0,eps);
same &= VerifyVectorIdentity(c_norm,c_graph,0,eps);
if (same)
cout << " All solutions are the same within tolerance of " << eps << endl;
else
cout << " Some solutions differ more than the allowed tolerance of " << eps << endl;
}

Definition in file solveLinear.C.

kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:100
e
#define e(i)
Definition: RSha256.hxx:103
TVectorD.h
TDecompSVD
Single Value Decomposition class.
Definition: TDecompSVD.h:24
TMatrixTSym< Double_t >::kAtA
@ kAtA
Definition: TMatrixTSym.h:51
Int_t
int Int_t
Definition: RtypesCore.h:45
x
Double_t x[n]
Definition: legend1.C:17
TMatrixTSym< Double_t >
ROOT::Math::Cephes::A
static double A[]
Definition: SpecFuncCephes.cxx:170
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Definition: TSystem.cxx:1854
TDecompSVD.h
TMatrixT< Double_t >
Bool_t
bool Bool_t
Definition: RtypesCore.h:63
TDecompChol.h
NormalEqn
TVectorD NormalEqn(const TMatrixD &A, const TVectorD &b)
Solve min {(A .
Definition: TDecompChol.cxx:404
TMatrixDRow
TMatrixTRow< Double_t > TMatrixDRow
Definition: TMatrixDUtilsfwd.h:62
TMatrixT::Invert
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1397
gr
TGraphErrors * gr
Definition: legend1.C:25
y
Double_t y[n]
Definition: legend1.C:17
TMatrixDColumn
TMatrixTColumn< Double_t > TMatrixDColumn
Definition: TMatrixDUtilsfwd.h:63
gSystem
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
TVectorT< Double_t >
f1
TF1 * f1
Definition: legend1.C:11
Double_t
double Double_t
Definition: RtypesCore.h:59
TGraphErrors
A TGraphErrors is a TGraph with error bars.
Definition: TGraphErrors.h:26
TF1.h
TGraph::Fit
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Axis_t xmin=0, Axis_t xmax=0)
Fit this graph with function with name fname.
Definition: TGraph.cxx:1073
Bool_t VerifyVectorIdentity(const TVectorT< Element > &m1, const TVectorT< Element > &m2, Int_t verbose, Element maxDevAllow)
Verify that elements of the two vectors are equal within maxDevAllow .
Definition: TVectorT.cxx:2297
TGraphErrors.h
TF1
1-Dim function class
Definition: TF1.h:213
TMatrixD.h
TF1::GetParameter
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:512
Riostream.h
TMatrixT::T
TMatrixT< Element > & T()
Definition: TMatrixT.h:150