Algorithm to get inverse parabola fitting

Question

I am coming from Stack overflow, where they closed my question because of being "mathemathical and nothing of programming", please, I know that in this forums usually mathemathical questions are more complex and elaborated, but I have nowhere else to go!

I am trying to fit some points to an inverse parabola of the form of $F(x)=1/(ax^2+bx+c)$.

My objective is to program a function in c++ that would take a set of 10-30 points and fit them to the inverse parabola.

I started trying to get an analytical expresion using Least squares, but I can't reach to get a result. I tried by hand (little crazy) and then I tried to solve analytically the expressions for a,b and c, but mupad doesn't give me a result (I am pretty new to Matlab's mupad so maybe i am not doing it correctly). i don't know anymore how to approach the problem.

Can I get a anallytical expresion for this specific problem? I have also seen algorithms for general least squares fitting but I don't need a so complicated algorithm, I just need it for this equation. If not, how would StackOverflow people approach the problem?

If needed I can post the equations, and the small Mupad code I've tried, but I think is unnecessary.

EDIT: some example

Im sorry the image is a little bit messy but it is the thing i need. The data is in blue (this data is particularly noisy). I need to use only the data that is between the vertical lines (a bunch of data in the left and another one in the right).

The result of the fit is in de red line.

all this has been made with matlab, but I need to make it in c++.

I'll try to post some data...

Edit 2: I actually did the fitting in Matlab as follows (not the actual code):

 create linear system Ax = b, with 
 A = [x²  x  1]
 x = [a; b; c]
 b = 1/y;

It should work, shouldn it? I can solve then using Moore-Penrose pseudoinv calculated with SVD. Isn't it?

Answer 1

If you have some data $y$ and want to fit it to a curve of the form $1/(ax^2+bx+c)$, then what we can do is fit the data $1/y$ to the curve $ax^2+bx+c$. To do this in the least squares sense, we construct the following system of equations:

$$ \begin{pmatrix} x_1^2 & x_1 & 1 \\ x_2^2 & x_2 & 1 \\ \vdots & \vdots & \vdots \\ x_n^2 & x_n & 1\end{pmatrix} \begin{pmatrix} a \\ b \\ c\end{pmatrix} = \begin{pmatrix} 1/y_1 \\ 1/y_2 \\ \vdots \\ 1/y_n \end{pmatrix} $$

The matix on the left-hand side is known as a Vandermonde matrix. Denote this matrix $A$

Since in general $n \neq 3$, to solve this we simply solve $A^TA \mathbf{p} = A^T\mathbf{y}$, where $\mathbf{p} = (a,b,c)^T$ is the parameter vector.

Thus, $$\mathbf{p} = (A^TA)^{-1}A^T\mathbf{y},$$ which may be computed in many ways (note that $A^TA$ is symmetric, so the LU factorization should be easy to compute).

EDIT: Regarding largeness of $\mathbf{y}$ (which should not matter for the data in your picture).

Note that we are solving $A^TA\mathbf{p} = A^T\mathbf{y}$. Look at $A^T\mathbf{y}$:

$$A^T \mathbf{y} = \begin{pmatrix} \sum_{i=1}^n \frac{x_i^2}{y_1} \\ \sum_{i=1}^n \frac{x_i}{y_1} \\ \sum_{i=1}^n \frac{1}{y_1} \end{pmatrix}$$

If $y_i$ is so small that you are concerned about numerical stability, solve the following problem instead:

$$10^m \hat{a} x^2 + 10^m \hat{b} x + 10^m \hat{c} = 1/y$$

for some $m > 0$. Then, you condition your Vandermonde matrix in such a way that $A^T\mathbf{y}$ becomes:

$$A^T \mathbf{y} = \begin{pmatrix} \sum_{i=1}^n \frac{10^m x_i^2}{y_1} \\ \sum_{i=1}^n \frac{10^m x_i}{y_1} \\ \sum_{i=1}^n \frac{10^m}{y_1} \end{pmatrix}$$

and the smallness is gone. Compute $\hat{a},\hat{b},\hat{c}$ in the traditional sense, then $a = 10^m\hat{a}$, and so forth.

Answer 2

As you yourself mentioned, this method can be unstable especially with the "small" values of y you are using (what we see in the picture). But this method is definitely easier and faster than nonlinear least squares fit since you want to write your own implementation. You can try this and see how sensitive this method is to your typical y-values (with typical errors). But keep in mind that if y is small, then 1/y will be large so a small change in y can lead to a huge change in 1/y so your answers can change radically.