Nonlinear Feature Engineering

Nonlinear Feature Engineering#

Learning Objectives

By the end of this chapter, you will be able to:

  • Explain why nonlinear feature transformations are necessary when linear combinations cannot improve a linear model’s performance.

  • Construct polynomial feature matrices with sklearn.preprocessing.PolynomialFeatures and estimate the combinatorial growth of feature count with polynomial order.

  • Apply LASSO regularization to a polynomial feature matrix to select a sparse, interpretable set of nonlinear features.

  • Use autofeat to generate and select physically meaningful nonlinear features from engineering data, with and without dimensional constraints.

  • Articulate the practical trade-offs of symbolic regression: interpretability and efficiency at test time vs. high training cost and sensitivity to train/test splits.

  • Describe how automated model discovery methods extend feature generation to searching for the functional form of a model, and identify representative tools used in science and in chemical engineering.

%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, Lasso
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler

plt.style.use('../settings/plot_style.mplstyle')

clrs = np.array(['#003057', '#EAAA00', '#4B8B9B', '#B3A369', '#377117',
                 '#1879DB', '#8E8B76', '#F5D580', '#002233'])

# Load Dow dataset — same preprocessing as Topic 2.4/2.5
df = pd.read_excel('data/impurity_dataset-training.xlsx')

def is_real_and_finite(x):
    try:
        val = float(x)
        return np.isfinite(val)
    except (TypeError, ValueError):
        return False

nondate_cols = df.columns[1:]
numeric_map = df[nondate_cols].apply(lambda col: col.map(is_real_and_finite))
real_rows = numeric_map.all(axis=1).values

all_data = df[nondate_cols].values
dow_feature_names = list(nondate_cols)

X_dow = np.array(all_data[real_rows, :-5], dtype='float')
y_dow = np.array(all_data[real_rows, -3],  dtype='float')

# Standardize features
X_dow_scaled = (X_dow - X_dow.mean(axis=0)) / X_dow.std(axis=0)
print(f'X shape: {X_dow_scaled.shape},  y shape: {y_dow.shape}')

X shape: (10297, 40),  y shape: (10297,)

This is the same preparation used in Topics 2.4 and 2.5, repeated here so the chapter runs on its own. The is_real_and_finite helper is needed because the raw spreadsheet contains text placeholders in some cells, so real_rows keeps only the rows where every column is a usable number. The column slicing selects the process variables as inputs (:-5) and the impurity as the target (-3), and the features are standardized to zero mean and unit variance. Standardization matters more than usual in this chapter: we are about to multiply features together, and a product of two variables with very different magnitudes would otherwise dominate the feature matrix for reasons of units rather than physics.

Motivation: When Linear Combinations Are Not Enough#

In Topic 2.5 we showed that linear combinations — including PCA and PLS — can reduce dimensionality and sometimes improve model performance. But they share a fundamental limitation: a linear combination of features fed into a linear model is still a linear model. No matter how we rotate or re-weight the input columns, we cannot capture interactions like \(x_i \cdot x_j\) or nonlinear effects like \(x_i^2\).

This matters for engineering systems, where the physics is rarely linear. A distillation column impurity may depend on the product of two flow rates (a multiplicative interaction), or on the square of a temperature difference. To capture such effects with a linear regression model, we must create the nonlinear features explicitly — add columns like \(x_i^2\) or \(x_i x_j\) to the feature matrix before fitting.

The challenge is that the number of such features grows rapidly, so we need both a systematic way to generate them and a regularization strategy to select the useful ones.

Exercise 42

Verify the linear-model ceiling on the Dow dataset.

  1. Using the scaled features X_dow_scaled and target y_dow from above, fit a LinearRegression() and record the train \(r^2\).

  2. Fit a Ridge(alpha=1.0) and compare the train and test \(r^2\) (use a 70/30 split with random_state=0). Is there evidence of overfitting?

  3. Compute the Pearson correlation between each individual feature and the target. Report the top-5 most correlated features. Are any correlations very strong (\(|r| > 0.5\))?

  4. Based on your results, argue briefly whether a linear model is likely to be the best achievable model for this dataset.

Polynomial Features#

Generating Polynomial Features with sklearn#

