{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Matrix properties\n", "This notebook aims to look at general matrix properties and operations, such as the differences between eigendecomposition, singular value decomposition, PCA and fourier analysis of general matrices and covariance matrices." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "ExecuteTime": { "end_time": "2018-11-30T00:05:11.859028Z", "start_time": "2018-11-30T00:05:11.639371Z" } }, "outputs": [], "source": [ "import numpy as np\n", "%matplotlib notebook\n", "import matplotlib.pyplot as plt\n", "import scipy\n", "np.set_printoptions(precision=4, linewidth=100, suppress=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Eigenvectors and Eigenvalues\n", "If we have a linear operator $T: V \\rightarrow V$ (T is a linear transformation from a vector space $V$ to itself), then applying $T$ to an eigenvector of $T$ (let us call it $\\mathbf{v}$) only scales $\\mathbf{v}$ and doesn't change its direction. I.e. \n", "$$\n", "T(\\mathbf{v}) = \\lambda \\mathbf{v}\n", "$$\n", "\n", "The mapping from $V\\rightarrow V$ is the bit that enforces eigendecomposition to work only on square matrices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the simplest case, $T$ is a matrix operation $A$, so $A\\mathbf{v} = \\lambda \\mathbf{v}$, but there are other linear operations like differentiation which has eigenvectors (or eigenfunctions). E.g. for differentiation $$\\frac{d}{dx}e^{\\lambda x} = \\lambda e^{\\lambda x}$$ So $e^{\\lambda x}$ is an eigenfunction of differentiation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One way to find the eigenvectors/eigenvalues is by solving the characteristic polynomial. I.e. \n", "$$\n", "|A-\\lambda I| = 0\n", "$$\n", "\n", "This is difficult to do for anything greater than a $2\\times 2$ matrix as it will involve finding the roots of a high order polynomial.\n", "\n", "Once we have found the eigenvectors and eigenvalues of $A$, then the eigendecomposition of A is:\n", "$$\n", "A = Q\\Lambda Q^{-1}\n", "$$\n", "\n", "where \n", "\n", "$$\n", "Q = \\left[\\mathbf{v}_1 | \\mathbf{v}_2 | \\cdots | \\mathbf{v}_n \\right]\n", "$$\n", "\n", "A nice way to interpret this is to consider $Q^{-1}$ as a change of basis matrix. Then we are equating $A$ and $\\Lambda$ to the same linear transformation, just that we need to change the basis by applying $Q^{-1}$, then apply the transformation $\\Lambda$, then change back to the original basis with $Q$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Eigendecomposition of Symmetric Matrices\n", "When you take the eigendecomposition of a symmetric positive semidefinite matrix, you get an orthogonal basis of eigenvectors. I.e. $Q^{-1} = Q^{T}$. In general this is not the case and you have to also invert the resulting $Q$ matrix." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Power Iteration method for finding the largest eigenvalue (and its vector)\n", "A very simple and often used method to find the eigenvector with the largest eigenvalue is called the power iteration method. If you start with a random vector $\\mathbf{b}_0$, we can apply the update step:\n", "$$\\mathbf{b}_{k+1} = \\frac{A\\mathbf{b}_k}{||A\\mathbf{b}_k||}$$\n", "\n", "Once we have the eigenvector, we can get its corresponding eigenvalue by using the Rayleigh quotient:\n", "$$\\lambda = \\frac{\\mathbf{b}_k^T A \\mathbf{b}_k}{\\mathbf{b}_k^T \\mathbf{b}_k}$$\n", "\n", "A modification of this - the Lanczos algorithm can be used to find several eigenvalues and eigenvectors of sparse matrices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us try do this in numpy:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "ExecuteTime": { "end_time": "2018-11-29T22:51:59.745689Z", "start_time": "2018-11-29T22:51:59.717424Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Result from power iteration method:\n", "Largest eigenvalue: 1.0000\n", "eigenvector: [0.7071 0.7071]\n", "\n", "Results from numpy linalg package:\n", "All eigenvectors:\n", "0.3000: [-0.9285 0.3714]\n", "1.0000: [-0.7071 -0.7071]\n" ] } ], "source": [ "def rayleigh_quotient(A, b):\n", " return (b.T @ A @ b)/(b.T @ b)\n", "\n", "def power_iteration(A, num_simulations):\n", " # Choose a random initial vector. We need to be careful as if the initial \n", " # vector is orthogonal to the eigenvector, it will not converge.\n", " b_k = np.random.rand(A.shape[1])\n", "\n", " for _ in range(num_simulations):\n", " # calculate the matrix-by-vector product Ab\n", " b_k1 = np.dot(A, b_k)\n", "\n", " # calculate the norm\n", " b_k1_norm = np.linalg.norm(b_k1)\n", "\n", " # re normalize the vector\n", " b_k = b_k1 / b_k1_norm\n", "\n", " λ_k = rayleigh_quotient(A, b_k)\n", " \n", " return b_k, λ_k\n", "\n", "# Declare an array\n", "A = np.array([[0.5, 0.5], [0.2, 0.8]])\n", "b, λ = power_iteration(A, 20)\n", "print('Result from power iteration method:')\n", "print('Largest eigenvalue: {:.4f}\\neigenvector: {}'.format(λ, b))\n", "\n", "print('\\nResults from numpy linalg package:')\n", "print('All eigenvectors:')\n", "w, v = np.linalg.eig(A)\n", "for i in range(len(w)):\n", " print('{:.4f}: {}'.format(w[i], v[:,i]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Rayleigh Quotient\n", "We introduced the Rayleigh quotient in the power method to find the eigenvalue of the corresponding eigenvector:\n", "$$\n", "R(M, x) = \\frac{x^T M x}{x^T x} = \\frac{\\sum_{i=1}^n \\lambda_i y_i^2}{\\sum_{i=1}^n y_i^2}\n", "$$\n", "\n", "where $(\\lambda _{i},v_{i})$ is the $i$th eigenpair after orthonormalization and $y_{i}=v_{i}^{*}x$ is the $i$th coordinate of x in the eigenbasis.\n", "\n", "There are a couple of other neat properties about the Rayleigh Quotient:\n", "\n", "- $R(M, x) \\in [\\lambda_{min}, \\lambda_{max}]$ \n", "- If $M$ is a symmetric matrix (such as a covariance matrix), then the vectors which maximise/minimise the Rayleigh quotient are their corresponding eigenvectors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Pseudo-Inverse\n", "I've used this a lot but have not really realized. Consider a general matrix $A\\in M(m,n;K)$ ($m$ rows, $n$ columns of values in the space $K$), the pseudo-inverse of A is defined as:\n", "$$\n", "A^+ \\in M(n,m; K)\n", "$$\n", "\n", "If $A$ has more rows than columns (i.e. a very tall matrix with linearly independent columns), then $A^*A \\in M(n,n;K)$ is invertible, and the psuedo-inverse is:\n", "$$\n", "A^+ = (A^*A)^{-1}A^*\n", "$$\n", "\n", "$$\n", "A^+A = I\n", "$$\n", "\n", "If $A$ has more columns than rows (i.e. a fat matrix with linearly independent rows), then $AA^* \\in M(m,m;K)$ is invertible and the pseudo-inverse is:\n", "\n", "$$\n", "A^+ = A^*(AA^*)^{-1}\n", "$$\n", "\n", "$$\n", "AA^+ = I\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# SVD\n", "This is a decomposition of an $m\\times n$ matrix \n", "$$\n", "M = U\\Sigma V^{*}\n", "$$\n", "\n", "where $U \\in \\mathbb{C}^{m\\times m}$, $V \\in \\mathbb{C}^{n \\times n}$ and $\\Sigma \\in \\mathbb{C}^{m\\times n}$. The diagonal entries $\\sigma_i$ of $\\Sigma$ are the singular values of $M$. The left singular values of $M$ ($U$) are a set of orthonormal eigenvectors of $MM^*$. The right singular vectors of $M$ ($V$) are a set of orthonormal eigenvectors of $M^*M$. The non-zero singular values of $M$ are the square roots of the non-zero eigenvalues of $M^*M$ and $MM^*$.\n", "\n", "A non-negative real number, $\\sigma$ is a singular value for $M$ iff there exist unit length vectors $u\\in \\mathbb{C}^{m\\times m}$ and $v \\in \\mathbb{C}^{n \\times n}$ where $Mv = \\sigma u$ and $M^*u = \\sigma v$.\n", "\n", "Also note that\n", "$$\n", "M^*M = V\\Sigma^*U^*U\\Sigma V^* = V\\Sigma^*\\Sigma V^*\n", "$$\n", "\n", "$$\n", "MM^*= U\\Sigma V* V\\Sigma^* U^* = U\\Sigma\\Sigma^* U^*\n", "$$\n", "\n", "So that the nonzero elements of $\\Sigma$ are the square roots of the non-zero eigenvalues of $M^*M$ or $MM*$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that if we have the SVD of a matrix, it is easier to calculate its psuedo-inverse, as \n", "\n", "$$\n", "M^+ = V\\Sigma^+U^*\n", "$$\n", "\n", "And $\\Sigma^+$ is simply formed by replacing every non-zero diagonal by its reciprocal and transposing the resulting matrix. This is used a lot in linear least squares problems." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Tikhonov Regularization\n", "Tikhonov regularization is the general term for the regularization we normally use, except that it allows for more complex terms. I.e. instead of minimizing\n", "$$\n", "||Ax-b||_2^2\n", "$$\n", "\n", "we minimize\n", "$$\n", "||Ax-b||^2_2 + ||\\Gamma x||^2_2\n", "$$\n", "\n", "for some suitably chosen Tikhonov matrix $\\Gamma$. In many cases $\\Gamma = \\alpha I$, in which case this collapses to $L2$ regularization, but it could be a more special matrix too.\n", "\n", "## Regularized Least Squares\n", "If we wanted to solve unregularized least squares, then for\n", "$$A = U\\Sigma V^T \\\\\n", "\\hat{x} = V\\Sigma^+U^Tb \\\\\n", "$$\n", "\n", "And with regularization:\n", "\n", "$$\n", "\\hat{x} = VDU^T b \\\\\n", "\\text{with} \\\\\n", "D_{ii} = \\frac{\\sigma_i}{\\sigma_i^2 + \\alpha^2}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Frobenius Norm and Singular Values\n", "There is a special case worth mentioning for the frobenius norm of a matrix:\n", "\n", "$$\n", "||A||_2^2 = \\sum_{i,j} A_{i,j}^2 = \\sum_i \\sigma_i^2\n", "$$\n", "\n", "Where $\\sigma_i$ is the $i$th singular value of $A$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Positive Semi-Definite Matrices and Diagonal Dominance\n", "If a matrix is diagonally dominant, i.e.\n", "$$\n", "h_{i,i} \\geq \\sum_{j \\neq i} |h_{ij}|\n", "$$\n", "\n", "and $h_{i,i} \\geq 0$, then we have a PSD matrix (this property is sufficient but not necessary). It is often why we sometimes add a small identity matrix to a covariance matrix to make sure we are satisfying its PSD conditions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# PCA\n", "As stated before, if a matrix is symmetric then the eigenvalues of it are orthogonal and $Q^{-1}=Q^{T}$. PCA finds the eigenvectors of the covariance matrix $C = (X-\\mu)^T (X-\\mu)$, and sorts them so that the eigenvalues are given in decreasing order. It is also equivalent to finding the singular values of $X-\\mu$." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "ExecuteTime": { "end_time": "2018-11-29T23:57:37.271695Z", "start_time": "2018-11-29T23:57:37.234278Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--------------------------------------------------\n", "PCA direct method\n", "--------------------------------------------------\n", "[[ 0.1739 -0.4942 0.8359 -0.1532 0.0565]\n", " [-0.5933 0.1449 0.1428 -0.0759 0.7752]\n", " [-0.2748 -0.4962 -0.0846 0.8189 -0.0219]\n", " [ 0.3267 0.6463 0.407 0.5461 0.1077]\n", " [ 0.6599 -0.2661 -0.3287 0.0429 0.6195]]\n", "\n", "--------------------------------------------------\n", "Eigendecomposition of sample covariance matrix\n", "--------------------------------------------------\n", "Eigenvalues in decreasing size:\n", "[ 9.1117e-01 5.3306e-01 1.0585e-01 2.5098e-03 -3.9139e-17]\n", "\n", "Corresponding eigenvectors\n", "[[-0.1739 0.4942 0.8359 0.1532 0.0565]\n", " [ 0.5933 -0.1449 0.1428 0.0759 0.7752]\n", " [ 0.2748 0.4962 -0.0846 -0.8189 -0.0219]\n", " [-0.3267 -0.6463 0.407 -0.5461 0.1077]\n", " [-0.6599 0.2661 -0.3287 -0.0429 0.6195]]\n", "\n", "--------------------------------------------------\n", "SVD of sample matrix\n", "--------------------------------------------------\n", "Square of the singular values\n", "[0.9112 0.5331 0.1058 0.0025 0. ]\n", "\n", "Right singular vectors:\n", "[[ 0.1739 -0.5933 -0.2748 0.3267 0.6599]\n", " [-0.4942 0.1449 -0.4962 0.6463 -0.2661]\n", " [ 0.8359 0.1428 -0.0846 0.407 -0.3287]\n", " [-0.1532 -0.0759 0.8189 0.5461 0.0429]\n", " [ 0.0565 0.7752 -0.0219 0.1077 0.6195]]\n" ] } ], "source": [ "from sklearn.decomposition import PCA\n", "\n", "print('-'*50)\n", "print('PCA direct method')\n", "print('-'*50)\n", "X = np.random.rand(5,5)\n", "pca = PCA()\n", "y = pca.fit(X)\n", "v1 = y.components_.T\n", "print(v1)\n", "\n", "print('\\n'+ '-'*50)\n", "print('Eigendecomposition of sample covariance matrix')\n", "print('-'*50)\n", "u = np.mean(X, axis=0)\n", "# Sort by eigenvalues\n", "w, v = np.linalg.eigh((X-u).T @ (X-u))\n", "k = np.argsort(w)[::-1]\n", "print('Eigenvalues in decreasing size:')\n", "print(w[k])\n", "print('\\nCorresponding eigenvectors')\n", "print(v[:,k])\n", "\n", "print('\\n'+ '-'*50)\n", "print('SVD of sample matrix')\n", "print('-'*50)\n", "u, s, v = np.linalg.svd(X-u)\n", "print('Square of the singular values')\n", "print(s**2)\n", "print('\\nRight singular vectors:')\n", "print(v)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Fourier Bases and PCA" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Olshausen and Field 97 - there is a comment that reads:\n", "\n", " Because the image statistics are roughly stationary, the eigenvectors of the covariance matrix will essentially be equivalent to the Fourier bases.\n", " \n", "This is something that is not entirely obvious to me, but empirically proveable. Let us take in a sample image and find the principle components of all $10\\times 10$ blocks of it." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "ExecuteTime": { "end_time": "2018-11-30T00:05:42.727906Z", "start_time": "2018-11-30T00:05:42.020302Z" } }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('');\n", " button.click(method_name, toolbar_event);\n", " button.mouseover(tooltip, toolbar_mouse_event);\n", " nav_element.append(button);\n", " }\n", "\n", " // Add the status bar.\n", " var status_bar = $('');\n", " nav_element.append(status_bar);\n", " this.message = status_bar[0];\n", "\n", " // Add the close button to the window.\n", " var buttongrp = $('
');\n", " var button = $('');\n", " button.click(function (evt) { fig.handle_close(fig, {}); } );\n", " button.mouseover('Stop Interaction', toolbar_mouse_event);\n", " buttongrp.append(button);\n", " var titlebar = this.root.find($('.ui-dialog-titlebar'));\n", " titlebar.prepend(buttongrp);\n", "}\n", "\n", "mpl.figure.prototype._root_extra_style = function(el){\n", " var fig = this\n", " el.on(\"remove\", function(){\n", "\tfig.close_ws(fig, {});\n", " });\n", "}\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(el){\n", " // this is important to make the div 'focusable\n", " el.attr('tabindex', 0)\n", " // reach out to IPython and tell the keyboard manager to turn it's self\n", " // off when our div gets focus\n", "\n", " // location in version 3\n", " if (IPython.notebook.keyboard_manager) {\n", " IPython.notebook.keyboard_manager.register_events(el);\n", " }\n", " else {\n", " // location in version 2\n", " IPython.keyboard_manager.register_events(el);\n", " }\n", "\n", "}\n", "\n", "mpl.figure.prototype._key_event_extra = function(event, name) {\n", " var manager = IPython.notebook.keyboard_manager;\n", " if (!manager)\n", " manager = IPython.keyboard_manager;\n", "\n", " // Check for shift+enter\n", " if (event.shiftKey && event.which == 13) {\n", " this.canvas_div.blur();\n", " event.shiftKey = false;\n", " // Send a \"J\" for go to next cell\n", " event.which = 74;\n", " event.keyCode = 74;\n", " manager.command_mode();\n", " manager.handle_keydown(event);\n", " }\n", "}\n", "\n", "mpl.figure.prototype.handle_save = function(fig, msg) {\n", " fig.ondownload(fig, null);\n", "}\n", "\n", "\n", "mpl.find_output_cell = function(html_output) {\n", " // Return the cell and output element which can be found *uniquely* in the notebook.\n", " // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n", " // IPython event is triggered only after the cells have been serialised, which for\n", " // our purposes (turning an active figure into a static one), is too late.\n", " var cells = IPython.notebook.get_cells();\n", " var ncells = cells.length;\n", " for (var i=0; i= 3 moved mimebundle to data attribute of output\n", " data = data.data;\n", " }\n", " if (data['text/html'] == html_output) {\n", " return [cell, data, j];\n", " }\n", " }\n", " }\n", " }\n", "}\n", "\n", "// Register the function which deals with the matplotlib target/channel.\n", "// The kernel may be null if the page has been refreshed.\n", "if (IPython.notebook.kernel != null) {\n", " IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n", "}\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "X = np.reshape(np.transpose(im, (2,0,1)), (-1,10*10))\n", "u = np.mean(X, axis=0)\n", "# Calculate all the covariance matrices\n", "C = 1/X.shape[0] * (X-u).T @ (X-u)\n", "\n", "# Get the eigenvectors of C\n", "w, v = np.linalg.eigh(C)\n", "k = np.argsort(w)[::-1]\n", "\n", "# Plot the result\n", "z = v.reshape((10,10,100))\n", "fig, ax = plt.subplots(1,2, figsize=(10,5))\n", "plotters.plot_activations(z[:,:,k], cols=10, ax=ax[0])\n", "ax[0].set_title('PCA Direct implementation');\n", "\n", "# Compare to Sklearn PCA\n", "pca = PCA(n_components=100)\n", "pca.fit(X)\n", "Z = pca.components_.reshape((100,10,10))\n", "plotters.plot_activations(Z.transpose((1,2,0)), cols=10, ax=ax[1]);\n", "ax[1].set_title('PCA Sklearn implementation');" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "ExecuteTime": { "end_time": "2018-11-30T00:05:50.617803Z", "start_time": "2018-11-30T00:05:50.375453Z" } }, "outputs": [ { "data": { "application/javascript": [ "/* Put everything inside the global mpl namespace */\n", "window.mpl = {};\n", "\n", "\n", "mpl.get_websocket_type = function() {\n", " if (typeof(WebSocket) !== 'undefined') {\n", " return WebSocket;\n", " } else if (typeof(MozWebSocket) !== 'undefined') {\n", " return MozWebSocket;\n", " } else {\n", " alert('Your browser does not have WebSocket support.' +\n", " 'Please try Chrome, Safari or Firefox ≥ 6. ' +\n", " 'Firefox 4 and 5 are also supported but you ' +\n", " 'have to enable WebSockets in about:config.');\n", " };\n", "}\n", "\n", "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n", " this.id = figure_id;\n", "\n", " this.ws = websocket;\n", "\n", " this.supports_binary = (this.ws.binaryType != undefined);\n", "\n", " if (!this.supports_binary) {\n", " var warnings = document.getElementById(\"mpl-warnings\");\n", " if (warnings) {\n", " warnings.style.display = 'block';\n", " warnings.textContent = (\n", " \"This browser does not support binary websocket messages. \" +\n", " \"Performance may be slow.\");\n", " }\n", " }\n", "\n", " this.imageObj = new Image();\n", "\n", " this.context = undefined;\n", " this.message = undefined;\n", " this.canvas = undefined;\n", " this.rubberband_canvas = undefined;\n", " this.rubberband_context = undefined;\n", " this.format_dropdown = undefined;\n", "\n", " this.image_mode = 'full';\n", "\n", " this.root = $('
');\n", " this._root_extra_style(this.root)\n", " this.root.attr('style', 'display: inline-block');\n", "\n", " $(parent_element).append(this.root);\n", "\n", " this._init_header(this);\n", " this._init_canvas(this);\n", " this._init_toolbar(this);\n", "\n", " var fig = this;\n", "\n", " this.waiting = false;\n", "\n", " this.ws.onopen = function () {\n", " fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n", " fig.send_message(\"send_image_mode\", {});\n", " if (mpl.ratio != 1) {\n", " fig.send_message(\"set_dpi_ratio\", {'dpi_ratio': mpl.ratio});\n", " }\n", " fig.send_message(\"refresh\", {});\n", " }\n", "\n", " this.imageObj.onload = function() {\n", " if (fig.image_mode == 'full') {\n", " // Full images could contain transparency (where diff images\n", " // almost always do), so we need to clear the canvas so that\n", " // there is no ghosting.\n", " fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n", " }\n", " fig.context.drawImage(fig.imageObj, 0, 0);\n", " };\n", "\n", " this.imageObj.onunload = function() {\n", " fig.ws.close();\n", " }\n", "\n", " this.ws.onmessage = this._make_on_message_function(this);\n", "\n", " this.ondownload = ondownload;\n", "}\n", "\n", "mpl.figure.prototype._init_header = function() {\n", " var titlebar = $(\n", " '
');\n", " var titletext = $(\n", " '
');\n", " titlebar.append(titletext)\n", " this.root.append(titlebar);\n", " this.header = titletext[0];\n", "}\n", "\n", "\n", "\n", "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "\n", "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n", "\n", "}\n", "\n", "mpl.figure.prototype._init_canvas = function() {\n", " var fig = this;\n", "\n", " var canvas_div = $('
');\n", "\n", " canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n", "\n", " function canvas_keyboard_event(event) {\n", " return fig.key_event(event, event['data']);\n", " }\n", "\n", " canvas_div.keydown('key_press', canvas_keyboard_event);\n", " canvas_div.keyup('key_release', canvas_keyboard_event);\n", " this.canvas_div = canvas_div\n", " this._canvas_extra_style(canvas_div)\n", " this.root.append(canvas_div);\n", "\n", " var canvas = $('');\n", " canvas.addClass('mpl-canvas');\n", " canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n", "\n", " this.canvas = canvas[0];\n", " this.context = canvas[0].getContext(\"2d\");\n", "\n", " var backingStore = this.context.backingStorePixelRatio ||\n", "\tthis.context.webkitBackingStorePixelRatio ||\n", "\tthis.context.mozBackingStorePixelRatio ||\n", "\tthis.context.msBackingStorePixelRatio ||\n", "\tthis.context.oBackingStorePixelRatio ||\n", "\tthis.context.backingStorePixelRatio || 1;\n", "\n", " mpl.ratio = (window.devicePixelRatio || 1) / backingStore;\n", "\n", " var rubberband = $('');\n", " rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n", "\n", " var pass_mouse_events = true;\n", "\n", " canvas_div.resizable({\n", " start: function(event, ui) {\n", " pass_mouse_events = false;\n", " },\n", " resize: function(event, ui) {\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " stop: function(event, ui) {\n", " pass_mouse_events = true;\n", " fig.request_resize(ui.size.width, ui.size.height);\n", " },\n", " });\n", "\n", " function mouse_event_fn(event) {\n", " if (pass_mouse_events)\n", " return fig.mouse_event(event, event['data']);\n", " }\n", "\n", " rubberband.mousedown('button_press', mouse_event_fn);\n", " rubberband.mouseup('button_release', mouse_event_fn);\n", " // Throttle sequential mouse events to 1 every 20ms.\n", " rubberband.mousemove('motion_notify', mouse_event_fn);\n", "\n", " rubberband.mouseenter('figure_enter', mouse_event_fn);\n", " rubberband.mouseleave('figure_leave', mouse_event_fn);\n", "\n", " canvas_div.on(\"wheel\", function (event) {\n", " event = event.originalEvent;\n", " event['data'] = 'scroll'\n", " if (event.deltaY < 0) {\n", " event.step = 1;\n", " } else {\n", " event.step = -1;\n", " }\n", " mouse_event_fn(event);\n", " });\n", "\n", " canvas_div.append(canvas);\n", " canvas_div.append(rubberband);\n", "\n", " this.rubberband = rubberband;\n", " this.rubberband_canvas = rubberband[0];\n", " this.rubberband_context = rubberband[0].getContext(\"2d\");\n", " this.rubberband_context.strokeStyle = \"#000000\";\n", "\n", " this._resize_canvas = function(width, height) {\n", " // Keep the size of the canvas, canvas container, and rubber band\n", " // canvas in synch.\n", " canvas_div.css('width', width)\n", " canvas_div.css('height', height)\n", "\n", " canvas.attr('width', width * mpl.ratio);\n", " canvas.attr('height', height * mpl.ratio);\n", " canvas.attr('style', 'width: ' + width + 'px; height: ' + height + 'px;');\n", "\n", " rubberband.attr('width', width);\n", " rubberband.attr('height', height);\n", " }\n", "\n", " // Set the figure to an initial 600x600px, this will subsequently be updated\n", " // upon first draw.\n", " this._resize_canvas(600, 600);\n", "\n", " // Disable right mouse context menu.\n", " $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n", " return false;\n", " });\n", "\n", " function set_focus () {\n", " canvas.focus();\n", " canvas_div.focus();\n", " }\n", "\n", " window.setTimeout(set_focus, 100);\n", "}\n", "\n", "mpl.figure.prototype._init_toolbar = function() {\n", " var fig = this;\n", "\n", " var nav_element = $('
')\n", " nav_element.attr('style', 'width: 100%');\n", " this.root.append(nav_element);\n", "\n", " // Define a callback function for later on.\n", " function toolbar_event(event) {\n", " return fig.toolbar_button_onclick(event['data']);\n", " }\n", " function toolbar_mouse_event(event) {\n", " return fig.toolbar_button_onmouseover(event['data']);\n", " }\n", "\n", " for(var toolbar_ind in mpl.toolbar_items) {\n", " var name = mpl.toolbar_items[toolbar_ind][0];\n", " var tooltip = mpl.toolbar_items[toolbar_ind][1];\n", " var image = mpl.toolbar_items[toolbar_ind][2];\n", " var method_name = mpl.toolbar_items[toolbar_ind][3];\n", "\n", " if (!name) {\n", " // put a spacer in here.\n", " continue;\n", " }\n", " var button = $('