# Algorithm to get inverse parabola fitting

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?

-
If you can do it in MATLAB, then you must understand the algorithm -- unless you're using a canned function from the Curve Fitting Toolbox or something similar. What exactly are you having trouble with? Also, have you considered using some of MATLAB's code-generation utilities? MATLAB to C is quite robust, if you have access to MATLAB Coder (fka Real Time Workshop). –  Arkamis Dec 17 '12 at 18:09
I did it in matlab as explained, creating a linear system but making b to be 1/y, but when somebody recommended me to do this they told me that it can be unestable. What do you think @EdGorcenski? –  Ander Biguri Dec 17 '12 at 18:16
I think I would need to see your MATLAB code. 1/y, when y is a vector, is not something that makes sense. –  Arkamis Dec 17 '12 at 18:20
well what i write is pseudocode, it will probably be 1/.y . Ill post it tomorrow, I can0t acces it now @EdGorcenski –  Ander Biguri Dec 17 '12 at 18:30
Ah, I see what you are doing now. I will post a solution. –  Arkamis Dec 17 '12 at 18:31

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.

-
But how stable is this method? If y's are very small (which we kind of see in his picture) then 1/y can have relatively large errors due to finite precision and so on. –  Fixed Point Dec 17 '12 at 18:47
The issue of $y_k$ being small is taken only in consideration with the product $A^T\mathbf{y}$. Generally, it doesn't matter if $b$ is large when solving $Ax = b$, unless you get to overflow range (which is very large indeed). You can easily compute the right hand side of the least-squares solution, note that $A^T\mathbf{y}$ is a vector, and each element is the dot product of $\mathbf{y}$ and one of the columns of the Vandermonde matrix. You can apply some scaling and shifting to the elements in the vandermonde matrix to cancel out smallness of $\mathbf{y}$. –  Arkamis Dec 17 '12 at 18:59
@PrinceAli I have edited my answer. –  Arkamis Dec 17 '12 at 19:05
Ahhhhhhhh, forgot about preconditioning. Another thing, with all this, wouldn't it be worth it to just do non-linear least squares fit? –  Fixed Point Dec 17 '12 at 19:45
It all depends on your needs. –  Arkamis Dec 17 '12 at 19:53

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.

-