sklearn.preprocessing.PolynomialFeatures creates all monomials up to a specified degree from an input matrix. For degree 2 with two input features \(x_0, x_1\), it produces \([1,\, x_0,\, x_1,\, x_0^2,\, x_0 x_1,\, x_1^2]\).

from sklearn.preprocessing import PolynomialFeatures

# Subsample to keep memory manageable (see note below)
X_sub = X_dow_scaled[::2]
y_sub = y_dow[::2]

X_train, X_test, y_train, y_test = train_test_split(
    X_sub, y_sub, test_size=0.5, random_state=0)

poly = PolynomialFeatures(degree=2, include_bias=False)
X_poly_train = poly.fit_transform(X_train)
X_poly_test  = poly.transform(X_test)

print(f'Original features:  {X_train.shape[1]}')
print(f'Degree-2 features:  {X_poly_train.shape[1]}')

Original features:  40
Degree-2 features:  860

Note

Why subsample? With 40 features and degree-2 polynomial expansion, the feature count is \(\binom{40 + 2}{2} - 1 = 860\). At degree 3 it exceeds 12,000, which is larger than many training sets. We subsample to avoid memory issues during the fitting step below. In practice you would use regularization (LASSO, ridge) rather than subsampling to handle the large-feature regime.

Combinatorial Explosion#

The exact count of degree-\(d\) polynomial features from \(p\) original features (without the intercept) is the multiset coefficient \(\binom{p + d}{d} - 1\):

from math import comb

p = 40  # original features
for d in range(1, 5):
    n_poly = comb(p + d, d) - 1
    print(f'Degree {d}: {n_poly:,} features')

Degree 1: 40 features
Degree 2: 860 features
Degree 3: 12,340 features
Degree 4: 135,750 features

Beyond degree 2 or 3, the number of features quickly exceeds the number of data points, making unregularized least squares completely ill-conditioned. This is one manifestation of the curse of dimensionality.

Regularized Polynomial Regression with LASSO#

With more features than samples, standard linear regression overfits badly. LASSO regularization shrinks many coefficients to exactly zero, effectively selecting a sparse subset of the polynomial features.

To see the effect clearly we fit three models on the same train/test split, changing one thing at a time. The first is a plain linear regression on the original 40 features, which serves as the baseline to beat:

# Baseline: plain linear regression on original features
linreg = LinearRegression()
linreg.fit(X_train, y_train)
r2_base_train = linreg.score(X_train, y_train)
r2_base_test  = linreg.score(X_test,  y_test)
print(f'Linear regression:  train r² = {r2_base_train:.3f},  test r² = {r2_base_test:.3f}')

Linear regression:  train r² = 0.723,  test r² = 0.660

The second uses the 860 polynomial features but no regularization. With that many features relative to the number of training samples, we expect the training \(r^2\) to improve substantially while the test \(r^2\) gets worse, which is the signature of overfitting from Topic 2.3:

# Unregularized: polynomial features + ordinary least squares — should overfit
linreg_poly = LinearRegression()
linreg_poly.fit(X_poly_train, y_train)
r2_poly_train = linreg_poly.score(X_poly_train, y_train)
r2_poly_test  = linreg_poly.score(X_poly_test,  y_test)
print(f'Poly + OLS:         train r² = {r2_poly_train:.3f},  test r² = {r2_poly_test:.3f}')

Poly + OLS:         train r² = 0.949,  test r² = -132958.239

The third uses the same polynomial features with a LASSO penalty, and also counts how many coefficients survive:

# LASSO on polynomial features — sparse selection
lasso = Lasso(alpha=1.0, max_iter=5000)
lasso.fit(X_poly_train, y_train)
r2_lasso_train = lasso.score(X_poly_train, y_train)
r2_lasso_test  = lasso.score(X_poly_test,  y_test)
n_nonzero = np.sum(lasso.coef_ != 0)
print(f'Poly + LASSO (α=1): train r² = {r2_lasso_train:.3f},  test r² = {r2_lasso_test:.3f}')
print(f'                    non-zero coefficients: {n_nonzero} of {len(lasso.coef_)}')

