{ "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" }, "colab": { "provenance": [] } }, "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Day 24: Matrix Calculus & Gradients\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "Differentiating with respect to a whole vector at once, and deriving the OLS, ridge, WLS and logit first-order conditions.\n", "\n", "**Learning outcomes**\n", "\n", "- Apply the core rules ∂(a′x)/∂x = a and ∂(x′Ax)/∂x = (A + A′)x.\n", "- Derive the gradient and Hessian of the sum of squared residuals and obtain the normal equations.\n", "- Derive the ridge and weighted-least-squares estimators.\n", "- Compute the gradient and Hessian of the logit log-likelihood and the sensitivity of OLS to one observation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### Layout\n", "\n", "For a scalar function $f(x)$ of a vector $x \\in \\mathbb{R}^k$, the **gradient** is the column vector $\\partial f/\\partial x = (\\partial f/\\partial x_1, \\dots, \\partial f/\\partial x_k)^{\\top}$, and the **Hessian** is the $k \\times k$ matrix $\\partial^2 f/\\partial x\\,\\partial x^{\\top}$. For a vector function $g(x) \\in \\mathbb{R}^m$, the **Jacobian** $\\partial g/\\partial x^{\\top}$ is $m \\times k$.\n", "\n", "### Four rules cover almost everything\n", "\n", "| Function | Gradient |\n", "|---|---|\n", "| $a^{\\top}x$ | $a$ |\n", "| $x^{\\top}Ax$ | $(A + A^{\\top})x$, which is $2Ax$ if $A$ is symmetric |\n", "| $Ax$ (vector) | Jacobian $A$ |\n", "| $f(g(x))$ | chain rule: $\\frac{\\partial g}{\\partial x^{\\top}}^{\\!\\top}\\nabla f(g(x))$ |\n", "\n", "### The least-squares objective\n", "\n", "$$\n", "S(\\beta) = (y - X\\beta)^{\\top}(y - X\\beta) = y^{\\top}y - 2\\beta^{\\top}X^{\\top}y + \\beta^{\\top}X^{\\top}X\\beta,\n", "$$\n", "\n", "$$\n", "\\nabla S(\\beta) = -2X^{\\top}y + 2X^{\\top}X\\beta = -2X^{\\top}(y - X\\beta), \\qquad \\nabla^2 S = 2X^{\\top}X.\n", "$$\n", "\n", "Setting the gradient to zero gives the **normal equations** $X^{\\top}X\\hat\\beta = X^{\\top}y$. The Hessian is PSD (Day 17), so $S$ is convex and the solution is a minimum (Day 25).\n", "\n", "### The same template, other estimators\n", "\n", "- **Ridge:** $S(\\beta) + \\lambda\\beta^{\\top}\\beta \\Rightarrow \\hat\\beta = (X^{\\top}X + \\lambda I)^{-1}X^{\\top}y$.\n", "- **Weighted least squares:** $(y - X\\beta)^{\\top}W(y - X\\beta) \\Rightarrow \\hat\\beta = (X^{\\top}WX)^{-1}X^{\\top}Wy$.\n", "- **Logit:** $\\nabla\\ell = X^{\\top}(y - p)$ and $\\nabla^2\\ell = -X^{\\top}VX$ with $V = \\operatorname{diag}(p_i(1 - p_i))$.\n", "\n", "> **Causal link.** Weighted least squares with weights $1/\\hat e(x)$ and $1/(1 - \\hat e(x))$ is how IPW regressions are run; ridge and lasso are the nuisance learners in double machine learning.\n", "\n", "## Part 1 · Warm-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · A linear function\n", "\n", "Find the gradient of $f(x) = a^{\\top}x$ with $a = (2, -1, 3)$.\n", "\n", "
Show solution\n", "\n", "$\\nabla f = a = (2, -1, 3)^{\\top}$, the same at every $x$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · A quadratic form\n", "\n", "For $A = \\begin{pmatrix}2 & 1\\\\1 & 3\\end{pmatrix}$ and $f(x) = x^{\\top}Ax$, find $\\nabla f$ at $x = (1, 1)$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $A$ is symmetric, so $\\nabla(x^{\\top}Ax) = 2Ax$.\n", "\n", "**Step 2.** $Ax = (2 + 1,\\; 1 + 3) = (3, 4)$ at $x = (1, 1)$.\n", "\n", "**Answer:** $\\nabla f = (6, 8)$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · Its Hessian\n", "\n", "Find the Hessian of $f$ in Exercise 2.\n", "\n", "
Show solution\n", "\n", "$2A = \\begin{pmatrix}4 & 2\\\\2 & 6\\end{pmatrix}$, constant because $f$ is quadratic.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · Component check\n", "\n", "Find the gradient of $f(x_1, x_2) = x_1^2 + x_1x_2 + 2x_2^2$ at $(1, 2)$ directly and via the matrix rule.\n", "\n", "
Show solution\n", "\n", "Directly: $(2x_1 + x_2,\\; x_1 + 4x_2) = (4, 9)$. Matrix: $A = \\begin{pmatrix}1 & 0.5\\\\0.5 & 2\\end{pmatrix}$, $2Ax = 2(1 + 1, 0.5 + 4) = (4, 9)$. ✓\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · A Jacobian\n", "\n", "Find the Jacobian of $g(x) = (x_1x_2,\\; x_1 + x_2^2)^{\\top}$.\n", "\n", "
Show solution\n", "\n", "$\\frac{\\partial g}{\\partial x^{\\top}} = \\begin{pmatrix}x_2 & x_1\\\\1 & 2x_2\\end{pmatrix}$.\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · The linear rule\n", "\n", "Prove $\\partial(a^{\\top}x)/\\partial x = a$.\n", "\n", "
Show solution\n", "\n", "$a^{\\top}x = \\sum_j a_jx_j$, so $\\partial/\\partial x_\\ell = a_\\ell$. Stacking gives $a$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · The quadratic rule\n", "\n", "Prove $\\partial(x^{\\top}Ax)/\\partial x = (A + A^{\\top})x$.\n", "\n", "
Show solution\n", "\n", "$x^{\\top}Ax = \\sum_i\\sum_j a_{ij}x_ix_j$. Differentiating with respect to $x_\\ell$ picks up terms with $i = \\ell$ and with $j = \\ell$: $\\sum_j a_{\\ell j}x_j + \\sum_i a_{i\\ell}x_i = (Ax)_\\ell + (A^{\\top}x)_\\ell$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · The normal equations\n", "\n", "Derive $\\nabla S(\\beta)$ and the normal equations.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Expand $S(\\beta) = y^{\\top}y - 2\\beta^{\\top}X^{\\top}y + \\beta^{\\top}X^{\\top}X\\beta$ (the two cross terms are equal scalars).\n", "\n", "**Step 2.** By the rules, $\\nabla S = -2X^{\\top}y + 2X^{\\top}X\\beta$.\n", "\n", "**Step 3.** Setting it to zero: $X^{\\top}X\\hat\\beta = X^{\\top}y$, so $\\hat\\beta = (X^{\\top}X)^{-1}X^{\\top}y$ under full rank.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · A unique minimum\n", "\n", "Show that the Hessian of $S$ is PD when $X$ has full column rank, so the solution is the unique minimiser.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\nabla^2S = 2X^{\\top}X$ and $v^{\\top}X^{\\top}Xv = \\|Xv\\|^2 > 0$ for $v \\ne 0$ under full rank (Day 18).\n", "\n", "**Step 2.** A function with a PD Hessian everywhere is strictly convex, so its stationary point is the unique global minimum (Day 25).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Ridge regression\n", "\n", "Derive the minimiser of $S(\\beta) + \\lambda\\|\\beta\\|^2$ and show it exists even when $X$ is rank-deficient (for $\\lambda > 0$).\n", "\n", "
Show solution\n", "\n", "**Step 1.** Gradient: $-2X^{\\top}y + 2X^{\\top}X\\beta + 2\\lambda\\beta = 0 \\Rightarrow (X^{\\top}X + \\lambda I)\\beta = X^{\\top}y$.\n", "\n", "**Step 2.** For $v \\ne 0$, $v^{\\top}(X^{\\top}X + \\lambda I)v = \\|Xv\\|^2 + \\lambda\\|v\\|^2 > 0$, so the matrix is PD and invertible whatever the rank of $X$.\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · A gradient from the farm data\n", "\n", "For the four farms of Day 17 (visits $x$, yields $y$), $X^{\\top}X = \\begin{pmatrix}4 & 10\\\\10 & 30\\end{pmatrix}$ and $X^{\\top}y = (16, 47)$. Compute $\\nabla S$ at $\\beta = (0, 1)$ (no intercept, 1 t/ha per visit) and at $\\hat\\beta = (0.5, 1.4)$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\nabla S(\\beta) = -2(X^{\\top}y - X^{\\top}X\\beta)$.\n", "\n", "**Step 2.** At $(0, 1)$: $X^{\\top}X\\beta = (10, 30)$, so $\\nabla S = -2\\big((16, 47) - (10, 30)\\big) = -2(6, 17) = (-12, -34)$.\n", "\n", "**Step 3.** Both entries are negative, so increasing either coefficient lowers the SSR.\n", "\n", "**Step 4.** At $\\hat\\beta$: $X^{\\top}X\\hat\\beta = (4(0.5) + 10(1.4),\\; 10(0.5) + 30(1.4)) = (16, 47)$, so $\\nabla S = 0$.\n", "\n", "**Answer:** $(-12, -34)$ at the guess, and zero at the OLS estimate, as the normal equations require.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · Weighted least squares\n", "\n", "Derive the WLS estimator, and explain how weights $w_i = 1/\\hat e(x_i)$ for treated and $1/(1 - \\hat e(x_i))$ for control units turn a regression of $y$ on $(1, T)$ into an IPW estimator.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Objective: $(y - X\\beta)^{\\top}W(y - X\\beta)$ with $W$ diagonal (one weight per unit).\n", "\n", "**Step 2.** Gradient: $-2X^{\\top}W(y - X\\beta)$. Setting it to zero gives $X^{\\top}WX\\beta = X^{\\top}Wy$, so $\\hat\\beta = (X^{\\top}WX)^{-1}X^{\\top}Wy$.\n", "\n", "**Step 3.** With $X = (\\mathbf 1, T)$, the fitted values are the weighted group means, so the slope is (weighted mean of treated outcomes) − (weighted mean of control outcomes).\n", "\n", "**Step 4.** With weights $1/\\hat e(x_i)$ for treated and $1/(1 - \\hat e(x_i))$ for controls, that is the Hájek IPW estimator.\n", "\n", "**Answer:** an IPW estimate of the ATE is a weighted regression of $y$ on a treatment dummy (Days 14 and 44).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · The Hessian of a propensity-score logit\n", "\n", "For a logit model of programme participation, $\\ell(\\beta) = \\sum_i[y_ix_i^{\\top}\\beta - \\ln(1 + e^{x_i^{\\top}\\beta})]$ ($y_i = 1$ if household $i$ participates). Derive the gradient and Hessian and show that $\\ell$ is concave.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Derivative of the second term: $\\dfrac{\\partial}{\\partial\\beta}\\ln(1 + e^{x^{\\top}\\beta}) = \\dfrac{e^{x^{\\top}\\beta}}{1 + e^{x^{\\top}\\beta}}x = p\\,x$.\n", "\n", "**Step 2.** Gradient: $\\nabla\\ell = \\sum_i(y_i - p_i)x_i = X^{\\top}(y - p)$.\n", "\n", "**Step 3.** Differentiate $-p_ix_i$ again using $\\partial p_i/\\partial\\beta^{\\top} = p_i(1 - p_i)x_i^{\\top}$: $\\nabla^2\\ell = -\\sum_ip_i(1 - p_i)x_ix_i^{\\top} = -X^{\\top}VX$ with $V = \\operatorname{diag}(p_i(1 - p_i))$.\n", "\n", "**Step 4.** For any $v$, $v^{\\top}X^{\\top}VXv = \\sum_ip_i(1 - p_i)(x_i^{\\top}v)^2 \\ge 0$, so the Hessian is negative semi-definite.\n", "\n", "**Answer:** $\\ell$ is concave, so any stationary point is a global maximum and Newton's method is reliable (Day 25).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · The MLE of the error variance in a wage regression\n", "\n", "For a log-wage regression with normal errors, $\\ell = -\\frac n2\\ln(2\\pi\\sigma^2) - \\frac{S(\\beta)}{2\\sigma^2}$, where $S$ is the sum of squared residuals. Find the MLE of $\\sigma^2$ and compare it with $s^2$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Differentiate with respect to $\\sigma^2$ (treat it as one variable): $\\dfrac{\\partial\\ell}{\\partial\\sigma^2} = -\\dfrac{n}{2\\sigma^2} + \\dfrac{S}{2\\sigma^4}$.\n", "\n", "**Step 2.** Set to zero and multiply by $2\\sigma^4$: $-n\\sigma^2 + S = 0$.\n", "\n", "**Step 3.** At $\\hat\\beta$: $\\hat\\sigma^2 = S(\\hat\\beta)/n = e^{\\top}e/n$.\n", "\n", "**Answer:** $e^{\\top}e/n$, which is biased downward; $s^2 = e^{\\top}e/(n - k)$ is unbiased (Day 27).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · How much does one farm move the estimates?\n", "\n", "$\\hat\\beta = (X^{\\top}X)^{-1}X^{\\top}y$ is linear in $y$, so $\\partial\\hat\\beta/\\partial y^{\\top} = (X^{\\top}X)^{-1}X^{\\top}$. For the four farms, by how much do the intercept and the yield-per-visit slope change if farm 4's yield is recorded 1 t/ha higher?\n", "\n", "
Show solution\n", "\n", "**Step 1.** $(X^{\\top}X)^{-1} = \\frac{1}{20}\\begin{pmatrix}30 & -10\\\\-10 & 4\\end{pmatrix}$ (determinant $4 \\cdot 30 - 10^2 = 20$).\n", "\n", "**Step 2.** Farm 4's column of $X^{\\top}$ is $(1, 4)^{\\top}$.\n", "\n", "**Step 3.** $\\frac{1}{20}\\begin{pmatrix}30 - 40\\\\-10 + 16\\end{pmatrix} = \\frac{1}{20}\\begin{pmatrix}-10\\\\6\\end{pmatrix} = \\begin{pmatrix}-0.5\\\\0.3\\end{pmatrix}$.\n", "\n", "**Answer:** the slope rises by 0.3 t/ha per visit and the intercept falls by 0.5: an end-point farm tilts the line. A data-entry error there matters a lot.\n", "\n", "
\n", "\n", "## Verify in Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "\n", "def num_grad(f, x, h=1e-6):\n", " return np.array([(f(x + h * e) - f(x - h * e)) / (2 * h) for e in np.eye(len(x))])\n", "\n", "A = np.array([[2, 1], [1, 3.]])\n", "print(\"grad x'Ax at (1,1):\", num_grad(lambda x: x @ A @ x, np.array([1, 1.])).round(4))\n", "\n", "x = np.array([1, 2, 3, 4.]); y = np.array([2, 3, 5, 6.]); X = np.column_stack([np.ones(4), x])\n", "S = lambda b: (y - X @ b) @ (y - X @ b)\n", "for b in (np.array([0, 1.]), np.array([0.5, 1.4])):\n", " print(\"grad S at\", b, \":\", num_grad(S, b).round(4), \" formula:\", (-2 * X.T @ (y - X @ b)).round(4))\n", "\n", "print(\"dbeta/dy4:\", (np.linalg.inv(X.T @ X) @ X.T)[:, 3])\n", "lam = 0.5\n", "print(\"ridge:\", np.linalg.solve(X.T @ X + lam * np.eye(2), X.T @ y).round(4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- $\\partial(a^{\\top}x)/\\partial x = a$ and $\\partial(x^{\\top}Ax)/\\partial x = 2Ax$ (symmetric $A$) do most of the work.\n", "- Setting the gradient of the SSR to zero gives the normal equations; its Hessian $2X^{\\top}X$ guarantees a minimum.\n", "- Ridge, WLS and logit follow the same recipe: write the objective, differentiate, solve." ] } ] }