problems 6-7
This commit is contained in:
Steven G. Johnson committed
1 parent
15d90574fa
commit
b70d326baa
1 file changed
+195
-4
+195
-4
@@ -278,16 +278,207 @@
|
|||||||
"id": "88ac965b",
|
"id": "88ac965b",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"source": [
|
"source": [
|
||||||
"## Problems 6, etc: coming soon"
|
"## 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",
|
"cell_type": "code",
|
||||||
"execution_count": null,
|
"execution_count": null,
|
||||||
"id": "e1fca19a",
|
"id": "7be7e0f5",
|
||||||
"metadata": {},
|
"metadata": {},
|
||||||
"outputs": [],
|
"outputs": [],
|
||||||
"source": []
|
"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": {
|
"metadata": {
|
||||||
@@ -300,7 +491,7 @@
|
|||||||
"file_extension": ".jl",
|
"file_extension": ".jl",
|
||||||
"mimetype": "application/julia",
|
"mimetype": "application/julia",
|
||||||
"name": "julia",
|
"name": "julia",
|
||||||
"version": "1.8.3"
|
"version": "1.8.5"
|
||||||
}
|
}
|
||||||
},
|
},
|
||||||
"nbformat": 4,
|
"nbformat": 4,
|
||||||
|
|||||||
Reference in new issue
Block a user