Fitting a low-degree polynomial to a function on a finite 1d grid - a combinatorial problem?

Question

I need to fit a low-degree polynomial $p$ (with $\text{deg}(p) \leq k$) to a function $f$ evaluated on the grid $\{0, 1, ... n-1\}$, so as to minimize the $L_\infty$ norm, i.e. minimize $\text{max}_{0 \leq i < n}{\left| f(i) - p(i) \right|}$.

Now, I know that for $L_2$, this can be solved using a system of linear equations (for getting a zero vector as the derivative of the error) - that is the "Least Squares" method. But - the $L_\infty$ error is not a derivable function. I also know that for a continuous interval rather than a finite grid, there's an iterative approximation algorithm - the Remez Algorithm (although I'm not entirely clear about the specifics - I just noticed it.) And I also had a vague hunch which I didn't manage to act on that Chebychev polynomials, each of which being the optimal solution to the problem with $f = 0$ and a fixed leading coefficient 1, might be relevant here somehow.

Anyway, what I'd like to know is:

• Is this problem is essentially combinatorial or essentially numero-analytic?
• Can it be solved in PTIME (assuming operations on real values are O(1))?
• If so, what's the numerical stability features of such a solution (since my n is large-ish)?
• If not, does it admit an efficient (PTIME) approximation scheme?
• If not, how do people deal with it in practice?

Answer 1

The problem can be solved in polynomial time using linear programming.

Write $p(x) = c_k x^k + \dots + c_1 x + c_0$, and think of $c_0,\dots,c_k$ as unknowns. Also introduce an unknown $d$, which will represent the $L_\infty$ norm. We will minimize $d$, subject to the following inequality:

$$-d \le f(i) - (c_k i^k + \dots + c_0) \le d.$$

Notice that, since $i$ and $f(i)$ are constants, this is a linear inequality in the unknowns. We obtain one inequality per value of $i$. Now the minimum achievable value of $d$ corresponds to the polynomial $p$ that minimize the $L_\infty$ norm, so you can solve this using any algorithm for linear programming. Since linear programming is in P, your problem is in P, too.

I expect that this algorithm will be efficient in practice, as linear programming over $k+2$ unknowns and $n$ inequalities should be pretty efficient for reasonable values of $k,n$.

I don't know what the numerical stability issues might be.

---

Another solution in practice might be to find a solution that minimizes the $L_p$ norm. You could start by minimizing the $L_2$ norm, using any standard method (e.g., least squares). Then you could start from that solution and find a solution that minimizes the $L_p$ norm, for some larger value of $p$, say using gradient descent. After increasing $p$ a few times I expect the resulting solution might often turn out to be a reasonable approximation, even if not optimal.

Answer 2

Here's my current idea. It's oriented towards a faster implementation - i.e., for me, $\Omega(n^2)$ is much too much. It also assumes a lower threshold value for distance (which is relevant for me, since the functions' domain in practice is values stored in O(1) space, e.g. 4-byte integers or floating point values).

Well, the idea is only a sketch, but in a nutshell:

• Look for feasible solutions instead of trying to optimize.
• Perform a binary search for the minimum feasible $L_\infty$ value - bewteen an initial value of some simplistic solution and a threshold value under which you don't care (or under which you can assume optimality).
• The feasibility constraint for for $L_\infty$-error $\epsilon$, at point i, is $$\left|\sum^k_{j = 0} a_j i^j - f(i)\right| \leq \epsilon$$ but this can easily be converted into two linear constraints: $$\sum^k_{j = 0} a_j i^j \leq f(i) + \epsilon \quad \text{ and } \quad \sum^k_{j = 0} a_j i^j \geq f(i) - \epsilon$$ _{(The variables here are the $a_j$'s, so it's not a problem to use powers of $i$.)}

... and this should be a reduction of the problem to checking the feasibility of a linear program, for a number of times which does not depend on $n$.

If that's the best I can do - although I have no reason to assume that is the case - then the answers to my questions would be:

• Combinatorial or numero-analytical? - Somewhere in the middle I guess you could say.
• Can it be solved in PTIME - ???
• Admits an efficient (PTIME) approximation scheme.
• How people deadl with it in practice - probably with this kind of approximation, perhaps with wiser guessing of initial solutions. For example, one might start with an $L_2$-optimal solution on a fixed-size sample (perhaps using over-defined least squares with a larger-than-k+1 sample?) to get the initial solution, reducing the number of iterations if $f$ is well-behaved.