pset 1 solutions
This commit is contained in:
Steven G. Johnson committed
1 parent
8a79abdb34
commit
261c79ed71
2 files changed
+731
-1
No files matched your search
@@ -115,7 +115,7 @@ The *nice* case of diagonalization is when you have **orthonormal eigenvectors**
|
||||
* Uᵣ and Vᵣ as the "right" bases for C(A) and C(Aᵀ). Avₖ=σₖvₖ
|
||||
* SVD and eigenvalues: U and V as eigenvectors of AAᵀ and AᵀA, respectively, and σₖ² as the positive eigenvalues.
|
||||
* Deriving the SVD from eigenvalues: the key step is showing that AAᵀ and AᵀA share the same positive eigenvalues, and λₖ=σₖ², with eigenvectors related by a factor of A. From there you can work backwards to the SVD.
|
||||
* pset 1 solutions: coming soon
|
||||
* [pset 1 solutions](psets/pset1sol.ipynb)
|
||||
* pset 2: coming soon
|
||||
|
||||
**Further reading**: [OCW lecture 6](https://ocw.mit.edu/courses/18-065-matrix-methods-in-data-analysis-signal-processing-and-machine-learning-spring-2018/resources/lecture-6-singular-value-decomposition-svd/) and textbook section I.8. The [Wikipedia SVD article](https://en.wikipedia.org/wiki/Singular_value_decomposition).
|
||||
@@ -0,0 +1,730 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "08cd437b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# 18.065 Pset 1 Solutions"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0d53b289",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 1 (4+4+4+4+4 points)\n",
|
||||
"\n",
|
||||
"Recall from class that multiplying an $m \\times p$ by a $p \\times n$ matrix costs $mnp$ scalar multiplications (and a similar number of additions) by the standard (practical) algorithms.\n",
|
||||
"\n",
|
||||
"Matrix multiplication is **not commutative** ($AB \\ne BA$ in general), but it **is associative**: $(AB)C=A(BC)$. It turns out that where you put the parentheses (i.e. in *what order* you do the multiplications) can make a *huge* difference in computational cost.\n",
|
||||
"\n",
|
||||
"**(a)** If $x \\in \\mathbb{R}^n$ and $A,B$ are $n \\times n$ matrices, compare the scalar multiplication counts of $(AB)x$ vs. $A(Bx)$, i.e. if we do the multiplications in the order indicated by the parentheses.\n",
|
||||
"\n",
|
||||
"**(b)** If $x, b\\in \\mathbb{R}^n$, **how many scalar multiplications** does the computation $$p = (I - (xx^T)/(x^T x)) b$$ take if we *do it in the order indicated by the parentheses*? (Note that dividing by a scalar $\\alpha$ is equivalent to multiplying by $\\alpha^{-1}$ at the negligible cost of one scalar division.)\n",
|
||||
"\n",
|
||||
"**(c)** Explain how to compute the *same* $p$ as in part (b) using as *few multiplications as possible*. Outline the sequence of computational steps, and give the count of multiplications.\n",
|
||||
"\n",
|
||||
"**(d)** $p^T x = $ what?\n",
|
||||
"\n",
|
||||
"**(e)** Implement your algorithm from (c) in Julia, filling in the code below, and time it for $n=1000$ using the `@btime` macro from the [BenchmarkTools package](https://github.com/JuliaCI/BenchmarkTools.jl), along with the algorithm from part (b), following the outline below. How does the ratio of the two times compare to your ratio of multiplication counts?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d5d43b59",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution:\n",
|
||||
"\n",
|
||||
"**(a)** $AB$ requires $n^3$ scalar multiplications ($n^2$ row×col dot products), and $Bx$ requires $n^2$ multiplications ($n$ dot products). So, $(AB)x$ requires $\\boxed{n^3 + n^2}$ multiplications, while $A(Bx)$ requires $\\boxed{2n^2}$ multiplications (a *lot* fewer for large $n$!).\n",
|
||||
"\n",
|
||||
"**(b)** $$p = (I - (xx^T)/(x^T x)) b$$, computed in this order, requires (1) $n$ multiplies for $x^T x$, (2) $n^2$ multiplies for $xx^T$, (3) $n^2$ multiplies to scale $xx^T$ by $(x^T x)^{-1}$ (+ 1 division), (4) no multiplies to subtract from $I$, and (5) $n^2$ multiplies to multiply by $b$. Total: $\\boxed{3n^2 + n}$ scalar multiplications.\n",
|
||||
"\n",
|
||||
"**(c)** If we re-arrange it as $$p = b - x((x^T b)/(x^T x))$$, then it requires only $\\boxed{3n+1}$ multiplications: $2n$ for the two dot products, 1 for the scalar division (counting it as a scalar multiplication by the inverse as instructed), and $n$ to scale $x$.\n",
|
||||
"\n",
|
||||
"**(d)** By remembering your linear algebra, you might recognize $P = (xx^T)/(x^T x)$ as **orthogonal projection** onto the direction of $x$, and hence $I-P$ is projection **perpendicular** to $x$. So, we must have $\\boxed{p^T x = 0}$. But we can also compute this explicitly:\n",
|
||||
"$$\n",
|
||||
"p^T x = x^T p = x^T \\left[ b - x((x^T b)/(x^T x)) \\right] = \\\\\n",
|
||||
"x^T b - (x^T x)(x^T b)/(x^T x) = x^T b - x^T b = 0 \\, .\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"**(e)** See the following code. On my laptop, my part (c) solution is about **3400× faster** than part (b), which is actually bigger than the ratio of $(3n^2 + n) / (3n+1) = n = 1000$ you might predict just by counting multiplications, but it is the right order of magnitude.\n",
|
||||
"\n",
|
||||
"Of course, computers are *much* more complicated than simple arithmetic counts than capture, and you will almost never be able to predict performance ratios simply by arithmetic counts. To get a sense of this complexity, it's interesting to plot $\\frac{1}{n}\\text{(time for b)}/\\text{(time for c)}$ versus $n$. Arithmetic counts would predict this ratio to be a constant, but in fact it is much more complicated. (You were *not* required to do this — I show it below only for interest.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "b149871c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"Hooray, part (c) and part (b) agree!\n",
|
||||
"\n",
|
||||
"part (b): \n",
|
||||
" 4.143 ms (7 allocations: 22.90 MiB)\n",
|
||||
"\n",
|
||||
"part (c): \n",
|
||||
" 1.279 μs (2 allocations: 15.88 KiB)\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using LinearAlgebra, BenchmarkTools\n",
|
||||
"\n",
|
||||
"# algorithm from part (b)\n",
|
||||
"function part_b(x, b)\n",
|
||||
" return (I - (x*x')*(x'*x)^-1) * b\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"# algorithm from part (c)\n",
|
||||
"function part_c(x, b)\n",
|
||||
" # CHANGE THIS:\n",
|
||||
" return b - x * ((x'*b)/(x'*x))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"# test and benchmark on random vectors:\n",
|
||||
"n = 1000\n",
|
||||
"x, b = rand(n), rand(n)\n",
|
||||
"\n",
|
||||
"# test it first — should give same answer up to roundoff error\n",
|
||||
"if part_c(x, b) ≈ part_b(x, b)\n",
|
||||
" println(\"Hooray, part (c) and part (b) agree!\")\n",
|
||||
"else\n",
|
||||
" error(\"You made a mistake: part (c) and part (b) do not agree!\")\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"# benchmark it:\n",
|
||||
"\n",
|
||||
"println(\"\\npart (b): \")\n",
|
||||
"@btime part_b($x, $b);\n",
|
||||
"\n",
|
||||
"println(\"\\npart (c): \")\n",
|
||||
"@btime part_c($x, $b);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "27552024",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"3450.351837372948"
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"4.413e-3 / 1.279e-6 # ratio of times"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"id": "bf2cee16",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"n = 10\n",
|
||||
"n = 14\n",
|
||||
"n = 21\n",
|
||||
"n = 30\n",
|
||||
"n = 43\n",
|
||||
"n = 62\n",
|
||||
"n = 89\n",
|
||||
"n = 127\n",
|
||||
"n = 183\n",
|
||||
"n = 264\n",
|
||||
"n = 379\n",
|
||||
"n = 546\n",
|
||||
"n = 785\n",
|
||||
"n = 1129\n",
|
||||
"n = 1624\n",
|
||||
"n = 2336\n",
|
||||
"n = 3360\n",
|
||||
"n = 4833\n",
|
||||
"n = 6952\n",
|
||||
"n = 10000\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"20-element Vector{Float64}:\n",
|
||||
" 0.4198663315479747\n",
|
||||
" 0.4024561123451314\n",
|
||||
" 0.5016266576985883\n",
|
||||
" 0.5217398342398343\n",
|
||||
" 0.5499016471424656\n",
|
||||
" 0.5297567217200534\n",
|
||||
" 0.7579014598586505\n",
|
||||
" 1.0874537900252532\n",
|
||||
" 1.0783124610718169\n",
|
||||
" 0.7767890148134123\n",
|
||||
" 0.9204295359856696\n",
|
||||
" 1.2168404617858037\n",
|
||||
" 2.2617932818717397\n",
|
||||
" 3.025106985400225\n",
|
||||
" 3.060177735357155\n",
|
||||
" 3.1275025975637285\n",
|
||||
" 3.0167103840430975\n",
|
||||
" 7.520125153516695\n",
|
||||
" 8.25342810247817\n",
|
||||
" 8.673908026056296"
|
||||
]
|
||||
},
|
||||
"execution_count": 7,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# for interest, let's plot (time for b) / (n*(time for c)) versus n,\n",
|
||||
"# and see how different it is from a constant:\n",
|
||||
"\n",
|
||||
"ns = round.(Int, 10 .^ range(1, 4, length=20)) # log-spaced points\n",
|
||||
"\n",
|
||||
"t_bc = map(ns) do n\n",
|
||||
" @show n\n",
|
||||
" x, b = rand(n), rand(n)\n",
|
||||
" @belapsed(part_b($x, $b)), @belapsed(part_c($x, $b))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"t_ratio = first.(t_bc) ./ (ns .* last.(t_bc))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 10,
|
||||
"id": "d6fa1641",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAHJCAYAAACMppPqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABbO0lEQVR4nO3deVxUVf8H8M+wLyKCgqyCqYiouJK5oJL7gibaY2nmUplmubboY2VqT2a2uISV9ZTZo6Ypmqllbii55JLSopYLKiqmogGKIgzn98f5DTAMywzMzJ3l83695sVw58693xkuzJdzvucclRBCgIiIiMgGOSgdABEREZGpMNEhIiIim8VEh4iIiGwWEx0iIiKyWUx0iIiIyGYx0SEiIiKbxUSHiIiIbBYTHSIiIrJZTHSIiIjIZjHRMbGcnBy8/PLL6NmzJ/z8/KBSqfDGG2+Uu39+fj4iIyPx9ttvG3yuv/76Cy4uLvjll1+qEXGx7777DvHx8ahbty5cXFzg6+uLbt26YeXKlcjPzzfKOZS2devWCn8eROZW+m9EcnIyVCoVkpOTTXK+K1eu4I033sDx48eNfuzly5dDpVLh/PnzFe5X0e+hSqXC888/b/TYbMWqVauwcOFCne3nz5+HSqXCu+++W+kx3njjDahUqiqdX/MzPnLkSJWebw5MdEwsMzMTy5YtQ15eHh555JFK91+6dClu3bqFF154weBzRUREYPjw4ZgyZUoVIi0mhMDo0aMxYMAAFBYW4v3338eOHTvw5ZdfokWLFnjuueewdOnSap3DUmzduhWzZ89WOgyicrVu3RoHDhxA69atTXL8K1euYPbs2SZJdPTF38OqKy/RMcTTTz+NAwcOGCcgC+SkdAC2LiwsDLdu3YJKpcKNGzfw2WeflbtvQUEBFixYgDFjxsDT07NK53v++efRtm1b7N+/Hx06dKjSMRYsWIDly5dj9uzZeP3117Uei4+Px8svv4wzZ85U6dj2Ijc3Fx4eHnZ3bmPIz8+HSqWCkxP/PAFAzZo18dBDDykdBtmwkJAQhISEKB2GybBFx8RUKpXeTYKbNm3C5cuXMWLECJ3HTp8+jWHDhsHf3x+urq5o0qQJEhMTdfZr06YNmjRpgo8//rhK8ebn52P+/PmIjIzEa6+9VuY+AQEB6NSpU9H3N2/exHPPPYfg4GC4uLjggQcewMyZM5GXl6f1PE0T9BdffIHGjRvD3d0dbdu2xcGDByGEwIIFC1C/fn3UqFEDDz/8sE4y1bVrVzRr1gwpKSl46KGH4O7ujuDgYLz22mtQq9VF+5XX1K9pyl2+fDkAYNSoUUXvoebnVLKZXQiBpUuXomXLlnB3d4ePjw+GDBmCc+fOlRnX3r170aFDB3h4eGDMmDEAgF27dqFr166oXbs23N3dUa9ePQwePBi5ubkV/hzWrFmDnj17IjAwEO7u7mjSpAmmT5+OO3fuaO03atQo1KhRA7/99ht69uwJLy8vdOvWDQBw//59vPnmm4iMjISrqyv8/PwwevRoXL9+vcJzlzzuqVOn0KtXL3h6eiIwMLCoS/XgwYPo1KkTPD09ERERgS+//FLnGL///jsGDhwIHx8fuLm5oWXLljr7aX5WX331FaZNm4bg4GC4uroW/ex37NiBbt26oWbNmvDw8EDHjh2xc+fOSuMHgH/++QfTpk3DAw88AFdXV/j7+6Nv3744depU0T6Wdu2Wpbzr+eeff0Z8fDxq164NNzc3NGjQAJMnTy56/MyZMxg9ejQaNWoEDw8PBAcHIz4+Hr/99pvWsWNiYgAAo0ePLvodKNmNdOTIEQwYMAC+vr5wc3NDq1atsHbtWp04Dx48iI4dO8LNzQ1BQUGYMWOGXl3clf0eanz11Vdo0qQJPDw80KJFC2zevFnrcX1eb8n3c/Xq1Zg5cyaCgoJQs2ZNdO/eHX/++WeFsaakpBQ9t7QVK1ZApVLh8OHDRds2bdqE9u3bw8PDA15eXujRo4dOy8moUaMQHh6uczx9upO6du2KLVu24MKFC1rvXWnvv/9+0fXZvn17HDx4sNJzhYeHo3///vjhhx/QunVruLu7IzIyEp9//nmFMQFARkYG2rRpg0aNGuH06dOV7m9ygszm+vXrAoCYNWtWmY+PGTNG+Pv762z/448/hLe3t2jevLlYsWKF+PHHH8W0adOEg4ODeOONN3T2Hz9+vKhTp44oLCws2rZ79+4Kz62xf/9+AUC88sorer2mu3fviujoaOHp6Sneffdd8eOPP4rXXntNODk5ib59+2rtC0CEhYWJDh06iKSkJLFhwwYREREhfH19xZQpU8TAgQPF5s2bxcqVK0XdunVFdHS01mvo0qWLqF27tggKChKLFy8W27ZtExMnThQAxIQJE3Re6+7du7XOn5aWJgCIL774QgghxJkzZ8SQIUMEAHHgwIGi271794QQQjzzzDPC2dlZTJs2Tfzwww9i1apVIjIyUtStW1dcvXpVKy5fX18RGhoqlixZInbv3i327Nkj0tLShJubm+jRo4fYuHGjSE5OFitXrhQjRowQt27dqvB9nTt3rvjggw/Eli1bRHJysvj4449F/fr1RVxcnNZ+I0eOFM7OziI8PFzMmzdP7Ny5U2zbtk2o1WrRu3dv4enpKWbPni22b98uPvvsMxEcHCyioqJEbm5uhecfOXKkcHFxEU2aNBGLFi0S27dvF6NHjxYAxIwZM0RERIT473//K7Zt2yb69+8vAIgjR44UPf/UqVPCy8tLNGjQQKxYsUJs2bJFPP744wKAmD9/vs7PKjg4WAwZMkRs2rRJbN68WWRmZoqvvvpKqFQq8cgjj4ikpCTx3Xffif79+wtHR0exY8eOCuPPzs4WTZs2FZ6enmLOnDli27ZtYv369WLSpEli165dQgjLvHY15yr5e1rW9fzDDz8IZ2dnER0dLZYvXy527dolPv/8c/HYY48V7bNnzx4xbdo0sW7dOrFnzx6xYcMG8cgjjwh3d3dx6tQpIYQQWVlZ4osvvhAAxKuvvlr0O5Ceni6EEGLXrl3CxcVFxMbGijVr1ogffvhBjBo1Suv3SAj5N8rDw0NERUWJ1atXi2+//Vb06tVL1KtXTwAQaWlp5f6sKvs9BCDCw8PFgw8+KNauXSu2bt0qunbtKpycnMTZs2cNer0l38/w8HAxfPhwsWXLFrF69WpRr1490ahRI1FQUFBurEII0apVK9GxY0ed7TExMSImJqbo+5UrVwoAomfPnmLjxo1izZo1ok2bNsLFxUWkpKQU7Tdy5EgRFhamc7xZs2aJyj6i//jjD9GxY0cREBCg9d4JUfz3Ljw8XPTu3Vts3LhRbNy4UTRv3lz4+PiIf/75p8JzhYWFiZCQEBEVFSVWrFghtm3bJh599FEBQOzZs6doP831c/jwYSGEEL/99psIDQ0V7du3F9evX68wfnNhomNGlSU6TZo0Eb1799bZ3qtXLxESEiKysrK0tj///PPCzc1N3Lx5U2v7p59+KgCIkydPFm1LTk4Wjo6OYvbs2RXG+PXXXwsA4uOPP9brNX388ccCgFi7dq3W9vnz5wsA4scffyzaBkAEBASI27dvF23buHGjACBatmyp9cGwcOFCAUD8+uuvRdu6dOkiAIhvv/1W61zPPPOMcHBwEBcuXBBC6J/oCCHEhAkTyvxjcuDAAQFAvPfee1rb09PThbu7u3j55Zd14tq5c6fWvuvWrRMAxPHjx3WOb4jCwkKRn58v9uzZIwCI1NTUosdGjhwpAIjPP/9c6zmrV68WAMT69eu1th8+fFgAEEuXLq3wnJrjlnx+fn6+8PPzEwDEL7/8UrQ9MzNTODo6iqlTpxZte+yxx4Srq6u4ePGi1nH79OkjPDw8iv7Ian5WnTt31trvzp07wtfXV8THx2ttV6vVokWLFuLBBx+sMP45c+YIAGL79u3l7mOJ167mXJUlOg0aNBANGjQQd+/erfB9KKmgoEDcv39fNGrUSEyZMqVou+aaKPl7oREZGSlatWol8vPztbb3799fBAYGCrVaLYQQYujQocLd3V3rH4CCggIRGRlZaaIjRPm/h0LI96Nu3boiOzu7aNvVq1eFg4ODmDdvnsGvV/N+lk5m165dW5RsVUTzwX7s2LGibYcOHRIAxJdffimEkNdpUFCQaN68edF7JIQQOTk5wt/fX3To0KFoW3USHSGE6NevX5nP1/y9a968uVbypol19erVFZ4rLCxMuLm5aV2bd+/eFb6+vuLZZ5/VeT8OHz4stm/fLmrWrCmGDBli0LVpauy6siBXrlyBv7+/1rZ79+5h586dGDRoEDw8PFBQUFB069u3L+7du6fTDKk5xuXLl4u2denSBQUFBTo1N9W1a9cueHp6YsiQIVrbR40aBQA63QxxcXFa9UdNmjQBAPTp00er6VSz/cKFC1rP9/LywoABA7S2DRs2DIWFhdi7d2/1XkwJmzdvhkqlwhNPPKH1ngcEBKBFixY63Qg+Pj54+OGHtba1bNkSLi4uGDt2LL788kudLq+KnDt3DsOGDUNAQAAcHR3h7OyMLl26AABOnjyps//gwYN14q9Vqxbi4+O14m/ZsiUCAgL0GsGjUqnQt2/fou+dnJzQsGFDBAYGolWrVkXbfX194e/vr/Wz2rVrF7p164bQ0FCtY44aNQq5ubk6zfel49+/fz9u3ryJkSNHasVfWFiI3r174/DhwzrdeCV9//33iIiIQPfu3cvdx1qv3b/++gtnz57FU089BTc3t3L3KygowFtvvYWoqCi4uLjAyckJLi4uOH36dJnXUGlnzpzBqVOnMHz48KLjlfzbk5GRUdTVs3v3bnTr1g1169Yter6joyOGDh2q9+uqSFxcHLy8vIq+r1u3rs41Z+jrLf2ziI6OBqD7cyvt8ccfh7+/v1bpwJIlS+Dn51f0ev/8809cuXIFI0aMgIND8cdsjRo1MHjwYBw8eLDS7mtj6devHxwdHYuLine truncated
|
||||
"text/plain": [
|
||||
"Figure(PyObject <Figure size 640x480 with 1 Axes>)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"PyObject Text(24.000000000000007, 0.5, '(time for b) / n×(time for c)')"
|
||||
]
|
||||
},
|
||||
"execution_count": 10,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using PyPlot\n",
|
||||
"loglog(ns, t_ratio, \"ro-\")\n",
|
||||
"title(\"1(e): Computers are more complicated than you think\")\n",
|
||||
"xlabel(L\"n\")\n",
|
||||
"ylabel(\"(time for b) / n×(time for c)\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "ef31bc89",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The ratio (the speedup from part c) actually **grows with n**, mainly because part (b) also uses more memory ($\\sim n^2$ rather than $\\sim n$) and [using more memory is slower](https://en.wikipedia.org/wiki/Memory_hierarchy).\n",
|
||||
"\n",
|
||||
"Interestingly, part (b) is actually *faster* than (c) for very small matrices, possibly because matrix–vector products are so heavily optimized in [linear-algebra software](https://www.openblas.net/). We could get a bit faster in (c) by optimizing `b - x * ((x'*b)/(x'*x))` in order to [fuse the loops](https://julialang.org/blog/2017/01/moredots/), or maybe by calling an optimized [BLAS](https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms) operation for [axpy](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.axpy!)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "bfa41bf1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 2 (8+4 points)\n",
|
||||
"\n",
|
||||
"**(a)** Describe the four fundamental subspaces of the **rank-1 matrix** $A = uv^T$ where $u \\in \\mathbb{R}^m$ and $n \\in \\mathbb{R}^n$.\n",
|
||||
"\n",
|
||||
"**(b)** For any column vectors $u,v \\in \\mathbb{R}^3$, the matrix $uv^T$ is rank 1, except when \\_\\_\\_\\_\\_\\_\\_\\_, in which case $uv^T$ has rank \\_\\_\\_\\_."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "54223d31",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution\n",
|
||||
"\n",
|
||||
"**(a)** Since the matrix is given to be \"rank-1\", we can assume $u, v \\ne 0$ (see part b). In that case, the four subspaces are:\n",
|
||||
"\n",
|
||||
"1. Column space $C(uv^T) \\subseteq \\mathbb{R}^m$ is 1-dimensional and is **spanned by u**.\n",
|
||||
"2. Row space $C(vu^T) \\subseteq \\mathbb{R}^n$ is 1-dimensional and is **spanned by v**.\n",
|
||||
"3. Nullspace $N(uv^T)$ is $n-1$ dimensional and is **every vector orthogonal to v**.\n",
|
||||
"3. Left nullspace $N(vu^T)$ is $m-1$ dimensional and is **every vector orthogonal to u**.\n",
|
||||
"\n",
|
||||
"**(b)** $uv^T$ is always rank 1 exept when **u and/or v is 0**, in which case it has **rank 0**."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "82b2c169",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 3 (5+4+4+4 points)\n",
|
||||
"\n",
|
||||
"**(a)** Pick the choices that makes this statement correct for arbitrary matrices A and B: $C(AB)$ (*contains / is contained in*) the column space of (*A / B*). Briefly justify your answer.\n",
|
||||
"\n",
|
||||
"**(b)** Suppose that $A$ is a $1000\\times 1000$ matrix of rank $< 10$. Suppose we multiply it by 10 random vectors $x_1, x_2, \\ldots, x_{10}$, e.g. generated by `randn(1000)`. How could we use the results to get a $10\\times 10$ matrix $C$ whose rank (almost certainly) matches $A$'s?\n",
|
||||
"\n",
|
||||
"**(c)** Suppose we instead make $1000\\times 10$ matrix $X$ whose columns are $x_1, x_2, \\ldots, x_{10}$. Give a formula for the *same* matrix $C$ in terms of matrix products involving $A$ and $X$.\n",
|
||||
"\n",
|
||||
"**(d)** Fill in the code for $C$ below, and compare the biggest 10 singular values of $A$ (chosen to be rank ≈ 4 in this case) to the corresponding 10 singular values of $C$. Does it match what you expect?"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4dc3782a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution\n",
|
||||
"\n",
|
||||
"**(a)** C(AB) **is contained in** the column space of **A**. The reason is that any $y \\in C(AB)$, we have $y = ABx = A(Bx)$, so $y \\in C(A)$. (The converse is not necessarily true, e.g. if $B = 0$.)\n",
|
||||
"\n",
|
||||
"**(b and c)** I will answer these together. The only reason for asking (c) separately was to make sure you know how to express your answer in matrix notation.\n",
|
||||
"\n",
|
||||
"First, realize that $AX$ *almost certainly* (probability 1) has the **same rank as A**, i.e. we have a $1000 \\times 10$ matrix with the same rank as A.\n",
|
||||
"* Proof: Suppose $C(AX)$ were lower-dimensional than $C(A)$, in which case consider the orthogonal complement of $C(AX)$ in $C(A)$: there must be some direction $y\\ne 0$ in $C(A) = N(A^T)^\\perp$ perpendicular to everything in $C(AX)$, i.e. $y^TAX = 0$. But that would require every one of our \"random\" vectors $x_k$ to be orthogonal to $A^T y$ (≠ 0 since $y \\notin N(A^T)$), which happens with probability 0.\n",
|
||||
"\n",
|
||||
"Now, how can we reduce this to a $10\\times 10$ matrix? If we are to use only elementary matrix operations (no fair using something like an SVD that you don't personally know how to compute), this suggests multipling $AX$ on the left by some $10 \\times 1000$ matrix. But which one? There are several valid choices:\n",
|
||||
"\n",
|
||||
"1. You could use $X^T AX$. (This could only reduce the rank if $X$ was $\\perp$ to some portion of $C(A)$, which happens with probability 0.)\n",
|
||||
"2. You could use $Y^T AX$, where $Y$ is another $1000 \\times 10$ random matrix like `randn(1000,10)`. (Again almost surely will have the same rank.)\n",
|
||||
"3. You could use $(AX)^T (AX) = X^T A^T A X$. In 18.06 and similar intro linear-algebra courses we prove that $B^T B$ has the same rank as $B$ for any $B$. (Proof by showing that $N(B) = N(B^T B)$: $Bx = 0 \\implies B^TBx = 0$, and $B^T Bx = 0 \\implies x^T B^T B x = 0 = \\Vert Bx \\Vert^2 \\implies Bx = 0$.)\n",
|
||||
"\n",
|
||||
"I think these are the easiest possibilities.\n",
|
||||
"\n",
|
||||
"**(d)** Let's try all three possibilities from above in Julia:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 11,
|
||||
"id": "e7b32901",
|
||||
"metadata": {
|
||||
"scrolled": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"(svdvals(A))[1:10] = [1069.3195877770895, 1032.123298970407, 998.5420983487572, 963.6168764929388, 1.197027719541634e-12, 9.397013502391205e-13, 8.816489090342173e-13, 8.368786422681453e-13, 7.163119064816855e-13, 6.558305069443912e-13]\n",
|
||||
"svdvals(X' * A * X) = [17355.334466001634, 10725.604519449384, 8360.035414898422, 2369.8158461371777, 7.68456896867752e-12, 5.9484402119558325e-12, 2.9567758084067988e-12, 2.0377094428213978e-12, 1.0361892803083682e-12, 1.7189758977956257e-13]\n",
|
||||
"svdvals(Y' * A * X) = [11903.189430493283, 11182.390116386241, 6413.062698992587, 3644.0437181289303, 4.420883286604323e-12, 3.670046234886016e-12, 3.053253554336238e-12, 2.327401156370041e-12, 2.2006586809078315e-12, 1.4777079446259352e-12]\n",
|
||||
"svdvals((A * X)' * A * X) = [1.8887124303351242e7, 1.2865209490083521e7, 7.204356096870981e6, 3.603491360462107e6, 4.0543542578357826e-9, 3.7280041982774055e-9, 2.130547526629226e-9, 1.2473174612561205e-9, 1.113501138973212e-9, 8.507451167986952e-12]\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"10-element Vector{Float64}:\n",
|
||||
" 1.8887124303351242e7\n",
|
||||
" 1.2865209490083521e7\n",
|
||||
" 7.204356096870981e6\n",
|
||||
" 3.603491360462107e6\n",
|
||||
" 4.0543542578357826e-9\n",
|
||||
" 3.7280041982774055e-9\n",
|
||||
" 2.130547526629226e-9\n",
|
||||
" 1.2473174612561205e-9\n",
|
||||
" 1.113501138973212e-9\n",
|
||||
" 8.507451167986952e-12"
|
||||
]
|
||||
},
|
||||
"execution_count": 11,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using LinearAlgebra\n",
|
||||
"\n",
|
||||
"# random 1000x1000 matrix of rank 4\n",
|
||||
"A = randn(1000, 4) * randn(4, 1000)\n",
|
||||
"@show svdvals(A)[1:10]\n",
|
||||
"\n",
|
||||
"# possibility 1: XᵀAX\n",
|
||||
"X = randn(1000, 10)\n",
|
||||
"@show svdvals(X'*A*X)\n",
|
||||
"\n",
|
||||
"# possibility 2: YᵀAX\n",
|
||||
"Y = randn(1000, 10)\n",
|
||||
"@show svdvals(Y'*A*X)\n",
|
||||
"\n",
|
||||
"# possibility 3: (AX)ᵀAX\n",
|
||||
"@show svdvals((A*X)'*A*X)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "626f50e3",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"A has rank 4, so it should only have 4 nonzero singular values in theory. In fact, *all* of its singular values are nonzero, but **only 4** of them (≈ 1069, 1032, 998, 963) are **non-negligible** — the remaining singular values are on the order of [machine precision 10⁻¹⁶](https://en.wikipedia.org/wiki/Machine_epsilon) times $\\sigma_1 \\approx 1000$, i.e. they are $\\sim 10^{-12}$. These are due to [roundoff errors](https://en.wikipedia.org/wiki/Round-off_error) in computer arithmetic.\n",
|
||||
"\n",
|
||||
"For all 3 of our $10\\times 10$ matrices, we indeed see the same thing: **4 non-negligible** singular values, and the rest are roundoff errors. Indeed, 4 singular values are much *larger* than those of $A$, which is a consequence of our lack of normalization of $X$ and $Y$. The biggest are the singular values of $(AX)^T AX = X^T A^T A X$, which shouldn't be surprising — this *squares* the singular values of $AX$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "ee67e16b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 4 (4+5+5 points)\n",
|
||||
"\n",
|
||||
"The famous Hadamard matrices are filled with ±1 and have orthogonal columns (orthonormal if we divide $H_n$ by $1/\\sqrt{n}$). The first few are:\n",
|
||||
"\\begin{align}\n",
|
||||
" H_1 &= \\begin{pmatrix} 1 \\end{pmatrix}, \\\\\n",
|
||||
" H_2 &= \\begin{pmatrix}\n",
|
||||
" 1 & 1 \\\\\n",
|
||||
" 1 & -1\n",
|
||||
" \\end{pmatrix}, \\\\\n",
|
||||
" H_4 &= \\begin{pmatrix} H_2 & H_2 \\\\ H_2 & -H_2 \\end{pmatrix} = \\begin{pmatrix}\n",
|
||||
" 1 & 1 & 1 & 1\\\\\n",
|
||||
" 1 & -1 & 1 & -1\\\\\n",
|
||||
" 1 & 1 & -1 & -1\\\\\n",
|
||||
" 1 & -1 & -1 & 1\n",
|
||||
" \\end{pmatrix} \\, .\n",
|
||||
"\\end{align}\n",
|
||||
"Notice that (for power-of-2 sizes), they are built up \"recursively\" out of smaller Hadamard matrices. Multiplying a vector by a Hadamard matrix requires no multiplications at all, only additions/subtractions.\n",
|
||||
"\n",
|
||||
"**(a)** If you multiply $H_4 x$ for some $x \\in \\mathbb{R}^4$ by the normal \"rows-times-columns\" method (*without* exploiting any special patterns), exactly how many scalar additions/subtractions are required?\n",
|
||||
"\n",
|
||||
"**(b)** Let's break $x$ into two blocks: $x = \\begin{pmatrix} x_1 \\\\ x_2 \\end{pmatrix}$ for $x_1, x_2 \\in \\mathbb{R}^2$. Write out $H_4 x$ in terms of a sequence of $2\\times 2$ block multiplications with $\\pm H_2$. You'll notice that some of these $2\\times 2$ multiplications are repeated. If we re-use these repeated multiplications rather than doing them twice, we can save a bunch of arithmetic — what is the new count of scalar additions/subtractions if you do this?\n",
|
||||
"\n",
|
||||
"**(c)** Similarly, the $8\\times 8$ Hadamard matrix $H_8 = \\begin{pmatrix} H_4 & H_4 \\\\ H_4 & -H_4 \\end{pmatrix}$ is made out of $H_4$ matrices. To multiply it by a vector $y \\in \\mathbb{R}^8$, the naive rows-times columns method would require \\_\\_\\_\\_ scalar additions/subtractions, whereas if you broke them up first into blocks of 4, used your solution from (b), and then re-used any repeated $H_4$ products, it would only require \\_\\_\\_\\_ scalar additions/subtractions."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "e419713b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solutions:\n",
|
||||
"\n",
|
||||
"**(a)** $H_4 x$ requires 4 dot products (with the rows of $H_4$), each of which requires 3 scalar additions or subtractions, so there are **12 additions/subtractions** overall.\n",
|
||||
"\n",
|
||||
"**(b)** We have:\n",
|
||||
"$$\n",
|
||||
"H_4 x = \\begin{pmatrix} H_2 x_1 + H_2 x_2 \\\\ H_2 x_1 - H_2 x_2 \\end{pmatrix}\n",
|
||||
"$$\n",
|
||||
"so we see that we need to compute $H_2 x_1$ and $H_2 x_2$ ($=2 \\times 2 \\times 1 = 4$ additions/subtractions) and then add/subtract the resulting 2-component vectors ($= 2\\times 2 = 4$ additions/subtractions), for **8 additions/subtractions overall**.\n",
|
||||
"\n",
|
||||
"**(c)** In the same way, a naive $H_8 y$ would require $8 \\times 7 = \\boxed{56}$ additions/subtractions, but breaking it into blocks would require 2 multiplications by $H_4$ (= 8 adds/subs each) and 2 additions/subtractions of 4-component vectors (= 8 add/subs), for $\\boxed{8 \\times 3 = 24}$ additions/subtractions.\n",
|
||||
"\n",
|
||||
"If you continue this process recursively, you obtain the [fast Walsh–Hadamard transform](https://en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform), which reduces the cost to multiply by $H_n$ from $\\Theta(n^2)$ to $\\Theta(n \\log n)$ for $n = 2^k$!\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "c2a2b8eb",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 5 (5+5 points)\n",
|
||||
"\n",
|
||||
"The famous \"discrete Fourier transform\" matrix $F$ has columns that are actually eigenvectors of the (unitary) permutation matrix:\n",
|
||||
"$$\n",
|
||||
"P = \\begin{pmatrix} & 1 & & \\\\ & &1 & \\\\ & & &1 \\\\ 1 & & & \\end{pmatrix}\n",
|
||||
"$$\n",
|
||||
"for the $4\\times 4$ case, and similarly for larger matrices.\n",
|
||||
"\n",
|
||||
"**(a)** One way of saying way Fourier transforms are practically important is that they *diagonalize* (are eigenvectors of) matrices that *commute* with $P$. If $A$ is a $4\\times 4$ matrix whose first row is $(a\\, b\\, c\\, d)$\n",
|
||||
"$$\n",
|
||||
"A = \\begin{pmatrix} a & b & c & d \\\\ ? & ? &? & ? \\\\ ? &? & ? &? \\\\ ? & ? & ? &? \\end{pmatrix}\n",
|
||||
"$$\n",
|
||||
"that commutes with $P$ (i.e. $AP=PA$), what must be true of the other (\"?\") entries of $A$?\n",
|
||||
"\n",
|
||||
"**(b)** Fill in the matrix `A` in Julia below and fill in and run the code to check that it commutes with $P$ and is diagonalized by $F$:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "14b558f8",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Solution:\n",
|
||||
"\n",
|
||||
"**(a)** By inspection:\n",
|
||||
"$$\n",
|
||||
"AP = \\begin{pmatrix} d & a & b & c \\\\ ? & ? &? & ? \\\\ ? &? & ? &? \\\\ ? & ? & ? &? \\end{pmatrix},\\; PA = \\begin{pmatrix} ? & ? &? & ? \\\\ ? &? & ? &? \\\\ ? & ? & ? &? \\\\ a & b & c & d\\end{pmatrix}\n",
|
||||
"$$\n",
|
||||
"i.e. $AP$ shifts the columns, and $PA$ shifts the rows. For these two to be equal, the second row of $A$ must be \"d a b c\". Continuing in this process, the third row of $A$ must be \"c d a b\", and the fourth row must be \"b c d a\", i.e. A must be of the form\n",
|
||||
"$$\n",
|
||||
"\\boxed{A = \\begin{pmatrix} a & b & c & d \\\\ d & a &b & c \\\\ c &d & a &b \\\\ b & c & d &a \\end{pmatrix}}\\, ,\n",
|
||||
"$$\n",
|
||||
"which is known as a [circulant matrix](https://en.wikipedia.org/wiki/Circulant_matrix) — it expresses a linear operation known as a [circular discrete convolution](https://en.wikipedia.org/wiki/Convolution#Circular_discrete_convolution).\n",
|
||||
"\n",
|
||||
"**(b)** See below:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 12,
|
||||
"id": "e7222137",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4×4 Matrix{Int64}:\n",
|
||||
" 1 7 3 2\n",
|
||||
" 2 1 7 3\n",
|
||||
" 3 2 1 7\n",
|
||||
" 7 3 2 1"
|
||||
]
|
||||
},
|
||||
"execution_count": 12,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"a, b, c, d = 1, 7, 3, 2 # 4 arbitrarily chosen values\n",
|
||||
"\n",
|
||||
"P = [0 1 0 0\n",
|
||||
" 0 0 1 0\n",
|
||||
" 0 0 0 1\n",
|
||||
" 1 0 0 0]\n",
|
||||
"\n",
|
||||
"F = im .^ ((0:3) .* (0:3)') # the 4×4 Fourier matrix\n",
|
||||
"\n",
|
||||
"# fill in:\n",
|
||||
"A = [a b c d\n",
|
||||
" d a b c\n",
|
||||
" c d a b\n",
|
||||
" b c d a]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 13,
|
||||
"id": "25dc7673",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"true"
|
||||
]
|
||||
},
|
||||
"execution_count": 13,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# check:\n",
|
||||
"P * A == A * P"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "b9344787",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Yay, it commutes!\n",
|
||||
"\n",
|
||||
"Now we want to check that $F$ \"diagonalizes\" $A$, which was defined (in the problem) as (columns of) $F$ being eigenvectors of $A$.\n",
|
||||
"\n",
|
||||
"There are several ways to check this, the easiest being to check that $A$ is diagonal in the **basis of F**, i.e. that $F^{-1} A F$ is diagonal."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 17,
|
||||
"id": "0ae85697",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4×4 Matrix{ComplexF64}:\n",
|
||||
" 13.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im\n",
|
||||
" 0.0+0.0im -2.0+5.0im 0.0+0.0im 0.0+0.0im\n",
|
||||
" 0.0+0.0im 0.0+0.0im -5.0+0.0im 0.0+0.0im\n",
|
||||
" 0.0+0.0im 0.0+0.0im 0.0+0.0im -2.0-5.0im"
|
||||
]
|
||||
},
|
||||
"execution_count": 17,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# A is diagonal in the basis of F:\n",
|
||||
"\n",
|
||||
"F^-1 * A * F # or F \\ A * F, or inv(F) * A * F, or F'/4 * A * F"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "897b25b0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"(From class, note also that $F^{-1} = \\overline{F^T} / 4$: F is unitary up to a scale factor.)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 16,
|
||||
"id": "e6492277",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"true"
|
||||
]
|
||||
},
|
||||
"execution_count": 16,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"inv(F) ≈ F'/4 # note that F' is the conjugate-transpose"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "714942d0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"We could also, of course, compute eigenvectors of $A$ and compare them to F, but this is a lot more work (both for Julia and for us). Julia will give us eigenvectors, but it might give them in a different order and will give them with a different normalization (normalized to unit length), and might give them with different phases or signs:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 19,
|
||||
"id": "dce982d1",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4×4 Matrix{ComplexF64}:\n",
|
||||
" 0.5+0.0im -3.67437e-17+0.5im … -0.5+0.0im\n",
|
||||
" -0.5+0.0im 0.5-0.0im -0.5+0.0im\n",
|
||||
" 0.5+0.0im 1.27272e-16-0.5im -0.5+0.0im\n",
|
||||
" -0.5+0.0im -0.5-4.19352e-16im -0.5+0.0im"
|
||||
]
|
||||
},
|
||||
"execution_count": 19,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using LinearAlgebra\n",
|
||||
"X = eigvecs(A)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "acebc0b5",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Also, it doesn't show the whole matrix to save space. To make it look more like $F$, we can divide the columns of $X$ so that their first row is \"1 1 1 1\", and then call a lower-level `show` function to make sure we show the whole matrix:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 21,
|
||||
"id": "c81825b9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"4×4 Matrix{ComplexF64}:\n",
|
||||
" 1.0+0.0im 1.0-0.0im 1.0+0.0im 1.0-0.0im\n",
|
||||
" -1.0+0.0im -7.34873e-17-1.0im -7.34873e-17+1.0im 1.0-0.0im\n",
|
||||
" 1.0+0.0im -1.0-1.81057e-16im -1.0+1.81057e-16im 1.0-0.0im\n",
|
||||
" -1.0+0.0im -7.65216e-16+1.0im -7.65216e-16-1.0im 1.0-0.0im"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"show(stdout, \"text/plain\", X ./ transpose(X[1,:]))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "72946b10",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"If you stare at this long enough, you'll see that they are indeed the columns of $F$, but in a permuted order:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 24,
|
||||
"id": "cb1467ac",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"true"
|
||||
]
|
||||
},
|
||||
"execution_count": 24,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"X ./ transpose(X[1,:]) ≈ F[:, [3,4,2,1]]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "dd09f649",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"But that is a lot of effort! It's a bit silly to check whether something is an eigenvector by anything other than multiplying it by $A$. (Analogously, if you wanted to check whether $x$ solves $Ax=b$, you wouldn't re-do Gaussian elimination and then compare the result to $x$, you would just check whether $Ax \\approx b$.)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.8.0",
|
||||
"language": "julia",
|
||||
"name": "julia-1.8"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.8.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Reference in new issue
Block a user