Classical Machine Learning Algorithms Questions

Direct answer

Add polynomial terms (x, x^2, x^3, ...) or interaction terms (x1*x2) as new columns and fit ordinary linear regression on the expanded feature set. The model stays linear in its parameters, so it keeps the same fitting machinery and interpretability, but as you raise the degree the parameter count grows and, with a fixed number of data points, you can eventually have as many parameters as data points, at which point the model can memorize the training data and generalizes badly.

Structured elaboration

Why it still works as linear regression: linear regression only requires linearity in the parameters, not the raw inputs. Feeding it phi(x) = [1, x, x^2, ...] lets it fit curves while remaining ordinary linear regression under the hood.

Parameter growth: a degree-d polynomial on one feature adds d terms; with p original features and all pairwise interactions, the interaction count grows like O(p^2), and a full degree-d expansion across p features grows combinatorially.

Effect on bias and variance: low degree underfits (high bias, misses real curvature) but has low variance. Higher degree lowers bias but variance grows quickly, because more parameters are being estimated from the same amount of data, and the added polynomial columns are highly correlated with each other (x and x^2 move together for x > 0), which is the same instability mechanism as ordinary multicollinearity.

Controls:

• Cross-validate the degree instead of picking it by training error.
• Standardize or use an orthogonal polynomial basis to reduce numerical correlation between terms.
• Apply ridge or lasso on top of the expanded features so unnecessary high-degree terms shrink toward zero.
• Prefer splines when only certain regions of x are nonlinear, since local basis functions don't oscillate far from the data the way a single high-degree global polynomial can. (This oscillation near the edges of the fitted range, worse at higher degree, is known as Runge's phenomenon.)

Worked example

Fit y from x with only n = 3 points: x = [0, 1, 2], y = [1, 3, 2].

A degree-1 fit (a line, 2 parameters) via least squares gives intercept 1.5 and slope 0.5:

$\hat y = 1.5 + 0.5x$

giving predictions [1.5, 2.0, 2.5] and a sum of squared errors of 1.5, an honest approximation that doesn't pass through every point.

A degree-2 fit (a parabola, 3 parameters, exactly matching the 3 data points) solves exactly:

$\hat y = 1 + 3.5x - 1.5x^2$

giving predictions [1.0, 3.0, 2.0], a perfect match, sum of squared errors = 0. That zero error is not evidence of a better model, it's evidence that zero degrees of freedom are left: with as many parameters as data points, the parabola threads exactly through every point, including whatever noise those points happen to contain, and there is no held-out data left to check whether it generalizes. In practice this parabola would swing wildly if evaluated between or beyond the three fitted points.

Trade-offs & pitfalls

