clarification
This commit is contained in:
Steven G. Johnson authored and 
GitHub committed
1 parent
139adf700b
commit
c3b34534b0
1 file changed
+499
-499
+499
-499
@@ -1,499 +1,499 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0f2a9965",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# 18.085 Problem Set 2\n",
|
||||
"\n",
|
||||
"Due Friday **March 3** at 1pm."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d9e8ae37",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 1\n",
|
||||
"\n",
|
||||
"**(a)** The eigenvalues of a real *anti-Hermitian* matrix $A = -A^T$ must be \\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_. Derive this by considering $\\overline{x^T} A x$ for an eigenvector $Ax = \\lambda x$, and demonstrate it numerically by constructing a random anti-Hermitian matrix `A = randn(n,n); A = A - A'` in Julia and looking at its eigenvalues `using LinearAlgebra; eigvals(A)` for `n=4` and `n=5`. (You might want to try the numerical experiment first if you aren't sure what the answer is.)\n",
|
||||
"\n",
|
||||
"**(b)** Suppose $A$ is a $3 \\times 3$ real-symmetric matrix with eigenvalues $\\lambda_1 = -2$, $\\lambda_2 = 3$, and $\\lambda_3 = -1$, with corresponding orthonormal eigenvectors $q_1, q_2, q_3$. In terms of these quantities, give the (full) SVD of $A$.\n",
|
||||
"\n",
|
||||
"**(c)** Construct a random $5 \\times 3$ matrix `A = randn(5, 3)` and form a related *real-symmetric* matrix `B = [ 0I A; A' 0I ]`, corresponding to\n",
|
||||
"$$\n",
|
||||
"B = \\begin{pmatrix} 0 & A \\\\ A^T & 0 \\end{pmatrix} .\n",
|
||||
"$$\n",
|
||||
"Compare the eigenvalues of $B$ (`eigvals(B)`) to the singular values of $A$ (`svdvals(A)`). What do you notice? Explain it by using the SVD $A = U \\Sigma V^T$ to construct eigenvalues and eigenvectors of $B$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8110f592",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 2\n",
|
||||
"\n",
|
||||
"**(a)** For any $m \\times n$ real matrix $A$ and any real unitary $m \\times m$ matrix $Q_1$ and any real unitary $n \\times n$ matrix $Q_2$, show that $\\Vert A \\Vert$ = $\\Vert Q_1 A Q_2 \\Vert$ for the norms:\n",
|
||||
"$$\n",
|
||||
"\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2}\\, , \\; \\; \\Vert A \\Vert_F = \\sqrt{\\text{tr}(A^T A) } \\, .\n",
|
||||
"$$\n",
|
||||
"Do *not* use the relationships of these norms to the singular values of $A$, from class; use only the definitions above. (Hint: a change of variables may be useful for the first norm, and the cyclic property of the trace for the second.)\n",
|
||||
"\n",
|
||||
"**(b)** Using the full SVD $A = U \\Sigma V^T$ and the unitary invariance from part (a), show that $\\Vert A \\Vert_2 = \\sigma_1$ and $\\Vert A \\Vert_F = \\sqrt{\\sum_k \\sigma_k^2}$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "32250444",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 3\n",
|
||||
"\n",
|
||||
"Find a closest-rank-1 matrix (in the Frobenius norm, for example) to:\n",
|
||||
"\n",
|
||||
"**(a)** $A = \\begin{pmatrix} 0 & 3 \\\\ 2 & 0 \\end{pmatrix}$\n",
|
||||
"\n",
|
||||
"**(b)** $A = \\begin{pmatrix} \\cos\\theta & -\\sin\\theta \\\\ \\sin\\theta & \\cos\\theta \\end{pmatrix}$ (where $\\theta$ is some real number).\n",
|
||||
"\n",
|
||||
"You should be able to do your calculations completely by hand (it's not too hard, honest!), but of course you may use Julia to check your answers if you wish."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d438a8e0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 4\n",
|
||||
"\n",
|
||||
"For an $m \\times m$ real-symmetric matrix $S = S^T$, we know that we have real eigenvalues $\\lambda_1,\\ldots,\\lambda_m$ and can find an orthonormal basis of eigenvectors $Q = \\begin{pmatrix} q_1 & \\cdots & q_m \\end{pmatrix}$. So, we can write any vector $x \\in \\mathbb{R}^m$ as $x = Qc = q_1 c_1 + \\cdots q_m c_m$ for some coefficient $c$.\n",
|
||||
"\n",
|
||||
"**(a)** Show that $x^T x = c_1^2 + \\cdots + c_m^2$ and $x^T S x = \\lambda_1 c_1^2 + \\cdots + \\lambda_m c_m^2$.\n",
|
||||
"\n",
|
||||
"**(b)** Show that the **Rayleigh quotient**\n",
|
||||
"$$\n",
|
||||
"R(x) = \\frac{x^T S x}{x^T x}\n",
|
||||
"$$\n",
|
||||
"is *maximized* (over *any* possible $x \\ne 0$, not just eigenvectors) by $R(q_1) = \\lambda_1$.\n",
|
||||
"\n",
|
||||
"**(c)** If $A$ is any $m \\times n$ real matrix, we know that the squared singular values $\\sigma_i^2$ are the nonzero eigenvalues of $A^T A$. Use this fact, combined with part (b), to give an alternative proof of why\n",
|
||||
"$$\n",
|
||||
"\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2} = \\sigma_1\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"attachments": {
|
||||
"image.png": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABd4AAAMwCAYAAAA+lnKaAAABYGlDQ1BJQ0MgUHJvZmlsZQAAKJF1kL9LQlEUx7/PDKEe9IOGNzS8oYjISlScWswkAomXJVm0PK+mwdMuzxfRFkHUH1CtbdHe4FKQQ0tTECQ0NdUaFDaUvM7V6mnROVy+H76cc+7hAC5Z59xwA8gXLDM+Pakml5ZVzxNcUCBjACGdFXlY02JUgm9tjeodJKG3o2LWflIaulZ6el9HdlfKx4dvf+tboiOdKTLSD3pBxk0LkHzE2qbFBW8T95m0FPGB4GyDTwWnGnxer1mIR4hviLtZTk8TPxB7U01+tonzxgb72kFsL2cKiXlShV4/phBFjFJFAgH4KQOYQ/SfnmC9J4J1cGzBxBqyyMGi7jA5HAYyxDMogGEMXmI/fGKmuPXvGzreThmY2CN4drzZR+DsQnzteIPjQBd5VxWum/rPZaWqu7ga8De4swS0H9n2yyLgGQZqFdt+L9l27QRouwcuq59VKGSmb1CEgwAAAIplWElmTU0AKgAAAAgABAEaAAUAAAABAAAAPgEbAAUAAAABAAAARgEoAAMAAAABAAIAAIdpAAQAAAABAAAATgAAAAAAAACQAAAAAQAAAJAAAAABAAOShgAHAAAAEgAAAHigAgAEAAAAAQAABd6gAwAEAAAAAQAAAzAAAAAAQVNDSUkAAABTY3JlZW5zaG90lHbXhwAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAAddpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDYuMC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iPgogICAgICAgICA8ZXhpZjpQaXhlbFlEaW1lbnNpb24+ODE2PC9leGlmOlBpeGVsWURpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6UGl4ZWxYRGltZW5zaW9uPjE1MDI8L2V4aWY6UGl4ZWxYRGltZW5zaW9uPgogICAgICAgICA8ZXhpZjpVc2VyQ29tbWVudD5TY3JlZW5zaG90PC9leGlmOlVzZXJDb21tZW50PgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4KWjx6AQAAABxpRE9UAAAAAgAAAAAAAAGYAAAAKAAAAZgAAAGYAAHovQmC9HEAAEAASURBVHgB7J0HmD019f4DiGIvgBURu2LHClhAsCNW7IjYwI4VFbGABQFRVAQU+w97QYooWLCDYOGPCIL6VcFCUVEQK84/n+h7PZvNzJ22u3d3z3me3Zk7k0l5c3JycnKSrFVFCk6OgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AIzAKAmu54X0UHD0SR8ARcAQcAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHICHghndnBEfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYixnVP/7xj/DVr341bL755uFqV7vaYibtaTkCjoAj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AINCDghvcGcGb51W677RYOPPDA8O53vzs84xnPmOWset4cAUfAEXAEViAC5513XvjWt74Vzj///PTH76222irssMMOK7C0K6tIf//738Nb3/rWVH+nn3564Pd1rnOdcPe73z284AUvCLe61a3CJZdcEu52t7uFd77zneG+973vygLASzMHgUsvvTQcf/zx4YILLggXXnhh+ttoo40CuqaTI7CYCHznO99J8simefnLXz7JpnXWWcc+nnN/6qmnhosvvjist9564YpXvGJ697e//S389a9/TQ5Kt7/97eeE9x+OgCPgCDgCjsBqQeBXv/pV+O53vztHz3vQgx4UHvjAB64WCJa8nAtmeD/22GPDVa5ylXCve92rcyFRkn7yk5+Es88+O6y77rrhxje+cfpbSM/uc889N5x55pmJGa9//euHTTbZJNzwhjcMa6+9duf8L/QH4HKnO90p/OUvfwmvf/3rwx577LHQSXr8jkAtAuecc05affHoRz86XOlKV6oNtxJe/Pvf/w6fP+bzYYMNNwj3uMc9VkKRVkQZvF6Wphqf+tSnhve///1zEn/uc58b3vGOd8x5ph+zJCtmKS/CZ7GuJ598cthpp53CGWeckQzrW2+9dbjmNa8Zfvazn4Wjjjoq6RZvectbwmmnnZbq8pOf/GRAvjutXAQwWt7xjnecU8Bb3/rW4cc//vGcZ/6jPQK//OUvw6te9arw61//OrziFa8I97vf/dp/vEpDlvhQULzpTW8KL3/5y/Vz3nXDDTdME0bzXsQHN7/5zcNZZ51VeuXPHAFHYMYQWO2yc7WXf1bYcaXVw7777ht23333OfDivIsTr9MiIVCNTMcdd1y12WabVTH7VfQ4qKJBpHUKJ514UhVnXaro2ZC+Jw77F43N1WGHHVbFbVZax9kUMHr4VG984xuraGifk47SvMY1rlE9/elPr+LgtCmaRX+33XbbTfL74he/eNHTJ0Hq+agjj1qStD3R2UEgDiirq171qokf99prr9nJWMectOXnXXbZJZV1/fXXr6KHaMdUPPhCIeD1slDINsf7pz/9qYrG2er5z3/+pE/61Kc+VfxoMWRF9NatDjnkkCp64BfzoIeLkRelNWvX3/72t1U0UKX62n///edlL3qIVo973OMm9Yk+FA3v88L5g+WFAP1VdNaozTTvf/CDH1Qf/ehHqzgJk+o/Gt5rwy/Wi2n5Xqx89EnHtqONN964TxSr8puf/vSn1Q9/+MMqGgjmyKHoENU4plyzZk11zDHHTMag0Vmriis2qqOPPrr6zW9+s+qwpH++7LLLVl25vcCzi0BbnlypsnO1l392ObOcs5XGh3/4wx+qU045pXr1q1896Vuj4b1c+JGetuX5kZKb+WjCkBweccQR1Uc+8pHqzW9+c7XjjjtWN7rRjSYVKeN1WyP5Bz/4wepyl7tc+v6mN71pFb0aqoMPPrh60YteNBkkKs7oNTLY+M4AJHr3pPSYIIjeXNV+++1XMRDddttt55QDw+IXvvCFIVCN9m1cCjwnbzvvvPNocbeNiHqnLqJ3c0WDclq9CHzsYx+b8CPtZzlSF36+zW1uk8rLZJ0b3mentr1elrYuXvjCF6Z2cYUrXKG66KKLiplZDFmhSf+4VUoxD3q4GHlRWrN21cT9Yx/72NqsYax56EMfOpHtbnivhWrZvMBJI65AbZVf6cBxu6FW4RcyUJd8L2Q++sS9/fbbT9oQ4xqnbghQ94w1GKNp/Fc3sWtjZjxJ+A9/+MP28aq6Z+yNE1s8C2xVldsLO7sIdOHJlSg7V3v5Z5cz63O2EvlQpZXT5EIa3rvwvPK10q+DDO+qNClEpSugT6O4lLWK29JMFKuPf/zjcz7BEz6P+3nPe96cMF1/2FksBuu5J9D973//OWle/epXr373u991TWbU8KweiHsUzsnXIx7xiFHTaBMZkyzUB15R7k3RBrGVG4b6f+1rX1vts88+FZ6Sy5G68DMTX3hhIZOcZgcBr5elrYu73vWuqU+4wx3uUJuRhZYVeLFLT9h7771r88GLhc5LY+JL+JKJcq0onGaQQZ7HPd8Tpm54X8JKGylpdN62hveHPOQhqd5nweO9S75Hgmq0aL7//e9XN7nJTaob3OAG1fve977R4l0tEcXzJSrGONb4cZ/73Gdq8eO2PtX1rne9Ru/4qZEs8wB4/9MfTpPzy7yYnv1lhEAXnlyJsnO1l38ZseokqyuRD1W4DTbYIPURC2l478LzytdKvw4yvL/kJS+pnvKUp1Rc8XxnOZ8Gvrq2Mbyz3YvCc2V5YG7Ei/vzzQnD776ELine truncated
|
||||
}
|
||||
},
|
||||
"cell_type": "markdown",
|
||||
"id": "b945636a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 5\n",
|
||||
"\n",
|
||||
"The following problem is adapted from the [VMLS book](https://web.stanford.edu/~boyd/vmls/), additional [exercise 13.3](https://web.stanford.edu/~boyd/vmls/vmls-additional-exercises.pdf):\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"**(a)** Describe the **size** and **entries** of the matrix $A$ and the vector $b$ for which $J = \\Vert A\\beta - b\\Vert^2 / (N-M)$. This allows you to find the coefficients $\\beta$ that minimize $J$, i.e. the auto-regressive predictor that minimizes the mean-square prediction error, simply by solving a standard least-square problem (`β = A \\ b` in Julia).\n",
|
||||
"\n",
|
||||
"**(b)** Implement an autoregressive model to predict German residential power demand with data from [this tutorial](https://vitalflux.com/autoregressive-ar-models-with-python-examples/#:~:text=Autoregressive%20models%20are%20based%20on,going%20up%20in%20the%20future.) via the code below.\n",
|
||||
"\n",
|
||||
"* Fill in the `???` in the `autoregressive_model` to form `A, b` according to (a) for given data `x` and `M`.\n",
|
||||
"* Apply it for `M=8` to find the prediction coefficients $\\beta$ from \"training\" data consisting of the power-demand data excluding the last 100 days.\n",
|
||||
"* Use the last 100 points as \"test\" data: compare the predictions of your fitted model by plotting the predictions for the last 100 days along with the actual data (fill in the ???)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "722752b9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"\"/var/folders/w_/rf8p8vtd5jb1d8rglz61jsz40000gp/T/jl_8RI3CLu1qv\""
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# download the raw data (CSV file)\n",
|
||||
"import Downloads\n",
|
||||
"datafile = Downloads.download(\"https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "3efb4edd",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4384×5 Matrix{Any}:\n",
|
||||
" \"Date\" \"Consumption\" \"Wind\" \"Solar\" \"Wind+Solar\"\n",
|
||||
" \"2006-01-01\" 1069.18 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-02\" 1380.52 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-03\" 1442.53 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-04\" 1457.22 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-05\" 1477.13 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-06\" 1403.43 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-07\" 1300.29 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-08\" 1207.98 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-09\" 1529.32 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-10\" 1576.91 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-11\" 1577.18 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-12\" 1553.28 \"\" \"\" \"\"\n",
|
||||
" ⋮ \n",
|
||||
" \"2017-12-20\" 1559.69 90.88 8.799 99.679\n",
|
||||
" \"2017-12-21\" 1520.37 259.039 7.313 266.352\n",
|
||||
" \"2017-12-22\" 1423.24 228.773 10.065 238.838\n",
|
||||
" \"2017-12-23\" 1272.17 748.074 8.45 756.524\n",
|
||||
" \"2017-12-24\" 1141.76 812.422 9.949 822.371\n",
|
||||
" \"2017-12-25\" 1111.28 587.81 15.765 603.575\n",
|
||||
" \"2017-12-26\" 1130.12 717.453 30.923 748.376\n",
|
||||
" \"2017-12-27\" 1263.94 394.507 16.53 411.037\n",
|
||||
" \"2017-12-28\" 1299.86 506.424 14.162 520.586\n",
|
||||
" \"2017-12-29\" 1295.09 584.277 29.854 614.131\n",
|
||||
" \"2017-12-30\" 1215.45 721.247 7.467 728.714\n",
|
||||
" \"2017-12-31\" 1107.11 721.176 19.98 741.156"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"summary(days) = \"4383-element Vector{Int64}\"\n",
|
||||
"summary(consumption) = \"4383-element Vector{Float64}\"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAPvCAYAAADTa0gLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd5jUVPfHvzPD0mHpZRkEBEQUbIg0FVQUxQVfG6CCoiugqIjS7f5UqoK9wVpeUaqCL+oLIiAWerOigi8rGzoISxEpO/f3R8xOJpNk0nMzcz7PM89uMpnMmeTm3nPPPSXEGGMgCIIgCIIgCEKXsN8CEARBEARBEEQQIMWZIAiCIAiCIAxAijNBEARBEARBGIAUZ4IgCIIgCIIwACnOBEEQBEEQBGEAUpwJgiAIgiAIwgCkOBMEQRAEQRCEAUhxJgiCIAiCIAgDkOJMEARBEARBEAYgxZkg0pjvv/8eeXl5aNy4McqVK4dy5cqhadOmGDBgANasWeO3eATBDa+++ireeeedpP0FBQUIhUKq7/lF37590bBhQ0uf1fqdBEEYI0QltwkiPXnjjTdw7733olmzZujfvz/OPPNMhEIhbNy4EdOmTcO3336LzZs3o3Hjxn6LShC+06JFC9SoUQNffvllwv5jx45h/fr1aNy4MWrWrOmPcAr69u2LL7/8EgUFBaY/q/U7CYIwRim/BSAIwnm+/fZbDBw4EFdffTVmz56N0qVLl7x36aWX4p577sGsWbNQrlw5R77v6NGjKFu2LEKhkCPnI4xx9OhRx+4hoU6ZMmXQtm1bv8UgCIITyFWDINKQ0aNHIxKJ4I033khQmuXceOONyMnJSdi3Zs0adO/eHdWqVUPZsmVx7rnnYubMmQnHvPPOOwiFQvj8889xxx13oGbNmihfvjyOHTuGTp06oUWLFli+fDnat2+PcuXKoWHDhnj77bcBAJ9++inOO+88lC9fHi1btsT8+fMTzr1582bcfvvtaNq0KcqXL4969eqhW7du+OGHHxKO+/LLLxEKhTBt2jQ8/PDDyMnJQeXKldG5c2f8+uuvKa/PE088gVAohPXr1+O6665D5cqVkZ2djd69e2PPnj0Jx8ZiMYwfPx6nn346ypQpg1q1auHWW2+FIAglx7zyyisIh8PYvXt3yb7nnnsOoVAI99xzT8K5qlatiiFDhpTsO378OJ5++umS89esWRO33357khwNGzZEbm4uPvroI5x77rkoW7YsnnzySd3fOX/+fFx22WXIzs5G+fLl0bx5c4wZMybhmP/85z9o164dypcvj0qVKuHyyy/H8uXLVa/XTz/9hJtuugnZ2dmoXbs27rjjDhQVFSUcO2vWLLRp06bkO0899VTccccdJe9L7UdpLZXuqdwSarc9Gb3PDRs2xE8//YSlS5ciFAohFAqVuEJouWp88803uOyyy1CpUiWUL18e7du3x6effppwjPRblyxZgrvvvhs1atRA9erVcd1112H79u3qN03BO++8g2bNmqFMmTJo3rw5/v3vf6se9+STT6JNmzaoVq0aKleujPPOOw/5+fmQLyrr/c6///4bQ4YMwTnnnIPs7GxUq1YN7dq1w8cff2xIToLIGBhBEGnFyZMnWbly5Vi7du1MfW7x4sWsdOnS7KKLLmIzZsxg8+fPZ3379mUA2Ntvv11y3Ntvv80AsHr16rH+/fuz//73v2z27Nns5MmTrGPHjqx69eqsWbNmLD8/ny1YsIDl5uYyAOzJJ59kLVu2ZNOmTWOfffYZa9u2LStTpgzbtm1bybmXLl3KhgwZwmbPns2WLl3K5syZw/71r3+xcuXKsV9++aXkuCVLljAArGHDhuyWW25hn376KZs2bRo75ZRTWNOmTdnJkyd1f+vjjz/OALAGDRqwYcOGsQULFrCJEyeyChUqsHPPPZcdP3685Nj+/fszAOzee+9l8+fPZ6+//jqrWbMmq1+/PtuzZw9jjLFffvmFAWAffPBByeeuvPJKVq5cOda0adOSfStXrmQA2GeffcYYY6y4uJhdeeWVrEKFCuzJJ59kCxcuZFOmTGH16tVjZ5xxBvvrr79KPtugQQNWt25dduqpp7K33nqLLVmyhK1atUrzN06ZMoWFQiHWqVMn9sEHH7AvvviCvfrqq2zgwIElx7z//vsMALviiivY3Llz2YwZM1irVq1Y6dKl2ddff510vZo1a8Yee+wxtnDhQjZx4kRWpkwZdvvtt5cct2zZMhYKhVivXr3YZ599xhYvXszefvtt1qdPn6T2s2XLlgR5pXu6ZMmSkn1225PR+7xu3Tp26qmnsnPPPZctX76cLV++nK1bt44xxtiWLVuSnoEvv/ySZWVlsVatWrEZM2awuXPnsiuuuIKFQiE2ffr0pN966qmnsvvuu48tWLCATZkyhVWtWpVdcsklmvdO+flrrrmGzZs3j02dOpU1adKE1a9fnzVo0CDh2L59+7L8/Hy2cOFCtnDhQvbUU0+xcuXKsSeffLLkGL3feeDAAda3b1/23nvvscWLF7P58+ezoUOHsnA4zN59992UshJEpkCKM0GkGTt37mQAWK9evZLeO3nyJDtx4kTJKxaLlbx3+umns3PPPZedOHEi4TO5ubmsbt26rLi4mDEWH8xvvfXWpPN37NiRAWBr1qwp2bdv3z4WiURYuXLlEpSaDRs2MADsxRdf1PwtJ0+eZMePH2dNmzZlDzzwQMl+Scnq2rVrwvEzZ85kANjy5cs1z8lYXKGSn5OxuCI5depUxhhjGzduZAASlE3G4grwQw89VLIvGo2yO+64gzHG2LFjx1iFChXYiBEjGAD2xx9/MMYYe+aZZ1hWVhY7fPgwY4yxadOmMQDsww8/TDj/6tWrGQD26quvluxr0KABi0Qi7Ndff9X9bYwxdujQIVa5cmV24YUXJtxjOcXFxSwnJ4e1bNmy5N5Kn61VqxZr37590vUaP358wjkGDhzIypYtW/Idzz77LAPADhw4oCmbWcXZTnsyep8ZY+zMM89kHTt2TJJXTXFu27Ytq1WrFjt06FDJvpMnT7IWLVqwaDRacj2k36psP+PHj2cA2I4dOzSvk3R/zjvvvIR7WFBQwLKyspIUZ+VnT5w4wf7v//6PVa9ePeHzWr9TidRX5OXlsXPPPTfl8QSRKZCrBkFkEK1atUJWVlbJ67nnngMgukj88ssvuOWWWwAAJ0+eLHl17doVO3bsSHKBuP7661W/o27dumjVqlXJdrVq1VCrVi2cc845Ca4hzZs3BwD88ccfJftOnjyJ0aNH44wzzkDp0qVRqlQplC5dGps2bcLGjRuTvqt79+4J22eddVbSOfWQfq9Ejx49UKpUKSxZsgQASv727ds34bgLLrgAzZs3x6JFi0r2XXbZZfjiiy8AAMuWLcNff/2FBx98EDVq1MDChQsBAF988QXatWuHChUqAAA++eQTVKlSBd26dUu45ueccw7q1KmTFMB11lln4bTTTkv5u5YtW4aDBw9i4MCBmn7nv/76K7Zv344+ffogHI4PBRUrVsT111+PFStW4K+//kr4jNr1/vvvv0tcVFq3bg1AvI4zZ87Etm3bUsqaCjvtSSLVfTbDkSNHsHLlStxwww2oWLFiyf5IJII+ffpAEISkZ8VKO5Xuz80335xwDxs0aID27dsnHb948WJ07twZ2dnZiEQiyMrKwmOPPYZ9+/YluBDpMWvWLHTo0AEVK1ZEqVKlkJWVhfz8fNVnjyAyFVKcCSLNqFGjBsqVK6c6KH/wwQdYvXo1/vOf/yTs37VrFwBg6NChCYp1VlYWBg4cCADYu3dvwmfq1q2r+v3VqlVL2le6dOmk/ZLv9d9//12y78EHH8Sjjz6Kf/3rX5g3bx5WrlyJ1atX4+yzz8bRo0eTzlu9evWE7TJlygCA6rFq1KlTJ2G7VKlSqF69Ovbt2wcAJX/VfmtOTk7J+wDQuXNnbN26FZs2bcIXX3yBc889F7Vq1cKll16KL774AkePHsWyZcvQuXPnks/s2rULBw4cQOnSpZOu+86dOw1fcyWS/240GtU8JtVvi8Vi2L9/f8L+VNf74osvxty5c3Hy5EnceuutiEajaNGiBaZNm2ZIbjXstCeJVPfZDPv37wdjTPO6AUg6r5V2Kp1DKbvavlWrVuGKK64AAEyePBnffvstVq9ejYcffjjl90h89NFH6NGjB+rVq4epU6di+fLlWL16Ne644w7Va0oQmQpl1SCINCMSieDSSy/F559/jh07diQM8GeccQYAJAVm1ahRAwAwatQoXHfddarnbdasWcK2Gxk0pk6diltvvRWjR49O2L93715UqVLF8e/buXMn6tWrV7J98uRJ7Nu3r0TRkf7u2LEjSQndvn17yXUDRIszIFqVFy5ciMsvv7xk/yOPPIKvvvoKx44dS1CcpWAxZVCbRKVKlRK2jV5zKW2aPIBRify3Kdm+fTvC4TCqVq1q6PvkXHPNNbjmmmtw7NgxrFixAmPGjMHNN9+Mhg0bol27dihbtiwAMc2bHOUkwUlS3WczVK1aFeFwWPO6AUhoF1aRZNu5c2fSe8p906dPR1ZWFj755JOS6wsAc+fONfx9U6dORaNGjTBjxoyEdqa8TwSR6ZDFmSDSkFGjRqG4uBh33XUXTpw4kfL4Zs2aoWnTpvjuu+9w/vnnq76USpwbhEKhEmucxKeffurIkr8a77//fsL2zJkzcfLkSXTq1AmAmLoPEJUKOatXr8bGjRtLlGVAtNyeccYZ+PDDD7F27doSxfnyyy/Hnj17MHHiRFSuXLnEnQEAcnNzsW/fPhQXF6tec+VkxSjt27dHdnY2Xn/99YSsCnKaNWuGevXq4YMPPkg45siRI/jwww9LMm1YpUyZMujYsSPGjRsHAFi/fj0AlGRx+P777xOOV66COEmq+yzJa8QyW6FCBbRp0wYfffRRwvGxWAxTp05FNBo15E6TimbNmqFu3bqYNm1awv35448Line truncated
|
||||
"text/plain": [
|
||||
"Figure(PyObject <Figure size 800x1200 with 2 Axes>)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(1.0, 200.0)"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# extract the data and plot it:\n",
|
||||
"\n",
|
||||
"using DelimitedFiles, Dates\n",
|
||||
"data = readdlm(datafile, ',')\n",
|
||||
"display(data)\n",
|
||||
"dates = parse.(Date, data[2:end, 1])\n",
|
||||
"days = Dates.value.(dates .- dates[1])\n",
|
||||
"consumption = Float64.(data[2:end, 2])\n",
|
||||
"@show summary(days)\n",
|
||||
"@show summary(consumption)\n",
|
||||
"\n",
|
||||
"using PyPlot\n",
|
||||
"\n",
|
||||
"figure(figsize=(8,12))\n",
|
||||
"subplot(2,1,1)\n",
|
||||
"title(\"German power consumption data\")\n",
|
||||
"plot(days, consumption, \"r.-\")\n",
|
||||
"xlabel(\"days (from 2006-01-01)\")\n",
|
||||
"ylabel(\"power consumption (MWh)\")\n",
|
||||
"\n",
|
||||
"subplot(2,1,2)\n",
|
||||
"title(\"German power consumption data (first 200 days)\")\n",
|
||||
"plot(days, consumption, \"r.-\")\n",
|
||||
"xlabel(\"days (from 2006-01-01)\")\n",
|
||||
"ylabel(\"power consumption (MWh)\")\n",
|
||||
"xlim(1,200)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e12d9551",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# fill in the ??? to get A and b from x and m\n",
|
||||
"function autoregressive_model(x, M)\n",
|
||||
" N = length(x)\n",
|
||||
" A = similar(x, ???, ???) # ??? x ??? matrix of same type as x\n",
|
||||
" for i = 1:size(A,1) # fill in A row-by row\n",
|
||||
" t = M-1+i\n",
|
||||
" A[i,:] = x[???:???] # hint: each row of A comes from a slice of x\n",
|
||||
" end\n",
|
||||
" b = ???\n",
|
||||
" return A, b\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"training_data = consumption[1:end-100]\n",
|
||||
"\n",
|
||||
"# construct your A and b from M=8 and find β\n",
|
||||
"M = 8\n",
|
||||
"A, b = autoregressive_model(training_data[1:end-100], M)\n",
|
||||
"β = A \\ b # least-square fit"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e066bd80",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# test data from last 100+M points:\n",
|
||||
"test_data = consumption[end-99-M:end]\n",
|
||||
"\n",
|
||||
"# plot the predictions for the last 100 points along with test data\n",
|
||||
"predictions = ???\n",
|
||||
"plot(predictions, \"ro-\")\n",
|
||||
"plot(test_data[end-99:end], \"b*-\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "88ac965b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 6\n",
|
||||
"\n",
|
||||
"Show that the pseudo-inverse $A^+$ obeys $AA^+A=A$ and $A^+ A A^+ = A^+$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"attachments": {
|
||||
"image.png": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABs4AAAMCCAYAAADJaDPVAAAMPmlDQ1BJQ0MgUHJvZmlsZQAASImVVwdYU8kWnluSkEBoAQSkhN4EkRpASggt9I4gKiEJEEqMgaBiRxcVXLuIgA1dFVGwAmJBETuLYu+LBRVlXSzYlTcpoOu+8r2Tb+7988+Z/5w5d24ZANSOc0SiXFQdgDxhgTg22J8+LjmFTnoKEKANfxiw43DzRczo6HAAbej8d3t3HXpDu2Iv1fpn/381DR4/nwsAEg1xOi+fmwfxAQDwaq5IXAAAUcqbTS0QSTFsQEsME4R4kRRnynG1FKfL8R6ZT3wsC+J2AJRUOBxxJgCqlyBPL+RmQg3VfogdhTyBEAA1OsQ+eXmTeRCnQWwNfUQQS/UZ6T/oZP5NM31Yk8PJHMbyuchMKUCQL8rlTP8/y/G/LS9XMhTDEjaVLHFIrHTOsG43cyaHSbEKxH3C9MgoiDUh/iDgyfwhRilZkpAEuT9qwM1nwZoBHYgdeZyAMIgNIA4S5kaGK/j0DEEQG2K4QtBpggJ2PMS6EC/i5wfGKXw2iSfHKmKhjRliFlPBn+WIZXGlse5LchKYCv3XWXy2Qh9TLcqKT4KYArF5oSAxEmJViB3yc+LCFD5ji7JYkUM+YkmsNH9ziGP5wmB/uT5WmCEOilX4l+blD80X25QlYEcq8L6CrPgQeX2wdi5Hlj+cC3aJL2QmDOnw88eFD82Fxw8IlM8de8YXJsQpdD6ICvxj5WNxiig3WuGPm/Jzg6W8KcQu+YVxirF4YgFckHJ9PENUEB0vzxMvyuaERsvzwZeDcMACAYAOJLClg8kgGwg6+5r64D95TxDgADHIBHxgr2CGRiTJeoTwGAeKwJ8Q8UH+8Dh/WS8fFEL+6zArP9qDDFlvoWxEDngCcR4IA7nwv0Q2SjgcLRE8hozgH9E5sHFhvrmwSfv/PT/EfmeYkAlXMJKhiHS1IU9iIDGAGEIMItrg+rgP7oWHw6MfbE44A/cYmsd3f8ITQhfhIeEaoZtwa5KgWPxTlhGgG+oHKWqR/mMtcEuo6Yr7495QHSrjOrg+sMddYBwm7gsju0KWpchbWhX6T9p/m8EPV0PhR3Yko+QRZD+y9c8jVW1VXYdVpLX+sT7yXNOH680a7vk5PuuH6vPgOexnT2wRth87g53AzmFHsCZAx1qxZqwDOyrFw6vrsWx1DUWLleWTA3UE/4g3dGWllcx3rHPsdfwi7yvgT5M+owFrsmi6WJCZVUBnwjcCn84Wch1G0Z0cnZwBkL5f5I+vNzGy9wai0/Gdm/8HAN6tg4ODh79zoa0A7HWHt/+h75w1A746lAE4e4grERfKOVx6IMCnhBq80/SAETAD1nA+TsANeAE/EAhCQRSIB8lgIsw+C65zMZgKZoJ5oASUgeVgDagEG8EWsAPsBvtAEzgCToDT4AK4BK6BO3D19IAXoB+8A58RBCEhVISG6CHGiAVihzghDMQHCUTCkVgkGUlDMhEhIkFmIvORMmQlUolsRmqRvcgh5ARyDulCbiEPkF7kNfIJxVAVVAs1RC3R0SgDZaJhaDw6Ac1Ep6BF6AJ0KVqB1qC70Eb0BHoBvYZ2oy/QAQxgypgOZoLZYwyMhUVhKVgGJsZmY6VYOVaD1WMt8DpfwbqxPuwjTsRpOB23hys4BE/AufgUfDa+BK/Ed+CNeDt+BX+A9+PfCFSCAcGO4ElgE8YRMglTCSWEcsI2wkHCKXgv9RDeEYlEHaIV0R3ei8nEbOIM4hLiemID8Tixi/iIOEAikfRIdiRvUhSJQyoglZDWkXaRWkmXST2kD0rKSsZKTkpBSilKQqVipXKlnUrHlC4rPVX6TFYnW5A9yVFkHnk6eRl5K7mFfJHcQ/5M0aBYUbwp8ZRsyjxKBaWecopyl/JGWVnZVNlDOUZZoDxXuUJ5j/JZ5QfKH1U0VWxVWCqpKhKVpSrbVY6r3FJ5Q6VSLal+1BRqAXUptZZ6knqf+kGVpuqgylblqc5RrVJtVL2s+lKNrGahxlSbqFakVq62X+2iWp86Wd1SnaXOUZ+tXqV+SP2G+oAGTWOMRpRGnsYSjZ0a5zSeaZI0LTUDNXmaCzS3aJ7UfETDaGY0Fo1Lm0/bSjtF69EiallpsbWytcq0dmt1avVra2q7aCdqT9Ou0j6q3a2D6VjqsHVydZbp7NO5rvNphOEI5gj+iMUj6kdcHvFed6Suny5ft1S3Qfea7ic9ul6gXo7eCr0mvXv6uL6tfoz+VP0N+qf0+0ZqjfQayR1ZOnLfyNsGqIGtQazBDIMtBh0GA4ZGhsGGIsN1hicN+4x0jPyMso1WGx0z6jWmGfsYC4xXG7caP6dr05n0XHoFvZ3eb2JgEmIiMdls0mny2dTKNMG02LTB9J4ZxYxhlmG22qzNrN/c2DzCfKZ5nfltC7IFwyLLYq3FGYv3llaWSZYLLZssn1npWrGtiqzqrO5aU619radY11hftSHaMGxybNbbXLJFbV1ts2yrbC/aoXZudgK79XZdowijPEYJR9WMumGvYs+0L7Svs3/goOMQ7lDs0OTwcrT56JTRK0afGf3N0dUx13Gr450xmmNCxxSPaRnz2snWietU5XTVmeoc5DzHudn5lYudC99lg8tNV5prhOtC1zbXr27ubmK3erded3P3NPdq9xsMLUY0YwnjrAfBw99jjscRj4+ebp4Fnvs8//Ky98rx2un1bKzVWP7YrWMfeZt6c7w3e3f70H3SfDb5dPua+HJ8a3wf+pn58fy2+T1l2jCzmbuYL/0d/cX+B/3fszxZs1jHA7CA4IDSgM5AzcCEwMrA+0GmQZlBdUH9wa7BM4KPhxBCwkJWhNxgG7K57Fp2f6h76KzQ9jCVsLiwyrCH4bbh4vCWCDQiNGJVxN1Ii0hhZFMUiGJHrYq6F20VPSX6cAwxJjqmKuZJ7JjYmbFn4mhxk+J2xr2L949fFn8nwTpBktCWqJaYmlib+D4pIGllUve40eNmjbuQrJ8sSG5OIaUkpmxLGRgfOH7N+J5U19SS1OsTrCZMm3Buov7E3IlHJ6lN4kzan0ZIS0rbmfaFE8Wp4Qyks9Or0/u5LO5a7gueH281r5fvzV/Jf5rhnbEy41mmd+aqzN4s36zyrD4BS1ApeJUdkr0x+31OVM72nMHcpNyGPKW8tLxDQk1hjrB9stHkaZO7RHaiElH3FM8pa6b0i8PE2/KR/An5zQVa8EO+Q2It+UXyoNCnsKrww9TEqfunaUwTTuuYbjt98fSnRUFFv83AZ3BntM00mTlv5oNZzFmbZyOz02e3zTGbs2BOz9zguTvmUeblzPu92LF4ZfHb+UnzWxYYLpi74NEvwb/UlaiWiEtuLPRauHERvkiwqHOx8+J1i7+V8krPlzmWlZd9WcJdcv7XMb9W/Dq4NGNp5zK3ZRuWE5cLl19f4btix0qNlUUrH62KWNW4mr66dPXbNZPWnCt3Kd+4lrJWsra7IryieZ35uuXrvlRmVV6r8q9qqDaoXlz9fj1v/eUNfhvqNxpuLNv4aZNg083NwZsbayxryrcQtxRuebI1ceuZ3xi/1W7T31a27et24fbuHbE72mvda2t3GuxcVofWSep6d6XuurQ7YHdzvX395gadhrI9YI9kz/O9aXuv7wvb17afsb/+gMWB6oO0g6WNSOP0xv6mrKbu5uTmrkOhh9pavFoOHnY4vP2IyZGqo9pHlx2jHFtwbLC1qHXguOh434nME4/aJrXdOTnu5NX2mPbOU2Gnzp4OOn3yDPNM61nvs0fOeZ47dJ5xvumC24XGDteOg7+7/n6w062z8aL7xeZLHpdausZ2Hbvse/nElYArp6+yr164Fnmt63rC9Zs3Um903+TdfHYr99ar24W3P9+Ze5dwt/Se+r3y+wb3a/6w+aOh26376IOABx0P4x7eecR99OJx/uMvPQueUJ+UPzV+WvvM6dmR3qDeS8/HP+95IXrxua/kT40/q19avzzwl99fHf3j+nteiV8Nvl7yRu/N9rcub9sGogfuv8t79/l96Qe9Dzs+Mj6e+ZT06ennqV9IXyq+2nxt+Rb27e5g3uCgiCPmyD4FMNjQjAwAXm8HgJoMAA3uzyjj5fs/mSHyPasMgf+E5XtEmbkBUA+/32P64NfNDQD2bIXbL6ivlgpANBWAeA+AOjsPt6G9mmxfKTUi3AdsivyanpcO/o3J95w/5P3zGUhVXcDP538B4cx8YTeMD0YAAACKZVhJZk1NACoAAAAIAAQBGgAFAAAAAQAAAD4BGwAFAAAAAQAAAEYBKAADAAAAAQACAACHaQAEAAAAAQAAAE4AAAAAAAAAkAAAAAEAAACQAAAAAQADkoYABwAAABIAAAB4oAIABAAAAAEAAAbOoAMABAAAAAEAAAMCAAAAAEFTQ0lJAAAAU2NyZWVuc2hvdFFtsQ8AAAAJcEhZcwAAFiUAABYlAUlSJPAAAAHXaVRYdFhNTDpjb20uYWRvYmUueG1wAAAAAAA8eDp4bXBtZXRhIHhtbG5zOng9ImFkb2JlOm5zOm1ldGEvIiB4OnhtcHRrPSJYTVAgQ29yZSA2LjAuMCI+CiAgIDxyZGY6UkRGIHhtbG5zOnJkZj0iaHR0cDovL3d3dy53My5vcmcvMTk5OS8wMi8yMi1yZGYtc3ludGF4LW5zIyI+CiAgICAgIDxyZGY6RGVzY3JpcHRpb24gcmRmOmFib3V0PSIiCiAgICAgICAgICAgIHhtbG5zOmV4aWY9Imh0dHA6Ly9ucy5hZG9iZS5jb20vZXhpZi8xLjAvIj4KICAgICAgICAgPGV4aWY6UGl4ZWxZRGltZW5zaW9uPjc3MDwvZXhpZjpQaXhlbFlEaW1lbnNpb24+CiAgICAgICAgIDxleGlmOlBpeGVsWERpbWVuc2lvbj4xNzQyPC9leGlmOlBpeGVsWERpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6Line truncated
|
||||
}
|
||||
},
|
||||
"cell_type": "markdown",
|
||||
"id": "c918611f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 7\n",
|
||||
"\n",
|
||||
"In this problem you will use **principal components analysis (PCA)** and the SVD to classify a images.\n",
|
||||
"\n",
|
||||
"The code below downloads a [collection of 80 monochrome/grayscale drawings](https://math.mit.edu/~stevenj/snow.zip) of [snowpeople](https://en.wikipedia.org/wiki/Snowman) and [snowflakes](https://en.wikipedia.org/wiki/Snowflake) found by a google image search in Spring 2023, shown in miniature here:\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Each image is resized to a common size of $156 \\times 128$, and then reshaped into a vector of $156\\times 156 = 19968$ numbers (the grayscale intensity of each pixel), and stored in a matrix $X$ of \"measurements\" for 80 images:\n",
|
||||
"$$\n",
|
||||
"X = 19968 \\times 80 \\mbox{ matrix} = \\mbox{(data per image)} \\times \\mbox{(# images)}\n",
|
||||
"$$\n",
|
||||
"Think of this as 80 points in a **19968-dimensional space**, i.e. each image is a point in $\\mathbb{R}^{19968}$!!\n",
|
||||
"\n",
|
||||
"PCA allows us to figure out which combinations of these 19968 variables explain most of the variations, and allows us to *project* from $\\mathbb{R}^{19968}$ to a *much* lower-dimensional space in order to help us *classify* the data and identify *clusters* of similar images. In this problem (parts a and b, below), you will perform the PCA clustering analysis yourself.\n",
|
||||
"\n",
|
||||
"You will need to install some Julia packages to help extract the image data. \n",
|
||||
"```jl\n",
|
||||
"import Pkg\n",
|
||||
"Pkg.add(\"ZipFile\")\n",
|
||||
"Pkg.add(\"Images\")\n",
|
||||
"Pkg.add(\"FileIO\")\n",
|
||||
"```\n",
|
||||
"(or type `] add ZipFile Images FileIO`).\n",
|
||||
"\n",
|
||||
"You can then run the following code to download the data and reshape it into our $X$ matrix:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "5f33349c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"\"/var/folders/q9/2vfl1wt51558xwpl4nv4wt580000gn/T/jl_aH2TS1Y9Mx\""
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using ZipFile, Images, FileIO, Statistics, LinearAlgebra, PyPlot\n",
|
||||
"\n",
|
||||
"# download the files as a zip archive:\n",
|
||||
"imagezip = download(\"https://math.mit.edu/~stevenj/snow.zip\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "7be7e0f5",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"rows, cols = 156,128 # common size to resize all of the images to\n",
|
||||
"\n",
|
||||
"# extract the images from the zip file, load them, resize them, and stick them into a big array Xdata\n",
|
||||
"Xdata = Float64[]\n",
|
||||
"r = ZipFile.Reader(imagezip)\n",
|
||||
"for f in r.files\n",
|
||||
" # load image and rescale to rows x cols\n",
|
||||
" image = imresize(load(IOBuffer(read(f))), rows, cols)\n",
|
||||
" imagedata = Float64.(channelview(image))\n",
|
||||
" if length(imagedata) == rows*cols # grayscale\n",
|
||||
" # convert to RGB\n",
|
||||
" bw = vec(imagedata)\n",
|
||||
" imagedata = permutedims([bw bw bw])\n",
|
||||
" elseif size(imagedata, 1) == 4 # RGBA\n",
|
||||
" α = reshape(imagedata[4, :, :], 1, rows, cols) # alpha (transparency)\n",
|
||||
" # convert to RGB, with transparent pixels converted to white:\n",
|
||||
" imagedata = imagedata[1:3, :, :] .* α .+ (1 .- α)\n",
|
||||
" end\n",
|
||||
" # convert back to grayscale (average RGB channels) and collect in Xdata\n",
|
||||
" append!(Xdata, reshape(mean(imagedata, dims=1), rows, cols))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"# array of file names, and indices of snowflakes vs snowpeople\n",
|
||||
"names = [f.name for f in r.files]\n",
|
||||
"snowperson_indices = findall(startswith(\"s\"), names)\n",
|
||||
"snowflake_indices = findall(startswith(\"f\"), names)\n",
|
||||
"\n",
|
||||
"# reshape into 19968x80 matrix X\n",
|
||||
"X = reshape(Xdata, :, length(r.files)) # reshape to (imagesize) x (numimages) matrix\n",
|
||||
"@show size(X)\n",
|
||||
"\n",
|
||||
"# a useful function to reshape a vector into a BW image and plot it\n",
|
||||
"function plotimg(imgdata, rows=rows, cols=cols)\n",
|
||||
" u = copy(imgdata)\n",
|
||||
" u .-= minimum(u)\n",
|
||||
" u ./= maximum(u)\n",
|
||||
" imshow(reshape(u, rows, cols), cmap=\"gray\")\n",
|
||||
" axis(\"off\")\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "14b8590f",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# first, let's plot all of the images to make sure they look okay:\n",
|
||||
"figure(figsize=(6.3, 3.15))\n",
|
||||
"r, c = 7, 14\n",
|
||||
"for i = 1:r, j=1:c\n",
|
||||
" k = (i-1)*c + j\n",
|
||||
" k > size(X,2) && break\n",
|
||||
" subplot(r, c, k)\n",
|
||||
" plotimg(X[:,k])\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4917563f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(a)** Fill in the code below to compute the matrix `A` with the *mean of each pixel* (over all images) subtracted, and take its SVD to plot the singular values.\n",
|
||||
"\n",
|
||||
"(Be careful: some treatments of PCA put the different samples in the rows of $X$ and $A$, and some put the samples in the columns. Which are we doing here?)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "776c2100",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"A = ??? # subtract means from X\n",
|
||||
"\n",
|
||||
"U, σ, V = svd(A);\n",
|
||||
"\n",
|
||||
"plot(σ, \"b.-\")\n",
|
||||
"title(\"Problem 7(a): Singular values of A\");\n",
|
||||
"xlabel(L\"k\");\n",
|
||||
"legend([L\"\\sigma_k\"]);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "52307e61",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(b)** Compute the projections `c1[i]` and `c2[i]` of each image\n",
|
||||
"onto the dominant two principal components (vectors u or v in this case?). That is, reduce the 19968-dimensional data to **2 dimensions**.\n",
|
||||
"\n",
|
||||
"(Does it matter much whether you do the projection to the images in $A$ or $X$? Why or why not?)\n",
|
||||
"\n",
|
||||
"Run the code below to plot `c2` vs `c1` for every image, coloring the points corresponding to snowpeople blue and the points corresponding to snowflakes red. Hopefully, in this 2d projection, you will see\n",
|
||||
"different types of images separated into mostly distinct coordinate regions, with similar images clustered together."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "4f260943",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"c1 = ???\n",
|
||||
"c2 = ???"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "16f9fdc0",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"figure(figsize=(20,15))\n",
|
||||
"plot(c1[snowperson_indices], c2[snowperson_indices], \"b.\")\n",
|
||||
"plot(c1[snowflake_indices], c2[snowflake_indices], \"r.\")\n",
|
||||
"text.(c1, c2, string.(1:length(c1)), fontsize=12);\n",
|
||||
"\n",
|
||||
"title(\"Problem 7(b): Compressed data, visualized on a plane\", fontsize=20);\n",
|
||||
"xlabel(L\"$c_1$\", fontsize=20);\n",
|
||||
"ylabel(L\"$c_2$\", fontsize=20);"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.8.3",
|
||||
"language": "julia",
|
||||
"name": "julia-1.8"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "0f2a9965",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# 18.085 Problem Set 2\n",
|
||||
"\n",
|
||||
"Due Friday **March 3** at 1pm."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d9e8ae37",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 1\n",
|
||||
"\n",
|
||||
"**(a)** The eigenvalues of a real *anti-Hermitian* matrix $A = -A^T$ must be \\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_\\_. Derive this by considering $\\overline{x^T} A x$ for an eigenvector $Ax = \\lambda x$, and demonstrate it numerically by constructing a random anti-Hermitian matrix `A = randn(n,n); A = A - A'` in Julia and looking at its eigenvalues `using LinearAlgebra; eigvals(A)` for `n=4` and `n=5`. (You might want to try the numerical experiment first if you aren't sure what the answer is.)\n",
|
||||
"\n",
|
||||
"**(b)** Suppose $A$ is a $3 \\times 3$ real-symmetric matrix with eigenvalues $\\lambda_1 = -2$, $\\lambda_2 = 3$, and $\\lambda_3 = -1$, with corresponding orthonormal eigenvectors $q_1, q_2, q_3$. In terms of these quantities, give the (full) SVD of $A$.\n",
|
||||
"\n",
|
||||
"**(c)** Construct a random $5 \\times 3$ matrix `A = randn(5, 3)` and form a related *real-symmetric* matrix `B = [ 0I A; A' 0I ]`, corresponding to\n",
|
||||
"$$\n",
|
||||
"B = \\begin{pmatrix} 0 & A \\\\ A^T & 0 \\end{pmatrix} .\n",
|
||||
"$$\n",
|
||||
"Compare the eigenvalues of $B$ (`eigvals(B)`) to the singular values of $A$ (`svdvals(A)`). What do you notice? Explain it by using the SVD $A = U \\Sigma V^T$ to construct eigenvalues and eigenvectors of $B$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "8110f592",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 2\n",
|
||||
"\n",
|
||||
"**(a)** For any $m \\times n$ real matrix $A$ and any real unitary $m \\times m$ matrix $Q_1$ and any real unitary $n \\times n$ matrix $Q_2$, show that $\\Vert A \\Vert$ = $\\Vert Q_1 A Q_2 \\Vert$ for the norms:\n",
|
||||
"$$\n",
|
||||
"\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2}\\, , \\; \\; \\Vert A \\Vert_F = \\sqrt{\\text{tr}(A^T A) } \\, .\n",
|
||||
"$$\n",
|
||||
"Do *not* use the relationships of these norms to the singular values of $A$, from class; use only the definitions above. (Hint: a change of variables may be useful for the first norm, and the cyclic property of the trace for the second.)\n",
|
||||
"\n",
|
||||
"**(b)** Using the full SVD $A = U \\Sigma V^T$ and the unitary invariance from part (a), show that $\\Vert A \\Vert_2 = \\sigma_1$ and $\\Vert A \\Vert_F = \\sqrt{\\sum_k \\sigma_k^2}$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "32250444",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 3\n",
|
||||
"\n",
|
||||
"Find a closest-rank-1 matrix (in the Frobenius norm, for example) to:\n",
|
||||
"\n",
|
||||
"**(a)** $A = \\begin{pmatrix} 0 & 3 \\\\ 2 & 0 \\end{pmatrix}$\n",
|
||||
"\n",
|
||||
"**(b)** $A = \\begin{pmatrix} \\cos\\theta & -\\sin\\theta \\\\ \\sin\\theta & \\cos\\theta \\end{pmatrix}$ (where $\\theta$ is some real number).\n",
|
||||
"\n",
|
||||
"You should be able to do your calculations completely by hand (it's not too hard, honest!), but of course you may use Julia to check your answers if you wish."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d438a8e0",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 4\n",
|
||||
"\n",
|
||||
"For an $m \\times m$ real-symmetric matrix $S = S^T$, we know that we have real eigenvalues $\\lambda_1 \\ge \\lambda_2 \\ge \\cdots \\ge \\lambda_m$ and can find an orthonormal basis of eigenvectors $Q = \\begin{pmatrix} q_1 & \\cdots & q_m \\end{pmatrix}$. So, we can write any vector $x \\in \\mathbb{R}^m$ as $x = Qc = q_1 c_1 + \\cdots q_m c_m$ for some coefficient $c$.\n",
|
||||
"\n",
|
||||
"**(a)** Show that $x^T x = c_1^2 + \\cdots + c_m^2$ and $x^T S x = \\lambda_1 c_1^2 + \\cdots + \\lambda_m c_m^2$.\n",
|
||||
"\n",
|
||||
"**(b)** Show that the **Rayleigh quotient**\n",
|
||||
"$$\n",
|
||||
"R(x) = \\frac{x^T S x}{x^T x}\n",
|
||||
"$$\n",
|
||||
"is *maximized* (over *any* possible $x \\ne 0$, not just eigenvectors) by $R(q_1) = \\lambda_1$.\n",
|
||||
"\n",
|
||||
"**(c)** If $A$ is any $m \\times n$ real matrix, we know that the squared singular values $\\sigma_i^2$ are the nonzero eigenvalues of $A^T A$. Use this fact, combined with part (b), to give an alternative proof of why\n",
|
||||
"$$\n",
|
||||
"\\Vert A \\Vert_2 = \\max_{x\\ne 0} \\frac{\\Vert A x \\Vert_2}{\\Vert x \\Vert_2} = \\sigma_1\n",
|
||||
"$$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"attachments": {
|
||||
"image.png": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABd4AAAMwCAYAAAA+lnKaAAABYGlDQ1BJQ0MgUHJvZmlsZQAAKJF1kL9LQlEUx7/PDKEe9IOGNzS8oYjISlScWswkAomXJVm0PK+mwdMuzxfRFkHUH1CtbdHe4FKQQ0tTECQ0NdUaFDaUvM7V6mnROVy+H76cc+7hAC5Z59xwA8gXLDM+Pakml5ZVzxNcUCBjACGdFXlY02JUgm9tjeodJKG3o2LWflIaulZ6el9HdlfKx4dvf+tboiOdKTLSD3pBxk0LkHzE2qbFBW8T95m0FPGB4GyDTwWnGnxer1mIR4hviLtZTk8TPxB7U01+tonzxgb72kFsL2cKiXlShV4/phBFjFJFAgH4KQOYQ/SfnmC9J4J1cGzBxBqyyMGi7jA5HAYyxDMogGEMXmI/fGKmuPXvGzreThmY2CN4drzZR+DsQnzteIPjQBd5VxWum/rPZaWqu7ga8De4swS0H9n2yyLgGQZqFdt+L9l27QRouwcuq59VKGSmb1CEgwAAAIplWElmTU0AKgAAAAgABAEaAAUAAAABAAAAPgEbAAUAAAABAAAARgEoAAMAAAABAAIAAIdpAAQAAAABAAAATgAAAAAAAACQAAAAAQAAAJAAAAABAAOShgAHAAAAEgAAAHigAgAEAAAAAQAABd6gAwAEAAAAAQAAAzAAAAAAQVNDSUkAAABTY3JlZW5zaG90lHbXhwAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAAddpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDYuMC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iPgogICAgICAgICA8ZXhpZjpQaXhlbFlEaW1lbnNpb24+ODE2PC9leGlmOlBpeGVsWURpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6UGl4ZWxYRGltZW5zaW9uPjE1MDI8L2V4aWY6UGl4ZWxYRGltZW5zaW9uPgogICAgICAgICA8ZXhpZjpVc2VyQ29tbWVudD5TY3JlZW5zaG90PC9leGlmOlVzZXJDb21tZW50PgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4KWjx6AQAAABxpRE9UAAAAAgAAAAAAAAGYAAAAKAAAAZgAAAGYAAHovQmC9HEAAEAASURBVHgB7J0HmD019f4DiGIvgBURu2LHClhAsCNW7IjYwI4VFbGABQFRVAQU+w97QYooWLCDYOGPCIL6VcFCUVEQK84/n+h7PZvNzJ22u3d3z3me3Zk7k0l5c3JycnKSrFVFCk6OgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AIzAKAmu54X0UHD0SR8ARcAQcAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHICHghndnBEfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYHpUj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AIOAKOgCPghnfnAUfAEXAEHAFHwBFwBBwBR8ARcAQcAUfAEXAEHAFHwBFwBBwBR2BEBNzwPiKYixnVP/7xj/DVr341bL755uFqV7vaYibtaTkCjoAj4Ag4Ao6AI+AIOAKOgCPgCDgCjoAj4Ag4Ao6AI+AINCDghvcGcGb51W677RYOPPDA8O53vzs84xnPmOWset4cAUfAEXAEViAC5513XvjWt74Vzj///PTH76222irssMMOK7C0K6tIf//738Nb3/rWVH+nn3564Pd1rnOdcPe73z284AUvCLe61a3CJZdcEu52t7uFd77zneG+973vygLASzMHgUsvvTQcf/zx4YILLggXXnhh+ttoo40CuqaTI7CYCHznO99J8simefnLXz7JpnXWWcc+nnN/6qmnhosvvjist9564YpXvGJ697e//S389a9/TQ5Kt7/97eeE9x+OgCPgCDgCjsBqQeBXv/pV+O53vztHz3vQgx4UHvjAB64WCJa8nAtmeD/22GPDVa5ylXCve92rcyFRkn7yk5+Es88+O6y77rrhxje+cfpbSM/uc889N5x55pmJGa9//euHTTbZJNzwhjcMa6+9duf8L/QH4HKnO90p/OUvfwmvf/3rwx577LHQSXr8jkAtAuecc05affHoRz86XOlKV6oNtxJe/Pvf/w6fP+bzYYMNNwj3uMc9VkKRVkQZvF6Wphqf+tSnhve///1zEn/uc58b3vGOd8x5ph+zJCtmKS/CZ7GuJ598cthpp53CGWeckQzrW2+9dbjmNa8Zfvazn4Wjjjoq6RZvectbwmmnnZbq8pOf/GRAvjutXAQwWt7xjnecU8Bb3/rW4cc//vGcZ/6jPQK//OUvw6te9arw61//OrziFa8I97vf/dp/vEpDlvhQULzpTW8KL3/5y/Vz3nXDDTdME0bzXsQHN7/5zcNZZ51VeuXPHAFHYMYQWO2yc7WXf1bYcaXVw7777ht23333OfDivIsTr9MiIVCNTMcdd1y12WabVTH7VfQ4qKJBpHUKJ514UhVnXaro2ZC+Jw77F43N1WGHHVbFbVZax9kUMHr4VG984xuraGifk47SvMY1rlE9/elPr+LgtCmaRX+33XbbTfL74he/eNHTJ0Hq+agjj1qStD3R2UEgDiirq171qokf99prr9nJWMectOXnXXbZJZV1/fXXr6KHaMdUPPhCIeD1slDINsf7pz/9qYrG2er5z3/+pE/61Kc+VfxoMWRF9NatDjnkkCp64BfzoIeLkRelNWvX3/72t1U0UKX62n///edlL3qIVo973OMm9Yk+FA3v88L5g+WFAP1VdNaozTTvf/CDH1Qf/ehHqzgJk+o/Gt5rwy/Wi2n5Xqx89EnHtqONN964TxSr8puf/vSn1Q9/+MMqGgjmyKHoENU4plyzZk11zDHHTMag0Vmriis2qqOPPrr6zW9+s+qwpH++7LLLVl25vcCzi0BbnlypsnO1l392ObOcs5XGh3/4wx+qU045pXr1q1896Vuj4b1c+JGetuX5kZKb+WjCkBweccQR1Uc+8pHqzW9+c7XjjjtWN7rRjSYVKeN1WyP5Bz/4wepyl7tc+v6mN71pFb0aqoMPPrh60YteNBkkKs7oNTLY+M4AJHr3pPSYIIjeXNV+++1XMRDddttt55QDw+IXvvCFIVCN9m1cCjwnbzvvvPNocbeNiHqnLqJ3c0WDclq9CHzsYx+b8CPtZzlSF36+zW1uk8rLZJ0b3mentr1elrYuXvjCF6Z2cYUrXKG66KKLiplZDFmhSf+4VUoxD3q4GHlRWrN21cT9Yx/72NqsYax56EMfOpHtbnivhWrZvMBJI65AbZVf6cBxu6FW4RcyUJd8L2Q++sS9/fbbT9oQ4xqnbghQ94w1GKNp/Fc3sWtjZjxJ+A9/+MP28aq6Z+yNE1s8C2xVldsLO7sIdOHJlSg7V3v5Z5cz63O2EvlQpZXT5EIa3rvwvPK10q+DDO+qNClEpSugT6O4lLWK29JMFKuPf/zjcz7BEz6P+3nPe96cMF1/2FksBuu5J9D973//OWle/epXr373u991TWbU8KweiHsUzsnXIx7xiFHTaBMZkyzUB15R7k3RBrGVG4b6f+1rX1vts88+FZ6Sy5G68DMTX3hhIZOcZgcBr5elrYu73vWuqU+4wx3uUJuRhZYVeLFLT9h7771r88GLhc5LY+JL+JKJcq0onGaQQZ7HPd8Tpm54X8JKGylpdN62hveHPOQhqd5nweO9S75Hgmq0aL7//e9XN7nJTaob3OAG1fve977R4l0tEcXzJSrGONb4cZ/73Gdq8eO2PtX1rne9Ru/4qZEs8wB4/9MfTpPzy7yYnv1lhEAXnlyJsnO1l38ZseokqyuRD1W4DTbYIPURC2l478LzytdKvw4yvL/kJS+pnvKUp1Rc8XxnOZ8Gvrq2Mbyz3YvCc2V5YG7Ei/vzzQnD776ELine truncated
|
||||
}
|
||||
},
|
||||
"cell_type": "markdown",
|
||||
"id": "b945636a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 5\n",
|
||||
"\n",
|
||||
"The following problem is adapted from the [VMLS book](https://web.stanford.edu/~boyd/vmls/), additional [exercise 13.3](https://web.stanford.edu/~boyd/vmls/vmls-additional-exercises.pdf):\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"**(a)** Describe the **size** and **entries** of the matrix $A$ and the vector $b$ for which $J = \\Vert A\\beta - b\\Vert^2 / (N-M)$. This allows you to find the coefficients $\\beta$ that minimize $J$, i.e. the auto-regressive predictor that minimizes the mean-square prediction error, simply by solving a standard least-square problem (`β = A \\ b` in Julia).\n",
|
||||
"\n",
|
||||
"**(b)** Implement an autoregressive model to predict German residential power demand with data from [this tutorial](https://vitalflux.com/autoregressive-ar-models-with-python-examples/#:~:text=Autoregressive%20models%20are%20based%20on,going%20up%20in%20the%20future.) via the code below.\n",
|
||||
"\n",
|
||||
"* Fill in the `???` in the `autoregressive_model` to form `A, b` according to (a) for given data `x` and `M`.\n",
|
||||
"* Apply it for `M=8` to find the prediction coefficients $\\beta$ from \"training\" data consisting of the power-demand data excluding the last 100 days.\n",
|
||||
"* Use the last 100 points as \"test\" data: compare the predictions of your fitted model by plotting the predictions for the last 100 days along with the actual data (fill in the ???)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 2,
|
||||
"id": "722752b9",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"\"/var/folders/w_/rf8p8vtd5jb1d8rglz61jsz40000gp/T/jl_8RI3CLu1qv\""
|
||||
]
|
||||
},
|
||||
"execution_count": 2,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# download the raw data (CSV file)\n",
|
||||
"import Downloads\n",
|
||||
"datafile = Downloads.download(\"https://raw.githubusercontent.com/jenfly/opsd/master/opsd_germany_daily.csv\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"id": "3efb4edd",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"4384×5 Matrix{Any}:\n",
|
||||
" \"Date\" \"Consumption\" \"Wind\" \"Solar\" \"Wind+Solar\"\n",
|
||||
" \"2006-01-01\" 1069.18 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-02\" 1380.52 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-03\" 1442.53 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-04\" 1457.22 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-05\" 1477.13 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-06\" 1403.43 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-07\" 1300.29 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-08\" 1207.98 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-09\" 1529.32 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-10\" 1576.91 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-11\" 1577.18 \"\" \"\" \"\"\n",
|
||||
" \"2006-01-12\" 1553.28 \"\" \"\" \"\"\n",
|
||||
" ⋮ \n",
|
||||
" \"2017-12-20\" 1559.69 90.88 8.799 99.679\n",
|
||||
" \"2017-12-21\" 1520.37 259.039 7.313 266.352\n",
|
||||
" \"2017-12-22\" 1423.24 228.773 10.065 238.838\n",
|
||||
" \"2017-12-23\" 1272.17 748.074 8.45 756.524\n",
|
||||
" \"2017-12-24\" 1141.76 812.422 9.949 822.371\n",
|
||||
" \"2017-12-25\" 1111.28 587.81 15.765 603.575\n",
|
||||
" \"2017-12-26\" 1130.12 717.453 30.923 748.376\n",
|
||||
" \"2017-12-27\" 1263.94 394.507 16.53 411.037\n",
|
||||
" \"2017-12-28\" 1299.86 506.424 14.162 520.586\n",
|
||||
" \"2017-12-29\" 1295.09 584.277 29.854 614.131\n",
|
||||
" \"2017-12-30\" 1215.45 721.247 7.467 728.714\n",
|
||||
" \"2017-12-31\" 1107.11 721.176 19.98 741.156"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"summary(days) = \"4383-element Vector{Int64}\"\n",
|
||||
"summary(consumption) = \"4383-element Vector{Float64}\"\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAs4AAAPvCAYAAADTa0gLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAEAAElEQVR4nOydd5jUVPfHvzPD0mHpZRkEBEQUbIg0FVQUxQVfG6CCoiugqIjS7f5UqoK9wVpeUaqCL+oLIiAWerOigi8rGzoISxEpO/f3R8xOJpNk0nMzcz7PM89uMpnMmeTm3nPPPSXEGGMgCIIgCIIgCEKXsN8CEARBEARBEEQQIMWZIAiCIAiCIAxAijNBEARBEARBGIAUZ4IgCIIgCIIwACnOBEEQBEEQBGEAUpwJgiAIgiAIwgCkOBMEQRAEQRCEAUhxJgiCIAiCIAgDkOJMEARBEARBEAYgxZkg0pjvv/8eeXl5aNy4McqVK4dy5cqhadOmGDBgANasWeO3eATBDa+++ireeeedpP0FBQUIhUKq7/lF37590bBhQ0uf1fqdBEEYI0QltwkiPXnjjTdw7733olmzZujfvz/OPPNMhEIhbNy4EdOmTcO3336LzZs3o3Hjxn6LShC+06JFC9SoUQNffvllwv5jx45h/fr1aNy4MWrWrOmPcAr69u2LL7/8EgUFBaY/q/U7CYIwRim/BSAIwnm+/fZbDBw4EFdffTVmz56N0qVLl7x36aWX4p577sGsWbNQrlw5R77v6NGjKFu2LEKhkCPnI4xx9OhRx+4hoU6ZMmXQtm1bv8UgCIITyFWDINKQ0aNHIxKJ4I033khQmuXceOONyMnJSdi3Zs0adO/eHdWqVUPZsmVx7rnnYubMmQnHvPPOOwiFQvj8889xxx13oGbNmihfvjyOHTuGTp06oUWLFli+fDnat2+PcuXKoWHDhnj77bcBAJ9++inOO+88lC9fHi1btsT8+fMTzr1582bcfvvtaNq0KcqXL4969eqhW7du+OGHHxKO+/LLLxEKhTBt2jQ8/PDDyMnJQeXKldG5c2f8+uuvKa/PE088gVAohPXr1+O6665D5cqVkZ2djd69e2PPnj0Jx8ZiMYwfPx6nn346ypQpg1q1auHWW2+FIAglx7zyyisIh8PYvXt3yb7nnnsOoVAI99xzT8K5qlatiiFDhpTsO378OJ5++umS89esWRO33357khwNGzZEbm4uPvroI5x77rkoW7YsnnzySd3fOX/+fFx22WXIzs5G+fLl0bx5c4wZMybhmP/85z9o164dypcvj0qVKuHyyy/H8uXLVa/XTz/9hJtuugnZ2dmoXbs27rjjDhQVFSUcO2vWLLRp06bkO0899VTccccdJe9L7UdpLZXuqdwSarc9Gb3PDRs2xE8//YSlS5ciFAohFAqVuEJouWp88803uOyyy1CpUiWUL18e7du3x6effppwjPRblyxZgrvvvhs1atRA9erVcd1112H79u3qN03BO++8g2bNmqFMmTJo3rw5/v3vf6se9+STT6JNmzaoVq0aKleujPPOOw/5+fmQLyrr/c6///4bQ4YMwTnnnIPs7GxUq1YN7dq1w8cff2xIToLIGBhBEGnFyZMnWbly5Vi7du1MfW7x4sWsdOnS7KKLLmIzZsxg8+fPZ3379mUA2Ntvv11y3Ntvv80AsHr16rH+/fuz//73v2z27Nns5MmTrGPHjqx69eqsWbNmLD8/ny1YsIDl5uYyAOzJJ59kLVu2ZNOmTWOfffYZa9u2LStTpgzbtm1bybmXLl3KhgwZwmbPns2WLl3K5syZw/71r3+xcuXKsV9++aXkuCVLljAArGHDhuyWW25hn376KZs2bRo75ZRTWNOmTdnJkyd1f+vjjz/OALAGDRqwYcOGsQULFrCJEyeyChUqsHPPPZcdP3685Nj+/fszAOzee+9l8+fPZ6+//jqrWbMmq1+/PtuzZw9jjLFffvmFAWAffPBByeeuvPJKVq5cOda0adOSfStXrmQA2GeffcYYY6y4uJhdeeWVrEKFCuzJJ59kCxcuZFOmTGH16tVjZ5xxBvvrr79KPtugQQNWt25dduqpp7K33nqLLVmyhK1atUrzN06ZMoWFQiHWqVMn9sEHH7AvvviCvfrqq2zgwIElx7z//vsMALviiivY3Llz2YwZM1irVq1Y6dKl2ddff510vZo1a8Yee+wxtnDhQjZx4kRWpkwZdvvtt5cct2zZMhYKhVivXr3YZ599xhYvXszefvtt1qdPn6T2s2XLlgR5pXu6ZMmSkn1225PR+7xu3Tp26qmnsnPPPZctX76cLV++nK1bt44xxtiWLVuSnoEvv/ySZWVlsVatWrEZM2awuXPnsiuuuIKFQiE2ffr0pN966qmnsvvuu48tWLCATZkyhVWtWpVdcsklmvdO+flrrrmGzZs3j02dOpU1adKE1a9fnzVo0CDh2L59+7L8/Hy2cOFCtnDhQvbUU0+xcuXKsSeffLLkGL3feeDAAda3b1/23nvvscWLF7P58+ezoUOHsnA4zN59992UshJEpkCKM0GkGTt37mQAWK9evZLeO3nyJDtx4kTJKxaLlbx3+umns3PPPZedOHEi4TO5ubmsbt26rLi4mDEWH8xvvfXWpPN37NiRAWBr1qwp2bdv3z4WiURYuXLlEpSaDRs2MADsxRdf1PwtJ0+eZMePH2dNmzZlDzzwQMl+Scnq2rVrwvEzZ85kANjy5cs1z8lYXKGSn5OxuCI5depUxhhjGzduZAASlE3G4grwQw89VLIvGo2yO+64gzHG2LFjx1iFChXYiBEjGAD2xx9/MMYYe+aZZ1hWVhY7fPgwY4yxadOmMQDsww8/TDj/6tWrGQD26quvluxr0KABi0Qi7Ndff9X9bYwxdujQIVa5cmV24YUXJtxjOcXFxSwnJ4e1bNmy5N5Kn61VqxZr37590vUaP358wjkGDhzIypYtW/Idzz77LAPADhw4oCmbWcXZTnsyep8ZY+zMM89kHTt2TJJXTXFu27Ytq1WrFjt06FDJvpMnT7IWLVqwaDRacj2k36psP+PHj2cA2I4dOzSvk3R/zjvvvIR7WFBQwLKyspIUZ+VnT5w4wf7v//6PVa9ePeHzWr9TidRX5OXlsXPPPTfl8QSRKZCrBkFkEK1atUJWVlbJ67nnngMgukj88ssvuOWWWwAAJ0+eLHl17doVO3bsSHKBuP7661W/o27dumjVqlXJdrVq1VCrVi2cc845Ca4hzZs3BwD88ccfJftOnjyJ0aNH44wzzkDp0qVRqlQplC5dGps2bcLGjRuTvqt79+4J22eddVbSOfWQfq9Ejx49UKpUKSxZsgQASv727ds34bgLLrgAzZs3x6JFi0r2XXbZZfjiiy8AAMuWLcNff/2FBx98EDVq1MDChQsBAF988QXatWuHChUqAAA++eQTVKlSBd26dUu45ueccw7q1KmTFMB11lln4bTTTkv5u5YtW4aDBw9i4MCBmn7nv/76K7Zv344+ffogHI4PBRUrVsT111+PFStW4K+//kr4jNr1/vvvv0tcVFq3bg1AvI4zZ87Etm3bUsqaCjvtSSLVfTbDkSNHsHLlStxwww2oWLFiyf5IJII+ffpAEISkZ8VKO5Xuz80335xwDxs0aID27dsnHb948WJ07twZ2dnZiEQiyMrKwmOPPYZ9+/YluBDpMWvWLHTo0AEVK1ZEqVKlkJWVhfz8fNVnjyAyFVKcCSLNqFGjBsqVK6c6KH/wwQdYvXo1/vOf/yTs37VrFwBg6NChCYp1VlYWBg4cCADYu3dvwmfq1q2r+v3VqlVL2le6dOmk/ZLv9d9//12y78EHH8Sjjz6Kf/3rX5g3bx5WrlyJ1atX4+yzz8bRo0eTzlu9evWE7TJlygCA6rFq1KlTJ2G7VKlSqF69Ovbt2wcAJX/VfmtOTk7J+wDQuXNnbN26FZs2bcIXX3yBc889F7Vq1cKll16KL774AkePHsWyZcvQuXPnks/s2rULBw4cQOnSpZOu+86dOw1fcyWS/240GtU8JtVvi8Vi2L9/f8L+VNf74osvxty5c3Hy5EnceuutiEajaNGiBaZNm2ZIbjXstCeJVPfZDPv37wdjTPO6AUg6r5V2Kp1DKbvavlWrVuGKK64AAEyePBnffvstVq9ejYcffjjl90h89NFH6NGjB+rVq4epU6di+fLlWL16Ne644w7Va0oQmQpl1SCINCMSieDSSy/F559/jh07diQM8GeccQYAJAVm1ahRAwAwatQoXHfddarnbdasWcK2Gxk0pk6diltvvRWjR49O2L93715UqVLF8e/buXMn6tWrV7J98uRJ7Nu3r0TRkf7u2LEjSQndvn17yXUDRIszIFqVFy5ciMsvv7xk/yOPPIKvvvoKx44dS1CcpWAxZVCbRKVKlRK2jV5zKW2aPIBRify3Kdm+fTvC4TCqVq1q6PvkXHPNNbjmmmtw7NgxrFixAmPGjMHNN9+Mhg0bol27dihbtiwAMc2bHOUkwUlS3WczVK1aFeFwWPO6AUhoF1aRZNu5c2fSe8p906dPR1ZWFj755JOS6wsAc+fONfx9U6dORaNGjTBjxoyEdqa8TwSR6ZDFmSDSkFGjRqG4uBh33XUXTpw4kfL4Zs2aoWnTpvjuu+9w/vnnq76USpwbhEKhEmucxKeffurIkr8a77//fsL2zJkzcfLkSXTq1AmAmLoPEJUKOatXr8bGjRtLlGVAtNyeccYZ+PDDD7F27doSxfnyyy/Hnj17MHHiRFSuXLnEnQEAcnNzsW/fPhQXF6tec+VkxSjt27dHdnY2Xn/99YSsCnKaNWuGevXq4YMPPkg45siRI/jwww9LMm1YpUyZMujYsSPGjRsHAFi/fj0AlGRx+P777xOOV66COEmq+yzJa8QyW6FCBbRp0wYfffRRwvGxWAxTp05FNBo15E6TimbNmqFu3bqYNm1awv35448Line truncated
|
||||
"text/plain": [
|
||||
"Figure(PyObject <Figure size 800x1200 with 2 Axes>)"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"(1.0, 200.0)"
|
||||
]
|
||||
},
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# extract the data and plot it:\n",
|
||||
"\n",
|
||||
"using DelimitedFiles, Dates\n",
|
||||
"data = readdlm(datafile, ',')\n",
|
||||
"display(data)\n",
|
||||
"dates = parse.(Date, data[2:end, 1])\n",
|
||||
"days = Dates.value.(dates .- dates[1])\n",
|
||||
"consumption = Float64.(data[2:end, 2])\n",
|
||||
"@show summary(days)\n",
|
||||
"@show summary(consumption)\n",
|
||||
"\n",
|
||||
"using PyPlot\n",
|
||||
"\n",
|
||||
"figure(figsize=(8,12))\n",
|
||||
"subplot(2,1,1)\n",
|
||||
"title(\"German power consumption data\")\n",
|
||||
"plot(days, consumption, \"r.-\")\n",
|
||||
"xlabel(\"days (from 2006-01-01)\")\n",
|
||||
"ylabel(\"power consumption (MWh)\")\n",
|
||||
"\n",
|
||||
"subplot(2,1,2)\n",
|
||||
"title(\"German power consumption data (first 200 days)\")\n",
|
||||
"plot(days, consumption, \"r.-\")\n",
|
||||
"xlabel(\"days (from 2006-01-01)\")\n",
|
||||
"ylabel(\"power consumption (MWh)\")\n",
|
||||
"xlim(1,200)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e12d9551",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# fill in the ??? to get A and b from x and m\n",
|
||||
"function autoregressive_model(x, M)\n",
|
||||
" N = length(x)\n",
|
||||
" A = similar(x, ???, ???) # ??? x ??? matrix of same type as x\n",
|
||||
" for i = 1:size(A,1) # fill in A row-by row\n",
|
||||
" t = M-1+i\n",
|
||||
" A[i,:] = x[???:???] # hint: each row of A comes from a slice of x\n",
|
||||
" end\n",
|
||||
" b = ???\n",
|
||||
" return A, b\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"training_data = consumption[1:end-100]\n",
|
||||
"\n",
|
||||
"# construct your A and b from M=8 and find β\n",
|
||||
"M = 8\n",
|
||||
"A, b = autoregressive_model(training_data[1:end-100], M)\n",
|
||||
"β = A \\ b # least-square fit"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "e066bd80",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# test data from last 100+M points:\n",
|
||||
"test_data = consumption[end-99-M:end]\n",
|
||||
"\n",
|
||||
"# plot the predictions for the last 100 points along with test data\n",
|
||||
"predictions = ???\n",
|
||||
"plot(predictions, \"ro-\")\n",
|
||||
"plot(test_data[end-99:end], \"b*-\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "88ac965b",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 6\n",
|
||||
"\n",
|
||||
"Show that the pseudo-inverse $A^+$ obeys $AA^+A=A$ and $A^+ A A^+ = A^+$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"attachments": {
|
||||
"image.png": {
|
||||
"image/png": "iVBORw0KGgoAAAANSUhEUgAABs4AAAMCCAYAAADJaDPVAAAMPmlDQ1BJQ0MgUHJvZmlsZQAASImVVwdYU8kWnluSkEBoAQSkhN4EkRpASggt9I4gKiEJEEqMgaBiRxcVXLuIgA1dFVGwAmJBETuLYu+LBRVlXSzYlTcpoOu+8r2Tb+7988+Z/5w5d24ZANSOc0SiXFQdgDxhgTg22J8+LjmFTnoKEKANfxiw43DzRczo6HAAbej8d3t3HXpDu2Iv1fpn/381DR4/nwsAEg1xOi+fmwfxAQDwaq5IXAAAUcqbTS0QSTFsQEsME4R4kRRnynG1FKfL8R6ZT3wsC+J2AJRUOBxxJgCqlyBPL+RmQg3VfogdhTyBEAA1OsQ+eXmTeRCnQWwNfUQQS/UZ6T/oZP5NM31Yk8PJHMbyuchMKUCQL8rlTP8/y/G/LS9XMhTDEjaVLHFIrHTOsG43cyaHSbEKxH3C9MgoiDUh/iDgyfwhRilZkpAEuT9qwM1nwZoBHYgdeZyAMIgNIA4S5kaGK/j0DEEQG2K4QtBpggJ2PMS6EC/i5wfGKXw2iSfHKmKhjRliFlPBn+WIZXGlse5LchKYCv3XWXy2Qh9TLcqKT4KYArF5oSAxEmJViB3yc+LCFD5ji7JYkUM+YkmsNH9ziGP5wmB/uT5WmCEOilX4l+blD80X25QlYEcq8L6CrPgQeX2wdi5Hlj+cC3aJL2QmDOnw88eFD82Fxw8IlM8de8YXJsQpdD6ICvxj5WNxiig3WuGPm/Jzg6W8KcQu+YVxirF4YgFckHJ9PENUEB0vzxMvyuaERsvzwZeDcMACAYAOJLClg8kgGwg6+5r64D95TxDgADHIBHxgr2CGRiTJeoTwGAeKwJ8Q8UH+8Dh/WS8fFEL+6zArP9qDDFlvoWxEDngCcR4IA7nwv0Q2SjgcLRE8hozgH9E5sHFhvrmwSfv/PT/EfmeYkAlXMJKhiHS1IU9iIDGAGEIMItrg+rgP7oWHw6MfbE44A/cYmsd3f8ITQhfhIeEaoZtwa5KgWPxTlhGgG+oHKWqR/mMtcEuo6Yr7495QHSrjOrg+sMddYBwm7gsju0KWpchbWhX6T9p/m8EPV0PhR3Yko+QRZD+y9c8jVW1VXYdVpLX+sT7yXNOH680a7vk5PuuH6vPgOexnT2wRth87g53AzmFHsCZAx1qxZqwDOyrFw6vrsWx1DUWLleWTA3UE/4g3dGWllcx3rHPsdfwi7yvgT5M+owFrsmi6WJCZVUBnwjcCn84Wch1G0Z0cnZwBkL5f5I+vNzGy9wai0/Gdm/8HAN6tg4ODh79zoa0A7HWHt/+h75w1A746lAE4e4grERfKOVx6IMCnhBq80/SAETAD1nA+TsANeAE/EAhCQRSIB8lgIsw+C65zMZgKZoJ5oASUgeVgDagEG8EWsAPsBvtAEzgCToDT4AK4BK6BO3D19IAXoB+8A58RBCEhVISG6CHGiAVihzghDMQHCUTCkVgkGUlDMhEhIkFmIvORMmQlUolsRmqRvcgh5ARyDulCbiEPkF7kNfIJxVAVVAs1RC3R0SgDZaJhaDw6Ac1Ep6BF6AJ0KVqB1qC70Eb0BHoBvYZ2oy/QAQxgypgOZoLZYwyMhUVhKVgGJsZmY6VYOVaD1WMt8DpfwbqxPuwjTsRpOB23hys4BE/AufgUfDa+BK/Ed+CNeDt+BX+A9+PfCFSCAcGO4ElgE8YRMglTCSWEcsI2wkHCKXgv9RDeEYlEHaIV0R3ei8nEbOIM4hLiemID8Tixi/iIOEAikfRIdiRvUhSJQyoglZDWkXaRWkmXST2kD0rKSsZKTkpBSilKQqVipXKlnUrHlC4rPVX6TFYnW5A9yVFkHnk6eRl5K7mFfJHcQ/5M0aBYUbwp8ZRsyjxKBaWecopyl/JGWVnZVNlDOUZZoDxXuUJ5j/JZ5QfKH1U0VWxVWCqpKhKVpSrbVY6r3FJ5Q6VSLal+1BRqAXUptZZ6knqf+kGVpuqgylblqc5RrVJtVL2s+lKNrGahxlSbqFakVq62X+2iWp86Wd1SnaXOUZ+tXqV+SP2G+oAGTWOMRpRGnsYSjZ0a5zSeaZI0LTUDNXmaCzS3aJ7UfETDaGY0Fo1Lm0/bSjtF69EiallpsbWytcq0dmt1avVra2q7aCdqT9Ou0j6q3a2D6VjqsHVydZbp7NO5rvNphOEI5gj+iMUj6kdcHvFed6Suny5ft1S3Qfea7ic9ul6gXo7eCr0mvXv6uL6tfoz+VP0N+qf0+0ZqjfQayR1ZOnLfyNsGqIGtQazBDIMtBh0GA4ZGhsGGIsN1hicN+4x0jPyMso1WGx0z6jWmGfsYC4xXG7caP6dr05n0XHoFvZ3eb2JgEmIiMdls0mny2dTKNMG02LTB9J4ZxYxhlmG22qzNrN/c2DzCfKZ5nfltC7IFwyLLYq3FGYv3llaWSZYLLZssn1npWrGtiqzqrO5aU619radY11hftSHaMGxybNbbXLJFbV1ts2yrbC/aoXZudgK79XZdowijPEYJR9WMumGvYs+0L7Svs3/goOMQ7lDs0OTwcrT56JTRK0afGf3N0dUx13Gr450xmmNCxxSPaRnz2snWietU5XTVmeoc5DzHudn5lYudC99lg8tNV5prhOtC1zbXr27ubmK3erded3P3NPdq9xsMLUY0YwnjrAfBw99jjscRj4+ebp4Fnvs8//Ky98rx2un1bKzVWP7YrWMfeZt6c7w3e3f70H3SfDb5dPua+HJ8a3wf+pn58fy2+T1l2jCzmbuYL/0d/cX+B/3fszxZs1jHA7CA4IDSgM5AzcCEwMrA+0GmQZlBdUH9wa7BM4KPhxBCwkJWhNxgG7K57Fp2f6h76KzQ9jCVsLiwyrCH4bbh4vCWCDQiNGJVxN1Ii0hhZFMUiGJHrYq6F20VPSX6cAwxJjqmKuZJ7JjYmbFn4mhxk+J2xr2L949fFn8nwTpBktCWqJaYmlib+D4pIGllUve40eNmjbuQrJ8sSG5OIaUkpmxLGRgfOH7N+J5U19SS1OsTrCZMm3Buov7E3IlHJ6lN4kzan0ZIS0rbmfaFE8Wp4Qyks9Or0/u5LO5a7gueH281r5fvzV/Jf5rhnbEy41mmd+aqzN4s36zyrD4BS1ApeJUdkr0x+31OVM72nMHcpNyGPKW8tLxDQk1hjrB9stHkaZO7RHaiElH3FM8pa6b0i8PE2/KR/An5zQVa8EO+Q2It+UXyoNCnsKrww9TEqfunaUwTTuuYbjt98fSnRUFFv83AZ3BntM00mTlv5oNZzFmbZyOz02e3zTGbs2BOz9zguTvmUeblzPu92LF4ZfHb+UnzWxYYLpi74NEvwb/UlaiWiEtuLPRauHERvkiwqHOx8+J1i7+V8krPlzmWlZd9WcJdcv7XMb9W/Dq4NGNp5zK3ZRuWE5cLl19f4btix0qNlUUrH62KWNW4mr66dPXbNZPWnCt3Kd+4lrJWsra7IryieZ35uuXrvlRmVV6r8q9qqDaoXlz9fj1v/eUNfhvqNxpuLNv4aZNg083NwZsbayxryrcQtxRuebI1ceuZ3xi/1W7T31a27et24fbuHbE72mvda2t3GuxcVofWSep6d6XuurQ7YHdzvX395gadhrI9YI9kz/O9aXuv7wvb17afsb/+gMWB6oO0g6WNSOP0xv6mrKbu5uTmrkOhh9pavFoOHnY4vP2IyZGqo9pHlx2jHFtwbLC1qHXguOh434nME4/aJrXdOTnu5NX2mPbOU2Gnzp4OOn3yDPNM61nvs0fOeZ47dJ5xvumC24XGDteOg7+7/n6w062z8aL7xeZLHpdausZ2Hbvse/nElYArp6+yr164Fnmt63rC9Zs3Um903+TdfHYr99ar24W3P9+Ze5dwt/Se+r3y+wb3a/6w+aOh26376IOABx0P4x7eecR99OJx/uMvPQueUJ+UPzV+WvvM6dmR3qDeS8/HP+95IXrxua/kT40/q19avzzwl99fHf3j+nteiV8Nvl7yRu/N9rcub9sGogfuv8t79/l96Qe9Dzs+Mj6e+ZT06ennqV9IXyq+2nxt+Rb27e5g3uCgiCPmyD4FMNjQjAwAXm8HgJoMAA3uzyjj5fs/mSHyPasMgf+E5XtEmbkBUA+/32P64NfNDQD2bIXbL6ivlgpANBWAeA+AOjsPt6G9mmxfKTUi3AdsivyanpcO/o3J95w/5P3zGUhVXcDP538B4cx8YTeMD0YAAACKZVhJZk1NACoAAAAIAAQBGgAFAAAAAQAAAD4BGwAFAAAAAQAAAEYBKAADAAAAAQACAACHaQAEAAAAAQAAAE4AAAAAAAAAkAAAAAEAAACQAAAAAQADkoYABwAAABIAAAB4oAIABAAAAAEAAAbOoAMABAAAAAEAAAMCAAAAAEFTQ0lJAAAAU2NyZWVuc2hvdFFtsQ8AAAAJcEhZcwAAFiUAABYlAUlSJPAAAAHXaVRYdFhNTDpjb20uYWRvYmUueG1wAAAAAAA8eDp4bXBtZXRhIHhtbG5zOng9ImFkb2JlOm5zOm1ldGEvIiB4OnhtcHRrPSJYTVAgQ29yZSA2LjAuMCI+CiAgIDxyZGY6UkRGIHhtbG5zOnJkZj0iaHR0cDovL3d3dy53My5vcmcvMTk5OS8wMi8yMi1yZGYtc3ludGF4LW5zIyI+CiAgICAgIDxyZGY6RGVzY3JpcHRpb24gcmRmOmFib3V0PSIiCiAgICAgICAgICAgIHhtbG5zOmV4aWY9Imh0dHA6Ly9ucy5hZG9iZS5jb20vZXhpZi8xLjAvIj4KICAgICAgICAgPGV4aWY6UGl4ZWxZRGltZW5zaW9uPjc3MDwvZXhpZjpQaXhlbFlEaW1lbnNpb24+CiAgICAgICAgIDxleGlmOlBpeGVsWERpbWVuc2lvbj4xNzQyPC9leGlmOlBpeGVsWERpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6Line truncated
|
||||
}
|
||||
},
|
||||
"cell_type": "markdown",
|
||||
"id": "c918611f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Problem 7\n",
|
||||
"\n",
|
||||
"In this problem you will use **principal components analysis (PCA)** and the SVD to classify a images.\n",
|
||||
"\n",
|
||||
"The code below downloads a [collection of 80 monochrome/grayscale drawings](https://math.mit.edu/~stevenj/snow.zip) of [snowpeople](https://en.wikipedia.org/wiki/Snowman) and [snowflakes](https://en.wikipedia.org/wiki/Snowflake) found by a google image search in Spring 2023, shown in miniature here:\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Each image is resized to a common size of $156 \\times 128$, and then reshaped into a vector of $156\\times 156 = 19968$ numbers (the grayscale intensity of each pixel), and stored in a matrix $X$ of \"measurements\" for 80 images:\n",
|
||||
"$$\n",
|
||||
"X = 19968 \\times 80 \\mbox{ matrix} = \\mbox{(data per image)} \\times \\mbox{(# images)}\n",
|
||||
"$$\n",
|
||||
"Think of this as 80 points in a **19968-dimensional space**, i.e. each image is a point in $\\mathbb{R}^{19968}$!!\n",
|
||||
"\n",
|
||||
"PCA allows us to figure out which combinations of these 19968 variables explain most of the variations, and allows us to *project* from $\\mathbb{R}^{19968}$ to a *much* lower-dimensional space in order to help us *classify* the data and identify *clusters* of similar images. In this problem (parts a and b, below), you will perform the PCA clustering analysis yourself.\n",
|
||||
"\n",
|
||||
"You will need to install some Julia packages to help extract the image data. \n",
|
||||
"```jl\n",
|
||||
"import Pkg\n",
|
||||
"Pkg.add(\"ZipFile\")\n",
|
||||
"Pkg.add(\"Images\")\n",
|
||||
"Pkg.add(\"FileIO\")\n",
|
||||
"```\n",
|
||||
"(or type `] add ZipFile Images FileIO`).\n",
|
||||
"\n",
|
||||
"You can then run the following code to download the data and reshape it into our $X$ matrix:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 1,
|
||||
"id": "5f33349c",
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"\"/var/folders/q9/2vfl1wt51558xwpl4nv4wt580000gn/T/jl_aH2TS1Y9Mx\""
|
||||
]
|
||||
},
|
||||
"execution_count": 1,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"using ZipFile, Images, FileIO, Statistics, LinearAlgebra, PyPlot\n",
|
||||
"\n",
|
||||
"# download the files as a zip archive:\n",
|
||||
"imagezip = download(\"https://math.mit.edu/~stevenj/snow.zip\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "7be7e0f5",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"rows, cols = 156,128 # common size to resize all of the images to\n",
|
||||
"\n",
|
||||
"# extract the images from the zip file, load them, resize them, and stick them into a big array Xdata\n",
|
||||
"Xdata = Float64[]\n",
|
||||
"r = ZipFile.Reader(imagezip)\n",
|
||||
"for f in r.files\n",
|
||||
" # load image and rescale to rows x cols\n",
|
||||
" image = imresize(load(IOBuffer(read(f))), rows, cols)\n",
|
||||
" imagedata = Float64.(channelview(image))\n",
|
||||
" if length(imagedata) == rows*cols # grayscale\n",
|
||||
" # convert to RGB\n",
|
||||
" bw = vec(imagedata)\n",
|
||||
" imagedata = permutedims([bw bw bw])\n",
|
||||
" elseif size(imagedata, 1) == 4 # RGBA\n",
|
||||
" α = reshape(imagedata[4, :, :], 1, rows, cols) # alpha (transparency)\n",
|
||||
" # convert to RGB, with transparent pixels converted to white:\n",
|
||||
" imagedata = imagedata[1:3, :, :] .* α .+ (1 .- α)\n",
|
||||
" end\n",
|
||||
" # convert back to grayscale (average RGB channels) and collect in Xdata\n",
|
||||
" append!(Xdata, reshape(mean(imagedata, dims=1), rows, cols))\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"# array of file names, and indices of snowflakes vs snowpeople\n",
|
||||
"names = [f.name for f in r.files]\n",
|
||||
"snowperson_indices = findall(startswith(\"s\"), names)\n",
|
||||
"snowflake_indices = findall(startswith(\"f\"), names)\n",
|
||||
"\n",
|
||||
"# reshape into 19968x80 matrix X\n",
|
||||
"X = reshape(Xdata, :, length(r.files)) # reshape to (imagesize) x (numimages) matrix\n",
|
||||
"@show size(X)\n",
|
||||
"\n",
|
||||
"# a useful function to reshape a vector into a BW image and plot it\n",
|
||||
"function plotimg(imgdata, rows=rows, cols=cols)\n",
|
||||
" u = copy(imgdata)\n",
|
||||
" u .-= minimum(u)\n",
|
||||
" u ./= maximum(u)\n",
|
||||
" imshow(reshape(u, rows, cols), cmap=\"gray\")\n",
|
||||
" axis(\"off\")\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "14b8590f",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# first, let's plot all of the images to make sure they look okay:\n",
|
||||
"figure(figsize=(6.3, 3.15))\n",
|
||||
"r, c = 7, 14\n",
|
||||
"for i = 1:r, j=1:c\n",
|
||||
" k = (i-1)*c + j\n",
|
||||
" k > size(X,2) && break\n",
|
||||
" subplot(r, c, k)\n",
|
||||
" plotimg(X[:,k])\n",
|
||||
"end"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "4917563f",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(a)** Fill in the code below to compute the matrix `A` with the *mean of each pixel* (over all images) subtracted, and take its SVD to plot the singular values.\n",
|
||||
"\n",
|
||||
"(Be careful: some treatments of PCA put the different samples in the rows of $X$ and $A$, and some put the samples in the columns. Which are we doing here?)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "776c2100",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"A = ??? # subtract means from X\n",
|
||||
"\n",
|
||||
"U, σ, V = svd(A);\n",
|
||||
"\n",
|
||||
"plot(σ, \"b.-\")\n",
|
||||
"title(\"Problem 7(a): Singular values of A\");\n",
|
||||
"xlabel(L\"k\");\n",
|
||||
"legend([L\"\\sigma_k\"]);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "52307e61",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"**(b)** Compute the projections `c1[i]` and `c2[i]` of each image\n",
|
||||
"onto the dominant two principal components (vectors u or v in this case?). That is, reduce the 19968-dimensional data to **2 dimensions**.\n",
|
||||
"\n",
|
||||
"(Does it matter much whether you do the projection to the images in $A$ or $X$? Why or why not?)\n",
|
||||
"\n",
|
||||
"Run the code below to plot `c2` vs `c1` for every image, coloring the points corresponding to snowpeople blue and the points corresponding to snowflakes red. Hopefully, in this 2d projection, you will see\n",
|
||||
"different types of images separated into mostly distinct coordinate regions, with similar images clustered together."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "4f260943",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"c1 = ???\n",
|
||||
"c2 = ???"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "16f9fdc0",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"figure(figsize=(20,15))\n",
|
||||
"plot(c1[snowperson_indices], c2[snowperson_indices], \"b.\")\n",
|
||||
"plot(c1[snowflake_indices], c2[snowflake_indices], \"r.\")\n",
|
||||
"text.(c1, c2, string.(1:length(c1)), fontsize=12);\n",
|
||||
"\n",
|
||||
"title(\"Problem 7(b): Compressed data, visualized on a plane\", fontsize=20);\n",
|
||||
"xlabel(L\"$c_1$\", fontsize=20);\n",
|
||||
"ylabel(L\"$c_2$\", fontsize=20);"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Julia 1.8.3",
|
||||
"language": "julia",
|
||||
"name": "julia-1.8"
|
||||
},
|
||||
"language_info": {
|
||||
"file_extension": ".jl",
|
||||
"mimetype": "application/julia",
|
||||
"name": "julia",
|
||||
"version": "1.8.5"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
Reference in new issue
Block a user