Poly + LASSO (α=1): train r² = 0.052,  test r² = 0.052
                    non-zero coefficients: 4 of 860

The LASSO reduces training accuracy compared to unregularized polynomial regression but improves test accuracy significantly — and with far fewer nonzero coefficients.

Interpreting the Selected Features#

A key advantage over black-box approaches is that PolynomialFeatures can name every output column, so a sparse model can be read as an equation. The cell below asks for those names, finds the coefficients that LASSO did not set to zero, and sorts them by magnitude so that the most influential terms appear first:

feature_names = poly.get_feature_names_out(
    [f'x{i}' for i in range(X_train.shape[1])])

coef = lasso.coef_
nonzero_idx = np.where(coef != 0)[0]
sorted_idx = nonzero_idx[np.argsort(np.abs(coef[nonzero_idx]))[::-1]]

print('Top 10 selected polynomial features (by |coefficient|):')
for i in sorted_idx[:10]:
    print(f'  {feature_names[i]:30s}  coef = {coef[i]:+.4f}')

Top 10 selected polynomial features (by |coefficient|):
  x21^2                           coef = -0.0099
  x35^2                           coef = -0.0077
  x7^2                            coef = -0.0050
  x34^2                           coef = -0.0001

The feature names reveal which interactions and squared terms the model found useful. In engineering problems, a feature like x5^2 or x3 x7 may have a clear physical interpretation (e.g., a quadratic temperature effect or a cross-flow interaction).

Exercise 43

Compare the test \(r^2\) of LASSO polynomial regression against plain linear regression as a function of the LASSO regularization strength \(\alpha\).

  1. Using the same 50/50 split above, sweep \(\alpha\) over np.logspace(-2, 2, 20).

  2. For each \(\alpha\), fit Lasso(alpha=alpha, max_iter=5000) on X_poly_train and record train and test \(r^2\) and the number of nonzero coefficients.

  3. Plot train \(r^2\), test \(r^2\), and nonzero coefficient count vs. \(\log_{10}(\alpha)\) on two subplots (shared x-axis).

  4. Identify the \(\alpha\) that maximizes test \(r^2\) and report the corresponding number of nonzero features. How does this compare to the baseline linear regression test \(r^2\)?

Symbolic Regression with autofeat#

Polynomial features restrict transformations to integer powers of original features. Symbolic regression extends this to arbitrary nonlinear combinations — ratios, exponentials, square roots, and products thereof — and uses a built-in feature selection step to keep only the combinations that improve prediction.

The autofeat library implements this workflow: it generates a large library of nonlinear feature candidates, then applies multi-stage feature selection to reduce to a compact, interpretable set.

Note

autofeat is not installed by default. Install it with: pip install autofeat or conda install -c conda-forge autofeat

Basic AutoFeatRegressor#

try:
    from autofeat import AutoFeatRegressor
    _autofeat_available = True
except ImportError:
    print('autofeat not installed — skipping symbolic regression cells.')
    print('Install with: pip install autofeat')
    _autofeat_available = False

if _autofeat_available:
    # Use the same subsampled/split data from above
    transforms = ['1/', 'exp', 'abs', 'sqrt', '^2', '^3']
    afreg = AutoFeatRegressor(
        verbose=1,
        feateng_steps=2,
        featsel_runs=1,
        transformations=transforms,
    )
    afreg.fit(X_train, y_train)

autofeat not installed — skipping symbolic regression cells.
Install with: pip install autofeat

The feateng_steps=2 parameter controls how many times transformations are nested (e.g., sqrt(x0 * exp(x1)) requires 2 steps). Higher values create exponentially more candidate features.

if _autofeat_available:
    # New nonlinear features selected as predictive
    print('New features generated:')
    print(afreg.new_feat_cols_)

These are the nonlinear combinations that autofeat constructed and judged predictive. It is worth also inspecting good_cols_, which is a different quantity: new_feat_cols_ lists only the newly generated features, while good_cols_ lists every feature the final model actually uses, including whichever of the original 40 variables survived selection:

if _autofeat_available:
    # All features (original + new) retained after final selection
    print('All selected features (used in model):')
    print(afreg.good_cols_)

