psets/pset1sol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "e0edcb1f",
   "metadata": {},
   "source": [
    "# Matrix Calculus (18.S096) Pset 1 Solutions\n",
    "\n",
    "This notebook contains computational solutions for problem set 1 of *18.S096: Matrix Calculus* at MIT in IAP 2024.  See also the solutions PDF for mathematical derivations and analytical results."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ab6299f",
   "metadata": {},
   "source": [
    "## Problem 3\n",
    "\n",
    "You were *not* required to computationally validate your solutions on this problem, but it is always a good idea to do a quick finite-difference check.\n",
    "\n",
    "We pick a random $x \\in \\mathbb{R}^5$ and $\\delta x \\in \\mathbb{R}^5$ on the order of $10^{-8}$ for use in all examples, and/or $A \\in \\mathbb{R}^{5 \\times 5}$ and $\\delta A \\in \\mathbb{R}^{5 \\times 5}$, and a function `relerr` to compute relative errors, which should be $\\ll 1$ compared to finite differences if our answers are correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2bcc633c",
   "metadata": {},
   "outputs": [],
   "source": [
    "using LinearAlgebra # for norm\n",
    "relerr(approx, exact) = norm(approx - exact) / norm(exact)\n",
    "\n",
    "x = randn(5)\n",
    "δx = randn(5)*1e-8\n",
    "A = randn(5,5)\n",
    "δA = randn(5,5) * 1e-8;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f641180",
   "metadata": {},
   "source": [
    "### 3.1\n",
    "\n",
    "$f(x) = \\Vert x \\Vert = \\sqrt{x^T x}$ gives $\\nabla f = x / \\Vert x \\Vert$.  Check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e48718a6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2336454571970445e-9"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_1(x) = norm(x)\n",
    "∇f_1(x) = x / norm(x)\n",
    "relerr(f_1(x + δx) - f_1(x), ∇f_1(x)'δx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad9ae91b",
   "metadata": {},
   "source": [
    "### 3.2\n",
    "\n",
    "$f(x) = \\frac{x^T (A + \\Vert x \\Vert^2 I) x}{x^T x}$ gives\n",
    "$$\n",
    "\\nabla f = \\frac{ \\left(A + A^T + 4 (x^T x) I\\right) x  - 2f(x) x}{x^T x}\n",
    "$$\n",
    "Check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "303f4e10",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.4905132535251805e-8"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_2(x) = (x' * (A + norm(x)^2 * I) * x) / (x'x)\n",
    "∇f_2(x) = ( (A + A' + 4(x'x) * I) * x - 2f_2(x) * x) / (x'x)\n",
    "relerr(f_2(x + δx) - f_2(x), ∇f_2(x)'δx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdffa6ea",
   "metadata": {},
   "source": [
    "### 3.3\n",
    "\n",
    "$f(A) = A^{-2}$ gives $f'(A)[dA] = - A^{-1}\\, dA \\, A^{-2} - A^{-2}\\, dA \\, A^{-1}$.  Check"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2ef40ff0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.2803102089311857e-8"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_3(A) = A^-2\n",
    "f_3′(A,dA) = -A^-1 * dA * A^-2 - A^-2 * dA * A^-1\n",
    "relerr(f_3(A + δA) - f_3(A), f_3′(A, δA))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36c59557",
   "metadata": {},
   "source": [
    "### 3.4\n",
    "\n",
    "$f(A) = (\\mathrm{trace} A)^9$ gives $f'(A)[dA] = 9(\\mathrm{trace} A)^8 * \\mathrm{trace} (dA)$.  Check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1b7ee4f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.217467729191711e-9"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_4(A) = tr(A)^9\n",
    "f_4′(A, dA) = 9tr(A)^8 * tr(dA)\n",
    "relerr(f_4(A + δA) - f_4(A), f_4′(A, δA))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2054b8e2",
   "metadata": {},
   "source": [
    "### 3.5\n",
    "\n",
    "$f_5(x) = A (x .* x)$ gives a Jacobian $f_5 = A .* x^T$.  Check:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "0194dc4f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.557718433805591e-9"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_5(x) = A * (x .* x)\n",
    "f_5′(x) = 2A .* x'\n",
    "relerr(f_5(x + δx) - f_5(x), f_5′(x) * δx)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88da72ce",
   "metadata": {},
   "source": [
    "## Problem 4.1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1969fd0f",
   "metadata": {},
   "source": [
    "Let's check the formula for the derivative of a matrix exponential with a random $3\\times 3$ matrix $A$ and a small perturbation $\\delta A$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "fd78f833",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       "  1.61348   -0.744903  -0.165515\n",
       " -1.36977    2.39593    0.914722\n",
       "  0.298162  -2.59377    0.0465305"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "A = randn(3,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0d7e9cf8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " -1.93183e-8  -1.38694e-8  -6.97525e-9\n",
       " -1.41878e-8   7.05109e-9   5.91334e-9\n",
       "  1.45326e-9  -1.13527e-8  -7.58662e-9"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "δA = randn(3,3) * 1e-8"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87afd4c2",
   "metadata": {},
   "source": [
    "First, we'll compute the finite difference:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "f97c2904",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " -6.10435e-8  -3.3479e-8   -3.41554e-8\n",
       " -4.74015e-8   6.77471e-8   4.9814e-8\n",
       "  5.6351e-8   -7.08204e-8  -5.54767e-8"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "approx = exp(A + δA) - exp(A)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "683cc95c",
   "metadata": {},
   "source": [
    "Now, we'll use the analytical expression via $M = \\begin{pmatrix} A & \\delta A \\\\ & A \\end{pmatrix}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "796ffc14",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Matrix{Float64}:\n",
       " -6.10435e-8  -3.3479e-8   -3.41554e-8\n",
       " -4.74016e-8   6.77471e-8   4.9814e-8\n",
       "  5.6351e-8   -7.08204e-8  -5.54767e-8"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = [ A  δA\n",
    "     0I   A ]\n",
    "exact = exp(M)[1:3, 4:6]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e9e0b1a",
   "metadata": {},
   "source": [
    "By inspection, they match pretty well!  Let's be quantitative using the `relerr` function we defined for problem 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "65cde689",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.244412175629453e-7"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "relerr(approx, exact)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84ed6f80",
   "metadata": {},
   "source": [
    "They match to about 8 digits, which is limited by the accuracy of the finite-difference approximation.  Success!"
   ]
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "kernelspec": {
   "display_name": "Julia 1.10.0",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}