I am currently working on determining the errors associated with a least-squares fit for powder X-ray diffraction data analysis. I am utilizing the numpy.linalg.lstsq function to accomplish this task. The system of equations for my case is presented below. The goal of this analysis is to determine the lattice parameters.
$$ \begin{bmatrix} h_1^2 & k_1^2 & l_1^2 \\ h_2^2 & k_2^2 & l_2^2 \\ \vdots & \vdots & \vdots \\ h_N^2 & k_N^2 & l_N^2 \end{bmatrix} \cdot \begin{pmatrix} u \\ v \\ w \end{pmatrix} = \begin{pmatrix} q_1 \\ q_2 \\ \vdots \\ q_N \end{pmatrix} $$
I have successfully used the numpy.linalg.lstsq() function to obtain results that are consistent with literature for this compound (i.e. lattice parameters). I am now focusing on estimating the uncertainty (i.e. error or variance) of the fit parameters u, v, w. For reference, when using the scipy.optimize.curve_fit() function, it returns the covariance matrix. My understanding is that the square root of the diagonal elements of the covariance matrix are the variances for each parameter. However, numpy.linalg.lstsq() does not return the covariance matrix, which makes the estimation of the uncertainty of the fit parameters more challenging.
I have examined scipy.optimize.curve_fit()'s source code. It states that the covariance matrix is equivalent to the inverse of the Hessian. Is this statement accurate? How can I verify this claim to be true?
For now, I am calculating the Hessian, $H$, of my matrix, $A$, using the following equation:
\begin{equation} H = A^{\text{T}}A \end{equation}
Therefore, my covariance matrix $P_{cov}$ ought to be:
\begin{equation} P_{cov} = \left(A^{\text{T}}A\right)^{-1} \end{equation}
This returns very believable fit uncertainties.
But why? And why are the fit uncertainties not dependent on my values of $q_1, . . ., q_N$.
I have checked out threads like: Relationship between Hessian Matrix and Covariance Matrix or Basic question about Fisher Information matrix and relationship to Hessian and standard errors but I realize that my understanding of linear algebra is not quite there yet ... I'm a chemist.
Here is my current understanding of what a covariance matrix and a Hessian matrix are:
The covariance matrix is a matrix that describes the covariance between multiple variables. It is used to represent the level of correlation and dependence between different elements of a dataset. The elements of the matrix are calculated as the covariance between each pair of variables, with the main diagonal elements representing the variance of each individual variable.
A Hessian matrix is a square matrix of second-order partial derivatives of a scalar-valued function. ... But here I'm already confused. My matrix $A$ is full of numbers. And the derivative of any number is just zero. So how can $H$ be a second derivative of $A$?
Any help in trying to understand this topic is greatly appreciated!