The distinction matters when interpreting the model, and it also explains the overfitting comparison at the end of this section: scoring the fitted AutoFeatRegressor uses only good_cols_, whereas transforming the data ourselves returns the full generated matrix.

if _autofeat_available:
    r2_af_test = afreg.score(X_test, y_test)
    print(f'autofeat test r²: {r2_af_test:.3f}')

The selected features include ratios, square roots, and products that a polynomial expansion would not generate. Note that results vary across runs due to the random train/test split; do not over-interpret any specific feature unless it has clear physical meaning and is consistent across multiple runs.

Transforming the Feature Matrix Directly#

We can also access the full transformed feature matrix and build our own model:

if _autofeat_available:
    X_af_train = afreg.transform(X_train)
    X_af_test  = afreg.transform(X_test)

    linreg2 = LinearRegression()
    linreg2.fit(X_af_train, y_train)
    r2_af_linreg = linreg2.score(X_af_test, y_test)
    print(f'Linear model on autofeat features: test r² = {r2_af_linreg:.3f}')

If this is much lower than afreg.score(), it indicates overfitting because the full transformed matrix includes features beyond those in good_cols_.

Units-Aware Feature Generation#

One of autofeat’s most powerful capabilities is dimensional analysis: if you supply physical units for each input feature, it will only generate dimensionally consistent combinations. This dramatically reduces the number of candidates and biases the search toward physically meaningful features.

Supplying those units requires knowing what each column represents. The Dow column names end in words like “Flow”, “Level”, or “Temperature”, so the cell below assigns units by reading that suffix and marks anything unrecognized as dimensionless. This is a convenient shortcut rather than a rigorous one: a mislabeled column would silently constrain the search in the wrong way, so in your own work it is worth checking the unit assignments against the process documentation. autofeat also needs the feature names themselves, which is why the arrays are wrapped in DataFrames:

if _autofeat_available:
    # Assign physical units to the Dow features based on column name suffixes
    unit_dict = {}
    for var in dow_feature_names[:X_train.shape[1]]:
        last_word = var.split(' ')[-1]
        if last_word == 'Flow':
            unit_dict[var] = 'L/s'
        elif last_word == 'Level':
            unit_dict[var] = 'm'
        elif last_word == 'DP':
            unit_dict[var] = 'bar'
        elif last_word == 'Pressure':
            unit_dict[var] = 'bar'
        elif last_word == 'Temperature':
            unit_dict[var] = 'K'
        else:
            unit_dict[var] = ''   # dimensionless

    # Build DataFrames so autofeat knows which feature is which
    dow_input_names = dow_feature_names[:X_train.shape[1]]
    X_train_df = pd.DataFrame(X_train, columns=dow_input_names)
    X_test_df  = pd.DataFrame(X_test,  columns=dow_input_names)

With units supplied, we can afford a larger pool of transformations than before, because the dimensional constraint will discard most of the combinations they generate. The units= argument is the only substantive change from the earlier fit:

if _autofeat_available:
    transforms_extended = ['1/', 'exp', 'log', 'abs', 'sqrt', '^2', '^3', '1+', '1-', 'exp-']
    afreg_units = AutoFeatRegressor(
        verbose=1,
        feateng_steps=2,
        featsel_runs=1,
        transformations=transforms_extended,
        units=unit_dict,
    )
    afreg_units.fit(X_train_df, y_train)

if _autofeat_available:
    r2_af_units = afreg_units.score(X_test_df, y_test)
    print(f'autofeat (units-aware) test r²: {r2_af_units:.3f}')
    print('\nNew features (dimensionless combinations):')
    print(afreg_units.new_feat_cols_)

By constraining the search to dimensionless combinations, autofeat finds features that correspond to physically interpretable ratios and products — the kind of dimensionless groups (analogous to Reynolds or Damköhler numbers) that engineers use to characterize process behavior.

Practical Guidance#

When is symbolic regression worth using?

Situation

Recommendation

Dataset is small–medium (\(n < 10{,}000\)) and you want an interpretable model

autofeat is viable; budget 5–30 min of training time

Features have known physical units and you want dimensionless groups

Use autofeat with units= — the dimensional constraint keeps the feature count tractable

You need fast training (batch re-training, CI pipelines)

