Finding polynomials from huge sets of points

Question

Let's say we have 1000 points, \begin{align*} (a_1 &, b_1) \\ (a_2 &, b_2) \\ & \vdots \\ (a_{1000} &, b_{1000}) \\ \end{align*} and we want to find a 999 degree polynomial that fits these points. The method I have in mind is to have a polynomial $$c_{999}x^{999}+c_{998}x^{99}+\ldots+c_1x+c_0=y$$ We would then plug in each of the 1000 points to get a huge matrix: $$\begin{bmatrix} a_{1000}^{999} & a_{1000}^{998} & \dots & a_{1000}^{0}\\ a_{999}^{999} & a_{999}^{998}\\ \vdots & \vdots& \ddots\\ a_1^{999} & \ldots & \ldots &a_1^0 \end{bmatrix} \begin{bmatrix} c_{999} \\ c_{998} \\ \vdots \\ c_0 \end{bmatrix} = \begin{bmatrix} b_{1000} \\ b_{999} \\ \vdots \\ b_1 \end{bmatrix} $$ However, this matrix is probably really hard to solve, because of the size of the numbers ($a_{1000}^{999}$ is huge), and the size of the matrix. I am doing this in Java, and I don't even know how to store numbers that big.

Is there a better way of doing this than using matrices that is more efficient?

Answer 1

The standard way of thinking about this is the Lagrange interpolation:

$$f(x)=\sum_{i=1}^{1000}b_i\prod_{\substack{j=0\\j \ne i}}^{1000} \frac{x-a_j}{a_i-a_j} .$$

I think if you're sufficiently clever you can write code for this that will be slightly more efficient than a standard matrix division but at $n=1000$ it might not even be worth it.

Of course, if you're using scipy to calculate these, there's already a function for that and if you're using a different system like MATLAB you can likely use a ready-made function like that as well, which deals with edge cases.

Answer 2

Yes! This is a standard problem in numerical analysis, and two possible solutions include using Lagrange interpolating polynomials, and Newton's method for interpolation.