Fitting closed polynomial curves to given points in the plane

Question

Consider the equation

$$ x^2 + y^2 - R^2 = 0 \tag{1}\label{470162_1} $$

which describes a circle. For a given point (e.g. $P_1 = (0.8, 1.5)$) the parameter $R$ can be easily determined.

I am interested in the situation where I have $n$ distinct given points in the plane. I want to "fit" a closed algebraic curve given by $F_{\mathbf q}(x, y)=0$ through these points, where $F_{\mathbf q}(x, y)$ is a bivariate polynomial of order $N$ parameterized by the coefficient vector $\mathbf q$.

I suppose for $N>2$ this is hard or even impossible to solve analytically but a numerical approach would be sufficient for me.

$N$ can be choosen arbitrarily. In general $\mathbf q$ contains $(N + 2)(N+1)/2$ free parameters and the $n$ given points together pose $2n$ conditions. Thus I would suppose that $(N + 2)(N+1)/2 \geq 2n$ is a condition for solvability. But there might be additional constraints come into play to ensure that the curve is closed (may be better formulated: has at least one closed branch).

Question: Is there an established way to calculate $\mathbf q$ e.g. like reformulating the original problem (system of polynomial equations with the same structure) as a well posed optimization problem or cleverly using a Groebner Basis?

Answer 1

Alright, this is my attempt (probably not complete).

First, let's consider the problem of bivariate interpolation:

Given points $(x_1, y_1), \dots, (x_n, y_n) \in \mathbb{R}^2$ and interpolation values $z_1,\dots,z_n \in \mathbb{R}$, we wish to find a polynomial $f(x,y) \in \mathbb{R}[x, y]$ such that $f(x_t, y_t) = z_t$ for $1 \le t \le n$.

This is easily solved using bivariate Lagrange polynomials, if we first partition the points based on their $x$-coordinates. So, let $\alpha_1,\dots, \alpha_m$ with $m \le n$ be the distinct $x$-coordimates, and for every $\alpha_i$ let $\beta_{i, j}$ be the $j$-th $y$-coordinate among those that belong to $\alpha_i$. Also, let $\gamma_{i,j}$ be the desired interpolation value for the point $(\alpha_i, \beta_{i,j})$. Finally, let $\ell(i)$ be the number of $y$-coordinates for $\alpha_i$. Our polynomial can be written down explicitly like so:

$ f(x,y) = \sum_{1 \le i \le m} \Big( \prod_{i' \ne i} \frac{x - \alpha_{i'}}{\alpha_i - \alpha_{i'}} \Big) \sum_{1 \le j \le \ell(i)} \gamma_{i,j} \prod_{j' \ne j} \frac{y - \beta_{j'}}{\beta_j - \beta_{j'}} $.

Setting $\gamma_{i,j} = 0$ would create some (possibly not irreducible) curve that fits all of your points, however it might not be closed/bounded. To enforce this, we could try the following trick: Since $f(x,y)$ has degree at most $n$ in either variable, then $f(x,y) + (x^2 + y^2)^{n+1}$ should approach infinity when we go far from the origin. In particular, this new polynomial is not zero at infinity (so it's closed?). Of course, by adding $(x^2 + y^2)^{n+1}$ we broke all of the interpolation values, but this is easy to fix. First, evaluate $(x^2 + y^2)^{n+1}$ at your points, and then use the Lagrange formula to interpolate the negation of these evaluations. This way $f(x,y) + (x^2 + y^2)^{n+1}$ would evaluate to $0$ at your points, and would approach infinity as you go far from the origin.

So it looks to me that the only remaining issue is that the created curve could be reducible, which maybe you don't want? Since this was not part of the original question, I'll ignore this issue. Also, maybe it's possible to enforce that the curve is bounded using projective coordinates, but I'm not sure about this. Also also, this should work in any number of dimensions, right?

Also also also, you could construct a linear system where the unknowns are the coefficients of your polynomial. This way, you could respect a prescribed monomial support (like if you want a circle).