# Is there any julia package for multivariate polynomial regression?

**URL:** <https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785>\
**Category:** Machine Learning\
**Tags:** question, package\
**Created:** [January 14, 2024, 6:02am UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785 "2024-01-14T06:02:49Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![WuSiren](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wusiren/32/42529_2.png) [@WuSiren](https://discourse.julialang.org/u/WuSiren)\
**Post date:** [January 14, 2024, 6:02am UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/1 "2024-01-14T06:02:49Z")

</div>

For example, I want to build a multivariate polynomial model for data \{x\_{1i},x\_{2i},\cdots,x\_{mi},y\_i\}\_{i=1}^N with given degree n.

Thanks for your attention!

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 14, 2024, 1:17pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/2 "2024-01-14T13:17:38Z")

</div>

[GLM.jl](https://juliastats.org/GLM.jl/stable/) and it’s formulas should get you quite far:

```julia
using GLM

vars = [:x1, :x2, :x3, :x4]
n = 2
lhs = filter(t -> StatsModels.degree(t) <= n, reduce(*, Term.(vars)))
rhs = Term(:y)
poly = FormulaTerm(rhs, lhs)

# Fit your model
glm(poly, data, Normal())

```

Please also see the [StatsModels.jl](https://juliastats.org/StatsModels.jl/stable/) documentation if the construction of the formula is unclear.

---

<div class="post-metadata">

**Author:** ![WuSiren](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wusiren/32/42529_2.png) [@WuSiren](https://discourse.julialang.org/u/WuSiren)\
**Post date:** [January 14, 2024, 2:16pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/4 "2024-01-14T14:16:45Z")

</div>

```julia
vars = [:x1, :x2, :x3, :x4]

```

In your code the vars are defined one by one. How can I do this if I want to define a larger number m variables?

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 14, 2024, 2:32pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/5 "2024-01-14T14:32:20Z")

</div>

You could either get them from your data directly

```julia
using DataFrames

df = DataFrame(x1 = rand(10), x2 = rand(10), x3 = rand(10), x4 = rand(10), y = rand(10))
vars = setdiff(Symbol.(names(df)), (:y, )) # Don't forget to remove target column

```

or just create them programmatically

```julia
vars = [Symbol("x" * string(i)) for i in 1:10]

```

PS: If your number of variables is large, you probably need to adapt the above approach to only create the terms up to degree n in the first place, i.e., currently it creates all interactions up to order m (of which there are 2^{m}-1) and then filters out the ones with degree \> n. This will not work if m is about ten or larger – on the other hand, polynomial regression does not make sense in that case as well.

---

<div class="post-metadata">

**Author:** ![WuSiren](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wusiren/32/42529_2.png) [@WuSiren](https://discourse.julialang.org/u/WuSiren)\
**Post date:** [January 14, 2024, 3:07pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/6 "2024-01-14T15:07:13Z")

</div>

Thank you very much! 🤝 🤝

Actually, I implemented the polynomial regression function I expected in a more complex way, and was preparing to work on the polynomial regression problem for m=25, when I saw you say

> [@bertschi](#):
>
> This will not work if m is about ten or larger – on the other hand, polynomial regression does not make sense in that case as well.

🙃

By the way, when you say _polynomial regression does not make sense in that case as well_, do you mean that the accuracy is not reliable or that the computational complexity is too great?

---

<div class="post-metadata">

**Author:** ![Ajay\_Shah](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ajay_shah/32/10520_2.png) [@Ajay\_Shah](https://discourse.julialang.org/u/Ajay_Shah)\
**Post date:** [January 14, 2024, 3:22pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/7 "2024-01-14T15:22:51Z")

</div>

There is too much collinearity between the regressors, e.g. between x^6 and x^7. You need something like an orthogonal polynomial basis system to do this. But even this is generally not worth doing. Higher dimensional polynomial regression is seldom useful.

---

<div class="post-metadata">

**Author:** ![lxvm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lxvm/32/50010_2.png) [@lxvm](https://discourse.julialang.org/u/lxvm)\
**Post date:** [January 14, 2024, 4:05pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/8 "2024-01-14T16:05:36Z")

</div>

`chebregression` from [FastChebInterp.jl](https://github.com/JuliaMath/FastChebInterp.jl) may also be worth trying.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [January 14, 2024, 5:22pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/9 "2024-01-14T17:22:47Z")

</div>

> [@bertschi](#):
>
> `lhs = filter(t -> StatsModels.degree(t) <= n, reduce(*, Term.(vars)))`

can be

```julia
using Combinatorics
lhs = sum(sum(map(splat(&), subsets(Term.(vars),d))) for d in 1:n)

```

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 14, 2024, 6:08pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/10 "2024-01-14T18:08:15Z")

</div>

Thanks, that’s a very neat one-liner creating only the desired interactions (instead of filtering them out of a much larger set as my code did). @WuSiren This code solves the performance problem I mentioned and works also when m is large, but n ≪ m.

---

<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 14, 2024, 6:17pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/11 "2024-01-14T18:17:42Z")

</div>

> [@lxvm](#):
>
> `chebregression` from [FastChebInterp.jl](https://github.com/JuliaMath/FastChebInterp.jl) may also be worth trying.

> [@Ajay\_Shah](#):
>
> There is too much collinearity between the regressors, e.g. between x^6 and x^7. You need something like an orthogonal polynomial basis system to do this. But even this is generally not worth doing. Higher dimensional polynomial regression is seldom useful

FastChebInterp does exactly this — it does multi-variate polynomial regression in the Chebyshev-polynomial basis, so it avoids the ill-conditioning that arises for the monomial x^n basis.

It’s very useful for smooth functions in a few dimensions, but in more than 4–5 dimensions you run into the [curse of dimensionality](https://en.wikipedia.org/wiki/Curse_of_dimensionality) quickly. A variant that sometimes works better in higher dimensions is to do something based on sparse grids, and there are various packages in Julia for that. You can also try Surrogates.jl for various other kinds of regression bases, e.g. radial basis functions.

---

<div class="post-metadata">

**Author:** ![Mattriks](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mattriks/32/351_2.png) [@Mattriks](https://discourse.julialang.org/u/Mattriks)\
**Post date:** [January 14, 2024, 6:38pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/12 "2024-01-14T18:38:09Z")

</div>

[CoupledFields.jl](https://github.com/Mattriks/CoupledFields.jl) can also do high dimensional polynomial regression. It’s based on [kernel ridge regression](https://discourse.julialang.org/t/smooth-delta-functions-to-compute-jacobian-and-hessian-from-samples/45282/10), and finds the directions of the active subspace via a kernel matrix. Both X and Y can be multivariable.

---

<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 14, 2024, 7:02pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/13 "2024-01-14T19:02:37Z")

</div>

> [@WuSiren](#):
>
> preparing to work on the polynomial regression problem for m=25

In 25 dimensions, you really need to think more carefully about what you are doing.

Even with just linear (degree-1 in each variable) multivariate models, this has 2^{25} (≈ 34 million) parameters. Do you have hundreds of millions of data points? For degree 2 (quadratic in each variable) models you will have 3^{25} ≈ 850 billion parameters, so you will need trillions (or more) data points. (Even if you do have that much data, solving for the fit parameters is then a supercomputer-scale system of linear equations.)

Unless you mean polynomial models where _each term_ has at most `n` variables? That is much more manageable. e.g. if each term has at most 2 variables (and at most one of each) then you only have \begin{pmatrix} 25 \\ 0 \end{pmatrix} + \begin{pmatrix} 25 \\ 1 \end{pmatrix} + \begin{pmatrix} 25 \\ 2 \end{pmatrix}= 326 unknowns. (If you allow the same variable to appear more than once in each term, you get a slightly larger number 351: n-choose-k becomes [n-choose-k with replacement](https://www.johndcook.com/blog/select_with_replacement/).) But of course, by doing so you have greatly restricted your polynomials (in fact, it’s somewhat similar in spirit to sparse grids: you’ve restricted the mixed derivatives).

There are lots of high-dimensional regression models (from sparse grids to kernel methods to neural nets) that try to avoid the curse of dimensionality. But they all trade off some generality, and go beyond simple polynomial regression.

---

<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 14, 2024, 7:19pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/14 "2024-01-14T19:19:46Z")

</div>

> [@bertschi](#):
>
> `lhs = filter(t -> StatsModels.degree(t) <= n, reduce(*, Term.(vars)))`

Maybe I’m misunderstanding the code here, but doesn’t this produce a single term? `reduce(*, Term.(vars))` is equivalent to `prod(Term, vars)`, no? If you wanted to go this route, wouldn’t you use Combinatorics.jl or similar to get all possible combinations of terms?

---

<div class="post-metadata">

**Author:** ![WuSiren](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wusiren/32/42529_2.png) [@WuSiren](https://discourse.julialang.org/u/WuSiren)\
**Post date:** [January 15, 2024, 3:37am UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/15 "2024-01-15T03:37:03Z")

</div>

Thank you, @stevengj , for your patient and professional answers! 🤝🤝

In fact, I was referring to the second case you mentioned, which is essentially a least square problem with 351 parameters. I’m sorry I didn’t make it clear at the beginning.

* * *

> [@stevengj](#):
>
> There are lots of high-dimensional regression models (from sparse grids to kernel methods to neural nets) that try to avoid the curse of dimensionality. But they all trade off some generality, and go beyond simple polynomial regression.

You mentioned kernel methods and neural networks, and I am wondering if least square regression based on polynomial function mapping can also be considered a kernel method? Or can it be seen as a special kind of neural network? It would be nice to have some kind of literature on how they relate to each other. 🤔

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 15, 2024, 7:22pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/16 "2024-01-15T19:22:44Z")

</div>

Well, `*` on terms does expand into all combinations as the formula `y ~ x1 * x2 * x3` means all interactions up to that order in R. Thus,

```julia-repl
julia> @formula y ~ x1 * x2 * x3
FormulaTerm
Response:
  y(unknown)
Predictors:
  x1(unknown)
  x2(unknown)
  x3(unknown)
  x1(unknown) & x2(unknown)
  x1(unknown) & x3(unknown)
  x2(unknown) & x3(unknown)
  x1(unknown) & x2(unknown) & x3(unknown)

julia> vars = [:x1, :x2, :x3];

julia> prod(Term, vars)
x1(unknown)
x2(unknown)
x1(unknown) & x2(unknown)
x3(unknown)
x1(unknown) & x3(unknown)
x2(unknown) & x3(unknown)
x1(unknown) & x2(unknown) & x3(unknown)

# This is the interaction between x1, x2 and x3 alone (written as x1 : x2 : x3 in R formulas)
julia> reduce(&, Term.(vars))
x1(unknown) & x2(unknown) & x3(unknown)

```

---

<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 15, 2024, 7:24pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/17 "2024-01-15T19:24:15Z")

</div>

> [@bertschi](#):
>
> Well, `*` on terms does expand into all combinations as the formula `y ~ x1 * x2 * x3` means all interactions up to that order in R.

Ah, thanks — I didn’t realize that `*` meant something different from multiplication here.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 15, 2024, 7:29pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/18 "2024-01-15T19:29:43Z")

</div>

Yes, polynomial regression can be considered as a kernel method. The kernel is non-stationary with a variance that quickly grows in |\mathbf{x}| and a specific weighting of each interaction order (for the glorious details see chapter 4.2.2 of [GPML](http://gaussianprocess.org/gpml/chapters/RW.pdf)).  
Shallow neural networks give rise to a different kernel which is also non-stationary, but becomes constant for large |\mathbf{x}| (chapter 4.2.3 of the same book).

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [January 15, 2024, 9:39pm UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/19 "2024-01-15T21:39:11Z")

</div>

The syntax is a DSL called Wilkinson-Roger notation, see also [Use some character other than \* for "main effects and interactions" · Issue #99 · JuliaStats/StatsModels.jl · GitHub](https://github.com/JuliaStats/StatsModels.jl/issues/99#issuecomment-633037379) for more discussion

---

<div class="post-metadata">

**Author:** ![WuSiren](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/wusiren/32/42529_2.png) [@WuSiren](https://discourse.julialang.org/u/WuSiren)\
**Post date:** [January 16, 2024, 2:34am UTC](https://discourse.julialang.org/t/is-there-any-julia-package-for-multivariate-polynomial-regression/108785/20 "2024-01-16T02:34:33Z")

</div>

Thanks a lot!!! 🤝 🤝