# 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:** 1
**Showing post:** 7

<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.)

---

_[View the full topic](https://discourse.julialang.org/t/multivariate-polynomial-regression-of-discrete-data-in-l-infinity-norm/125369)._