# Multivariate polynomial regression of discrete data in L-infinity norm

**URL:** <https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369>\
**Category:** Numerics\
**Tags:** interpolations, approximation, chebychev\
**Created:** [January 30, 2025, 5:32am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369 "2025-01-30T05:32:34Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 5:32am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/1 "2025-01-30T05:32:34Z")

</div>

Hi everyone,

Does anyone know of a Julia package that implements polynomial approximations (preferably Chebyshev polynomials) given data on a multivariate function with a L-infinity norm?

Specifically, I have x, f(x) pairs and I aim to fit a polynomial that minimizes the maximum error at the data points. It needs to be a polynomial because I need its partial derivatives but it can either be an interpolation or a regression approach (I do not need the approximation to go through all data points).

[FastChebInterp.jl](https://github.com/JuliaMath/FastChebInterp.jl) does almost what I want but its interpolation command (`chebinterp()`) needs Chebychev points (which I do not have) and its regression command `chebregression()` uses a L-2 norm.

Any ideas are welcome!

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [January 30, 2025, 5:52am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/2 "2025-01-30T05:52:37Z")

</div>

Remez.jl does what you’re looking for (although it returns BigFloats in the monomial basis)

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 6:34am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/3 "2025-01-30T06:34:44Z")

</div>

Thanks! This is a great start but it seems to only accept one-dimensional functions, right?

---

<div class="post-metadata">

**Author:** ![dlfivefifty](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlfivefifty/32/1959_2.png) [@dlfivefifty](https://discourse.julialang.org/u/dlfivefifty)\
**Post date:** [January 30, 2025, 7:45am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/4 "2025-01-30T07:45:32Z")

</div>

Just wanted to point out Chebyshev can mean two different things in this context which might mean people misunderstand your question: (1) the Chebyshev polynomials `cos(k*acos(x))` and (2) Chebyshev approximation, which is the best L^∞ fit and the error equi-oscillates.

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [January 30, 2025, 9:05am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/5 "2025-01-30T09:05:33Z")

</div>

> [@CeterisPartybus](#):
>
> polynomial approximations (preferably Chebyshev polynomials) given data on a multivariate function with a L-infinity norm?

What dimensionality do you have in mind? This is a hard problem I think. Can you point to a package for some other language that does what you want?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 30, 2025, 2:06pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/6 "2025-01-30T14:06:08Z")

</div>

> [@CeterisPartybus](#):
>
> Does anyone know of a Julia package that implements polynomial approximations (preferably Chebyshev polynomials) given data on a multivariate function with a L-infinity norm? […] Specifically, I have x, f(x) pairs and I aim to fit a polynomial that minimizes the maximum error at the data points.

If I understand correctly what you are doing, what you want is a _lot_ easier than the question the other posters seem to think they are answering.

Packages like Remez.jl are solving the _continuous_ minimax problem in which one wants to minimize \Vert p(x) - f(x) \Vert\_\infty over a _continuous interval_ x \in [a,b], whereas it sounds like you only want to minimize the L\_\infty error over a **discrete set** of data points.

In particular, suppose you have a set of m data points \{ x\_k, y\_k \}\_{k = 1,\ldots, m} and you want to find a degree-n polynomial p(x) that minimizes the L\_\infty error _at those points_:

\min\_p \left[\max\_k |p(x\_k) - y\_k| \right]

This problem is straightforward to solve because it is equivalent to linear programming (LP). In particular, the “epigraph” reformulation into an LP works by adding an additional real parameter t \in \mathbb{R}

\min\_{p \in \text{polynomials}, t\in \mathbb{R}} t \\ \text{s.t. }\qquad t \ge p(x\_k) - y\_k, \; t \ge y\_k - p(x\_k) \qquad \text{for } k = 1,\ldots,m

(This is an LP because p(x\_k) is a linear function of the polynomial coefficients in whatever basis. That’s true regardless of the dimensionality of x.)

So, you can use any LP solver in one of Julia’s many convex-optimization packages to solve this.

PS. People may have been confused by original the title of your post, “Function Approximation in L-infinity Norm”, whereas you are really approximating _data_ (i.e. doing regression, i.e. discrete domain), not functions on a continuous domain. I’ve edited the thread title accordingly.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 30, 2025, 2:30pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/7 "2025-01-30T14:30:12Z")

</div>

More concisely, if the m \times (n+1)^d matrix V is the [Vandermonde matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix) of your polynomial basis (of degree n along d dimensions) for the points x\_k, then what you want to do is to solve for the coefficients c using:

\min\_{c \in \mathbb{R}^{(n+1)^d}} \Vert Vc - y \Vert\_\infty

which is a convex optimization problem that you should be able to give [almost directly](https://jump.dev/JuMP.jl/dev/tutorials/linear/tips_and_tricks/#Infinity-norm) to JuMP.jl or similar, something like:

```julia
using JuMP, MathOptInterface
function solve_for_coefs(V, y)
    N = size(V,2)
    @variable(model, coefs[1:N])
    @variable(model, t)
    @constraint(model, [t; V*coefs - y] in MathOptInterface.NormInfinityCone(1 + length(y)))
    @objective(model, Min, t)
    optimize!(model)
    return value.(coefs)
end

```

In the Chebyshev-polynomial basis, V can be computed by the function `FastChebInterp.chebvandermonde(x, lb, ub, order)` where `x` is a length-m array of d-dimensional `SVector`s `x[k]`, `lb` and `ub` are `SVector`s of the lower and upper bounds you want to use for your Chebyshev polynomials (a box containing the `x[k]`), and `order = (n,n,....)` is a tuple of the polynomial degrees along each dimension. From this, you can construct a `ChebPoly` object `p` from the coefficients, so that you can then evaluate `p(x)` (and its derivatives) at arbitrary points `x`. Something like:

```julia
V = FastChebInterp.chebvandermonde(x, lb, ub, order)
coefs = solve_for_coefs(V, y)
p = FastChebInterp.ChebPoly(reshape(coefs, order .+ 1), lb, ub)

```

(I haven’t tested this code, but it should give the general idea.)

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [January 30, 2025, 2:48pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/8 "2025-01-30T14:48:07Z")

</div>

There are some alternative, possibly more efficient, LP formulations in my package FindMinimaxPolynomial.jl. The package itself currently only supports univariate approximations, but the LP formulations could be useful as an inspiration.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 30, 2025, 3:41pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/9 "2025-01-30T15:41:21Z")

</div>

> [@nsajko](#):
>
> There are some alternative, possibly more efficient, LP formulations in my package [FindMinimaxPolynomial.jl](https://juliaregistries.github.io/General/packages/redirect_to_repo/FindMinimaxPolynomial). The package itself currently only supports univariate approximations,

It looks like that package is also for the continuous minimax problem, not the discrete one? I guess somewhere inside the code it is solving a sequence of discrete minimax problems at samples?

---

<div class="post-metadata">

**Author:** ![nsajko](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nsajko/32/221187_2.png) [@nsajko](https://discourse.julialang.org/u/nsajko)\
**Post date:** [January 30, 2025, 4:32pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/10 "2025-01-30T16:32:13Z")

</div>

Yes. The [`PolynomialPassingThroughIntervals`](https://gitlab.com/nsajko/FindMinimaxPolynomial.jl/-/blob/main/src/PolynomialPassingThroughIntervals.jl) module should have an implementation for the discrete case, using which the continuous case is then attacked. This is an iterative Remez-inspired process, calling a MathOptInterface.jl LP solver multiple times.

Regarding the linear programming formulations: I took a short look at the code now, and the only smart LP formulation is what I termed `ReplacedVariablesOption{:determinants}` in the code, mostly implemented in the [`OverdeterminedSystemConsistencyConstraints`](https://gitlab.com/nsajko/FindMinimaxPolynomial.jl/-/blob/main/src/OverdeterminedSystemConsistencyConstraints.jl) module. I forgot the details (although I have something written down on paper somewhere), but it seems the idea is to take the Vandermonde matrix, which might be very tall and thin, and transform it into an equivalent system of constraints with much less rows.

BTW this was/is just a toy hobby project of mine, in case it’s not clear.

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 6:43pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/11 "2025-01-30T18:43:56Z")

</div>

Thanks. I meant an approximation by a sum of (Chebyshev) Polynomials using the L-infinity norm as a goodness of fit measure.

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 8:57pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/12 "2025-01-30T20:57:22Z")

</div>

It should not get to more than 100 dimension. I did not find a package in another language but I also did not search around too much outside of Julia.

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 9:15pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/13 "2025-01-30T21:15:57Z")

</div>

Thanks for the answer and for adjusting the title. This is for a large part what I was looking for.

Just to clarify: I aim to approximate a function on a continuous domain [a,b]. Or in other words, I have a set of points x and the function value f(x) only at those points. I cannot get more evaluations at other points. However, I do care about the quality of the approximation on the whole interval (e.g. to do simulations on other data). I understand that there is no way to test the goodness of fit on other points than those I have data on but that is where I believe approximation theory helps me.

Now, what I really care about is the partial derivative of the function. Do you have an idea on how I could get an approximation of f(x) that also best approximates \frac{\partial f(x)}{\partial x} while only having data on x and f(x)? And if so, how an approximation error in f(x) affects the error in the approximation of \frac{\partial f(x)}{\partial x}?

Thanks! This helps already a lot!

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 9:17pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/14 "2025-01-30T21:17:48Z")

</div>

Awesome, thanks! I will check your package out. Maybe, I can use parts of it!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 30, 2025, 10:09pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/15 "2025-01-30T22:09:55Z")

</div>

> [@CeterisPartybus](#):
>
> Or in other words, I have a set of points x and the function value f(x) only at those points. I cannot get more evaluations at other points.

Where are the points coming from? Experiments? Do they contain measurement noise?

---

<div class="post-metadata">

**Author:** ![CeterisPartybus](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ceterispartybus/32/46868_2.png) [@CeterisPartybus](https://discourse.julialang.org/u/CeterisPartybus)\
**Post date:** [January 30, 2025, 11:52pm UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/16 "2025-01-30T23:52:02Z")

</div>

I am developing a method that should be applicable in different scenarios. So, the data could come from different sources (reports, interview, etc). It is likely that both the x and f(x) contain measurement noise. But if this complicates the approximation substantially, assuming no measurement noise would be okay to start with.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 31, 2025, 1:10am UTC](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369/17 "2025-01-31T01:10:55Z")

</div>

> [@CeterisPartybus](#):
>
> It is likely that both the x and f(x) contain measurement noise. But if this complicates the approximation substantially, assuming no measurement noise would be okay to start with.

My understanding is that L\_\infty regression is only better than L\_2 for certain types of noise? e.g. uniform noise within bounded support (see e.g. [Yi and Neykov, 2021](https://arxiv.org/abs/2108.07630)), as opposed to Gaussian noise that might be more common in experimental data? What led you to L\_\infty in your case?

Anyway, I’m not aware of any approximation methods that specifically target the derivative (as opposed to trying to approximate a function and getting the derivative as a side-effect), though maybe someone else on here knows of something.

Instead of doing a global polynomial fit, you might also consider some kind of local smoothing/interpolation, e.g. a smoothing spline. Depends on how much you know about your function, I guess.