Introduction

A multivariate polynomial fit models one output variable from two or more input variables using polynomial terms such as x^2, xy, or y^3. In Python, the easiest way to build this kind of best-fit surface is usually to expand the input features with polynomial terms and then fit a linear regression model on that expanded design matrix.

Why Polynomial Regression Still Uses Linear Regression

Polynomial regression sounds nonlinear, but once you create the polynomial features, the model is linear in its coefficients. For two variables x1 and x2, a degree-2 model might look like:

• constant term
• 'x1'
• 'x2'
• 'x1^2'
• 'x1 x2'
• 'x2^2'

You estimate the coefficients with an ordinary regression algorithm.

Use PolynomialFeatures with LinearRegression

Scikit-learn makes this workflow straightforward.

1import numpy as np
2from sklearn.preprocessing import PolynomialFeatures
3from sklearn.linear_model import LinearRegression
4
5X = np.array([
6    [1.0, 2.0],
7    [2.0, 1.0],
8    [3.0, 4.0],
9    [4.0, 3.0],
10    [5.0, 5.0],
11])
12
13y = np.array([5.0, 4.0, 15.0, 14.0, 25.0])
14
15poly = PolynomialFeatures(degree=2, include_bias=False)
16X_poly = poly.fit_transform(X)
17
18model = LinearRegression()
19model.fit(X_poly, y)
20
21print(poly.get_feature_names_out(["x1", "x2"]))
22print(model.coef_)
23print(model.intercept_)

X_poly now contains the original variables plus interaction and power terms.

Make Predictions on New Points

Once trained, predict like any other regression model.

1new_points = np.array([
2    [2.5, 2.5],
3    [6.0, 4.0],
4])
5
6predictions = model.predict(poly.transform(new_points))
7print(predictions)

This gives you the fitted polynomial surface evaluated at the new coordinates.

Build It as a Pipeline

A pipeline is cleaner because it keeps feature expansion and fitting together.

1from sklearn.pipeline import make_pipeline
2
3pipeline = make_pipeline(
4    PolynomialFeatures(degree=3, include_bias=False),
5    LinearRegression(),
6)
7
8pipeline.fit(X, y)
9print(pipeline.predict(new_points))

This is often the best form for real projects because training and prediction always apply the same transformation steps.

Choosing the Degree

A higher degree makes the model more flexible, but also raises the risk of overfitting. With multiple variables, the number of terms grows quickly. That means:

• Degree 2 is often a sensible starting point
• Higher degrees need more data
• Validation matters more than raw training accuracy

For noisy data, a lower-degree polynomial often generalizes better than an elaborate surface that chases every fluctuation.

Visualizing a Two-Variable Fit

If you only have two predictors, you can visualize the fitted surface by evaluating it on a grid.

1x1 = np.linspace(1, 5, 30)
2x2 = np.linspace(1, 5, 30)
3grid_x1, grid_x2 = np.meshgrid(x1, x2)
4
5grid = np.column_stack([grid_x1.ravel(), grid_x2.ravel()])
6grid_y = pipeline.predict(grid).reshape(grid_x1.shape)

You can then plot grid_y with Matplotlib as a surface or contour map. That helps you see whether the fitted polynomial shape is reasonable.

Common Pitfalls

A common mistake is using a very high polynomial degree with a small dataset. The fit may look excellent on training data while behaving wildly between points.

Another mistake is forgetting that feature count grows rapidly with more variables and higher degrees. The model can become unstable or slow sooner than expected.

A third mistake is calling it a “curve” when the problem is multivariate. With two or more inputs, the fitted object is usually a surface or hypersurface rather than a single 2D curve.

Summary

• Multivariate polynomial fitting is usually implemented by expanding features and then fitting linear regression.
• 'PolynomialFeatures and LinearRegression are the standard Python tools for this task.'
• Pipelines make the workflow cleaner and safer.
• Degree selection matters because feature counts grow quickly.
• Validate the model instead of assuming a higher-degree polynomial is automatically better.