Use LASSO polynomial instead; faster and more predictable

Dataset is large (\(n > 50{,}000\)) or you have \(>\)50 input features

Symbolic regression becomes very slow; prefer tree-based models or neural networks

You need to deploy a model that other engineers can inspect

Symbolic/polynomial regression produces plain algebraic expressions — ideal for documentation

The key caution is train/test sensitivity: because many feature combinations give similar \(r^2\) values on the training set, different random splits may yield very different selected features. Always run a few different splits and check whether the same features appear consistently before drawing conclusions about physical importance.

Exercise 44

Apply units-aware AutoFeatRegressor to the Dow dataset.

  1. Verify that the features in afreg_units.new_feat_cols_ are dimensionless (every generated feature should be a dimensionless combination of input variables, given the unit assignments above).

  2. Compare the test \(r^2\) of the units-aware model with the non-units-aware model and the LASSO polynomial model from the earlier exercise. Which performs best?

  3. Pick one of the selected features and write a one-sentence physical interpretation of what it represents (e.g., a ratio of flow rates, a product of pressures).

Beyond Feature Selection: Automated Model Discovery#

It is worth stepping back to see what autofeat does and does not do. It generates a large library of candidate nonlinear terms and then selects a useful subset, but the model itself is still a linear combination of those terms. The functional form is therefore fixed by the transformations we allowed in the search: we chose the pool, and the algorithm chose from it. The same is true of the LASSO polynomial models earlier in this chapter, and of the regularized models in Complexity and Optimization.

Model discovery (often called symbolic regression in its more general form) removes that restriction and searches over the structure of the equation itself, with the goal of returning a compact expression that a person can read, interpret, and check against physical intuition. The difficulty is that the space of possible equations is astronomically large, so every practical method imposes strong priors on what a reasonable model looks like. This is the same principle as regularization, applied to the structure of the model rather than to the size of its coefficients. Three examples are worth knowing:

  • AI-Feynman (Udrescu & Tegmark, 2020, Science Advances) exploits properties that physical laws tend to possess: dimensional consistency, symmetry, separability into independent parts, and smoothness. A neural network is fit to the data first, not as the final model, but as a probe that can be queried to test for these properties, which then break the problem into simpler sub-problems recursively. The method recovered all 100 equations from the Feynman Lectures on Physics, where the best previously available software found 71. The code is openly available. Note that the dimensional-consistency idea is exactly the units= constraint we used with autofeat above, applied to a much larger search.

  • AI-DARWIN (Chakraborty, Sivaram & Venkatasubramanian, 2021, Computers & Chemical Engineering) is aimed specifically at chemical engineering problems, where data are typically limited and noisy rather than abundant. It uses a genetic algorithm to evolve candidate nonlinear terms, statistical testing to retain only the terms the data actually support, and nonlinear regression to fit the remaining parameters. The result is a mechanistic-looking model rather than a black box, which is the motivation behind the broader argument that hybrid models are often more appropriate than purely data-driven ones in our field.

  • HyMech (Rossi, Bezzo & Barolo, 2026, Computers & Chemical Engineering) builds directly on the AI-DARWIN engine and adds prior process knowledge to the search through physics-informed pools of allowed functions and variables. It addresses a situation that is common in practice: you already have a first-principles model, it does not quite match the plant data, and you want to know where the model structure is deficient. Conventional hybrid modeling patches such a mismatch with a black-box correction term, which improves predictions without explaining anything; HyMech instead searches for an interpretable equation for the correction itself.

These methods are directly relevant to industrial practice, and in fact to the dataset used throughout this course. The Dow impurity data was originally posed by researchers at Dow as a “data challenge” problem: given more than 40 process variables from an integrated multi-column process, identify which of them actually drive the impurity measured at the primary column outlet, and build an inferential sensor to predict it. Identifying the relevant variables and producing a model that an engineer can inspect and trust is precisely the goal of the methods above, and it is the same problem we have been working on in this chapter with polynomial features and autofeat. Work in this direction is ongoing, including methods developed jointly with industry such as SyMANTIC (Muthyala et al., 2025, Ind. Eng. Chem. Res.), which combines mutual-information-based feature screening with sparse regression to keep the search tractable.

