{
 "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 17: Matrix Multiplication, Transpose & the Design Matrix\n",
    "\n",
    "*Beginner to Ninja: Math for Causal Inference*\n",
    "\n",
    "Multiplying matrices, transposing products, and the cross-product matrices X′X and X′y that every regression is built from.\n",
    "\n",
    "**Learning outcomes**\n",
    "\n",
    "- Multiply conformable matrices and read Ax as a linear combination of the columns of A.\n",
    "- Prove (AB)′ = B′A′ and tr(AB) = tr(BA).\n",
    "- Build X′X and X′y from data and interpret their entries as sums of squares and cross-products.\n",
    "- Use the centring matrix M = I − (1/n)11′ to compute deviations from the mean."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory\n",
    "\n",
    "### The product\n",
    "\n",
    "If $A$ is $m \\times n$ and $B$ is $n \\times p$, the product $C = AB$ is $m \\times p$ with\n",
    "\n",
    "$$\n",
    "c_{ij} = \\sum_{\\ell=1}^n a_{i\\ell}b_{\\ell j} = (\\text{row } i \\text{ of } A) \\cdot (\\text{column } j \\text{ of } B).\n",
    "$$\n",
    "\n",
    "The inner dimensions must match (\"conformable\"). Matrix multiplication is **associative** and **distributive** but **not commutative**: in general $AB \\ne BA$, and one of them may not even exist.\n",
    "\n",
    "### Two ways to read Ax\n",
    "\n",
    "$$\n",
    "Ax = \\sum_{j} x_j\\,a_{\\cdot j}\n",
    "$$\n",
    "\n",
    "is a **linear combination of the columns** of $A$ with weights $x_j$. So $X\\beta$, the vector of fitted values, is \"$\\beta_0$ times the intercept column plus $\\beta_1$ times the education column plus …\". Alternatively, entry $i$ of $Ax$ is row $i$ of $A$ dotted with $x$: unit $i$'s prediction.\n",
    "\n",
    "A matrix also acts as a **linear transformation**: it maps vectors to vectors, preserving sums and scalar multiples. In 2D you can watch it stretch, rotate, shear or flip the plane (Simulation Lab).\n",
    "\n",
    "### Transpose and trace\n",
    "\n",
    "$$\n",
    "(AB)^{\\top} = B^{\\top}A^{\\top}, \\qquad \\operatorname{tr}(A) = \\sum_i a_{ii}, \\qquad \\operatorname{tr}(AB) = \\operatorname{tr}(BA).\n",
    "$$\n",
    "\n",
    "### The cross-product matrices\n",
    "\n",
    "For a design matrix $X$ ($n \\times k$) and outcome $y$ ($n \\times 1$):\n",
    "\n",
    "$$\n",
    "X^{\\top}X = \\Big[\\sum_i x_{ij}x_{i\\ell}\\Big]_{j,\\ell}, \\qquad X^{\\top}y = \\Big[\\sum_i x_{ij}y_i\\Big]_{j}.\n",
    "$$\n",
    "\n",
    "With an intercept and one regressor,\n",
    "\n",
    "$$\n",
    "X^{\\top}X = \\begin{pmatrix} n & \\sum x_i \\\\ \\sum x_i & \\sum x_i^2 \\end{pmatrix}, \\qquad\n",
    "X^{\\top}y = \\begin{pmatrix} \\sum y_i \\\\ \\sum x_iy_i \\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "On Day 27 the OLS estimator will be $\\hat\\beta = (X^{\\top}X)^{-1}X^{\\top}y$: everything it needs is in these two objects.\n",
    "\n",
    "### The centring matrix\n",
    "\n",
    "$$\n",
    "M_{\\mathbf 1} = I_n - \\tfrac1n\\mathbf{1}\\mathbf{1}^{\\top}, \\qquad M_{\\mathbf 1}y = y - \\bar y\\mathbf{1}.\n",
    "$$\n",
    "\n",
    "Multiplying by $M_{\\mathbf 1}$ subtracts the mean, the matrix version of \"demeaning\", which becomes fixed-effects estimation on Day 28.\n",
    "\n",
    "## Part 1 · Warm-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1 · Order matters\n",
    "\n",
    "Let $A = \\begin{pmatrix}1 & 2\\\\3 & 4\\end{pmatrix}$ and $B = \\begin{pmatrix}0 & 1\\\\1 & 0\\end{pmatrix}$. Compute $AB$ and $BA$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$AB = \\begin{pmatrix}2 & 1\\\\4 & 3\\end{pmatrix}$ (columns of $A$ swapped), $BA = \\begin{pmatrix}3 & 4\\\\1 & 2\\end{pmatrix}$ (rows swapped). They differ.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2 · Conformability\n",
    "\n",
    "$A$ is $2 \\times 3$ and $B$ is $3 \\times 4$. Which of $AB$, $BA$, $A^{\\top}A$, $B^{\\top}A^{\\top}$ exist, and what shapes are they?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $AB$: $2 \\times 4$.\n",
    "\n",
    "**Step 2.** $BA$: not defined ($4 \\ne 2$).\n",
    "\n",
    "**Step 3.** $A^{\\top}A$: $3 \\times 3$.\n",
    "\n",
    "**Step 4.** $B^{\\top}A^{\\top}$: $4 \\times 2$ (it equals $(AB)^{\\top}$).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3 · Columns, combined\n",
    "\n",
    "$A = \\begin{pmatrix}1 & 2\\\\3 & 1\\end{pmatrix}$, $x = (1, 1)^{\\top}$. Compute $Ax$ both row-by-row and as a combination of columns.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "Rows: $(1 + 2,\\; 3 + 1) = (3, 4)$. Columns: $1\\cdot(1, 3) + 1\\cdot(2, 1) = (3, 4)$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4 · Transpose of a product\n",
    "\n",
    "Check $(AB)^{\\top} = B^{\\top}A^{\\top}$ for the matrices in Exercise 1.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$(AB)^{\\top} = \\begin{pmatrix}2 & 4\\\\1 & 3\\end{pmatrix}$. $B^{\\top}A^{\\top} = \\begin{pmatrix}0 & 1\\\\1 & 0\\end{pmatrix}\\begin{pmatrix}1 & 3\\\\2 & 4\\end{pmatrix} = \\begin{pmatrix}2 & 4\\\\1 & 3\\end{pmatrix}$. ✓\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5 · Trace\n",
    "\n",
    "Compute $\\operatorname{tr}(A)$, $\\operatorname{tr}(AB)$ and $\\operatorname{tr}(BA)$ for Exercise 1.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$\\operatorname{tr}(A) = 5$, $\\operatorname{tr}(AB) = 2 + 3 = 5$, $\\operatorname{tr}(BA) = 3 + 2 = 5$. Equal even though $AB \\ne BA$.\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 2 · Proofs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6 · Transpose of a product\n",
    "\n",
    "Prove $(AB)^{\\top} = B^{\\top}A^{\\top}$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$[(AB)^{\\top}]_{ij} = (AB)_{ji} = \\sum_\\ell a_{j\\ell}b_{\\ell i} = \\sum_\\ell (B^{\\top})_{i\\ell}(A^{\\top})_{\\ell j} = (B^{\\top}A^{\\top})_{ij}$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7 · Cyclic trace\n",
    "\n",
    "Prove $\\operatorname{tr}(AB) = \\operatorname{tr}(BA)$ for $A$ ($m \\times n$) and $B$ ($n \\times m$).\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$\\operatorname{tr}(AB) = \\sum_i\\sum_\\ell a_{i\\ell}b_{\\ell i} = \\sum_\\ell\\sum_i b_{\\ell i}a_{i\\ell} = \\operatorname{tr}(BA)$. Day 20 uses this to show $\\operatorname{tr}(P) = k$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 8 · The structure of X′X\n",
    "\n",
    "Show that $X^{\\top}X$ is symmetric, that its diagonal holds the sums of squares of the columns, and derive its form with an intercept and one regressor.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $(X^{\\top}X)^{\\top} = X^{\\top}(X^{\\top})^{\\top} = X^{\\top}X$.\n",
    "\n",
    "**Step 2.** Entry $(j, \\ell)$ is column $j$ dotted with column $\\ell$, so the diagonal is $\\sum_i x_{ij}^2$.\n",
    "\n",
    "**Step 3.** With columns $\\mathbf 1$ and $x$: $\\mathbf 1^{\\top}\\mathbf 1 = n$, $\\mathbf 1^{\\top}x = \\sum x_i$, $x^{\\top}x = \\sum x_i^2$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 9 · Quadratic forms\n",
    "\n",
    "Show that $u^{\\top}Au = \\sum_i\\sum_j a_{ij}u_iu_j$, and write $u^{\\top}Au$ for $A = \\begin{pmatrix}2 & 1\\\\1 & 3\\end{pmatrix}$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$(Au)_i = \\sum_j a_{ij}u_j$, so $u^{\\top}Au = \\sum_i u_i\\sum_j a_{ij}u_j$. Here $u^{\\top}Au = 2u_1^2 + 2u_1u_2 + 3u_2^2$ (the off-diagonal 1 appears twice).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 10 · X′X is positive semi-definite\n",
    "\n",
    "Show that $v^{\\top}X^{\\top}Xv \\ge 0$ for every $v$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $v^{\\top}X^{\\top}Xv = (Xv)^{\\top}(Xv) = \\|Xv\\|^2 \\ge 0$.\n",
    "\n",
    "**Step 2.** It is zero only if $Xv = 0$, i.e. if some combination of the columns vanishes: the collinearity problem of Day 18.\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 3 · Applied"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 11 · Cross-products from data\n",
    "\n",
    "An agricultural extension service records, for four farms, the number of advisory visits $x = (1, 2, 3, 4)$ and maize yield $y = (2, 3, 5, 6)$ tonnes per hectare. For a regression of yield on an intercept and visits, build $X^{\\top}X$ and $X^{\\top}y$. (Days 20 and 24 reuse these four farms.)\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** The design matrix has a column of ones and the visits column, so we need $n$, $\\sum x$, $\\sum x^2$, $\\sum y$ and $\\sum xy$.\n",
    "\n",
    "**Step 2.** $n = 4$, $\\sum x = 10$, $\\sum x^2 = 1 + 4 + 9 + 16 = 30$.\n",
    "\n",
    "**Step 3.** $\\sum y = 16$, $\\sum xy = 1(2) + 2(3) + 3(5) + 4(6) = 2 + 6 + 15 + 24 = 47$.\n",
    "\n",
    "**Answer:**\n",
    "\n",
    "$$\n",
    "X^{\\top}X = \\begin{pmatrix}4 & 10\\\\10 & 30\\end{pmatrix}, \\qquad X^{\\top}y = \\begin{pmatrix}16\\\\47\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 12 · Fitted yields\n",
    "\n",
    "For the four farms, use $\\beta = (0.5, 1.4)^{\\top}$ (intercept 0.5 t/ha, 1.4 t/ha per visit). Compute fitted yields $X\\beta$ and the residuals $y - X\\beta$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Fitted yield $= 0.5 + 1.4x$: $1.9,\\; 3.3,\\; 4.7,\\; 6.1$.\n",
    "\n",
    "**Step 2.** Residuals $= y - \\hat y$: $2 - 1.9 = 0.1$, $3 - 3.3 = -0.3$, $5 - 4.7 = 0.3$, $6 - 6.1 = -0.1$.\n",
    "\n",
    "**Step 3.** Checks: $\\sum e = 0$ and $\\sum x_ie_i = 0.1 - 0.6 + 0.9 - 0.4 = 0$.\n",
    "\n",
    "**Answer:** $\\hat y = (1.9, 3.3, 4.7, 6.1)$ and $e = (0.1, -0.3, 0.3, -0.1)$. The zero sums are no accident: $(0.5, 1.4)$ is the OLS solution (Day 27).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 13 · A treatment dummy\n",
    "\n",
    "A sample has $n = 10$ units of which $n_1 = 4$ are treated. Write $X^{\\top}X$ for a regression on an intercept and the treatment dummy $T$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Entries: $n = 10$; $\\sum T_i = 4$; $\\sum T_i^2 = 4$ (because $T^2 = T$ for a 0/1 variable).\n",
    "\n",
    "**Answer:**\n",
    "\n",
    "$$\n",
    "X^{\\top}X = \\begin{pmatrix}10 & 4\\\\4 & 4\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 14 · Demeaning household incomes\n",
    "\n",
    "Three households have monthly incomes $y = (1, 2, 6)$ (thousand dollars). Write the centring matrix $M_{\\mathbf 1}$ for $n = 3$ and compute $M_{\\mathbf 1}y$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $M_{\\mathbf 1} = I - \\frac13\\mathbf 1\\mathbf 1^{\\top}$: diagonal entries $1 - \\frac13 = \\frac23$, off-diagonal $-\\frac13$:\n",
    "\n",
    "$$\n",
    "M_{\\mathbf 1} = \\frac13\\begin{pmatrix}2 & -1 & -1\\\\-1 & 2 & -1\\\\-1 & -1 & 2\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "**Step 2.** First entry: $\\frac13(2 \\cdot 1 - 2 - 6) = -2$; similarly the others.\n",
    "\n",
    "**Step 3.** Equivalently, subtract the mean $\\bar y = 3$: $(1 - 3, 2 - 3, 6 - 3)$.\n",
    "\n",
    "**Answer:** $M_{\\mathbf 1}y = (-2, -1, 3)$: each household's income relative to the average.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 15 · Computing X′y for a large survey\n",
    "\n",
    "A survey has $n = 1000$ households and $k = 3$ regressors, so $X^{\\top}$ is $3 \\times 1000$, $X$ is $1000 \\times 3$ and $y$ is $1000 \\times 1$. Software must compute $X(X^{\\top}y)$ (needed for fitted values). Compare the multiplications needed for $(XX^{\\top})y$ and $X(X^{\\top}y)$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** $(XX^{\\top})y$: forming $XX^{\\top}$ ($1000 \\times 1000$) costs $1000 \\times 1000 \\times 3 = 3{,}000{,}000$ multiplications.\n",
    "\n",
    "**Step 2.** Multiplying that by $y$ costs another $1000 \\times 1000 = 1{,}000{,}000$: total about **4 million**.\n",
    "\n",
    "**Step 3.** $X(X^{\\top}y)$: $X^{\\top}y$ ($3 \\times 1$) costs $3 \\times 1000 = 3{,}000$.\n",
    "\n",
    "**Step 4.** Then $X$ times a $3 \\times 1$ vector costs $1000 \\times 3 = 3{,}000$: total **6,000**.\n",
    "\n",
    "**Answer:** same result, about 700 times cheaper, which is why software never builds an $n \\times n$ matrix when it can avoid it.\n",
    "\n",
    "</details>\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",
    "A = np.array([[1, 2], [3, 4]]); B = np.array([[0, 1], [1, 0]])\n",
    "print(\"AB =\\n\", A @ B, \"\\nBA =\\n\", B @ A)\n",
    "print(\"(AB)' == B'A':\", np.array_equal((A @ B).T, B.T @ A.T), \" tr(AB), tr(BA):\", np.trace(A @ B), np.trace(B @ A))\n",
    "\n",
    "x = np.array([1, 2, 3, 4.]); y = np.array([2, 3, 5, 6.])\n",
    "X = np.column_stack([np.ones(4), x])\n",
    "print(\"X'X =\\n\", X.T @ X, \"\\nX'y =\", X.T @ y)\n",
    "beta = np.array([0.5, 1.4]); e = y - X @ beta\n",
    "print(\"residuals:\", e.round(2), \" X'e =\", (X.T @ e).round(10))\n",
    "\n",
    "n = 3; M = np.eye(n) - np.ones((n, n)) / n\n",
    "print(\"M y =\", M @ np.array([1, 2, 6.]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key takeaways\n",
    "\n",
    "- $Ax$ is a weighted combination of the columns of $A$; $X\\beta$ is the vector of fitted values.\n",
    "- $(AB)^{\\top} = B^{\\top}A^{\\top}$ and $\\operatorname{tr}(AB) = \\operatorname{tr}(BA)$ are used constantly in proofs.\n",
    "- $X^{\\top}X$ and $X^{\\top}y$ hold all the sums of squares and cross-products that regression needs."
   ]
  }
 ]
}