Hi Estelle,

I see your “brute force” solution and want to draw your attention

to an efficient implementation in the matrix package.

Have a look at the non-member function

https://root.cern/doc/master/TDecompChol_8cxx.html#a6f8408a544cb8b5ac2a24db5840ab8ea

```
////////////////////////////////////////////////////////////////////////////////
/// Solve min {(A . x - b)^T W (A . x - b)} for vector x where
/// A : (m x n) matrix, m >= n
/// b : (m) vector
/// x : (n) vector
/// W : (m x m) weight matrix with W(i,j) = 1/std(i)^2 for i == j
/// = 0 for i != j
TVectorD NormalEqn(const TMatrixD &A,const TVectorD &b,const TVectorD &std)
{
if (!AreCompatible(b,std)) {
::Error("NormalEqn","vectors b and std are not compatible");
return TVectorD();
}
TMatrixD mAw = A;
TVectorD mBw = b;
for (Int_t irow = 0; irow < A.GetNrows(); irow++) {
TMatrixDRow(mAw,irow) *= 1/std(irow);
mBw(irow) /= std(irow);
}
TDecompChol ch(TMatrixDSym(TMatrixDSym::kAtA,mAw));
Bool_t ok;
return ch.Solve(TMatrixD(TMatrixD::kTransposed,mAw)*mBw,ok);
}
```

Here a diagonal matrix is assumed but trivially replaced by a matrix.

The tutorial solveLinear.C might also give some inspiration:

https://root.cern/doc/master/solveLinear_8C.html

-Eddy