16. Matrix calculus for least squares and logistic regression

Every tool this module has built — gradients, convexity, the Hessian, the chain rule — converges here, on the two most fundamental models in machine learning. Linear regression turns out to have a closed-form exact solution, derivable by setting a gradient to zero exactly like §10.6's critical points; logistic regression does not, but its gradient takes a remarkably clean form that makes gradient descent (§14.0) straightforward to run.

Prerequisite note. This is the most matrix-heavy lesson in the series, and it uses linear algebra the calculus series does not develop: the transpose X^T, matrix multiplication, the matrix inverse (X^TX)^{-1}, and positive semidefiniteness. The calculus below — take a gradient, set it to zero — is exactly §10.6's method; the matrix notation is what is new. If §13.8 and §14.2 went smoothly, so will this — and the companion linear-algebra course's §20.5 (least squares and orthogonal projections) derives this same (X^TX)^{-1}X^Ty formula from scratch, via orthogonality rather than calculus, as a second, independent route to it.

Linear regression, in matrix form

Given n data points, each with p features, stack them into a design matrix X (one row per example, one column per feature, typically with a constant column of 1s for the bias/intercept) and a target vector \vec y. A linear model predicts \hat y=Xw for weight vector w, and the standard least-squares loss is the sum of squared errors:

L(w)=|X w-\vec y|^2=(Xw-\vec y)^T(Xw-\vec y)

— a direct matrix generalization of §5.6's "sum of squared deviations" idea from statistics, now applied to a vector of residuals instead of a single scalar spread.

The gradient of the least-squares loss

