linear_algebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "57e2472d",
   "metadata": {},
   "source": [
    "\n",
    "<a id='linear-algebra'></a>\n",
    "<div id=\"qe-notebook-header\" align=\"right\" style=\"text-align:right;\">\n",
    "        <a href=\"https://quantecon.org/\" title=\"quantecon.org\">\n",
    "                <img style=\"width:250px;display:inline;\" width=\"250px\" src=\"https://assets.quantecon.org/img/qe-menubar-logo.svg\" alt=\"QuantEcon\">\n",
    "        </a>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1acfd898",
   "metadata": {},
   "source": [
    "# Linear Algebra\n",
    "\n",
    "\n",
    "<a id='index-0'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fd9b513",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Linear Algebra](#Linear-Algebra)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Vectors](#Vectors)  \n",
    "  - [Matrices](#Matrices)  \n",
    "  - [Solving Systems of Equations](#Solving-Systems-of-Equations)  \n",
    "  - [Eigenvalues and Eigenvectors](#Eigenvalues-and-Eigenvectors)  \n",
    "  - [Further Topics](#Further-Topics)  \n",
    "  - [Exercises](#Exercises)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b77d5df",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "Linear algebra is one of the most useful branches of applied mathematics for economists to invest in.\n",
    "\n",
    "For example, many applied problems in economics and finance require the solution of a linear system of equations, such as\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    y_1 = a x_1 + b x_2 \\\\\n",
    "    y_2 = c x_1 + d x_2\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "or, more generally,\n",
    "\n",
    "\n",
    "<a id='equation-la-se'></a>\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    y_1 = a_{11} x_1 + a_{12} x_2 + \\cdots + a_{1k} x_k \\\\\n",
    "    \\vdots  \\\\\n",
    "    y_n = a_{n1} x_1 + a_{n2} x_2 + \\cdots + a_{nk} x_k\n",
    "\\end{aligned} \\tag{2.1}\n",
    "$$\n",
    "\n",
    "The objective here is to solve for the “unknowns” $ x_1, \\ldots, x_k $ given $ a_{11}, \\ldots, a_{nk} $ and $ y_1, \\ldots, y_n $.\n",
    "\n",
    "When considering such problems, it is essential that we first consider at least some of the following questions\n",
    "\n",
    "- Does a solution actually exist?  \n",
    "- Are there in fact many solutions, and if so how should we interpret them?  \n",
    "- If no solution exists, is there a best “approximate” solution?  \n",
    "- If a solution exists, how should we compute it?  \n",
    "\n",
    "\n",
    "These are the kinds of topics addressed by linear algebra.\n",
    "\n",
    "In this lecture we will cover the basics of linear and matrix algebra, treating both theory and computation.\n",
    "\n",
    "We admit some overlap with [this lecture](https://python-programming.quantecon.org/numpy.html), where operations on NumPy arrays were first explained.\n",
    "\n",
    "Note that this lecture is more theoretical than most, and contains background\n",
    "material that will be used in applications as we go along.\n",
    "\n",
    "Let’s start with some imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3185f84",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.rcParams[\"figure.figsize\"] = (11, 5)  #set default figure size\n",
    "import numpy as np\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from scipy.linalg import inv, solve, det, eig"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d977698",
   "metadata": {},
   "source": [
    "## Vectors\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "A *vector* of length $ n $ is just a sequence (or array, or tuple) of $ n $ numbers, which we write as $ x = (x_1, \\ldots, x_n) $ or  $ x = [x_1, \\ldots, x_n] $.\n",
    "\n",
    "We will write these sequences either horizontally or vertically as we please.\n",
    "\n",
    "(Later, when we wish to perform certain matrix operations, it will become necessary to distinguish between the two)\n",
    "\n",
    "The set of all $ n $-vectors is denoted by $ \\mathbb R^n $.\n",
    "\n",
    "For example, $ \\mathbb R ^2 $ is the plane, and a vector in $ \\mathbb R^2 $ is just a point in the plane.\n",
    "\n",
    "Traditionally, vectors are represented visually as arrows from the origin to\n",
    "the point.\n",
    "\n",
    "The following figure represents three vectors in this manner"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "429a7b26",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "# Set the axes through the origin\n",
    "for spine in ['left', 'bottom']:\n",
    "    ax.spines[spine].set_position('zero')\n",
    "for spine in ['right', 'top']:\n",
    "    ax.spines[spine].set_color('none')\n",
    "\n",
    "ax.set(xlim=(-5, 5), ylim=(-5, 5))\n",
    "ax.grid()\n",
    "vecs = ((2, 4), (-3, 3), (-4, -3.5))\n",
    "for v in vecs:\n",
    "    ax.annotate('', xy=v, xytext=(0, 0),\n",
    "                arrowprops=dict(facecolor='blue',\n",
    "                shrink=0,\n",
    "                alpha=0.7,\n",
    "                width=0.5))\n",
    "    ax.text(1.1 * v[0], 1.1 * v[1], str(v))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41c1d687",
   "metadata": {},
   "source": [
    "### Vector Operations\n",
    "\n",
    "\n",
    "<a id='index-3'></a>\n",
    "The two most common operators for vectors are addition and scalar multiplication, which we now describe.\n",
    "\n",
    "As a matter of definition, when we add two vectors, we add them element-by-element\n",
    "\n",
    "$$\n",
    "x + y =\n",
    "\\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    x_n\n",
    "\\end{bmatrix} +\n",
    "\\begin{bmatrix}\n",
    "     y_1 \\\\\n",
    "     y_2 \\\\\n",
    "    \\vdots \\\\\n",
    "     y_n\n",
    "\\end{bmatrix} :=\n",
    "\\begin{bmatrix}\n",
    "    x_1 + y_1 \\\\\n",
    "    x_2 + y_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    x_n + y_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Scalar multiplication is an operation that takes a number $ \\gamma $ and a\n",
    "vector $ x $ and produces\n",
    "\n",
    "$$\n",
    "\\gamma x :=\n",
    "\\begin{bmatrix}\n",
    "    \\gamma x_1 \\\\\n",
    "    \\gamma x_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    \\gamma x_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Scalar multiplication is illustrated in the next figure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e007c4a",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "# Set the axes through the origin\n",
    "for spine in ['left', 'bottom']:\n",
    "    ax.spines[spine].set_position('zero')\n",
    "for spine in ['right', 'top']:\n",
    "    ax.spines[spine].set_color('none')\n",
    "\n",
    "ax.set(xlim=(-5, 5), ylim=(-5, 5))\n",
    "x = (2, 2)\n",
    "ax.annotate('', xy=x, xytext=(0, 0),\n",
    "            arrowprops=dict(facecolor='blue',\n",
    "            shrink=0,\n",
    "            alpha=1,\n",
    "            width=0.5))\n",
    "ax.text(x[0] + 0.4, x[1] - 0.2, '$x$', fontsize='16')\n",
    "\n",
    "\n",
    "scalars = (-2, 2)\n",
    "x = np.array(x)\n",
    "\n",
    "for s in scalars:\n",
    "    v = s * x\n",
    "    ax.annotate('', xy=v, xytext=(0, 0),\n",
    "                arrowprops=dict(facecolor='red',\n",
    "                shrink=0,\n",
    "                alpha=0.5,\n",
    "                width=0.5))\n",
    "    ax.text(v[0] + 0.4, v[1] - 0.2, f'${s} x$', fontsize='16')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2a100d9",
   "metadata": {},
   "source": [
    "In Python, a vector can be represented as a list or tuple, such as `x = (2, 4, 6)`, but is more commonly\n",
    "represented as a [NumPy array](https://python-programming.quantecon.org/numpy.html#numpy-arrays).\n",
    "\n",
    "One advantage of NumPy arrays is that scalar multiplication and addition have very natural syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9628f8bb",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x = np.ones(3)            # Vector of three ones\n",
    "y = np.array((2, 4, 6))   # Converts tuple (2, 4, 6) into array\n",
    "x + y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27852495",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "4 * x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d9d98a3",
   "metadata": {},
   "source": [
    "### Inner Product and Norm\n",
    "\n",
    "\n",
    "<a id='index-5'></a>\n",
    "The *inner product* of vectors $ x,y \\in \\mathbb R ^n $ is defined as\n",
    "\n",
    "$$\n",
    "x' y := \\sum_{i=1}^n x_i y_i\n",
    "$$\n",
    "\n",
    "Two vectors are called *orthogonal* if their inner product is zero.\n",
    "\n",
    "The *norm* of a vector $ x $ represents its “length” (i.e., its distance from the zero vector) and is defined as\n",
    "\n",
    "$$\n",
    "\\| x \\| := \\sqrt{x' x} := \\left( \\sum_{i=1}^n x_i^2 \\right)^{1/2}\n",
    "$$\n",
    "\n",
    "The expression $ \\| x - y\\| $ is thought of as the distance between $ x $ and $ y $.\n",
    "\n",
    "Continuing on from the previous example, the inner product and norm can be computed as\n",
    "follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6c553e1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.sum(x * y)          # Inner product of x and y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2fb8cd3",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.sqrt(np.sum(x**2))  # Norm of x, take one"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be2b9d83",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "np.linalg.norm(x)      # Norm of x, take two"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43f1a801",
   "metadata": {},
   "source": [
    "### Span\n",
    "\n",
    "\n",
    "<a id='index-6'></a>\n",
    "Given a set of vectors $ A := \\{a_1, \\ldots, a_k\\} $ in $ \\mathbb R ^n $, it’s natural to think about the new vectors we can create by performing linear operations.\n",
    "\n",
    "New vectors created in this manner are called *linear combinations* of $ A $.\n",
    "\n",
    "In particular, $ y \\in \\mathbb R ^n $ is a linear combination of $ A := \\{a_1, \\ldots, a_k\\} $ if\n",
    "\n",
    "$$\n",
    "y = \\beta_1 a_1 + \\cdots + \\beta_k a_k\n",
    "\\text{ for some scalars } \\beta_1, \\ldots, \\beta_k\n",
    "$$\n",
    "\n",
    "In this context, the values $ \\beta_1, \\ldots, \\beta_k $ are called the *coefficients* of the linear combination.\n",
    "\n",
    "The set of linear combinations of $ A $ is called the *span* of $ A $.\n",
    "\n",
    "The next figure shows the span of $ A = \\{a_1, a_2\\} $ in $ \\mathbb R ^3 $.\n",
    "\n",
    "The span is a two-dimensional plane passing through these two points and the origin.\n",
    "\n",
    "\n",
    "<a id='la-3dvec'></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9179576",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "ax = plt.figure(figsize=(10, 8)).add_subplot(projection='3d')\n",
    "\n",
    "x_min, x_max = -5, 5\n",
    "y_min, y_max = -5, 5\n",
    "\n",
    "α, β = 0.2, 0.1\n",
    "\n",
    "ax.set(xlim=(x_min, x_max), ylim=(x_min, x_max), zlim=(x_min, x_max),\n",
    "       xticks=(0,), yticks=(0,), zticks=(0,))\n",
    "\n",
    "gs = 3\n",
    "z = np.linspace(x_min, x_max, gs)\n",
    "x = np.zeros(gs)\n",
    "y = np.zeros(gs)\n",
    "ax.plot(x, y, z, 'k-', lw=2, alpha=0.5)\n",
    "ax.plot(z, x, y, 'k-', lw=2, alpha=0.5)\n",
    "ax.plot(y, z, x, 'k-', lw=2, alpha=0.5)\n",
    "\n",
    "\n",
    "# Fixed linear function, to generate a plane\n",
    "def f(x, y):\n",
    "    return α * x + β * y\n",
    "\n",
    "# Vector locations, by coordinate\n",
    "x_coords = np.array((3, 3))\n",
    "y_coords = np.array((4, -4))\n",
    "z = f(x_coords, y_coords)\n",
    "for i in (0, 1):\n",
    "    ax.text(x_coords[i], y_coords[i], z[i], f'$a_{i+1}$', fontsize=14)\n",
    "\n",
    "# Lines to vectors\n",
    "for i in (0, 1):\n",
    "    x = (0, x_coords[i])\n",
    "    y = (0, y_coords[i])\n",
    "    z = (0, f(x_coords[i], y_coords[i]))\n",
    "    ax.plot(x, y, z, 'b-', lw=1.5, alpha=0.6)\n",
    "\n",
    "\n",
    "# Draw the plane\n",
    "grid_size = 20\n",
    "xr2 = np.linspace(x_min, x_max, grid_size)\n",
    "yr2 = np.linspace(y_min, y_max, grid_size)\n",
    "x2, y2 = np.meshgrid(xr2, yr2)\n",
    "z2 = f(x2, y2)\n",
    "ax.plot_surface(x2, y2, z2, rstride=1, cstride=1, cmap=cm.jet,\n",
    "                linewidth=0, antialiased=True, alpha=0.2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d68624c1",
   "metadata": {},
   "source": [
    "#### Examples\n",
    "\n",
    "If $ A $ contains only one vector $ a_1 \\in \\mathbb R ^2 $, then its\n",
    "span is just the scalar multiples of $ a_1 $, which is the unique line passing through both $ a_1 $ and the origin.\n",
    "\n",
    "If $ A = \\{e_1, e_2, e_3\\} $ consists  of the *canonical basis vectors* of $ \\mathbb R ^3 $, that is\n",
    "\n",
    "$$\n",
    "e_1 :=\n",
    "\\begin{bmatrix}\n",
    "     1 \\\\\n",
    "     0 \\\\\n",
    "     0\n",
    "\\end{bmatrix}\n",
    ", \\quad\n",
    "e_2 :=\n",
    "\\begin{bmatrix}\n",
    "     0 \\\\\n",
    "     1 \\\\\n",
    "     0\n",
    "\\end{bmatrix}\n",
    ", \\quad\n",
    "e_3 :=\n",
    "\\begin{bmatrix}\n",
    "     0 \\\\\n",
    "     0 \\\\\n",
    "     1\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "then the span of $ A $ is all of $ \\mathbb R ^3 $, because, for any\n",
    "$ x = (x_1, x_2, x_3) \\in \\mathbb R ^3 $, we can write\n",
    "\n",
    "$$\n",
    "x = x_1 e_1 + x_2 e_2 + x_3 e_3\n",
    "$$\n",
    "\n",
    "Now consider $ A_0 = \\{e_1, e_2, e_1 + e_2\\} $.\n",
    "\n",
    "If $ y = (y_1, y_2, y_3) $ is any linear combination of these vectors, then $ y_3 = 0 $ (check it).\n",
    "\n",
    "Hence $ A_0 $ fails to span all of $ \\mathbb R ^3 $.\n",
    "\n",
    "\n",
    "<a id='la-li'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7690375b",
   "metadata": {},
   "source": [
    "### Linear Independence\n",
    "\n",
    "\n",
    "<a id='index-7'></a>\n",
    "As we’ll see, it’s often desirable to find families of vectors with relatively large span, so that many vectors can be described by linear operators on a few vectors.\n",
    "\n",
    "The condition we need for a set of vectors to have a large span is what’s called linear independence.\n",
    "\n",
    "In particular, a collection of vectors $ A := \\{a_1, \\ldots, a_k\\} $ in $ \\mathbb R ^n $ is said to be\n",
    "\n",
    "- *linearly dependent* if some strict subset of $ A $ has the same span as $ A $.  \n",
    "- *linearly independent* if it is not linearly dependent.  \n",
    "\n",
    "\n",
    "Put differently, a set of vectors is linearly independent if no vector is redundant to the span and linearly dependent otherwise.\n",
    "\n",
    "To illustrate the idea, recall [the figure](#la-3dvec) that showed the span of vectors $ \\{a_1, a_2\\} $ in $ \\mathbb R ^3 $ as a plane through the origin.\n",
    "\n",
    "If we take a third vector $ a_3 $ and form the set $ \\{a_1, a_2, a_3\\} $, this set will be\n",
    "\n",
    "- linearly dependent if $ a_3 $ lies in the plane  \n",
    "- linearly independent otherwise  \n",
    "\n",
    "\n",
    "As another illustration of the concept, since $ \\mathbb R ^n $ can be spanned by $ n $ vectors\n",
    "(see the discussion of canonical basis vectors above), any collection of\n",
    "$ m > n $ vectors in $ \\mathbb R ^n $ must be linearly dependent.\n",
    "\n",
    "The following statements are equivalent to linear independence of $ A := \\{a_1, \\ldots, a_k\\} \\subset \\mathbb R ^n $\n",
    "\n",
    "1. No vector in $ A $ can be formed as a linear combination of the other elements.  \n",
    "1. If $ \\beta_1 a_1 + \\cdots \\beta_k a_k = 0 $ for scalars $ \\beta_1, \\ldots, \\beta_k $, then $ \\beta_1 = \\cdots = \\beta_k = 0 $.  \n",
    "\n",
    "\n",
    "(The zero in the first expression is the origin of $ \\mathbb R ^n $)\n",
    "\n",
    "\n",
    "<a id='la-unique-reps'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f55f659f",
   "metadata": {},
   "source": [
    "### Unique Representations\n",
    "\n",
    "Another nice thing about sets of linearly independent vectors is that each element in the span has a unique representation as a linear combination of these vectors.\n",
    "\n",
    "In other words, if $ A := \\{a_1, \\ldots, a_k\\} \\subset \\mathbb R ^n $ is\n",
    "linearly independent and\n",
    "\n",
    "$$\n",
    "y = \\beta_1 a_1 + \\cdots \\beta_k a_k\n",
    "$$\n",
    "\n",
    "then no other coefficient sequence $ \\gamma_1, \\ldots, \\gamma_k $ will produce\n",
    "the same vector $ y $.\n",
    "\n",
    "Indeed, if we also have $ y = \\gamma_1 a_1 + \\cdots \\gamma_k a_k $,\n",
    "then\n",
    "\n",
    "$$\n",
    "(\\beta_1 - \\gamma_1) a_1 + \\cdots + (\\beta_k - \\gamma_k) a_k = 0\n",
    "$$\n",
    "\n",
    "Linear independence now implies $ \\gamma_i = \\beta_i $ for all $ i $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e083e858",
   "metadata": {},
   "source": [
    "## Matrices\n",
    "\n",
    "\n",
    "<a id='index-8'></a>\n",
    "Matrices are a neat way of organizing data for use in linear operations.\n",
    "\n",
    "An $ n \\times k $ matrix is a rectangular array $ A $ of numbers with $ n $ rows and $ k $ columns:\n",
    "\n",
    "$$\n",
    "A =\n",
    "\\begin{bmatrix}\n",
    "    a_{11} & a_{12} & \\cdots & a_{1k} \\\\\n",
    "    a_{21} & a_{22} & \\cdots & a_{2k} \\\\\n",
    "    \\vdots & \\vdots &  & \\vdots \\\\\n",
    "    a_{n1} & a_{n2} & \\cdots & a_{nk}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Often, the numbers in the matrix represent coefficients in a system of linear equations, as discussed at the start of this lecture.\n",
    "\n",
    "For obvious reasons, the matrix $ A $ is also called a vector if either $ n = 1 $ or $ k = 1 $.\n",
    "\n",
    "In the former case, $ A $ is called a *row vector*, while in the latter it is called a *column vector*.\n",
    "\n",
    "If $ n = k $, then $ A $ is called *square*.\n",
    "\n",
    "The matrix formed by replacing $ a_{ij} $ by $ a_{ji} $ for every $ i $ and $ j $ is called the *transpose* of $ A $ and denoted $ A' $ or $ A^{\\top} $.\n",
    "\n",
    "If $ A = A' $, then $ A $ is called *symmetric*.\n",
    "\n",
    "For a square matrix $ A $, the $ i $ elements of the form $ a_{ii} $ for $ i=1,\\ldots,n $ are called the *principal diagonal*.\n",
    "\n",
    "$ A $ is called *diagonal* if the only nonzero entries are on the principal diagonal.\n",
    "\n",
    "If, in addition to being diagonal, each element along the principal diagonal is equal to 1, then $ A $ is called the *identity matrix* and denoted by $ I $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdbc1264",
   "metadata": {},
   "source": [
    "### Matrix Operations\n",
    "\n",
    "\n",
    "<a id='index-9'></a>\n",
    "Just as was the case for vectors, a number of algebraic operations are defined for matrices.\n",
    "\n",
    "Scalar multiplication and addition are immediate generalizations of the vector case:\n",
    "\n",
    "$$\n",
    "\\gamma A =\n",
    "\\gamma\n",
    "\\begin{bmatrix}\n",
    "    a_{11} &  \\cdots & a_{1k} \\\\\n",
    "    \\vdots & \\vdots  & \\vdots \\\\\n",
    "    a_{n1} &  \\cdots & a_{nk}\n",
    "\\end{bmatrix} :=\n",
    "\\begin{bmatrix}\n",
    "    \\gamma a_{11} & \\cdots & \\gamma a_{1k} \\\\\n",
    "    \\vdots & \\vdots & \\vdots \\\\\n",
    "    \\gamma a_{n1} & \\cdots & \\gamma a_{nk}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "and\n",
    "\n",
    "$$\n",
    "A + B =\n",
    "\\begin{bmatrix}\n",
    "    a_{11} & \\cdots & a_{1k} \\\\\n",
    "    \\vdots & \\vdots & \\vdots \\\\\n",
    "    a_{n1} & \\cdots & a_{nk}\n",
    "\\end{bmatrix} +\n",
    "\\begin{bmatrix}\n",
    "    b_{11} & \\cdots & b_{1k} \\\\\n",
    "    \\vdots & \\vdots & \\vdots \\\\\n",
    "    b_{n1} & \\cdots & b_{nk}\n",
    "\\end{bmatrix} :=\n",
    "\\begin{bmatrix}\n",
    "    a_{11} + b_{11} &  \\cdots & a_{1k} + b_{1k} \\\\\n",
    "    \\vdots & \\vdots & \\vdots \\\\\n",
    "    a_{n1} + b_{n1} &  \\cdots & a_{nk} + b_{nk}\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "In the latter case, the matrices must have the same shape in order for the definition to make sense.\n",
    "\n",
    "We also have a convention for *multiplying* two matrices.\n",
    "\n",
    "The rule for matrix multiplication generalizes the idea of inner products discussed above\n",
    "and is designed to make multiplication play well with basic linear operations.\n",
    "\n",
    "If $ A $ and $ B $ are two matrices, then their product $ A B $ is formed by taking as its\n",
    "$ i,j $-th element the inner product of the $ i $-th row of $ A $ and the\n",
    "$ j $-th column of $ B $.\n",
    "\n",
    "There are many tutorials to help you visualize this operation, such as [this one](https://www.mathsisfun.com/algebra/matrix-multiplying.html), or the discussion on the [Wikipedia page](https://en.wikipedia.org/wiki/Matrix_multiplication).\n",
    "\n",
    "If $ A $ is $ n \\times k $ and $ B $ is $ j \\times m $, then\n",
    "to multiply $ A $ and $ B $ we require $ k = j $, and the\n",
    "resulting matrix $ A B $ is $ n \\times m $.\n",
    "\n",
    "As perhaps the most important special case, consider multiplying $ n \\times k $ matrix $ A $ and $ k \\times 1 $ column vector $ x $.\n",
    "\n",
    "According to the preceding rule, this gives us an $ n \\times 1 $ column vector\n",
    "\n",
    "\n",
    "<a id='equation-la-atx'></a>\n",
    "$$\n",
    "A x =\n",
    "\\begin{bmatrix}\n",
    "    a_{11} &  \\cdots & a_{1k} \\\\\n",
    "    \\vdots & \\vdots  & \\vdots \\\\\n",
    "    a_{n1} &  \\cdots & a_{nk}\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "    x_{1}  \\\\\n",
    "    \\vdots  \\\\\n",
    "    x_{k}\n",
    "\\end{bmatrix} :=\n",
    "\\begin{bmatrix}\n",
    "    a_{11} x_1 + \\cdots + a_{1k} x_k \\\\\n",
    "    \\vdots \\\\\n",
    "    a_{n1} x_1 + \\cdots + a_{nk} x_k\n",
    "\\end{bmatrix} \\tag{2.2}\n",
    "$$\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">$ A B $ and $ B A $ are not generally the same thing.\n",
    "\n",
    "Another important special case is the identity matrix.\n",
    "\n",
    "You should check that if $ A $ is $ n \\times k $ and $ I $ is the $ k \\times k $ identity matrix, then $ AI = A $.\n",
    "\n",
    "If $ I $ is the $ n \\times n $ identity matrix, then $ IA = A $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45601955",
   "metadata": {},
   "source": [
    "### Matrices in NumPy\n",
    "\n",
    "\n",
    "<a id='index-10'></a>\n",
    "NumPy arrays are also used as matrices, and have fast, efficient functions and methods for all the standard matrix operations <sup><a href=#fn-mdt id=fn-mdt-link>[1]</a></sup>.\n",
    "\n",
    "You can create them manually from tuples of tuples (or lists of lists) as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d55f95d7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = ((1, 2),\n",
    "     (3, 4))\n",
    "\n",
    "type(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e633d59f",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = np.array(A)\n",
    "\n",
    "type(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "177d3f27",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "256bc68a",
   "metadata": {},
   "source": [
    "The `shape` attribute is a tuple giving the number of rows and columns —\n",
    "see [here](https://python-programming.quantecon.org/numpy.html#shape-and-dimension)\n",
    "for more discussion.\n",
    "\n",
    "To get the transpose of `A`, use `A.transpose()` or, more simply, `A.T`.\n",
    "\n",
    "There are many convenient functions for creating common matrices (matrices of zeros,\n",
    "ones, etc.) — see [here](https://python-programming.quantecon.org/numpy.html#creating-arrays).\n",
    "\n",
    "Since operations are performed elementwise by default, scalar multiplication and addition have very natural syntax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a074ce8b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = np.identity(3)\n",
    "B = np.ones((3, 3))\n",
    "2 * A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29a300a1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A + B"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b92ff63f",
   "metadata": {},
   "source": [
    "To multiply matrices we use the `@` symbol.\n",
    "\n",
    "In particular, `A @ B` is matrix multiplication, whereas `A * B` is element-by-element multiplication.\n",
    "\n",
    "See [here](https://python-programming.quantecon.org/numpy.html#matrix-multiplication) for more discussion.\n",
    "\n",
    "\n",
    "<a id='la-linear-map'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "526de3e6",
   "metadata": {},
   "source": [
    "### Matrices as Maps\n",
    "\n",
    "\n",
    "<a id='index-11'></a>\n",
    "Each $ n \\times k $ matrix $ A $ can be identified with a function $ f(x) = Ax $ that maps $ x \\in \\mathbb R ^k $ into $ y = Ax \\in \\mathbb R ^n $.\n",
    "\n",
    "These kinds of functions have a special property: they are *linear*.\n",
    "\n",
    "A function $ f \\colon \\mathbb R ^k \\to \\mathbb R ^n $ is called *linear* if, for all $ x, y \\in \\mathbb R ^k $ and all scalars $ \\alpha, \\beta $, we have\n",
    "\n",
    "$$\n",
    "f(\\alpha x + \\beta y) = \\alpha f(x) + \\beta f(y)\n",
    "$$\n",
    "\n",
    "You can check that this holds for the function $ f(x) = A x + b $ when $ b $ is the zero vector and fails when $ b $ is nonzero.\n",
    "\n",
    "In fact, it’s [known](https://en.wikipedia.org/wiki/Linear_map#Matrices) that $ f $ is linear if and *only if* there exists a matrix $ A $ such that $ f(x) = Ax $ for all $ x $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22efbf0e",
   "metadata": {},
   "source": [
    "## Solving Systems of Equations\n",
    "\n",
    "\n",
    "<a id='index-12'></a>\n",
    "Recall again the system of equations [(2.1)](#equation-la-se).\n",
    "\n",
    "If we compare [(2.1)](#equation-la-se) and [(2.2)](#equation-la-atx), we see that [(2.1)](#equation-la-se) can now be\n",
    "written more conveniently as\n",
    "\n",
    "\n",
    "<a id='equation-la-se2'></a>\n",
    "$$\n",
    "y = Ax \\tag{2.3}\n",
    "$$\n",
    "\n",
    "The problem we face is to determine a vector $ x \\in \\mathbb R ^k $ that solves [(2.3)](#equation-la-se2), taking $ y $ and $ A $ as given.\n",
    "\n",
    "This is a special case of a more general problem: Find an $ x $ such that $ y = f(x) $.\n",
    "\n",
    "Given an arbitrary function $ f $ and a $ y $, is there always an $ x $ such that $ y = f(x) $?\n",
    "\n",
    "If so, is it always unique?\n",
    "\n",
    "The answer to both these questions is negative, as the next figure shows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fc6b299",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    return 0.6 * np.cos(4 * x) + 1.4\n",
    "\n",
    "\n",
    "xmin, xmax = -1, 1\n",
    "x = np.linspace(xmin, xmax, 160)\n",
    "y = f(x)\n",
    "ya, yb = np.min(y), np.max(y)\n",
    "\n",
    "fig, axes = plt.subplots(2, 1, figsize=(10, 10))\n",
    "\n",
    "for ax in axes:\n",
    "    # Set the axes through the origin\n",
    "    for spine in ['left', 'bottom']:\n",
    "        ax.spines[spine].set_position('zero')\n",
    "    for spine in ['right', 'top']:\n",
    "        ax.spines[spine].set_color('none')\n",
    "\n",
    "    ax.set(ylim=(-0.6, 3.2), xlim=(xmin, xmax),\n",
    "           yticks=(), xticks=())\n",
    "\n",
    "    ax.plot(x, y, 'k-', lw=2, label='$f$')\n",
    "    ax.fill_between(x, ya, yb, facecolor='blue', alpha=0.05)\n",
    "    ax.vlines([0], ya, yb, lw=3, color='blue', label='range of $f$')\n",
    "    ax.text(0.04, -0.3, '$0$', fontsize=16)\n",
    "\n",
    "ax = axes[0]\n",
    "\n",
    "ax.legend(loc='upper right', frameon=False)\n",
    "ybar = 1.5\n",
    "ax.plot(x, x * 0 + ybar, 'k--', alpha=0.5)\n",
    "ax.text(0.05, 0.8 * ybar, '$y$', fontsize=16)\n",
    "for i, z in enumerate((-0.35, 0.35)):\n",
    "    ax.vlines(z, 0, f(z), linestyle='--', alpha=0.5)\n",
    "    ax.text(z, -0.2, f'$x_{i}$', fontsize=16)\n",
    "\n",
    "ax = axes[1]\n",
    "\n",
    "ybar = 2.6\n",
    "ax.plot(x, x * 0 + ybar, 'k--', alpha=0.5)\n",
    "ax.text(0.04, 0.91 * ybar, '$y$', fontsize=16)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92ee93a4",
   "metadata": {},
   "source": [
    "In the first plot, there are multiple solutions, as the function is not one-to-one, while\n",
    "in the second there are no solutions, since $ y $ lies outside the range of $ f $.\n",
    "\n",
    "Can we impose conditions on $ A $ in [(2.3)](#equation-la-se2) that rule out these problems?\n",
    "\n",
    "In this context, the most important thing to recognize about the expression\n",
    "$ Ax $ is that it corresponds to a linear combination of the columns of $ A $.\n",
    "\n",
    "In particular, if $ a_1, \\ldots, a_k $ are the columns of $ A $, then\n",
    "\n",
    "$$\n",
    "Ax = x_1 a_1 + \\cdots + x_k a_k\n",
    "$$\n",
    "\n",
    "Hence the range of $ f(x) = Ax $ is exactly the span of the columns of $ A $.\n",
    "\n",
    "We want the range to be large so that it contains arbitrary $ y $.\n",
    "\n",
    "As you might recall, the condition that we want for the span to be large is [linear independence](#la-li).\n",
    "\n",
    "A happy fact is that linear independence of the columns of $ A $ also gives us uniqueness.\n",
    "\n",
    "Indeed, it follows from our [earlier discussion](#la-unique-reps) that if $ \\{a_1, \\ldots, a_k\\} $ are linearly independent and $ y = Ax = x_1 a_1 + \\cdots + x_k a_k $, then no $ z \\not= x $ satisfies $ y = Az $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5e51c2b",
   "metadata": {},
   "source": [
    "### The Square Matrix Case\n",
    "\n",
    "Let’s discuss some more details, starting with the case where $ A $ is $ n \\times n $.\n",
    "\n",
    "This is the familiar case where the number of unknowns equals the number of equations.\n",
    "\n",
    "For arbitrary $ y \\in \\mathbb R ^n $, we hope to find a unique $ x \\in \\mathbb R ^n $ such that $ y = Ax $.\n",
    "\n",
    "In view of the observations immediately above, if the columns of $ A $ are\n",
    "linearly independent, then their span, and hence the range of $ f(x) =\n",
    "Ax $, is all of $ \\mathbb R ^n $.\n",
    "\n",
    "Hence there always exists an $ x $ such that $ y = Ax $.\n",
    "\n",
    "Moreover, the solution is unique.\n",
    "\n",
    "In particular, the following are equivalent\n",
    "\n",
    "1. The columns of $ A $ are linearly independent.  \n",
    "1. For any $ y \\in \\mathbb R ^n $, the equation $ y = Ax $ has a unique solution.  \n",
    "\n",
    "\n",
    "The property of having linearly independent columns is sometimes expressed as having *full column rank*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42037c84",
   "metadata": {},
   "source": [
    "#### Inverse Matrices\n",
    "\n",
    "\n",
    "<a id='index-13'></a>\n",
    "Can we give some sort of expression for the solution?\n",
    "\n",
    "If $ y $ and $ A $ are scalar with $ A \\not= 0 $, then the\n",
    "solution is $ x = A^{-1} y $.\n",
    "\n",
    "A similar expression is available in the matrix case.\n",
    "\n",
    "In particular, if square matrix $ A $ has full column rank, then it possesses a multiplicative\n",
    "*inverse matrix* $ A^{-1} $, with the property that $ A A^{-1} = A^{-1} A = I $.\n",
    "\n",
    "As a consequence, if we pre-multiply both sides of $ y = Ax $ by $ A^{-1} $, we get $ x = A^{-1} y $.\n",
    "\n",
    "This is the solution that we’re looking for."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9a6cd63",
   "metadata": {},
   "source": [
    "#### Determinants\n",
    "\n",
    "\n",
    "<a id='index-14'></a>\n",
    "Another quick comment about square matrices is that to every such matrix we\n",
    "assign a unique number called the *determinant* of the matrix — you can find\n",
    "the expression for it [here](https://en.wikipedia.org/wiki/Determinant).\n",
    "\n",
    "If the determinant of $ A $ is not zero, then we say that $ A $ is\n",
    "*nonsingular*.\n",
    "\n",
    "Perhaps the most important fact about determinants is that $ A $ is nonsingular if and only if $ A $ is of full column rank.\n",
    "\n",
    "This gives us a useful one-number summary of whether or not a square matrix can be\n",
    "inverted."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf22b7e7",
   "metadata": {},
   "source": [
    "### More Rows than Columns\n",
    "\n",
    "This is the $ n \\times k $ case with $ n > k $.\n",
    "\n",
    "This case is very important in many settings, not least in the setting of linear regression (where $ n $ is the number of observations, and $ k $ is the number of explanatory variables).\n",
    "\n",
    "Given arbitrary $ y \\in \\mathbb R ^n $, we seek an $ x \\in \\mathbb R ^k $ such that $ y = Ax $.\n",
    "\n",
    "In this setting, the existence of a solution is highly unlikely.\n",
    "\n",
    "Without much loss of generality, let’s go over the intuition focusing on the case where the columns of\n",
    "$ A $ are linearly independent.\n",
    "\n",
    "It follows that the span of the columns of $ A $ is a $ k $-dimensional subspace of $ \\mathbb R ^n $.\n",
    "\n",
    "This span is very “unlikely” to contain arbitrary $ y \\in \\mathbb R ^n $.\n",
    "\n",
    "To see why, recall the [figure above](#la-3dvec), where $ k=2 $ and $ n=3 $.\n",
    "\n",
    "Imagine an arbitrarily chosen $ y \\in \\mathbb R ^3 $, located somewhere in that three-dimensional space.\n",
    "\n",
    "What’s the likelihood that $ y $ lies in the span of $ \\{a_1, a_2\\} $ (i.e., the two dimensional plane through these points)?\n",
    "\n",
    "In a sense, it must be very small, since this plane has zero “thickness”.\n",
    "\n",
    "As a result, in the $ n > k $ case we usually give up on existence.\n",
    "\n",
    "However, we can still seek the best approximation, for example, an\n",
    "$ x $ that makes the distance $ \\| y - Ax\\| $ as small as possible.\n",
    "\n",
    "To solve this problem, one can use either calculus or the theory of orthogonal\n",
    "projections.\n",
    "\n",
    "The solution is known to be $ \\hat x = (A'A)^{-1}A'y $ — see for example chapter 3 of [these notes](https://python.quantecon.org/_static/lecture_specific/linear_algebra/course_notes.pdf)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8d4ead0",
   "metadata": {},
   "source": [
    "### More Columns than Rows\n",
    "\n",
    "This is the $ n \\times k $ case with $ n < k $, so there are fewer\n",
    "equations than unknowns.\n",
    "\n",
    "In this case there are either no solutions or infinitely many — in other words, uniqueness never holds.\n",
    "\n",
    "For example, consider the case where $ k=3 $ and $ n=2 $.\n",
    "\n",
    "Thus, the columns of $ A $ consists of 3 vectors in $ \\mathbb R ^2 $.\n",
    "\n",
    "This set can never be linearly independent, since it is possible to find two vectors that span\n",
    "$ \\mathbb R ^2 $.\n",
    "\n",
    "(For example, use the canonical basis vectors)\n",
    "\n",
    "It follows that one column is a linear combination of the other two.\n",
    "\n",
    "For example, let’s say that $ a_1 = \\alpha a_2 + \\beta a_3 $.\n",
    "\n",
    "Then if $ y = Ax = x_1 a_1 + x_2 a_2 + x_3 a_3 $, we can also write\n",
    "\n",
    "$$\n",
    "y\n",
    "= x_1 (\\alpha a_2 + \\beta a_3) + x_2 a_2 + x_3 a_3\n",
    "= (x_1 \\alpha + x_2) a_2 + (x_1 \\beta + x_3) a_3\n",
    "$$\n",
    "\n",
    "In other words, uniqueness fails."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ec462fdd",
   "metadata": {},
   "source": [
    "### Linear Equations with SciPy\n",
    "\n",
    "\n",
    "<a id='index-15'></a>\n",
    "Here’s an illustration of how to solve linear equations with SciPy’s `linalg` submodule.\n",
    "\n",
    "All of these routines are Python front ends to time-tested and highly optimized FORTRAN code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b07d483",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = ((1, 2), (3, 4))\n",
    "A = np.array(A)\n",
    "y = np.ones((2, 1))  # Column vector\n",
    "det(A)  # Check that A is nonsingular, and hence invertible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e05cfa07",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A_inv = inv(A)  # Compute the inverse\n",
    "A_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7afee3e7",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "x = A_inv @ y  # Solution\n",
    "A @ x          # Should equal y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05e2a9dd",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "solve(A, y)  # Produces the same solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ccf3869",
   "metadata": {},
   "source": [
    "Observe how we can solve for $ x = A^{-1} y $ by either via `inv(A) @ y`, or using `solve(A, y)`.\n",
    "\n",
    "The latter method uses a different algorithm (LU decomposition) that is numerically more stable, and hence should almost always be preferred.\n",
    "\n",
    "To obtain the least-squares solution $ \\hat x = (A'A)^{-1}A'y $, use `scipy.linalg.lstsq(A, y)`.\n",
    "\n",
    "\n",
    "<a id='la-eigen'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35008e78",
   "metadata": {},
   "source": [
    "## Eigenvalues and Eigenvectors\n",
    "\n",
    "\n",
    "<a id='index-19'></a>\n",
    "Let $ A $ be an $ n \\times n $ square matrix.\n",
    "\n",
    "If $ \\lambda $ is scalar and $ v $ is a non-zero vector in $ \\mathbb R ^n $ such that\n",
    "\n",
    "$$\n",
    "A v = \\lambda v\n",
    "$$\n",
    "\n",
    "then we say that $ \\lambda $ is an *eigenvalue* of $ A $, and\n",
    "$ v $ is an *eigenvector*.\n",
    "\n",
    "Thus, an eigenvector of $ A $ is a vector such that when the map $ f(x) = Ax $ is applied, $ v $ is merely scaled.\n",
    "\n",
    "The next figure shows two eigenvectors (blue arrows) and their images under $ A $ (red arrows).\n",
    "\n",
    "As expected, the image $ Av $ of each $ v $ is just a scaled version of the original"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a90f04cc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = ((1, 2),\n",
    "     (2, 1))\n",
    "A = np.array(A)\n",
    "evals, evecs = eig(A)\n",
    "evecs = evecs[:, 0], evecs[:, 1]\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 8))\n",
    "# Set the axes through the origin\n",
    "for spine in ['left', 'bottom']:\n",
    "    ax.spines[spine].set_position('zero')\n",
    "for spine in ['right', 'top']:\n",
    "    ax.spines[spine].set_color('none')\n",
    "ax.grid(alpha=0.4)\n",
    "\n",
    "xmin, xmax = -3, 3\n",
    "ymin, ymax = -3, 3\n",
    "ax.set(xlim=(xmin, xmax), ylim=(ymin, ymax))\n",
    "\n",
    "# Plot each eigenvector\n",
    "for v in evecs:\n",
    "    ax.annotate('', xy=v, xytext=(0, 0),\n",
    "                arrowprops=dict(facecolor='blue',\n",
    "                shrink=0,\n",
    "                alpha=0.6,\n",
    "                width=0.5))\n",
    "\n",
    "# Plot the image of each eigenvector\n",
    "for v in evecs:\n",
    "    v = A @ v\n",
    "    ax.annotate('', xy=v, xytext=(0, 0),\n",
    "                arrowprops=dict(facecolor='red',\n",
    "                shrink=0,\n",
    "                alpha=0.6,\n",
    "                width=0.5))\n",
    "\n",
    "# Plot the lines they run through\n",
    "x = np.linspace(xmin, xmax, 3)\n",
    "for v in evecs:\n",
    "    a = v[1] / v[0]\n",
    "    ax.plot(x, a * x, 'b-', lw=0.4)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2a7851d3",
   "metadata": {},
   "source": [
    "The eigenvalue equation is equivalent to $ (A - \\lambda I) v = 0 $, and\n",
    "this has a nonzero solution $ v $ only when the columns of $ A -\n",
    "\\lambda I $ are linearly dependent.\n",
    "\n",
    "This in turn is equivalent to stating that the determinant is zero.\n",
    "\n",
    "Hence to find all eigenvalues, we can look for $ \\lambda $ such that the\n",
    "determinant of $ A - \\lambda I $ is zero.\n",
    "\n",
    "This problem can be expressed as one of solving for the roots of a polynomial\n",
    "in $ \\lambda $ of degree $ n $.\n",
    "\n",
    "This in turn implies the existence of $ n $ solutions in the complex\n",
    "plane, although some might be repeated.\n",
    "\n",
    "Some nice facts about the eigenvalues of a square matrix $ A $ are as follows\n",
    "\n",
    "1. The determinant of $ A $ equals  the product of the eigenvalues.  \n",
    "1. The trace of $ A $ (the sum of the elements on the principal diagonal) equals the sum of the eigenvalues.  \n",
    "1. If $ A $ is symmetric, then all of its eigenvalues are real.  \n",
    "1. If $ A $ is invertible and $ \\lambda_1, \\ldots, \\lambda_n $ are its eigenvalues, then the eigenvalues of $ A^{-1} $ are $ 1/\\lambda_1, \\ldots, 1/\\lambda_n $.  \n",
    "\n",
    "\n",
    "A corollary of the first statement is that a matrix is invertible if and only if all its eigenvalues are nonzero.\n",
    "\n",
    "Using SciPy, we can solve for the eigenvalues and eigenvectors of a matrix as\n",
    "follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be75200b",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "A = ((1, 2),\n",
    "     (2, 1))\n",
    "\n",
    "A = np.array(A)\n",
    "evals, evecs = eig(A)\n",
    "evals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb656a4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "evecs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b857e5f8",
   "metadata": {},
   "source": [
    "Note that the *columns* of `evecs` are the eigenvectors.\n",
    "\n",
    "Since any scalar multiple of an eigenvector is an eigenvector with the same\n",
    "eigenvalue (check it), the eig routine normalizes the length of each eigenvector\n",
    "to one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad15e2e6",
   "metadata": {},
   "source": [
    "### Generalized Eigenvalues\n",
    "\n",
    "It is sometimes useful to consider the *generalized eigenvalue problem*, which, for given\n",
    "matrices $ A $ and $ B $, seeks generalized eigenvalues\n",
    "$ \\lambda $ and eigenvectors $ v $ such that\n",
    "\n",
    "$$\n",
    "A v = \\lambda B v\n",
    "$$\n",
    "\n",
    "This can be solved in SciPy via `scipy.linalg.eig(A, B)`.\n",
    "\n",
    "Of course, if $ B $ is square and invertible, then we can treat the\n",
    "generalized eigenvalue problem as an ordinary eigenvalue problem $ B^{-1}\n",
    "A v = \\lambda v $, but this is not always the case."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ace79f42",
   "metadata": {},
   "source": [
    "## Further Topics\n",
    "\n",
    "We round out our discussion by briefly mentioning several other important\n",
    "topics."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4a6d30b2",
   "metadata": {},
   "source": [
    "### Series Expansions\n",
    "\n",
    "\n",
    "<a id='index-20'></a>\n",
    "Recall the usual summation formula for a geometric progression, which states\n",
    "that if $ |a| < 1 $, then $ \\sum_{k=0}^{\\infty} a^k = (1 - a)^{-1} $.\n",
    "\n",
    "A generalization of this idea exists in the matrix setting.\n",
    "\n",
    "\n",
    "<a id='la-mn'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "729f99ca",
   "metadata": {},
   "source": [
    "#### Matrix Norms\n",
    "\n",
    "\n",
    "<a id='index-21'></a>\n",
    "Let $ A $ be a square matrix, and let\n",
    "\n",
    "$$\n",
    "\\| A \\| := \\max_{\\| x \\| = 1} \\| A x \\|\n",
    "$$\n",
    "\n",
    "The norms on the right-hand side are ordinary vector norms, while the norm on\n",
    "the left-hand side is a *matrix norm* — in this case, the so-called\n",
    "*spectral norm*.\n",
    "\n",
    "For example, for a square matrix $ S $, the condition $ \\| S \\| < 1 $ means that $ S $ is *contractive*, in the sense that it pulls all vectors towards the origin <sup><a href=#cfn id=cfn-link>[2]</a></sup>.\n",
    "\n",
    "\n",
    "<a id='la-neumann'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d1aa063",
   "metadata": {},
   "source": [
    "#### Neumann’s Theorem\n",
    "\n",
    "\n",
    "<a id='index-23'></a>\n",
    "Let $ A $ be a square matrix and let $ A^k := A A^{k-1} $ with $ A^1 := A $.\n",
    "\n",
    "In other words, $ A^k $ is the $ k $-th power of $ A $.\n",
    "\n",
    "Neumann’s theorem states the following: If $ \\| A^k \\| < 1 $ for some\n",
    "$ k \\in \\mathbb{N} $, then $ I - A $ is invertible, and\n",
    "\n",
    "\n",
    "<a id='equation-la-neumann'></a>\n",
    "$$\n",
    "(I - A)^{-1} = \\sum_{k=0}^{\\infty} A^k \\tag{2.4}\n",
    "$$\n",
    "\n",
    "\n",
    "<a id='la-neumann-remarks'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd3bd420",
   "metadata": {},
   "source": [
    "#### Spectral Radius\n",
    "\n",
    "\n",
    "<a id='index-25'></a>\n",
    "A result known as Gelfand’s formula tells us that, for any square matrix $ A $,\n",
    "\n",
    "$$\n",
    "\\rho(A) = \\lim_{k \\to \\infty} \\| A^k \\|^{1/k}\n",
    "$$\n",
    "\n",
    "Here $ \\rho(A) $ is the *spectral radius*, defined as $ \\max_i |\\lambda_i| $, where $ \\{\\lambda_i\\}_i $ is the set of eigenvalues of $ A $.\n",
    "\n",
    "As a consequence of Gelfand’s formula, if all eigenvalues are strictly less than one in modulus,\n",
    "there exists a $ k $ with $ \\| A^k \\| < 1 $.\n",
    "\n",
    "In which case [(2.4)](#equation-la-neumann) is valid."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32a33ff4",
   "metadata": {},
   "source": [
    "### Positive Definite Matrices\n",
    "\n",
    "\n",
    "<a id='index-27'></a>\n",
    "Let $ A $ be a symmetric $ n \\times n $ matrix.\n",
    "\n",
    "We say that $ A $ is\n",
    "\n",
    "1. *positive definite* if $ x' A x > 0 $ for every $ x \\in \\mathbb R ^n \\setminus \\{0\\} $  \n",
    "1. *positive semi-definite* or *nonnegative definite* if $ x' A x \\geq 0 $ for every $ x \\in \\mathbb R ^n $  \n",
    "\n",
    "\n",
    "Analogous definitions exist for negative definite and negative semi-definite matrices.\n",
    "\n",
    "It is notable that if $ A $ is positive definite, then all of its eigenvalues\n",
    "are strictly positive, and hence $ A $ is invertible (with positive\n",
    "definite inverse).\n",
    "\n",
    "\n",
    "<a id='la-mcalc'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4938aee",
   "metadata": {},
   "source": [
    "### Differentiating Linear and Quadratic Forms\n",
    "\n",
    "\n",
    "<a id='index-28'></a>\n",
    "The following formulas are useful in many economic contexts.  Let\n",
    "\n",
    "- $ z, x $ and $ a $ all be $ n \\times 1 $ vectors  \n",
    "- $ A $ be an $ n  \\times n $ matrix  \n",
    "- $ B $ be an $ m \\times n $ matrix and $ y $ be an $ m  \\times 1 $ vector  \n",
    "\n",
    "\n",
    "Then\n",
    "\n",
    "1. $ \\frac{\\partial a' x}{\\partial x} = a $  \n",
    "1. $ \\frac{\\partial A x}{\\partial x} = A' $  \n",
    "1. $ \\frac{\\partial x'A x}{\\partial x} = (A + A') x $  \n",
    "1. $ \\frac{\\partial y'B z}{\\partial y} = B z $  \n",
    "1. $ \\frac{\\partial y'B z}{\\partial B} = y z' $  \n",
    "\n",
    "\n",
    "Exercise 2.1 below asks you to apply these formulas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6977f73d",
   "metadata": {},
   "source": [
    "### Further Reading\n",
    "\n",
    "The documentation of the `scipy.linalg` submodule can be found [here](https://docs.scipy.org/doc/scipy/reference/linalg.html).\n",
    "\n",
    "Chapters 2 and 3 of the [Econometric Theory](https://johnstachurski.net/emet.html) contains\n",
    "a discussion of linear algebra along the same lines as above, with solved exercises.\n",
    "\n",
    "If you don’t mind a slightly abstract approach, a nice intermediate-level text on linear algebra\n",
    "is [[Jänich, 1994](https://python.quantecon.org/zreferences.html#id178)]."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a752442c",
   "metadata": {},
   "source": [
    "## Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "abcbc5fd",
   "metadata": {},
   "source": [
    "## Exercise 2.1\n",
    "\n",
    "Let $ x $ be a given $ n \\times 1 $ vector and consider the problem\n",
    "\n",
    "$$\n",
    "v(x) =  \\max_{y,u} \\left\\{ - y'P y - u' Q u \\right\\}\n",
    "$$\n",
    "\n",
    "subject to the linear constraint\n",
    "\n",
    "$$\n",
    "y = A x + B u\n",
    "$$\n",
    "\n",
    "Here\n",
    "\n",
    "- $ P $ is an $ n \\times n $ matrix and $ Q $ is an $ m \\times m $ matrix  \n",
    "- $ A $ is an $ n \\times n $ matrix and $ B $ is an $ n \\times m $ matrix  \n",
    "- both $ P $ and $ Q $ are symmetric and positive semidefinite  \n",
    "\n",
    "\n",
    "(What must the dimensions of $ y $ and $ u $ be to make this a well-posed problem?)\n",
    "\n",
    "One way to solve the problem is to form the Lagrangian\n",
    "\n",
    "$$\n",
    "\\mathcal L = - y' P y - u' Q u + \\lambda' \\left[A x + B u - y\\right]\n",
    "$$\n",
    "\n",
    "where $ \\lambda $ is an $ n \\times 1 $ vector of Lagrange multipliers.\n",
    "\n",
    "Try applying the formulas given above for differentiating quadratic and linear forms to obtain the first-order conditions for maximizing $ \\mathcal L $ with respect to $ y, u $ and minimizing it with respect to $ \\lambda $.\n",
    "\n",
    "Show that these conditions imply that\n",
    "\n",
    "1. $ \\lambda = - 2 P y $.  \n",
    "1. The optimizing choice of $ u $ satisfies $ u = - (Q + B' P B)^{-1} B' P A x $.  \n",
    "1. The function $ v $ satisfies $ v(x) = - x' \\tilde P x $ where $ \\tilde P = A' P A - A'P B (Q + B'P B)^{-1} B' P A $.  \n",
    "\n",
    "\n",
    "As we will see, in economic contexts Lagrange multipliers often are shadow prices.\n",
    "\n",
    ">**Note**\n",
    ">\n",
    ">If we don’t care about the Lagrange multipliers, we can substitute the constraint into the objective function, and then just maximize $ -(Ax + Bu)'P (Ax + Bu) - u' Q u $ with respect to $ u $.  You can verify that this leads to the same maximizer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f04d0224",
   "metadata": {},
   "source": [
    "## Solution to[ Exercise 2.1](https://python.quantecon.org/#la_ex1)\n",
    "\n",
    "We have an optimization problem:\n",
    "\n",
    "$$\n",
    "v(x) = \\max_{y,u} \\{ -y'Py - u'Qu \\}\n",
    "$$\n",
    "\n",
    "s.t.\n",
    "\n",
    "$$\n",
    "y = Ax + Bu\n",
    "$$\n",
    "\n",
    "with primitives\n",
    "\n",
    "- $ P $ be a symmetric and positive semidefinite $ n \\times n $\n",
    "  matrix  \n",
    "- $ Q $ be a symmetric and positive semidefinite $ m \\times m $\n",
    "  matrix  \n",
    "- $ A $ an $ n \\times n $ matrix  \n",
    "- $ B $ an $ n \\times m $ matrix  \n",
    "\n",
    "\n",
    "The associated Lagrangian is:\n",
    "\n",
    "$$\n",
    "L = -y'Py - u'Qu + \\lambda' \\lbrack Ax + Bu - y \\rbrack\n",
    "$$\n",
    "\n",
    "**Step 1.**\n",
    "\n",
    "Differentiating Lagrangian equation w.r.t y and setting its derivative\n",
    "equal to zero yields\n",
    "\n",
    "$$\n",
    "\\frac{ \\partial L}{\\partial y} = - (P + P') y - \\lambda = - 2 P y - \\lambda = 0 \\:,\n",
    "$$\n",
    "\n",
    "since P is symmetric.\n",
    "\n",
    "Accordingly, the first-order condition for maximizing L w.r.t. y implies\n",
    "\n",
    "$$\n",
    "\\lambda = -2 Py \\:\n",
    "$$\n",
    "\n",
    "**Step 2.**\n",
    "\n",
    "Differentiating Lagrangian equation w.r.t. u and setting its derivative\n",
    "equal to zero yields\n",
    "\n",
    "$$\n",
    "\\frac{ \\partial L}{\\partial u} = - (Q + Q') u - B'\\lambda = - 2Qu + B'\\lambda = 0 \\:\n",
    "$$\n",
    "\n",
    "Substituting $ \\lambda = -2 P y $ gives\n",
    "\n",
    "$$\n",
    "Qu + B'Py = 0 \\:\n",
    "$$\n",
    "\n",
    "Substituting the linear constraint $ y = Ax + Bu $ into above\n",
    "equation gives\n",
    "\n",
    "$$\n",
    "Qu + B'P(Ax + Bu) = 0\n",
    "$$\n",
    "\n",
    "$$\n",
    "(Q + B'PB)u + B'PAx = 0\n",
    "$$\n",
    "\n",
    "which is the first-order condition for maximizing $ L $ w.r.t. $ u $.\n",
    "\n",
    "Thus, the optimal choice of u must satisfy\n",
    "\n",
    "$$\n",
    "u = -(Q + B'PB)^{-1}B'PAx \\:,\n",
    "$$\n",
    "\n",
    "which follows from the definition of the first-order conditions for\n",
    "Lagrangian equation.\n",
    "\n",
    "**Step 3.**\n",
    "\n",
    "Rewriting our problem by substituting the constraint into the objective\n",
    "function, we get\n",
    "\n",
    "$$\n",
    "v(x) = \\max_{u} \\{ -(Ax+ Bu)'P(Ax+Bu) - u'Qu \\} \\:\n",
    "$$\n",
    "\n",
    "Since we know the optimal choice of u satisfies $ u = -(Q +\n",
    "B'PB)^{-1}B'PAx $, then\n",
    "\n",
    "$$\n",
    "v(x) =  -(Ax+ B u)'P(Ax+B u) - u'Q u  \\,\\,\\,\\, with \\,\\,\\,\\, u = -(Q + B'PB)^{-1}B'PAx\n",
    "$$\n",
    "\n",
    "To evaluate the function\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "v(x) &=  -(Ax+ B u)'P(Ax+Bu) - u'Q u \\\\\n",
    "&= -(x'A' + u'B')P(Ax+Bu) - u'Q u \\\\\n",
    "&= - x'A'PAx - u'B'PAx - x'A'PBu - u'B'PBu - u'Qu \\\\\n",
    "&= - x'A'PAx - 2u'B'PAx - u'(Q + B'PB) u\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "For simplicity, denote by $ S := (Q + B'PB)^{-1} B'PA $, then $ u = -Sx $.\n",
    "\n",
    "Regarding the second term $ - 2u'B'PAx $,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "-2u'B'PAx &= -2 x'S'B'PAx  \\\\\n",
    "& = 2 x'A'PB( Q + B'PB)^{-1} B'PAx\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Notice that the term $ (Q + B'PB)^{-1} $ is symmetric as both P and Q\n",
    "are symmetric.\n",
    "\n",
    "Regarding the third term $ - u'(Q + B'PB) u $,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "-u'(Q + B'PB) u &= - x'S' (Q + B'PB)Sx \\\\\n",
    "&= -x'A'PB(Q + B'PB)^{-1}B'PAx\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Hence, the summation of second and third terms is\n",
    "$ x'A'PB(Q + B'PB)^{-1}B'PAx $.\n",
    "\n",
    "This implies that\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    " v(x) &= - x'A'PAx - 2u'B'PAx - u'(Q + B'PB) u\\\\\n",
    " &= - x'A'PAx + x'A'PB(Q + B'PB)^{-1}B'PAx \\\\\n",
    " &= -x'[A'PA - A'PB(Q + B'PB)^{-1}B'PA] x\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Therefore, the solution to the optimization problem\n",
    "$ v(x) = -x' \\tilde{P}x $ follows the above result by denoting\n",
    "$ \\tilde{P} := A'PA - A'PB(Q + B'PB)^{-1}B'PA $\n",
    "\n",
    "<p><a id=fn-mdt href=#fn-mdt-link><strong>[1]</strong></a> Although there is a specialized matrix data type defined in NumPy, it’s more standard to work with ordinary NumPy arrays.\n",
    "See [this discussion](https://python-programming.quantecon.org/numpy.html#matrix-multiplication).\n",
    "\n",
    "<p><a id=cfn href=#cfn-link><strong>[2]</strong></a> Suppose that $ \\|S \\| < 1 $. Take any nonzero vector $ x $, and let $ r := \\|x\\| $. We have $ \\| Sx \\| = r \\| S (x/r) \\| \\leq r \\| S \\| < r = \\| x\\| $. Hence every point is pulled towards the origin."
   ]
  }
 ],
 "metadata": {
  "date": 1736552528.595347,
  "filename": "linear_algebra.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Linear Algebra"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}