Linear Fitting
Over the last three years I have slowly been writing a linear fitting algorithm for hyperspy. It has turned out to be a much, much bigger project than I imagined. The following text documents some of the hurdles I've gone through, partly with regards to hyperspy, partly with regards to linear fitting in general.
The basic idea is to solve the linear algebra problem Ax=b, where A is our component data, b is the data we are fitting against, and x is the fit coefficient (or "how much do we need to scale our components by in order to fit the data").
For the purposes of highlighting my point, I'm only considering the Expression component in this discussion. A model will consist of several components. These are all added together in a linear combination - that is, you might have a model: background + gaussian1 + gaussian2. However, it turns out that the various components themselves can be a linear combination of "pseudo"-components. The easiest example is a polynomial - a*x**2 + b*x + c. Writing the code would be easier if the parts were separated into three separate components a*x**2, b*x and c, but that is not reality, and would be rather unfriendly to the user. In the code I handle this by breaking up the components into pseudocomponents. However, I only handle this for Expression components. Other "custom", non-Expression components raise an error.
Most people are used to using linear regression to fit a straight line to data, and so, fitting a gaussian to raw data might seem a bit unsettling. It is easiest to visualise linear regression as the act of simply scaling the component along the y-direction. Considering the expression a*x**2 without the rest of the Polynomial, one can then see that by multiplying the expression by a number n, it doesn't move left or right, nor stretch the expression horizontally. This process is applied to all parts of the Polynomial, yielding n1, n2 and n3.
Hyperspy has the ability to fix parameters. In the code, I use Sympy to find all pseudocomponents that are independent of any free parameters, sum them all, and then subtract this value from the signal in the pixel.
One might notice that the main fitting algorithm (when the model is not fitted with bounded=True) actually doesn't use a scipy or numpy fitting function. That's because numpy and scipy's algorithms are actually slower and potentially more memory intensive than the simple operations we need to do.
We multiply the component array A by its transpose, take the inverse, multiply that again by the component array. Finally we dot the "raw" signal data with the transpose of that. The reference to matmul is basically just elementwise multiplication, but it supports the case where the matrices are just a stack of matrices, which is the case for multiple pixels.
Hyperspy allows twinning of parameters. This is particularly applicable for EDS data, where the Gaussian of a weak peak might have a quantum-mechanical yield dependency as a function of a stronger peak (aka, the area of G2 is always 0.6*G1). Hyperspy also allows chaining of this, aka G3 = 0.5*G2. This means that we need to calculate G3 in terms of G1. Even though these components are not free, they will change as G1 changes, and we do want to take them into account when fitting data, in order to deal with peak overlap. Therefore, the pseudocomponent representing G1 must be the linear combination of the computed function for G1, G2 and G3.