Expanding L(w)=w^TX^TXw-2w^TX^T\vec y+\vec y^T\vec y and differentiating with respect to the vector w (using matrix-calculus rules that are themselves just §10.5's gradient, applied component by component — each entry of \nabla_wL is an ordinary partial derivative $\frac{\partial L}{\partial w_i}$):

\nabla_wL=2X^TXw-2X^T\vec y=2X^T(Xw-\vec y)

Setting this gradient to zero (§10.6's critical-point condition, applied to a vector variable) gives the normal equations:

X^TXw=X^T\vec y\ \Longrightarrow\ \boxed{w=(X^TX)^{-1}X^T\vec y}

A closed-form exact solution — no gradient descent required. This is possible precisely because L(w) is a quadratic function of w, whose Hessian is the constant matrix 2X^TX — always positive semidefinite (§14.2), so L is convex (§14.1) and this single critical point is guaranteed to be the global minimum, exactly the payoff convexity promised two lessons ago, here delivering an explicit formula rather than just a convergence guarantee.

Logistic regression: a clean gradient, no closed form

For binary classification, replace the linear prediction with a sigmoid (§14.3's activation function): \hat y=\sigma(Xw), and replace squared-error loss with cross-entropy (log-loss):

L(w)=-\sum_{i=1}^n\Big[y_i\ln\hat y_i+(1-y_i)\ln(1-\hat y_i)\Big]

Differentiating (via §14.3's exact chain-rule mechanics — the sigmoid derivative \sigma'(z)=\sigma(z)(1-\sigma(z)) cancels beautifully against terms from differentiating the logarithms) produces a gradient with a strikingly simple final form:

\nabla_wL=X^T\big(\sigma(Xw)-\vec y\big)

This is exactly the same "prediction minus target" pattern that appeared as the shared factor in §14.3's single-neuron backpropagation example, now written for an entire dataset at once via matrix multiplication. Unlike linear regression, setting this gradient to zero has no closed-form algebraic solution — \sigma is nonlinear, so X^T(\sigma(Xw)-\vec y)=0 can't be solved by simple matrix inversion. This is exactly why logistic regression is trained with gradient descent (§14.0) rather than a formula: the loss is still convex (its Hessian is positive semidefinite everywhere, guaranteeing gradient descent finds the global minimum, per §14.1), but no shortcut around the iteration exists.

Doing it in Python

Solving a small linear regression problem via the normal equations, and confirming the resulting gradient is genuinely zero:

import sympy as sp

X = sp.Matrix([[1, 1], [1, 2], [1, 3]])   # bias column + one feature, 3 data points
y = sp.Matrix([2, 3, 5])

w = (X.T * X).inv() * X.T * y
print(f"w (normal equations) = {w.T}")

gradient_at_solution = 2 * X.T * (X * w - y)
print(f"gradient at this w = {gradient_at_solution.T}   (should be zero)")

Confirming the normal-equations solution matches directly solving \partial L/\partial w_0=0, \partial L/\partial w_1=0 by hand — an independent check using ordinary scalar calculus:

import sympy as sp

w0, w1 = sp.symbols('w0 w1')
xs, ys = [1, 2, 3], [2, 3, 5]

L = sum((w0 + w1*xi - yi)**2 for xi, yi in zip(xs, ys))
solution = sp.solve([sp.diff(L, w0), sp.diff(L, w1)], [w0, w1])
print(f"direct scalar-calculus solution: {solution}")

Confirming the logistic regression gradient formula X^T(\sigma(Xw)-y) matches direct differentiation of the cross-entropy loss, for a tiny 2-example dataset:

import sympy as sp

w0, w1 = sp.symbols('w0 w1')
x_data, y_data = [1, 2], [0, 1]

def sigmoid(z):
    return 1 / (1 + sp.exp(-z))

predictions = [sigmoid(w0 + w1*x) for x in x_data]
L = -sum(y*sp.log(p) + (1-y)*sp.log(1-p) for y, p in zip(y_data, predictions))

dL_dw0_direct = sp.simplify(sp.diff(L, w0))
formula_w0 = sum((p - y) for p, y in zip(predictions, y_data))   # X^T(sigma(Xw)-y), first row (bias column of 1s)

print(f"direct dL/dw0:  {dL_dw0_direct}")
print(f"formula match:  {sp.simplify(dL_dw0_direct - formula_w0) == 0}")

Worked example

Fit a linear model to the three points (1,2), (2,3), (3,5) using the normal equations.

X=\begin{pmatrix}1&1\\1&2\\1&3\end{pmatrix},\qquad\vec y=\begin{pmatrix}2\\3\\5\end{pmatrix}

X^TX=\begin{pmatrix}3&6\\6&14\end{pmatrix},\qquad X^T\vec y=\begin{pmatrix}10\\23\end{pmatrix}

(X^TX)^{-1}=\frac1{3(14)-6(6)}\begin{pmatrix}14&-6\\-6&3\end{pmatrix}=\frac1{6}\begin{pmatrix}14&-6\\-6&3\end{pmatrix}

w=(X^TX)^{-1}X^T\vec y=\frac16\begin{pmatrix}14&-6\\-6&3\end{pmatrix}\begin{pmatrix}10\\23\end{pmatrix}=\frac16\begin{pmatrix}140-138\\-60+69\end{pmatrix}=\frac16\begin{pmatrix}2\\9\end{pmatrix}

\boxed{w=\left\langle\frac13,\ \frac32\right\rangle}

— the fitted line is \hat y=\frac13+\frac32x.

Sanity check. At x=1: \hat y=\frac13+\frac32\approx1.83, close to the actual y=2. At x=3: \hat y=\frac13+\frac92\approx4.83, close to the actual y=5 — the fitted line tracks the three points reasonably, never matching exactly (since three points don't lie on a single line), which is exactly what a least-squares fit is supposed to do: minimize total squared error, not necessarily pass through any point exactly. The gradient 2X^T(Xw-\vec y) evaluates to the zero vector at this w, confirmed directly in the "Doing it in Python" section — the defining property of the critical point the normal equations were built to find. ✓

Your turn

1. For X=\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}, \vec y=\begin{pmatrix}1\\3\\5\end{pmatrix}, set up (don't necessarily compute by hand) the normal equations X^TXw=X^T\vec y.

2. Explain why the least-squares loss L(w)=|Xw-\vec y|^2 is guaranteed convex for any design matrix X, referencing its Hessian.

3. True or false: logistic regression's cross-entropy loss is convex, just like linear regression's squared-error loss.

Solutions

1.

X^TX=\begin{pmatrix}1&1&1\\0&1&2\end{pmatrix}\begin{pmatrix}1&0\\1&1\\1&2\end{pmatrix}=\begin{pmatrix}3&3\\3&5\end{pmatrix}

X^T\vec y=\begin{pmatrix}1&1&1\\0&1&2\end{pmatrix}\begin{pmatrix}1\\3\\5\end{pmatrix}=\begin{pmatrix}9\\13\end{pmatrix}

\boxed{\begin{pmatrix}3&3\\3&5\end{pmatrix}w=\begin{pmatrix}9\\13\end{pmatrix}}

2. The Hessian of L(w)=w^TX^TXw-2w^TX^T\vec y+\vec y^T\vec y is the constant matrix 2X^TX. For any matrix X, X^TX is positive semidefinite: for any vector v, v^T(X^TX)v=(Xv)^T(Xv)=|Xv|^2\ge0 — a sum of squares can never be negative. Since the Hessian is positive semidefinite everywhere (not just at one point), L is convex regardless of what the data actually is — a structural guarantee, not something that needs checking case by case.

3. True. Cross-entropy loss is convex in w (its Hessian, built from products of the form \sigma(z)(1-\sigma(z)) times outer products of data vectors, is also positive semidefinite everywhere, by an argument parallel to problem 2's) — which is exactly why gradient descent reliably finds logistic regression's global optimum despite the lack of a closed-form formula. Convexity and having a closed-form solution are genuinely separate properties: linear regression happens to have both, logistic regression has only the first.

Check yourself in code

Solve the normal equations for X=\begin{pmatrix}1&1\\1&2\\1&3\end{pmatrix}, \vec y=\begin{pmatrix}2\\3\\5\end{pmatrix}, and confirm the gradient at the solution is zero.

Print exactly this:

w = [1/3, 3/2]
gradient at solution = [0, 0]

import sympy as sp

X = sp.Matrix([[1, 1], [1, 2], [1, 3]])
y = sp.Matrix([2, 3, 5])

w = (X.T * X).inv() * X.T * y
print("w = ...")

gradient_at_solution = 2 * X.T * (X * w - y)
print("gradient at solution = ...")

import sympy as sp

X = sp.Matrix([[1, 1], [1, 2], [1, 3]])
y = sp.Matrix([2, 3, 5])

w = (X.T * X).inv() * X.T * y
print(f"w = {list(w)}")

gradient_at_solution = 2 * X.T * (X * w - y)
print(f"gradient at solution = {list(gradient_at_solution)}")

Linear regression's least-squares loss, L(w)=|Xw-\vec y|^2, has gradient 2X^T(Xw-\vec y) and Hessian 2X^TX — always positive semidefinite (§14.2), so §14.1's convexity guarantees the single critical point found by setting the gradient to zero, the normal equations w=(X^TX)^{-1}X^T\vec y, is the exact global minimum, no iteration required. Logistic regression's cross-entropy loss has an equally clean gradient, X^T(\sigma(Xw)-\vec y) — the same "prediction minus target" pattern from §14.3's backpropagation, generalized to a whole dataset — but the sigmoid's nonlinearity rules out a closed form, making gradient descent the only route to its otherwise-guaranteed convex optimum.

Next: extending optimization one final step — minimizing a loss subject to constraints, generalizing §10.7's Lagrange multipliers into the inequality-constraint framework (KKT conditions) that underlies support vector machines and much of modern optimization.