Gradient Descent for Linear Regression in C++
The direct solvers in the OLS and QR/SVD posts process the entire training set in one factorization. That is ideal when $n$ and $d$ are modest, but the $O(nd^2)$ factorization cost can dominate for large datasets. Gradient descent replaces that factorization with repeated $O(nd)$ passes and updates the parameter vector one batch at a time. The batch implementation below still stores an $n \times d$ standardized training copy, while each update needs only an $O(d)$ gradient buffer; the online method needs only $O(d)$ state per arriving sample. This post implements batch, mini-batch, stochastic (SGD), and online variants in standard-library C++. A train-only standardization layer is applied before any gradient step, and frozen constants convert the optimized weights back to the same raw-space fit / predict / evaluate contract used by every post in this sub-series. The final program fits the same synthetic two-feature dataset used in Parts 1 and 2 and reaches the same coefficients as the direct solvers.
The Normalized MSE Gradient
The training objective is the mean squared error over $n$ samples:
$$J(\beta) = \frac{1}{n} \sum_{i=1}^{n} \bigl(y_i - \hat y_i\bigr)^2, \qquad \hat y_i = \beta_0 + \sum_{j=1}^{d} \beta_j x_{ij}.$$Differentiating with respect to the intercept $\beta_0$ and each coefficient $\beta_j$:
$$\frac{\partial J}{\partial \beta_0} = \frac{-2}{n} \sum_{i=1}^{n} r_i, \qquad \frac{\partial J}{\partial \beta_j} = \frac{-2}{n} \sum_{i=1}^{n} r_i\, x_{ij}, \qquad r_i = y_i - \hat y_i.$$The gradient update rule subtracts a step in the gradient direction:
$$\beta_0 \leftarrow \beta_0 - \alpha \cdot \frac{\partial J}{\partial \beta_0}, \qquad \beta_j \leftarrow \beta_j - \alpha \cdot \frac{\partial J}{\partial \beta_j},$$where $\alpha > 0$ is the learning rate. Because $J$ is a convex quadratic in $\beta$, gradient descent converges to the global minimum for any sufficiently small $\alpha$. The division by $n$ (the normalization) is important for learning-rate stability: without it, the effective step size scales with the training set size, and a learning rate that works for $n = 100$ fails for $n = 10\,000$.
For a mini-batch of size $B$ rows $\{i_1, \dots, i_B\}$, the gradients above are estimated with the sum over those $B$ rows divided by $B$ instead of $n$. This is a consistent estimator of the true gradient; the variance decreases as $B$ increases toward $n$.
Train-Only Standardization
The convergence speed of gradient descent is sensitive to the scale of the features. If one feature ranges from $-10^3$ to $10^3$ while another ranges from $0$ to $1$, the loss surface is a very elongated ellipsoid and the gradient points mostly along the wide axis, making progress along the narrow axis slow. Standardizing features to have zero mean and unit variance converts the ellipsoid into a shape closer to a sphere and makes a single learning rate work reasonably well across all directions.
The standardization statistics are computed from the training set only. Let $\mu_j$ and $\sigma_j$ denote the empirical mean and standard deviation (population formula, dividing by $n$) of the $j$-th feature column over the training rows, and let $\mu_y$ and $\sigma_y$ denote the same quantities for the training targets. The standardized inputs and outputs are
$$\tilde x_{ij} = \frac{x_{ij} - \mu_j}{\sigma_j}, \qquad \tilde y_i = \frac{y_i - \mu_y}{\sigma_y}.$$Gradient descent then minimizes MSE in the standardized space, finding standardized weights $\tilde\beta_0$ and $\tilde\beta_1, \dots, \tilde\beta_d$. Crucially, $\mu_j$, $\sigma_j$, $\mu_y$, and $\sigma_y$ are computed once, stored as members of the regressor struct, and never updated again after fit returns. Any test point or online update uses the same frozen constants, preventing information leakage from test data into the normalization step.
A constant feature has $\sigma_j = 0$, so direct division is undefined. The implementation handles that case explicitly by assigning the stored scale $\sigma_j = 1$. Centering still maps every training value in that column to zero, its gradient remains zero from the zero initialization, and its raw coefficient remains zero. The intercept therefore absorbs the constant column without introducing a duplicate trainable direction. Small but nonconstant columns keep their measured scale rather than being mistaken for constants. The standard deviation calculation divides deviations by their maximum magnitude before squaring, which avoids overflow and underflow for unusually large or small finite scales. A constant target is different: no slope optimization is needed and standardized target MSE has no meaningful scale, so fit rejects an exactly zero target standard deviation.
Converting Standardized Weights to Raw Parameters
After optimization, the standardized weights $\tilde\beta_j$ represent the relationship between standardized inputs and standardized outputs. To produce raw-space predictions, the final weights must be converted back. For each feature coefficient:
$$\beta_j = \tilde\beta_j \cdot \frac{\sigma_y}{\sigma_j}, \qquad j = 1, \dots, d.$$The raw intercept follows from expanding the standardized model back into raw space. Substituting $\tilde x_j = (x_j - \mu_j)/\sigma_j$ and $\tilde y = (y - \mu_y)/\sigma_y$ into $\tilde y = \tilde\beta_0 + \sum_j \tilde\beta_j \tilde x_j$ and solving for $y$:
$$\beta_0 = \mu_y + \sigma_y\,\tilde\beta_0 - \sum_{j=1}^{d} \beta_j\, \mu_j.$$The $\sigma_y\,\tilde\beta_0$ term carries the contribution of the optimized standardized intercept. At the global minimum of MSE on zero-mean standardized data, the gradient forces $\tilde\beta_0 \to 0$ at convergence. However, during training after any finite number of gradient steps, $\tilde\beta_0$ is not exactly zero, and dropping the term would corrupt the converted intercept. The formula guarantees that the raw model $\hat y = \beta_0 + \sum_j \beta_j x_j$ is algebraically identical to the standardized model for any input at any point in training, not only at convergence.
In code the conversion is encapsulated in a private syncRawParams() method called when batch training completes. The online path computes the same conversion into temporary values after every update, validates them, and then commits the synchronized standardized and raw state together:
void syncRawParams() {
for (int j = 0; j < nFeatures; ++j)
coef[j] = wTilde[j] * yStd / xStd[j];
intercept = yMean + yStd * w0Tilde;
for (int j = 0; j < nFeatures; ++j)
intercept -= coef[j] * xMean[j];
}Raw to Standardized and Back
A round trip on one concrete point confirms that the standardized model and its converted raw form predict identically. The batch fit in the complete program freezes these training-set constants: $\mu = (0.007450, 0.114314)$, $\sigma = (1.628750, 1.778122)$, $\mu_y = 2.805188$, $\sigma_y = 4.061318$, with converged standardized weights $\tilde\beta_0 = 0$ and $\tilde\beta = (0.793143, -0.664297)$. The raw query point for this check is $x = (1.5, -2.0)$.
- Raw to standardized. $\tilde x_1 = (1.5 - 0.007450)/1.628750 = 0.916378$ and $\tilde x_2 = (-2.0 - 0.114314)/1.778122 = -1.189071$.
- Predict in standardized space. $\tilde y = \tilde\beta_0 + \tilde\beta_1 \tilde x_1 + \tilde\beta_2 \tilde x_2 = 0 + (0.793143)(0.916378) + (-0.664297)(-1.189071) = 1.516715$.
- Standardized back to raw. $\hat y = \mu_y + \sigma_y\,\tilde y = 2.805188 + 4.061318 \times 1.516715 = 8.965050$.
- Convert the weights, then predict directly. The coefficients are $\beta_1 = \tilde\beta_1 \sigma_y/\sigma_1 = 1.977717$ and $\beta_2 = \tilde\beta_2 \sigma_y/\sigma_2 = -1.517287$, and the intercept is $\beta_0 = \mu_y + \sigma_y\tilde\beta_0 - \beta_1\mu_1 - \beta_2\mu_2 = 2.963902$. Predicting in raw space gives $\hat y = 2.963902 + (1.977717)(1.5) + (-1.517287)(-2.0) = 8.965050$.
Both routes land on $8.965050$, differing only at the level of floating-point rounding. The converted $\beta_0 = 2.963902$ and $\beta = (1.977717, -1.517287)$ are exactly the values the complete program prints for the batch fit, which is the guarantee that the raw predict path never needs the standardization constants at inference time.
Batch, Mini-Batch, and SGD
The batch size $B$ controls the trade-off between gradient quality and update frequency:
| Variant | $B$ | Updates per epoch | Gradient noise | Typical use |
|---|---|---|---|---|
| Batch GD | $n$ | 1 | None (exact gradient) | Small datasets; stable convergence |
| Mini-batch GD | $16{-}256$ | $\lceil n/B \rceil$ | Low to moderate | Most practical cases; GPU-friendly |
| SGD | $1$ | $n$ | High | Online-like updates; slow wall-clock convergence |
For batch GD the rows are processed in the same order every epoch, so no shuffle is needed. For mini-batch and SGD the rows are shuffled once per epoch using a deterministic RNG seeded from the user-supplied GDConfig::seed. Deterministic shuffling ensures reproducibility: the same seed always produces the same sequence of updates and therefore the same converged parameters (to floating-point precision).
One subtlety about batch semantics: the gradient for a batch of size $B$ is divided by $B$, not by $n$. This makes the random mini-batch average an unbiased estimate of the full-data mean gradient, so the learning rate has the same scale across batch sizes modulo gradient noise. Dividing the mini-batch sum by $n$ instead would shrink its expectation by $B/n$ and would require a compensating batch-size-dependent learning rate.
One Update by Hand
A single mini-batch step, worked out on paper, shows exactly where the $2/B$ factor lands. The example mini-batch contains $B = 2$ already-standardized rows, so the intercept participates alongside the two feature weights. Row A has $\tilde x = (0.5, -1.0)$ and $\tilde y = 0.4$; row B has $\tilde x = (-0.5, 1.0)$ and $\tilde y = -0.6$. The weights start at the zero vector $(\tilde\beta_0, \tilde\beta_1, \tilde\beta_2) = (0, 0, 0)$ with learning rate $\alpha = 0.1$.
- Residuals. Both predictions start at zero, so $r_A = \hat y_A - \tilde y_A = 0 - 0.4 = -0.4$ and $r_B = 0 - (-0.6) = 0.6$. The residual is prediction minus target, matching the sign used in the code.
- Gradient sums. The intercept gradient sum is $g_0 = r_A + r_B = 0.2$. The feature gradients are $g_1 = r_A \tilde x_{A1} + r_B \tilde x_{B1} = (-0.4)(0.5) + (0.6)(-0.5) = -0.5$ and $g_2 = r_A \tilde x_{A2} + r_B \tilde x_{B2} = (-0.4)(-1.0) + (0.6)(1.0) = 1.0$.
- Scale by $2/B$. With $B = 2$ the factor is $2/B = 1.0$, so each raw sum is the actual gradient component here. For a larger batch the same sums would be divided by that larger $B$.
- One step. Subtract $\alpha \cdot (2/B) \cdot g$: $\tilde\beta_0 = 0 - 0.1(1.0)(0.2) = -0.02$, $\tilde\beta_1 = 0 - 0.1(1.0)(-0.5) = 0.05$, $\tilde\beta_2 = 0 - 0.1(1.0)(1.0) = -0.10$.
- Recomputed loss. Before the step the batch MSE is $\tfrac{1}{2}(0.4^2 + 0.6^2) = 0.26$. After the step the predictions become $0.105$ and $-0.145$, the residuals $-0.295$ and $0.455$, and the batch MSE is $\tfrac{1}{2}(0.295^2 + 0.455^2) \approx 0.1470$. One step nearly halves the batch loss, confirming the update moves downhill.
The animation below scales the same mechanics up to a full dataset and contrasts full-batch, two-row mini-batch, and per-sample paths in the standardized weight plane. Every plotted point is computed from the stated rows and the same $2/B$ gradient used by the C++ implementation.
▶ Batch, Mini-Batch, and SGD Paths on the Standardized Loss Contour
Dataset (raw): $x_1 = [1,2,3,4,5]$, $x_2 = [3,1,4,2,5]$, $y = [1,2,5,4,8]$, standardized to zero mean and unit population variance. With zero-mean data the intercept stays at zero, so this view fixes $\tilde\beta_0 = 0$ and descends over $(\tilde\beta_1, \tilde\beta_2)$. The exact objective is $J(w)=1-2q^\top w+w^\top Mw$, where $q=(0.9238,0.7506)$ and $M=\left[\begin{smallmatrix}1&0.5\\0.5&1\end{smallmatrix}\right]$. Full batch uses lr $0.15$; mini-batch uses lr $0.08$, $B=2$, and slices $(0,1),(2,3),(4)$; SGD uses lr $0.10$. Both noisy variants visit rows in the stated order each epoch. Each arrow marks the latest parameter update computed from those rows and the stated $2/B$ gradient.
Learning-Rate Decay
A fixed learning rate that is large enough to make rapid progress early in training may prevent tight convergence late in training, where the iterates are close to the minimum and need smaller steps. Step decay reduces the learning rate as a function of the epoch number:
$$\alpha_t = \frac{\alpha_0}{1 + \lambda t},$$where $\alpha_0$ is the initial learning rate and $\lambda \ge 0$ is the decay coefficient. Setting $\lambda = 0$ recovers a constant rate. With $\lambda > 0$ the rate decreases monotonically, falling to $\alpha_0 / 2$ at epoch $1/\lambda$. In practice the decay coefficient is tuned alongside the initial rate; the implementation uses $\lambda = 0.005$ in the examples below, which halves the rate around epoch $200$ and provides a clean convergence curve on the $240$-sample synthetic dataset.
Convergence and Divergence Checks
Convergence is declared when the relative change in epoch loss falls below a tolerance:
$$\frac{|J_t - J_{t-1}|}{|J_{t-1}| + 10^{-8}} < \tau,$$where $\tau$ is the user-supplied GDConfig::tol. The $10^{-8}$ in the denominator prevents division by zero when the loss reaches machine zero, which can happen for noise-free synthetic data. After convergence the loop exits early; the number of epochs actually run is stored in epochsRun for inspection.
Divergence (numerically unstable updates that drive the loss to infinity or NaN) is detected by std::isfinite checks on the epoch loss and final parameters. For normalized MSE, the gradient has Lipschitz constant $L = 2\sigma_{\max}^2(\tilde X_b)/n$, so plain batch gradient descent is stable for $0 < \alpha < 2/L = n/\sigma_{\max}^2(\tilde X_b)$. If a check fails, a std::runtime_error is thrown immediately with a message to reduce the learning rate. Common causes are a learning rate above that stability interval or inputs whose finite values overflow while standardization statistics are accumulated. Because the implementation standardizes before updating, an extremely large learning rate is the most likely trigger in practice.
if (!isfinite(loss))
throw runtime_error(
"GDRegressor::fit: loss became non-finite at epoch " + to_string(epoch) +
"; reduce learning rate or check input scale");The animation below makes the convergence criterion concrete. It traces the standardized MSE loss for the batch configuration (lr = 0.05, decay = 0.005) on the 192-row training split, stepping through measured checkpoints from the deterministic run until the relative-change test halts training at epoch 161. The side panel shows the raw-space parameter estimates that the standardized weights convert to at each checkpoint, including the standardized-intercept correction.
▶ Batch Gradient Descent Convergence (lr=0.05, decay=0.005)
MSE loss in standardized space vs epoch. Starting from zero weights, the loss falls from 1.000 (the mean-predictor baseline) toward the optimal ~0.022 as the parameters converge.
Online Updates
fitOnline applies a single gradient step to one new observation after the initial batch fit. The frozen standardization constants from training are used to standardize the new point before computing the gradient: the new $x$ is mapped to $\tilde x_j = (x_j - \mu_j) / \sigma_j$ and the new $y$ to $\tilde y = (y - \mu_y) / \sigma_y$ using the training-time statistics stored in xMean, xStd, yMean, and yStd. The gradient is then the single-sample gradient in standardized space. Candidate standardized weights and their raw-space conversion are checked for finiteness before both representations are committed, so the public parameters remain synchronized after every online update.
Critically, fitOnline does not recompute the standardization statistics from the new point. Re-estimating $\mu_j$ and $\sigma_j$ per arrival would require streaming algorithms (Welford updates) and would change the meaning of all previously stored raw parameters, since old predictions $\hat y = \beta_0 + \sum_j \beta_j x_j$ would become inconsistent with the new $\mu_j$ and $\sigma_j$. The frozen-standardizer approach avoids this: the new point is transformed into the same standardized coordinate system that the batch training used, the weights move a small step in that fixed system, and the raw parameters are updated by the same algebraic conversion as always. The practical consequence is that fitOnline should be used only when the new data comes from the same distribution as the training data; a distribution shift that changes the true $\mu_j$ or $\sigma_j$ calls for a full retrain.
GDRegressor Implementation
struct GDConfig {
int batchSize = -1; // -1: full batch; 1: SGD; B: mini-batch
double lr = 0.01; // initial learning rate
double lrDecay = 0.0; // step-decay: lr_t = lr / (1 + lrDecay * epoch)
int maxEpochs = 1000;
double tol = 1e-7; // relative-change convergence tolerance
int seed = 0; // RNG seed for deterministic shuffle
};
struct GDRegressor {
bool fitted = false;
int nFeatures = 0;
double intercept = 0.0;
vector<double> coef;
// Frozen training-set standardization constants.
vector<double> xMean, xStd;
double yMean = 0.0, yStd = 1.0;
// Diagnostic fields.
vector<double> lossHistory; // per-epoch MSE in standardized space
int epochsRun = 0;
private:
vector<double> wTilde; // standardized feature weights (length d)
double w0Tilde = 0.0; // standardized intercept
void syncRawParams() {
for (int j = 0; j < nFeatures; ++j)
coef[j] = wTilde[j] * yStd / xStd[j];
intercept = yMean + yStd * w0Tilde;
for (int j = 0; j < nFeatures; ++j)
intercept -= coef[j] * xMean[j];
}
public:
void fit(const Matrix& X, const vector<double>& y, GDConfig cfg = GDConfig{}) {
fitted = false;
// ... shape, finiteness, and configuration validation ...
int n = (int)X.size(), d = (int)X[0].size();
nFeatures = d;
coef.assign(d, 0.0);
// Compute train-only means and standard deviations.
xMean.assign(d, 0.0); xStd.assign(d, 0.0); yMean = 0.0; yStd = 0.0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < d; ++j) xMean[j] += X[i][j];
yMean += y[i];
}
for (int j = 0; j < d; ++j) xMean[j] /= n;
yMean /= n;
// Scale deviations before squaring to avoid overflow and underflow.
vector<double> xMaxDev(d, 0.0);
double yMaxDev = 0.0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < d; ++j)
xMaxDev[j] = max(xMaxDev[j], fabs(X[i][j] - xMean[j]));
yMaxDev = max(yMaxDev, fabs(y[i] - yMean));
}
for (int j = 0; j < d; ++j) {
if (xMaxDev[j] == 0.0) { xStd[j] = 1.0; continue; }
double scaled = 0.0;
for (int i = 0; i < n; ++i) {
double z = (X[i][j] - xMean[j]) / xMaxDev[j];
scaled += z * z;
}
xStd[j] = xMaxDev[j] * sqrt(scaled / n);
}
if (yMaxDev == 0.0)
throw invalid_argument("GDRegressor::fit: target is constant");
double yScaled = 0.0;
for (int i = 0; i < n; ++i) {
double z = (y[i] - yMean) / yMaxDev;
yScaled += z * z;
}
yStd = yMaxDev * sqrt(yScaled / n);
// Build standardized copies.
Matrix Xs(n, vector<double>(d)); vector<double> ys(n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < d; ++j) Xs[i][j] = (X[i][j]-xMean[j])/xStd[j];
ys[i] = (y[i] - yMean) / yStd;
}
wTilde.assign(d, 0.0); w0Tilde = 0.0;
lossHistory.clear(); epochsRun = 0;
int bsz = cfg.batchSize == -1 ? n : cfg.batchSize;
mt19937 rng((unsigned)cfg.seed);
vector<int> idx(n); iota(idx.begin(), idx.end(), 0);
double prevLoss = numeric_limits<double>::infinity();
for (int epoch = 0; epoch < cfg.maxEpochs; ++epoch) {
if (bsz < n) shuffle(idx.begin(), idx.end(), rng);
double decayDenom = 1.0 + cfg.lrDecay * epoch;
if (!isfinite(decayDenom))
throw runtime_error(
"GDRegressor::fit: learning-rate decay overflowed at epoch " +
to_string(epoch));
double lrE = cfg.lr / decayDenom;
for (int start = 0; start < n; start += bsz) {
int end = min(start + bsz, n), B = end - start;
double gw0 = 0.0; vector<double> gw(d, 0.0);
for (int ii = start; ii < end; ++ii) {
int i = (bsz < n) ? idx[ii] : ii;
double pred = w0Tilde;
for (int j = 0; j < d; ++j) pred += wTilde[j] * Xs[i][j];
double err = pred - ys[i]; // residual = prediction - target
gw0 += err;
for (int j = 0; j < d; ++j) gw[j] += err * Xs[i][j];
}
double scale = 2.0 / B;
w0Tilde -= lrE * scale * gw0;
for (int j = 0; j < d; ++j) wTilde[j] -= lrE * scale * gw[j];
}
// Epoch loss on full training set.
double loss = 0.0;
for (int i = 0; i < n; ++i) {
double pred = w0Tilde;
for (int j = 0; j < d; ++j) pred += wTilde[j] * Xs[i][j];
double e = pred - ys[i]; loss += e*e;
}
loss /= n;
if (!isfinite(loss))
throw runtime_error(
"GDRegressor::fit: loss became non-finite at epoch " + to_string(epoch) +
"; reduce learning rate or check input scale");
lossHistory.push_back(loss); ++epochsRun;
double rel = fabs(prevLoss - loss) / (fabs(prevLoss) + 1e-8);
if (rel < cfg.tol) break;
prevLoss = loss;
}
syncRawParams();
fitted = true;
}
void fitOnline(const vector<double>& x, double yVal, double lr = 0.01) {
if (!fitted) throw logic_error("GDRegressor::fitOnline: call fit() before fitOnline()");
// ... validate x size, values, target, and positive finite learning rate ...
// Standardize with frozen training constants.
vector<double> xs(nFeatures);
for (int j = 0; j < nFeatures; ++j) xs[j] = (x[j] - xMean[j]) / xStd[j];
double ys = (yVal - yMean) / yStd;
double pred = w0Tilde;
for (int j = 0; j < nFeatures; ++j) pred += wTilde[j] * xs[j];
double err = pred - ys;
double nextW0 = w0Tilde - lr * 2.0 * err;
vector<double> nextW = wTilde;
vector<double> nextCoef(nFeatures);
for (int j = 0; j < nFeatures; ++j) {
nextW[j] -= lr * 2.0 * err * xs[j];
nextCoef[j] = nextW[j] * yStd / xStd[j];
if (!isfinite(nextW[j]) || !isfinite(nextCoef[j]))
throw runtime_error(
"GDRegressor::fitOnline: update became non-finite at index " +
to_string(j));
}
double nextIntercept = yMean + yStd * nextW0;
for (int j = 0; j < nFeatures; ++j)
nextIntercept -= nextCoef[j] * xMean[j];
if (!isfinite(nextW0) || !isfinite(nextIntercept))
throw runtime_error(
"GDRegressor::fitOnline: intercept update became non-finite");
w0Tilde = nextW0;
wTilde = nextW;
coef = nextCoef;
intercept = nextIntercept;
}
double predict(const vector<double>& x) const {
if (!fitted) throw logic_error("GDRegressor::predict: model is not fitted");
if ((int)x.size() != nFeatures)
throw invalid_argument("GDRegressor::predict: expected " +
to_string(nFeatures) + " features, got " +
to_string(x.size()));
double out = intercept;
for (int j = 0; j < nFeatures; ++j) {
if (!isfinite(x[j]))
throw invalid_argument(
"GDRegressor::predict: non-finite feature at index " +
to_string(j));
out += coef[j] * x[j];
}
if (!isfinite(out))
throw runtime_error("GDRegressor::predict: prediction became non-finite");
return out;
}
vector<double> predict(const Matrix& X) const {
vector<double> out; out.reserve(X.size());
for (const auto& row : X) out.push_back(predict(row));
return out;
}
};Inside One Epoch
The order of operations inside the epoch loop is deliberate, and each step depends on the one before it. Reading the loop body top to bottom:
- Deterministic shuffle. When $B < n$ (mini-batch or SGD), the index vector is shuffled once at the start of the epoch with an RNG seeded from
GDConfig::seed. Full batch skips the shuffle because processing order does not change an exact gradient. Seeding from the config keeps every run reproducible. - Learning-rate schedule. The decayed rate $\alpha_t = \alpha_0/(1 + \lambda t)$ is computed once per epoch from the current epoch index, and the denominator is checked for overflow before use.
- Batch slices. The row range is walked in strides of $B$; the final slice may be shorter than $B$, so the code uses the actual slice length when scaling the gradient.
- Standardized update. For each slice the residuals are accumulated in standardized space, scaled by $2/B$ (using the slice length, not $n$), and subtracted from the standardized weights. The intercept and feature weights update together in the same pass.
- Full-training loss. After all slices, one extra pass over every training row computes the standardized MSE for the whole epoch. This is a measurement pass, not an update, so it always uses the complete training set regardless of batch size.
- Convergence and history. The epoch loss is pushed onto
lossHistory,epochsRunis incremented, and the relative-change test decides whether to stop. The loss is recorded before the early-exit check so the history always ends with the loss that triggered convergence.
One caveat follows from measuring the full-training loss after a shuffled, partially updated epoch: for mini-batch and SGD the per-epoch loss is not monotone. Each mini-batch step moves the weights along a noisy gradient estimate, so the whole-set loss can tick upward between epochs even while the trend falls. The relative-change convergence test can therefore trip early on a lucky flat pair of epochs or run to the epoch budget when the noise floor keeps the change above tolerance, which is exactly why the mini-batch and SGD runs in the complete program exhaust more epochs than full batch. The raw-space parameters are synchronized only once, after the loop exits, so intermediate raw coefficients are never observed mid-epoch during batch training.
Validation and Errors
The input validation in fit follows the same conventions as the rest of the series:
- Row count mismatch between
Xandy, ragged rows, empty input, or non-finite entries throwstd::invalid_argument. - Underdetermined data ($n < d{+}1$) is accepted. Zero initialization gives batch gradient descent an implicit minimum-norm bias in standardized coordinates, but it does not provide the raw-design rank diagnostics or controlled pseudoinverse tolerance of SVD.
- An exactly constant feature receives a scale of $1$, centers to zero, and keeps coefficient zero. Small nonconstant features retain their scale through a max-deviation variance calculation. An exactly constant target throws
std::invalid_argumentbecause standardized target MSE has no meaningful scale. - The learning rate and tolerance must be finite and positive, decay must be finite and non-negative, epochs must be positive, and batch size must be
-1or between $1$ and $n$. - A non-finite epoch loss throws
std::runtime_errorwith the epoch index; the message directs the caller to reduce the learning rate. fitOnlinethrowsstd::logic_errorwhen called beforefit, rejects invalid values or learning rates, and validates candidate parameters before committing the update.- Calling
predictbeforefitthrowsstd::logic_error; a feature-count mismatch or non-finite feature throwsstd::invalid_argument.
Complete Program and Output
The program below uses the same $240$-point two-feature synthetic dataset as the OLS and QR/SVD posts (seed $7$, $80/20$ split), then runs three gradient-descent configurations: full-batch with 500-epoch budget, mini-batch of size 32, and SGD. A fourth block demonstrates fitOnline by applying 20 single-sample updates from the test set. A fifth block shows the divergence detector catching a learning rate of $50$.
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;
using Matrix = vector<vector<double>>;
// ============================================================
// Metric helpers.
// ============================================================
static void validateMetricInputs(const vector<double>& yTrue,
const vector<double>& yPred,
const string& name) {
if (yTrue.size() != yPred.size())
throw invalid_argument(name + ": size mismatch");
if (yTrue.empty())
throw invalid_argument(name + ": vectors must be non-empty");
for (size_t i = 0; i < yTrue.size(); ++i)
if (!isfinite(yTrue[i]) || !isfinite(yPred[i]))
throw invalid_argument(name + ": non-finite value at index " + to_string(i));
}
static double mse(const vector<double>& yTrue, const vector<double>& yPred) {
validateMetricInputs(yTrue, yPred, "mse");
double s = 0.0;
for (size_t i = 0; i < yTrue.size(); ++i) { double e = yTrue[i]-yPred[i]; s += e*e; }
return s / yTrue.size();
}
static double rmse(const vector<double>& yTrue, const vector<double>& yPred) {
return sqrt(mse(yTrue, yPred));
}
static double r2(const vector<double>& yTrue, const vector<double>& yPred) {
validateMetricInputs(yTrue, yPred, "r2");
double meanY = accumulate(yTrue.begin(), yTrue.end(), 0.0) / yTrue.size();
double ssRes = 0.0, ssTot = 0.0;
for (size_t i = 0; i < yTrue.size(); ++i) {
ssRes += (yTrue[i]-yPred[i])*(yTrue[i]-yPred[i]);
ssTot += (yTrue[i]-meanY)*(yTrue[i]-meanY);
}
if (ssTot == 0.0)
throw invalid_argument("r2: target values are constant, so R^2 is undefined");
return 1.0 - ssRes / ssTot;
}
// ============================================================
// GDConfig: hyperparameters for gradient descent.
// ============================================================
struct GDConfig {
int batchSize = -1; // -1: full batch; 1: SGD; B: mini-batch
double lr = 0.01; // initial learning rate
double lrDecay = 0.0; // step-decay: lr_t = lr / (1 + lrDecay * epoch)
int maxEpochs = 1000;
double tol = 1e-7; // relative-change convergence tolerance
int seed = 0; // RNG seed for deterministic shuffle
};
// ============================================================
// GDRegressor: gradient descent for linear regression.
// ============================================================
struct GDRegressor {
bool fitted = false;
int nFeatures = 0;
double intercept = 0.0;
vector<double> coef;
// Frozen training-set standardization constants.
vector<double> xMean, xStd;
double yMean = 0.0, yStd = 1.0;
// Diagnostic fields.
vector<double> lossHistory; // per-epoch MSE in standardized space
int epochsRun = 0;
private:
vector<double> wTilde; // standardized feature weights (length d)
double w0Tilde = 0.0; // standardized intercept
// Convert standardized weights to raw-space parameters.
// Derivation: the standardized model predicts
// y = yMean + yStd * (w0Tilde + sum_j wTilde[j] * (x_j - xMean[j]) / xStd[j])
// Rearranging to the raw form y = intercept + sum_j coef[j] * x_j gives:
// coef[j] = wTilde[j] * yStd / xStd[j]
// intercept = yMean + yStd * w0Tilde - sum_j coef[j] * xMean[j]
void syncRawParams() {
for (int j = 0; j < nFeatures; ++j)
coef[j] = wTilde[j] * yStd / xStd[j];
intercept = yMean + yStd * w0Tilde;
for (int j = 0; j < nFeatures; ++j)
intercept -= coef[j] * xMean[j];
}
public:
void fit(const Matrix& X, const vector<double>& y, GDConfig cfg = GDConfig{}) {
fitted = false;
// Shape and finiteness validation.
if (X.size() != y.size())
throw invalid_argument("GDRegressor::fit: X has " + to_string(X.size()) +
" rows but y has " + to_string(y.size()) + " entries");
int n = (int)X.size();
if (n == 0) throw invalid_argument("GDRegressor::fit: empty training set");
int d = (int)X[0].size();
if (d == 0) throw invalid_argument("GDRegressor::fit: zero feature columns");
for (int i = 0; i < n; ++i) {
if ((int)X[i].size() != d)
throw invalid_argument("GDRegressor::fit: row " + to_string(i) +
" has " + to_string(X[i].size()) +
" columns, expected " + to_string(d));
for (int j = 0; j < d; ++j)
if (!isfinite(X[i][j]))
throw invalid_argument("GDRegressor::fit: non-finite feature in row " +
to_string(i));
if (!isfinite(y[i]))
throw invalid_argument("GDRegressor::fit: non-finite target at row " +
to_string(i));
}
if (!isfinite(cfg.lr) || cfg.lr <= 0.0)
throw invalid_argument("GDRegressor::fit: lr must be finite and positive");
if (!isfinite(cfg.lrDecay) || cfg.lrDecay < 0.0)
throw invalid_argument("GDRegressor::fit: lrDecay must be finite and non-negative");
if (!isfinite(cfg.tol) || cfg.tol <= 0.0)
throw invalid_argument("GDRegressor::fit: tol must be finite and positive");
if (cfg.maxEpochs < 1)
throw invalid_argument("GDRegressor::fit: maxEpochs must be >= 1");
if (cfg.batchSize != -1 && (cfg.batchSize < 1 || cfg.batchSize > n))
throw invalid_argument(
"GDRegressor::fit: batchSize must be -1 or between 1 and the row count");
nFeatures = d;
coef.assign(d, 0.0);
// Compute train-only means and standard deviations (population formula).
xMean.assign(d, 0.0); xStd.assign(d, 0.0); yMean = 0.0; yStd = 0.0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < d; ++j) xMean[j] += X[i][j];
yMean += y[i];
}
for (int j = 0; j < d; ++j) xMean[j] /= n;
yMean /= n;
for (int j = 0; j < d; ++j)
if (!isfinite(xMean[j]))
throw runtime_error(
"GDRegressor::fit: feature mean became non-finite; rescale the input");
if (!isfinite(yMean))
throw runtime_error(
"GDRegressor::fit: target mean became non-finite; rescale the target");
// Scale deviations before squaring to avoid overflow and underflow.
vector<double> xMaxDev(d, 0.0);
double yMaxDev = 0.0;
for (int i = 0; i < n; ++i) {
for (int j = 0; j < d; ++j) {
double dx = fabs(X[i][j] - xMean[j]);
if (!isfinite(dx))
throw runtime_error(
"GDRegressor::fit: feature deviation became non-finite; "
"rescale the input");
xMaxDev[j] = max(xMaxDev[j], dx);
}
double dy = fabs(y[i] - yMean);
if (!isfinite(dy))
throw runtime_error(
"GDRegressor::fit: target deviation became non-finite; "
"rescale the target");
yMaxDev = max(yMaxDev, dy);
}
for (int j = 0; j < d; ++j) {
if (xMaxDev[j] == 0.0) {
xStd[j] = 1.0; // constant column centers to all zeros
continue;
}
double scaled = 0.0;
for (int i = 0; i < n; ++i) {
double z = (X[i][j] - xMean[j]) / xMaxDev[j];
scaled += z * z;
}
xStd[j] = xMaxDev[j] * sqrt(scaled / n);
if (!isfinite(xStd[j]) || xStd[j] == 0.0)
throw runtime_error(
"GDRegressor::fit: feature scale became invalid at index " +
to_string(j));
}
if (yMaxDev == 0.0)
throw invalid_argument("GDRegressor::fit: target is constant");
double yScaled = 0.0;
for (int i = 0; i < n; ++i) {
double z = (y[i] - yMean) / yMaxDev;
yScaled += z * z;
}
yStd = yMaxDev * sqrt(yScaled / n);
if (!isfinite(yStd) || yStd == 0.0)
throw runtime_error("GDRegressor::fit: target scale became invalid");
// Build standardized copies.
Matrix Xs(n, vector<double>(d)); vector<double> ys(n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < d; ++j) {
Xs[i][j] = (X[i][j] - xMean[j]) / xStd[j];
if (!isfinite(Xs[i][j]))
throw runtime_error(
"GDRegressor::fit: standardized feature became non-finite");
}
ys[i] = (y[i] - yMean) / yStd;
if (!isfinite(ys[i]))
throw runtime_error("GDRegressor::fit: standardized target became non-finite");
}
wTilde.assign(d, 0.0); w0Tilde = 0.0;
lossHistory.clear(); epochsRun = 0;
int bsz = cfg.batchSize == -1 ? n : cfg.batchSize;
mt19937 rng((unsigned)cfg.seed);
vector<int> idx(n); iota(idx.begin(), idx.end(), 0);
double prevLoss = numeric_limits<double>::infinity();
for (int epoch = 0; epoch < cfg.maxEpochs; ++epoch) {
// Shuffle row indices once per epoch for mini-batch / SGD.
if (bsz < n) shuffle(idx.begin(), idx.end(), rng);
double decayDenom = 1.0 + cfg.lrDecay * epoch;
if (!isfinite(decayDenom))
throw runtime_error(
"GDRegressor::fit: learning-rate decay overflowed at epoch " +
to_string(epoch));
double lrE = cfg.lr / decayDenom;
for (int start = 0; start < n; start += bsz) {
int end = min(start + bsz, n), B = end - start;
double gw0 = 0.0; vector<double> gw(d, 0.0);
for (int ii = start; ii < end; ++ii) {
int i = (bsz < n) ? idx[ii] : ii;
double pred = w0Tilde;
for (int j = 0; j < d; ++j) pred += wTilde[j] * Xs[i][j];
double err = pred - ys[i]; // residual = prediction - target
gw0 += err;
for (int j = 0; j < d; ++j) gw[j] += err * Xs[i][j];
}
double scale = 2.0 / B;
w0Tilde -= lrE * scale * gw0;
for (int j = 0; j < d; ++j) wTilde[j] -= lrE * scale * gw[j];
}
// Epoch loss on full training set (standardized space).
double loss = 0.0;
for (int i = 0; i < n; ++i) {
double pred = w0Tilde;
for (int j = 0; j < d; ++j) pred += wTilde[j] * Xs[i][j];
double e = pred - ys[i]; loss += e * e;
}
loss /= n;
if (!isfinite(loss))
throw runtime_error(
"GDRegressor::fit: loss became non-finite at epoch " +
to_string(epoch) +
"; reduce learning rate or check input scale");
lossHistory.push_back(loss); ++epochsRun;
double rel = fabs(prevLoss - loss) / (fabs(prevLoss) + 1e-8);
if (rel < cfg.tol) break;
prevLoss = loss;
}
syncRawParams();
if (!isfinite(intercept))
throw runtime_error("GDRegressor::fit: raw intercept became non-finite");
for (int j = 0; j < nFeatures; ++j)
if (!isfinite(coef[j]))
throw runtime_error(
"GDRegressor::fit: raw coefficient became non-finite at index " +
to_string(j));
fitted = true;
}
// Apply a single-sample gradient step using frozen standardization constants.
void fitOnline(const vector<double>& x, double yVal, double lr = 0.01) {
if (!fitted)
throw logic_error("GDRegressor::fitOnline: call fit() before fitOnline()");
if ((int)x.size() != nFeatures)
throw invalid_argument("GDRegressor::fitOnline: expected " +
to_string(nFeatures) + " features, got " +
to_string(x.size()));
for (int j = 0; j < nFeatures; ++j)
if (!isfinite(x[j]))
throw invalid_argument("GDRegressor::fitOnline: non-finite feature at index " +
to_string(j));
if (!isfinite(yVal))
throw invalid_argument("GDRegressor::fitOnline: non-finite target");
if (!isfinite(lr) || lr <= 0.0)
throw invalid_argument(
"GDRegressor::fitOnline: learning rate must be finite and positive");
// Standardize with frozen training-set constants.
vector<double> xs(nFeatures);
for (int j = 0; j < nFeatures; ++j) {
xs[j] = (x[j] - xMean[j]) / xStd[j];
if (!isfinite(xs[j]))
throw runtime_error(
"GDRegressor::fitOnline: standardized feature became non-finite");
}
double ys = (yVal - yMean) / yStd;
if (!isfinite(ys))
throw runtime_error(
"GDRegressor::fitOnline: standardized target became non-finite");
double pred = w0Tilde;
for (int j = 0; j < nFeatures; ++j) pred += wTilde[j] * xs[j];
double err = pred - ys;
if (!isfinite(err))
throw runtime_error(
"GDRegressor::fitOnline: prediction error became non-finite");
double nextW0 = w0Tilde - lr * 2.0 * err;
vector<double> nextW = wTilde;
vector<double> nextCoef(nFeatures);
for (int j = 0; j < nFeatures; ++j) {
nextW[j] -= lr * 2.0 * err * xs[j];
nextCoef[j] = nextW[j] * yStd / xStd[j];
if (!isfinite(nextW[j]) || !isfinite(nextCoef[j]))
throw runtime_error(
"GDRegressor::fitOnline: update became non-finite at index " +
to_string(j));
}
double nextIntercept = yMean + yStd * nextW0;
for (int j = 0; j < nFeatures; ++j)
nextIntercept -= nextCoef[j] * xMean[j];
if (!isfinite(nextW0) || !isfinite(nextIntercept))
throw runtime_error(
"GDRegressor::fitOnline: intercept update became non-finite");
w0Tilde = nextW0;
wTilde = nextW;
coef = nextCoef;
intercept = nextIntercept;
}
double predict(const vector<double>& x) const {
if (!fitted)
throw logic_error("GDRegressor::predict: model is not fitted");
if ((int)x.size() != nFeatures)
throw invalid_argument("GDRegressor::predict: expected " +
to_string(nFeatures) + " features, got " +
to_string(x.size()));
double out = intercept;
for (int j = 0; j < nFeatures; ++j) {
if (!isfinite(x[j]))
throw invalid_argument(
"GDRegressor::predict: non-finite feature at index " +
to_string(j));
out += coef[j] * x[j];
}
if (!isfinite(out))
throw runtime_error("GDRegressor::predict: prediction became non-finite");
return out;
}
vector<double> predict(const Matrix& X) const {
vector<double> out; out.reserve(X.size());
for (const auto& row : X) out.push_back(predict(row));
return out;
}
};
// ============================================================
// main: five demonstration blocks
// ============================================================
int main() {
// Dataset: y = 3 + 2*x1 - 1.5*x2 + noise, 240 points, seed 7
mt19937 rng(7);
uniform_real_distribution<double> ux(-3.0, 3.0);
normal_distribution<double> noise(0.0, 0.6);
const int N = 240;
Matrix X; vector<double> y;
for (int i = 0; i < N; ++i) {
double x1 = ux(rng), x2 = ux(rng);
X.push_back({x1, x2});
y.push_back(3.0 + 2.0*x1 - 1.5*x2 + noise(rng));
}
// 80/20 split with same shuffle as OLS post
vector<int> order(N); iota(order.begin(), order.end(), 0);
shuffle(order.begin(), order.end(), rng);
Matrix Xtr, Xte; vector<double> ytr, yte;
int nTrain = (int)(0.8 * N);
for (int k = 0; k < N; ++k) {
int i = order[k];
if (k < nTrain) { Xtr.push_back(X[i]); ytr.push_back(y[i]); }
else { Xte.push_back(X[i]); yte.push_back(y[i]); }
}
// --- Batch GD ---
GDConfig cfgBatch; cfgBatch.batchSize = -1; cfgBatch.lr = 0.05;
cfgBatch.lrDecay = 0.005; cfgBatch.maxEpochs = 500; cfgBatch.tol = 1e-9;
GDRegressor mBatch; mBatch.fit(Xtr, ytr, cfgBatch);
auto pBTr = mBatch.predict(Xtr), pBTe = mBatch.predict(Xte);
printf("=== Batch GD (full batch, lr=0.05, decay=0.005) ===\n");
printf("intercept=%.4f coef=[%.4f, %.4f] epochs=%d\n",
mBatch.intercept, mBatch.coef[0], mBatch.coef[1], mBatch.epochsRun);
printf("train RMSE=%.4f R^2=%.4f\n", rmse(ytr, pBTr), r2(ytr, pBTr));
printf("test RMSE=%.4f R^2=%.4f\n", rmse(yte, pBTe), r2(yte, pBTe));
printf("first loss=%.6f final loss=%.6f\n",
mBatch.lossHistory.front(), mBatch.lossHistory.back());
// --- Mini-batch GD ---
GDConfig cfgMini; cfgMini.batchSize = 32; cfgMini.lr = 0.05;
cfgMini.lrDecay = 0.005; cfgMini.maxEpochs = 500; cfgMini.tol = 1e-9; cfgMini.seed = 42;
GDRegressor mMini; mMini.fit(Xtr, ytr, cfgMini);
auto pMTr = mMini.predict(Xtr), pMTe = mMini.predict(Xte);
printf("\n=== Mini-batch GD (batch=32, lr=0.05, decay=0.005) ===\n");
printf("intercept=%.4f coef=[%.4f, %.4f] epochs=%d\n",
mMini.intercept, mMini.coef[0], mMini.coef[1], mMini.epochsRun);
printf("train RMSE=%.4f R^2=%.4f\n", rmse(ytr, pMTr), r2(ytr, pMTr));
printf("test RMSE=%.4f R^2=%.4f\n", rmse(yte, pMTe), r2(yte, pMTe));
// --- SGD ---
GDConfig cfgSGD; cfgSGD.batchSize = 1; cfgSGD.lr = 0.01;
cfgSGD.lrDecay = 0.01; cfgSGD.maxEpochs = 300; cfgSGD.tol = 1e-9; cfgSGD.seed = 99;
GDRegressor mSGD; mSGD.fit(Xtr, ytr, cfgSGD);
auto pSTr = mSGD.predict(Xtr), pSTe = mSGD.predict(Xte);
printf("\n=== SGD (batch=1, lr=0.01, decay=0.01) ===\n");
printf("intercept=%.4f coef=[%.4f, %.4f] epochs=%d\n",
mSGD.intercept, mSGD.coef[0], mSGD.coef[1], mSGD.epochsRun);
printf("train RMSE=%.4f R^2=%.4f\n", rmse(ytr, pSTr), r2(ytr, pSTr));
printf("test RMSE=%.4f R^2=%.4f\n", rmse(yte, pSTe), r2(yte, pSTe));
// --- Online update demo ---
GDRegressor mOnline; mOnline.fit(Xtr, ytr, cfgBatch);
printf("\n=== Online update demo ===\n");
printf("Before: intercept=%.4f coef=[%.4f, %.4f]\n",
mOnline.intercept, mOnline.coef[0], mOnline.coef[1]);
for (int i = 0; i < 20 && i < (int)Xte.size(); ++i)
mOnline.fitOnline(Xte[i], yte[i], 0.001);
auto pOTe = mOnline.predict(Xte);
printf("After 20 online steps: intercept=%.4f coef=[%.4f, %.4f]\n",
mOnline.intercept, mOnline.coef[0], mOnline.coef[1]);
printf("test RMSE=%.4f R^2=%.4f\n", rmse(yte, pOTe), r2(yte, pOTe));
// --- Divergence detection demo ---
printf("\n=== Divergence detection (lr=50) ===\n");
try {
GDConfig cfgBig; cfgBig.lr = 50.0; cfgBig.maxEpochs = 100;
GDRegressor mBig; mBig.fit(Xtr, ytr, cfgBig);
} catch (const runtime_error& e) {
printf("caught: %s\n", e.what());
}
return 0;
}Compiled with g++ -O2 -std=c++17 -Wall -Wextra -pedantic, this prints:
Batch GD converges in 161 epochs to the same intercept and coefficients as the OLS and QR/SVD posts (intercept $\approx 2.964$, $\beta_1 \approx 1.978$, $\beta_2 \approx -1.517$), confirming that all three solvers find the same global minimum. The standardized loss falls from $0.830$ to $0.023$: the absolute minimum for this dataset is not zero because the injected noise ($\sigma = 0.6$) contributes irreducible error. The fixed syncRawParams includes the $\sigma_y\,\tilde\beta_0$ term; for batch GD on this dataset $\tilde\beta_0 \approx 0$ at convergence so the batch-GD intercept is unchanged from the old formula. Mini-batch GD runs all 500 epochs because per-epoch loss fluctuates from mini-batch noise, giving intercept $2.9633$ slightly below the batch result: the residual $\tilde\beta_0$ after 500 epochs is small but nonzero, and the corrected conversion carries it into the raw intercept. SGD hits the 300-epoch budget; the higher gradient noise leaves a larger $\tilde\beta_0$ residual, shifting the converted intercept more visibly ($2.9532$) while the coefficients and test $R^2$ differ by less than $0.002$. Online updates move the parameters by a tiny amount (step size $0.001$) and apply the same corrected conversion; test metrics are essentially unchanged. Divergence detection fires at epoch 75 for lr=50, where each update overshoots the minimum by a growing margin until the loss overflows a double.
Complexity and Memory
| Step | Time | Memory |
|---|---|---|
| Compute $\mu_j$, $\sigma_j$, $\mu_y$, $\sigma_y$ | $O(nd)$ | $O(d)$ for the four statistics |
| Build standardized copies $\tilde X$, $\tilde y$ | $O(nd)$ | $O(nd + n)$, the main allocation |
| One epoch, batch size $B$ | $O(nd)$ (one update pass plus one full-loss pass) | $O(d)$ for the accumulated gradient and $O(n)$ for shuffled indices |
| $E$ epochs total | $O(End)$ | $O(nd)$ dominant |
| Predict on $m$ new rows | $O(md)$ | $O(m)$ for output |
The dominant memory cost of batch training is the $O(nd)$ allocation for $\tilde X$; the gradient accumulator itself is only $O(d)$ for every batch size because rows are consumed one at a time. Consequently, this teaching implementation does not provide an out-of-core batch fit. The separate fitOnline method has $O(d)$ per-sample working memory once the initial model and frozen standardizer exist. The $E$ parameter is controllable: stopping early or tightening the convergence tolerance trades accuracy for runtime. Direct solvers (QR, SVD) have no such dial; they run in one fixed-cost factorization.
Tuning Guide
- Learning rate. A practical starting range for standardized data is $\alpha = 0.05{-}0.1$. Loss that oscillates or grows after the first few epochs suggests halving $\alpha$; smooth but slow descent suggests doubling it. For the normalized MSE used here, batch gradient descent is stable when $\alpha < n/\sigma_{\max}^2(\tilde X_b)$, but estimating that bound adds a spectral computation, so empirical tuning is simpler in practice.
- Decay coefficient. A value of $0.001{-}0.01$ works well for a few hundred epochs. Too large and the rate drops below a useful threshold before convergence; too small and the rate barely decreases and late-training oscillations persist.
- Batch size. Full batch is the most stable choice for small $n$. Mini-batches of $32{-}64$ are a good default when $n$ is large. SGD ($B=1$) is rarely preferable to a small mini-batch for linear regression; the gradient variance is very high without any benefit over a direct solver.
- Max epochs. A generous budget (e.g.\ $1000{-}5000$) allows the convergence check to exit early. The final
epochsRunvalue reveals whether the budget is the binding constraint. - Convergence tolerance. $10^{-7}$ to $10^{-9}$ is appropriate for double precision. Tighter tolerances simply run more epochs without improving the fit beyond what the data noise floor allows.
Practical Diagnostics
The regressor exposes enough state to diagnose a run without adding instrumentation. Each check below reads fields the implementation already maintains.
- Learning-rate health. The first few entries of
lossHistoryare the fastest signal. A healthy rate shows the loss dropping monotonically for full batch, as in the complete program where it falls from $0.830$ to $0.023$. A loss that grows across the opening epochs means the rate sits above the stability interval $\alpha < n/\sigma_{\max}^2(\tilde X_b)$, and the divergence guard eventually throws; the fix is to halve the rate. A loss that decreases by a negligible amount per epoch means the rate is too small and can be doubled. - Loss history shape. Plotting or scanning
lossHistoryseparates the two failure modes above from a healthy plateau. A curve that flattens well above the irreducible noise floor points to an over-aggressive decay coefficient that shrank the rate before convergence; a curve that never flattens within the budget meansmaxEpochsis the binding constraint, visible asepochsRunequal tomaxEpochs. - Gradient norm. The Euclidean norm of the standardized gradient is the natural stopping quantity, since it is zero exactly at the minimum. Accumulating $\lVert (g_0, g_1, \dots, g_d) \rVert_2$ for the full-batch gradient and watching it shrink toward zero confirms progress even on runs where the loss change is dominated by mini-batch noise. A gradient norm that stalls at a nonzero value while the loss oscillates indicates the noise floor of a stochastic run rather than a true optimum.
- Online distribution shift. Because
fitOnlinestandardizes each new sample with the frozen training constants, a shift in the incoming distribution shows up as standardized inputs $\tilde x_j = (x_j - \mu_j)/\sigma_j$ drifting away from zero mean and unit spread. Tracking the running mean and variance of the standardized new points, and comparing them against the expected $0$ and $1$, flags when the frozen standardizer no longer matches the live data and a full retrain is warranted.
Limitations and Comparison with Direct Solvers
- No closed-form convergence guarantee for fixed $\alpha$. Direct solvers find the exact (or min-norm) solution in one factorization. Gradient descent with a fixed learning rate converges asymptotically and the final iterate may differ from the true OLS solution by a small but nonzero amount that depends on $\alpha$, $\lambda$, and the number of epochs.
- Near-collinear features. Collinear features make the loss surface highly elongated, and gradient descent can take many epochs to traverse the narrow valley. A regularized variant (ridge) or a preconditioned optimizer is more effective for near-singular designs; the next post covers ridge regression, which avoids the elongated landscape by adding a penalty term.
- Rank-deficient designs. Zero-initialized full-batch gradient descent converges to the minimum-norm solution in the coordinates being optimized when the learning rate is stable and optimization reaches the exact-arithmetic limit. Here those coordinates are standardized, so unequal feature scales can select a different raw-space coefficient vector than the unscaled SVD pseudoinverse. Finite epoch budgets, stochastic updates, and floating-point roundoff add further differences, and gradient descent exposes neither singular values nor a numerical rank. The SVD solver is appropriate when raw-design rank detection and a controlled raw-space minimum-norm pseudoinverse are required.
- Scaling sensitivity. Even with standardization, gradient descent can struggle if the features have very different curvatures in the loss landscape (e.g.\ one feature with 100$\times$ more data variance than another in the original space). The standardization step largely neutralizes this for linear regression, but other models may need additional preconditioning.
- When direct solvers are better. For moderate dense datasets, QR or SVD often reaches a more accurate answer with less tuning. Gradient descent becomes attractive when repeated $O(nd)$ passes are cheaper than an $O(nd^2)$ factorization, when data arrives in a stream, or when the model is extended beyond linear regression. The batch implementation here still allocates a standardized copy, so only
fitOnlinedemonstrates true $O(d)$ incremental updates.
Next Step
All three solvers so far (OLS, QR/SVD, gradient descent) minimize the unpenalized MSE: they find the best-fitting linear model with no constraint on the size of $\beta$. When features are correlated, the OLS estimate has high variance even if it is unbiased; adding a penalty term on $\lVert \beta \rVert_2^2$ (ridge) or $\lVert \beta \rVert_1$ (lasso) shrinks the coefficients and can substantially reduce test error. The next post implements ridge as a small modification of the normal equations and lasso via coordinate descent, both sharing the same fit/predict contract.