Fitting data to a 3rd degree polynomial

Question

I'm currently writing a C++ program where I have vectors of independent and dependent data that I would like to fit to a cubic function. However, I'm having trouble generating a polynomial that can fit my data.

Part of the problem is that I can't use various numerical packages, such as GSL (long story); it's possible that it might even be overkill for my case. I don't need a very generalized solution for least squares fitting. I specifically want to fit my data to a cubic function. I do have access to Sony's vector library, which supports 4x4 matrices and can calculate their inverses, among other things.

While prototyping this in Scilab, I used a function like:

function p = polyfit(x, y, n)
    m = length(x);
    aa = zeros(m, n+1)
    aa(:,1) = ones(m,1)
    for k = 2:n+1
        aa(:,k) = x.^(k-1)
    end
    p = aa\y
endfunction

Unfortunately, this doesn't map well to my current environment. The above example needs to support a matrix of M x N+1 dimensions. In my case, that's M x 4, where M depends on how much sample data that I have. There's also the problem of left division. I would need a matrix library that supported the inverse of matrices of arbitrary dimensions.

Is there an algorithm for least squares where I can avoid having to calculate aa\y, or at least limit it to a 4x4 matrix? I suppose that I'm trying to rewrite the above algorithm into a simpler case that works for fitting to a cubic polynomial. I'm not looking for a code solution, but if someone can point me in the right direction, I'd appreciate it.

Answer 1

Here is the page I am working from, although that page itself doesn't address your question directly. The summary of my answer would be:

If you can't work with Nx4 matrices directly, then do those matrix computations "manually" until you have the problem down to something that has only 4x4 or smaller matrices. In this answer I'll outline how to do the specific matrix computations you need "manually."

--

Let's suppose you have a bunch of data points (x1,y1)...(xn,yn) and you are looking for the cubic equation y = ax^3 + bx^2 + cx + d that fits those points best.

Then following the link above, you'd write this equation:

I'll write A, x, and B for those matrices. Then following my link above, you'd like to multiply by the transpose of A, which will give you the 4x4 matrix A^T*A that you can invert. In equations, the following is the plan:

A * x = B .................... [what we started with]

(A^T * A) * x = A^T * B ..... [multiply by A^T]

x = (A^T * A)^-1 * A^T * B ... [multiply by the inverse of A^T * A]

You said you are happy with inverting 4x4 matrices, so if we can code a way to get at these matrices without actually using matrix objects, we should be okay.

So, here is a method, although it might be a little bit too much like making your own matrix library for your taste. :)

• Write an explicit equation for each of the 16 entries of the 4x4 matrix. The (i,j)th entry (I'm starting with (0,0)) is given by x₁ⁱ * x₁^j + x₂ⁱ * x₂^j + ... + x_Nⁱ * x_N^j.

• Invert that 4x4 matrix using your matrix library. That is (A^T * A)^-1.

• Now all we need is A^T * B, which is a 4x1 matrix. The ith entry of it is given by x₁ⁱ * y₁ + x₂ⁱ * y₂ + ... + x_Nⁱ * y_N.

• Multiply our hand-created 4x4 matrix (A^T * A)^-1 by our hand-created 4x1 matrix A^T * B to get the 4x1 matrix of least-squares coefficients for your cubic.

Good luck!