Full model discovery is beyond the scope of this course, and these tools require more care and computational effort than the methods we have used here. They are worth knowing about, however, because they represent the natural endpoint of the progression in this chapter: from choosing features by hand, to generating and selecting them automatically, to searching for the form of the model itself.

Exercise 45

Estimate why unrestricted equation search is impractical, and how much a physical prior helps.

  1. Consider building expressions as binary trees with \(n\) internal nodes. Each internal node is one of \(k\) binary operators (\(+, -, \times, \div\), so \(k = 4\)) and each of the \(n+1\) leaves is one of \(p\) input variables. The number of tree shapes is the Catalan number \(C_n = \binom{2n}{n}/(n+1)\), so the total count of expressions is \(C_n \, k^n \, p^{\,n+1}\). Using math.comb, compute and print this count for \(p = 40\) (the Dow features) and \(n = 1, 2, \ldots, 6\).

  2. Plot the count versus \(n\) on a log scale. At what \(n\) does the count exceed the roughly \(10^3\) candidate features that PolynomialFeatures produced at degree 2?

  3. Now suppose dimensional analysis rules out all but 1% of these expressions as dimensionally inconsistent. Re-plot the surviving count. Does the constraint change the rate of growth, or only the prefactor? What does your answer imply about relying on dimensional analysis alone to make the search tractable?

Summary#

  • Linear feature combinations (PCA, PLS, random projections) cannot improve a linear model beyond what the original features can achieve; nonlinear feature engineering is necessary to capture interactions and higher-order effects.

  • Polynomial features with PolynomialFeatures generate all monomials up to a specified degree. Feature count grows as \(\binom{p+d}{d}\), making unregularized fitting impossible beyond degree 2–3 for typical engineering datasets.

  • LASSO regularization on a polynomial feature matrix simultaneously performs feature selection (setting many coefficients to zero) and regularization (preventing overfitting), often yielding better test accuracy than ordinary linear regression.

  • Symbolic regression (autofeat) generates a broader library of nonlinear candidates (ratios, exponentials, roots) and selects the most predictive subset. Results are compact and interpretable but sensitive to train/test splits; physical units constraints (units=) improve robustness and chemical interpretability.

  • Automated model discovery goes one step further by searching over the functional form of the model rather than selecting from a fixed pool of generated terms. Tools such as AI-Feynman, AI-DARWIN, and HyMech make the search tractable by imposing physical priors (dimensional consistency, symmetry, or an existing first-principles model), which is regularization applied to model structure.

Additional Reading#

  • Kanter, J. M. & Veeramachaneni, K. (2015), “Deep Feature Synthesis: Towards Automating Data Science Endeavors” — foundational automated feature engineering

  • Orzechowski, P. et al. (2018), “Where are we now? A large benchmark study of recent symbolic regression methods” — survey of symbolic regression approaches

  • Tibshirani, R. (1996), “Regression Shrinkage and Selection via the Lasso,” Journal of the Royal Statistical Society B — original LASSO paper

  • Udrescu, S.-M. & Tegmark, M. (2020), “AI Feynman: A physics-inspired method for symbolic regression,” Science Advances 6(16), eaay2631 — doi:10.1126/sciadv.aay2631

  • Chakraborty, A., Sivaram, A. & Venkatasubramanian, V. (2021), “AI-DARWIN: A first principles-based model discovery engine using machine learning,” Computers & Chemical Engineering 154, 107470 — doi:10.1016/j.compchemeng.2021.107470

  • Rossi, L., Bezzo, F. & Barolo, M. (2026), “HyMech: AI-driven framework for physics-informed discovery of interpretable hybrid models,” Computers & Chemical Engineering 210, 109634 — doi:10.1016/j.compchemeng.2026.109634

  • Muthyala, M. R., Sorourifar, F., Peng, Y. & Paulson, J. A. (2025), “SyMANTIC: An efficient symbolic regression method for interpretable and parsimonious model discovery,” Industrial & Engineering Chemistry Research 64(6), 3354–3369 — doi:10.1021/acs.iecr.4c03503

  • scikit-learn User Guide: PolynomialFeatures, Lasso

Contents