Chasing training R^2 by raising the degree without cross-validation is the classic trap this example illustrates directly. Interaction terms without polynomial terms should still be checked for multicollinearity using VIF (variance inflation factor: how much a term's coefficient variance is inflated by its correlation with the other predictors, with values above roughly 5-10 flagged as concerning), since x1*x2 correlates with both x1 and x2 individually. Extrapolating a polynomial fit outside the training range is especially dangerous because polynomial curvature can diverge quickly. High-degree polynomial design matrices are also numerically unstable to compute directly (they resemble ill-conditioned Vandermonde matrices); centering and scaling x, or using an orthogonal polynomial basis, mitigates this on top of the statistical fix.

Direct answer

Nested cross-validation is two CV loops in one: an outer loop that estimates how well the whole modeling pipeline generalizes, and an inner loop, run entirely inside each outer training fold, that picks hyperparameters or features. You need it whenever the same data is used both to tune the model and to report its performance, because tuning on the same data you evaluate on leaks information and inflates the reported score.

Structured elaboration

Why plain k-fold CV is not enough here. If you run k-fold CV once to pick the best hyperparameters (by, say, taking the config with the highest mean CV score) and then report that same mean CV score as your generalization estimate, you have used the test folds to make a selection decision, so the reported number is optimistically biased. The gap grows with the size of the hyperparameter search space: the more configurations you try, the more likely one of them fits the validation folds' noise, not just signal.

Outer/inner structure.

• Outer loop ($K_\text{outer}$ folds): each outer fold is held out entirely and touched only once, at the very end, purely for scoring. It never influences any modeling decision.
• Inner loop ($K_\text{inner}$ folds, run on the outer-training portion only): performs the hyperparameter search or feature selection, using its own train/validation splits. Whatever it selects (say, the config with the best mean inner-validation score) is refit on the full outer-training set.
• The refit model is then scored once on the untouched outer-test fold. Averaging that score across all outer folds gives an unbiased estimate of how well "the pipeline, including its tuning procedure" generalizes.
• Any preprocessing that looks at the target (target encoding, feature selection by correlation with $y$, scaling parameters) must be fit only on the current inner-training data, never on the inner-validation or outer-test data, or the same leakage reappears one level down.

Computational cost. Every hyperparameter configuration gets fit $K_\text{inner}$ times per outer fold just to be scored, then the winner gets refit once more on the full outer-training set. Total model fits:

where $G$ is the number of hyperparameter configurations evaluated. This is roughly $K_\text{inner}$ times more expensive than a single k-fold search (precisely $K_\text{inner} + 1/G$), which is why nested CV is usually reserved for the final reported generalization number rather than for every exploratory tuning pass.

Worked example

Suppose $K_\text{outer} = 5$, $K_\text{inner} = 3$, and a grid search over $G = 10$ hyperparameter configurations.

Compare that to a single (non-nested) k-fold grid search used just to pick hyperparameters, at $G \times K = 10 \times 5 = 50$ fits, with no separately reported unbiased generalization estimate. Nested CV costs about $155 / 50 \approx 3.1\times$ more fits here, and the ratio grows directly with $K_\text{inner}$: doubling $K_\text{inner}$ to 6 gives $5 \times (10 \times 6 + 1) = 305$ fits, roughly double, because the inner search dominates the total.

Trade-offs & pitfalls

• Nested CV answers "how good is this tuning procedure, on average," not "what hyperparameters should I ship." The winning configuration can differ across outer folds; for deployment, refit once on the full dataset using the inner-loop procedure (or the single most frequently selected configuration) after nested CV has validated that the procedure is trustworthy.
• Under a tight compute budget, replace grid search with randomized or Bayesian search in the inner loop to shrink $G$ without shrinking the search space explored, or reduce $K_\text{inner}$ to 3 (a common compromise, since the inner loop only needs to rank configurations relatively, not report a final number).
• Skipping the inner loop and just using a single train/val split inside each outer fold is a cheaper approximation, but reintroduces some tuning variance into the outer score; it is a reasonable trade-off for very expensive models, not for cheap ones where full nested CV is affordable.
• The single most common implementation bug: fitting a preprocessing step (scaler, target encoder, feature selector) once on the whole dataset before either loop starts. That silently defeats the entire point of nesting.

Direct answer

For a 100k-row, 200-feature tabular dataset with high-cardinality categoricals, I'd default to gradient boosting (LightGBM/CatBoost) over a neural network with embeddings, because GBMs typically reach strong accuracy on data of this size with less tuning, handle missingness and mixed scales natively, and have native categorical handling that avoids the embedding-design and data-volume concerns a neural net brings. A neural net becomes the better choice mainly when the dataset is substantially larger, the categoricals carry rich learnable semantics worth sharing across tasks, or the model needs to plug into an end-to-end differentiable pipeline.

Structured elaboration

Why sample size at 100k matters. Neural nets, and embeddings for high-cardinality categoricals in particular, need enough examples per category to learn a meaningful representation. With 200 features and 100k rows, a high-cardinality categorical (say, 5,000 distinct values) may only have on the order of 20 examples per category on average, not much for a network to learn a useful embedding vector from scratch, whereas a tree-based split on that same categorical (using target statistics) is a much lower-variance way to extract signal from a similarly sparse per-category sample.

When the neural net does start to win. As data grows well past 100k, embeddings get enough examples per category to learn genuinely useful representations, sometimes capturing similarity structure between categories (e.g. two product SKUs behaving similarly) that a GBM's per-feature splits don't naturally exploit. A neural net is also the better choice when there's a reason to want a differentiable end-to-end pipeline (e.g. joint training with another neural component, or transfer of learned embeddings to a related task), independent of which model wins on raw accuracy.

Worked example

Suppose one of the 200 features is a "merchant ID" categorical with 5,000 distinct values across the 100k rows. With even distribution across categories, that's $100{,}000 / 5{,}000 = 20$ rows per merchant on average. A tree-based split using target statistics (e.g. the historical target rate for that merchant, computed with appropriate regularization/smoothing to avoid leakage) reduces this to a single informative numeric feature per split candidate, requiring only enough examples per category to estimate a stable rate, not enough to fit a multi-dimensional embedding vector. An embedding of size, say, 16 for that same categorical has 16 free parameters per category to learn, i.e. $5{,}000 \times 16 = 80{,}000$ parameters for that one feature alone, competing for signal against only 20 rows per category; that parameter count relative to available data per category is the concrete reason embeddings are data-hungry in exactly the regime this question describes.

Trade-offs and pitfalls

• GBMs can plateau in accuracy once real interaction structure gets very deep or the categorical semantics matter (e.g. two rarely co-occurring categories that behave similarly); a neural net with shared embeddings can sometimes capture that structure a tree-based split cannot.
• Ensembling a GBM and a neural net sometimes captures complementary error patterns and outperforms either alone, but that's a meaningfully larger engineering and maintenance surface for a modest accuracy gain, worth it only when the accuracy gain is validated and the maintenance cost is acceptable.
• A common wrong turn is reaching for embeddings because "neural nets are more powerful" in the abstract, without checking whether the per-category sample size actually supports learning a useful embedding; at 100k rows with high-cardinality categoricals, that assumption often doesn't hold.
• The right first step in almost every case is to fit the GBM as a fast baseline; if it saturates and the use case justifies the added complexity (larger expected future data volume, need for a differentiable pipeline), only then invest in the neural net path.

Direct answer

Naive softmax or sigmoid computation calls exp() on raw logits, and if any logit is large (a few hundred or more, easily reached with unnormalized network outputs or a diverging training run), exp() overflows to infinity in floating point; when that happens in both the numerator and denominator of softmax, you get an inf / inf = NaN, silently poisoning the loss and every gradient downstream. The log-sum-exp trick fixes this by subtracting the maximum logit before exponentiating, which provably gives the exact same mathematical result while keeping every intermediate value in a safe numeric range.

Structured elaboration

Why overflow happens. IEEE 754 double-precision floats overflow to inf above roughly $1.8\times10^{308}$, but exp(x) reaches that threshold once $x$ is only around 709; single precision (float32, the common default in ML) overflows around $x\approx 88$. Logits in the hundreds or thousands are entirely plausible from an untrained or diverging model, or simply from unnormalized linear-layer outputs before a softmax, so this isn't a contrived edge case.

Why it specifically breaks softmax. Softmax computes $\mathrm{softmax}(z)_k = \dfrac{e^{z_k}}{\sum_j e^{z_j}}$. If the logits are all large (say all near 1000), every $e^{z_j}$ overflows to inf, and the result is inf / inf, which IEEE 754 defines as NaN, not an error you can catch and not a value close to the correct answer. Symmetrically, if the logits are all very negative (say all near -1000), every $e^{z_j}$ underflows to exactly 0.0, and you get 0/0, again NaN.

The log-sum-exp identity that fixes it. For any constant $m$, $\log \sum_j e^{z_j} = m + \log \sum_j e^{z_j - m}$, and picking $m = \max_j z_j$ guarantees every exponent $z_j - m \le 0$, so every exp() call is now bounded between 0 and 1, never overflowing. Applying the same shift inside softmax itself, $\mathrm{softmax}(z)_k = \dfrac{e^{z_k-m}}{\sum_j e^{z_j-m}}$, gives the mathematically identical result (the $e^{-m}$ factor cancels between numerator and denominator) while keeping every intermediate value numerically safe.

A subtler failure mode: log-probabilities, not just softmax itself. Even after stabilizing softmax with the max-subtraction trick, if you then separately compute log(softmax_output) for a cross-entropy loss, a genuinely correct but very small probability (from a large negative logit) can underflow to an exact 0.0, and log(0.0) = -inf, a real, silent failure rather than an imprecise approximation. The fix is the same idea applied one level up: compute $\text{cross-entropy for class } k = \mathrm{logsumexp}(z) - z_k$ directly, never materializing the intermediate probability at all, which is both more accurate and cheaper (one less division and log per element).

Production guardrails. (1) Always call a library's built-in logsumexp / stable softmax (scipy.special.logsumexp, torch.logsumexp, torch.nn.functional.log_softmax) rather than hand-rolling exp and log calls; these are the actual guardrail, not a coding-style preference. (2) Assert or log when raw logits exceed a sanity threshold (say $|z| > 50$) before they reach any exponential, since that's usually itself a symptom of a diverging training run or an unnormalized upstream computation, worth catching early rather than papering over downstream. (3) Prefer computing losses via the identity above (never materializing an intermediate probability that could underflow to exact zero) rather than composing separately-correct softmax and log calls.

Worked example

import numpy as np

logits = np.array([1000.0, 1001.0, 999.0])

def softmax_naive(z):
    e = np.exp(z)
    return e / np.sum(e)

def softmax_stable(z):
    shifted = z - np.max(z)
    e = np.exp(shifted)
    return e / np.sum(e)

def logsumexp_naive(z):
    return np.log(np.sum(np.exp(z)))

def logsumexp_stable(z):
    m = np.max(z)
    return m + np.log(np.sum(np.exp(z - m)))

Running this (float64, verified): softmax_naive(logits) returns [nan, nan, nan], since every exp() call overflowed to inf first, giving inf/inf. softmax_stable(logits) returns [0.2447, 0.6652, 0.0900], which matches scipy.special.softmax(logits) exactly (max absolute difference 0.0). logsumexp_naive(logits) returns inf (the log of an already-overflowed infinity), while logsumexp_stable(logits) correctly returns 1001.408, matching scipy.special.logsumexp(logits) exactly.

The subtler failure is just as concrete and just as reproducible:

def sigmoid_naive(z):
    return 1.0 / (1.0 + np.exp(-z))

def log_sigmoid_stable(z):
    return np.minimum(z, 0) - np.log1p(np.exp(-np.abs(z)))

z = np.array([-1000.0])
s = sigmoid_naive(z)          # [0.0], the correct limiting value
naive_log = np.log(s)         # [-inf], a real divide-by-zero
stable_log = log_sigmoid_stable(z)  # [-1000.0], exact

sigmoid_naive(-1000) correctly underflows to exactly 0.0 (that's the right limiting value, not an error itself), but np.log of that 0.0 then returns -inf, a real RuntimeWarning: divide by zero, not just an imprecise number. log_sigmoid_stable(-1000) instead returns exactly -1000.0, matching the asymptotic identity $\log\sigma(z) \to z$ as $z \to -\infty$ exactly. Both forms agree to machine precision at moderate logit values like $z=-5,0,5$ (verified: naive and stable give the identical [-5.00672, -0.69315, -0.00672] there), so the stable form costs nothing in the well-behaved regime and only pays off exactly where it's needed.

Trade-offs & pitfalls

• The max-subtraction trick is mathematically exact, not an approximation; there's no accuracy-versus-safety trade-off to weigh, which is exactly why every production ML library implements it as the default, not an opt-in flag.
• Don't stabilize softmax alone and assume you're safe; a naive log() applied afterward to the (now-correct) softmax output can still underflow-then-log-of-zero if any resulting probability rounds to exactly 0.0 in floating point. Compute cross-entropy via the direct logsumexp(z) - z_k identity instead of composing separately-stabilized pieces.
• Sanity-checking raw logit magnitude is a cheap, high-value guardrail: a logit of 1000 is very rarely legitimate model output, it's almost always a symptom of an upstream bug (missing normalization, a learning-rate-driven divergence, or accidentally feeding pre-activation values where post-activation ones were expected), worth surfacing immediately rather than silently absorbing.
• These stability tricks matter more, not less, in mixed-precision or float16 training, where the safe range for exp() is far narrower than float64, and where accumulation of many "successfully" masked overflow/underflow events over a training run's many steps is exactly the kind of failure that's hard to trace back to its root cause after the fact.

Direct answer

The L1 penalty $\lambda|\beta|$ has a sharp corner (a "kink") at $\beta = 0$: its left-hand derivative there is $-1$ and its right-hand derivative is $+1$, so no single tangent line, and hence no ordinary gradient, exists at that point. Solvers like coordinate descent sidestep this by optimizing one coefficient at a time, where the one-dimensional L1 subproblem has a closed-form solution called soft-thresholding, so the algorithm never actually needs a gradient at the kink.

Structured elaboration

Why the kink breaks vanilla gradient descent. For a smooth penalty like L2 ($\lambda\beta^2$), the derivative $2\lambda\beta$ is continuous everywhere, including at $\beta=0$, so gradient descent works unmodified. For L1:

These disagree, so $|\beta|$ is not differentiable at exactly $\beta=0$, which is precisely the point L1 is designed to push many coefficients toward.

Coordinate descent's workaround. Instead of computing a joint gradient over all coefficients, coordinate descent fixes every coefficient except one, $\beta_j$, and solves that one-dimensional subproblem exactly. Because the subproblem is one-dimensional, it has a closed-form minimizer even though the objective isn't smooth, no gradient step is needed at all for that coordinate; the algorithm just evaluates the closed form and moves on to the next coordinate.

Deriving the univariate soft-thresholding update. Fix all coefficients except $\beta_j$ (drop the subscript $j$ for readability). Let $r$ be the partial residual (the target minus the prediction from all other fixed coefficients), $x$ the feature column for this coordinate, $s = \|x\|^2$, and $z = x^\top r$. The one-dimensional subproblem (dropping the constant term $\|r\|^2$ that doesn't depend on $\beta$) is:

Consider the two smooth cases separately, since $|\beta|=\beta$ for $\beta>0$ and $|\beta|=-\beta$ for $\beta<0$:

• For $\beta>0$: setting the derivative to zero, $s\beta - z + \lambda = 0 \Rightarrow \beta = (z-\lambda)/s$, valid only if this is actually positive, i.e. $z>\lambda$.
• For $\beta<0$: $s\beta - z - \lambda = 0 \Rightarrow \beta = (z+\lambda)/s$, valid only if $z<-\lambda$.
• If $-\lambda \le z \le \lambda$, neither case is feasible. At a kink like $\beta=0$ there's no single tangent slope, but there's a whole range of slopes between the left-hand and right-hand derivatives (here $[-1,1]$, called the subdifferential) that all count as valid generalized gradients (subgradients) at that point. The optimality condition for a non-smooth function is that zero must be reachable by some value in that range; since $z/\lambda$ falls inside $[-1,1]$ here, that condition holds and confirms $\beta=0$ is the minimizer.

Combining all three cases into one closed form (the soft-thresholding operator):

If the feature column is standardized so that $s=1$, this simplifies to $\beta^\star = S(z,\lambda)$: shrink $z$ toward zero by $\lambda$, and clip to exactly zero if $|z|\le\lambda$.

Worked example

Take $s = 1$ (standardized feature) and $\lambda = 0.3$. For three candidate partial-correlation values $z$:

• $z = 0.5$: $|z| - \lambda = 0.2 > 0$, so $\beta^\star = \operatorname{sign}(0.5)\times 0.2 = 0.2$.
• $z = -0.5$: $|z|-\lambda = 0.2>0$, so $\beta^\star = \operatorname{sign}(-0.5)\times 0.2 = -0.2$.
• $z = 0.2$: $|z|-\lambda = -0.1 < 0$, so the max clips to $0$, giving $\beta^\star = 0$.

That third case is the entire point of the exercise: any coordinate whose signal $z$ is smaller in magnitude than the penalty $\lambda$ gets set to exactly zero, not just shrunk close to zero, which is the mechanism by which L1 produces genuinely sparse solutions rather than merely small coefficients (which is what L2 does instead).

Trade-offs and pitfalls

• Coordinate descent's closed-form update is fast per step but the algorithm cycles through coordinates, so convergence depends on the number of features and how correlated they are; highly correlated features can slow convergence.
• Proximal-gradient methods, ISTA (Iterative Shrinkage-Thresholding Algorithm) and its accelerated version FISTA, are the alternative to coordinate descent: they apply the same soft-thresholding operator to a full gradient step on the smooth part of the objective, useful when a closed-form per-coordinate update isn't available (e.g. more complex smooth losses), at the cost of needing a step size tied to the Lipschitz constant of the gradient (a bound on how fast the gradient itself can change; it sets the largest step size the update can safely take before it overshoots and diverges).
• A common wrong turn is trying to "fix" the non-differentiability by smoothing the L1 penalty (e.g. approximating $|\beta|$ with $\sqrt{\beta^2+\epsilon}$) so that plain gradient descent works; this sacrifices the exact-zero property that makes L1 useful for sparsity and feature selection in the first place, defeating the purpose of choosing L1 over L2.
• The derivation above assumes the feature column's scale $s=\|x\|^2$ is factored in correctly; skipping standardization and using a shared $\lambda$ across differently-scaled features silently penalizes large-scale features more, a common and easy-to-miss bug in from-scratch coordinate descent implementations.