1473 lines
653 KiB
Plaintext
1473 lines
653 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 45,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "skip"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"plt.rcParams.update({\n",
|
||
|
|
" \"text.usetex\": True,\n",
|
||
|
|
" \"font.size\": 14,\n",
|
||
|
|
" \"font.serif\": [\"Palatino\"],\n",
|
||
|
|
"})"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Matrix rank "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Matrix rank: dimension of the column space $\\equiv$ dimension of the row space"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"- We have not discussed yet, how to compute the matrix rank\n",
|
||
|
|
"- One (the most reliable but somewhat expensive) way to compute the matrix rank is to perform its SVD decomposition \n",
|
||
|
|
"- Before we turn to this, lets play with the function matrix_rank from numpy"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 56,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"100"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 56,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# A random matrix has full rank\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import numpy.linalg as la\n",
|
||
|
|
"n = 100\n",
|
||
|
|
"A = np.random.normal(size=(n,n))\n",
|
||
|
|
"la.matrix_rank(A)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Lets create a rank-1 matrix "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 59,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"1"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 59,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"u = np.random.normal(size=(n, 1))\n",
|
||
|
|
"A = u @ u.T\n",
|
||
|
|
"la.matrix_rank(A)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"$A_{ij} = u_j u_j$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 60,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"# Consider a small perturbation of A\n",
|
||
|
|
"delta_A = 1e-12*np.random.normal(size=(n,n))\n",
|
||
|
|
"Ap = A + delta_A"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 61,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[ 0.40534717, -0.07139809],\n",
|
||
|
|
" [-0.07139809, 0.0125761 ]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 61,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# upper left corner of A\n",
|
||
|
|
"A[:2,:2]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 62,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[ 0.40534717, -0.07139809],\n",
|
||
|
|
" [-0.07139809, 0.0125761 ]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 62,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# upper left corner of Ap\n",
|
||
|
|
"Ap[:2,:2]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 63,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"84"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 63,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"la.matrix_rank(Ap)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"- Matrix rank is very sensitive to small perturbations of a matrix\n",
|
||
|
|
"- Dimensionality of a space, spanned by a set of vectors can be poorly defined\n",
|
||
|
|
"- Today we resolve this issue to come with a better understanding of related matrix properties"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## What does it imply if we know that the matrix $A$ has low rank?"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Consider generic $n\\times n$ matrix $A = \\bigg[a_1\\bigg|a_2\\bigg|...\\bigg|a_n\\bigg]$:\n",
|
||
|
|
"- It takes $\\mathcal{O}(n^2)$ elements to store $A$\n",
|
||
|
|
"- It takes $\\mathcal{O}(n^2)$ FLOPs to compute $u=Av$ for a given $v$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Consider rank-one matrix $A$\n",
|
||
|
|
"- It takes $\\mathcal{O}(n)$ elements to store $A$. Indeed, all columns are collinear: \n",
|
||
|
|
"it takes $n$ elements to store one column $a_1$ and $n-1$ elements to store the proportionality coefficients $b_i$ such that $$a_i = b_i a_1$$\n",
|
||
|
|
"As a result, the matrix is compressed:\n",
|
||
|
|
"$$\n",
|
||
|
|
"A_{ij} = a_i b_j\\quad\\textrm{or}\\quad A=a^T b\n",
|
||
|
|
"$$\n",
|
||
|
|
"- It takes $\\mathcal{O}(n)$ flops to compute $u=Av$ for a given $v$. Why?"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Low-rank approximation"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 65,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"u = np.random.normal(size=(n, 1))\n",
|
||
|
|
"A = u @ u.T\n",
|
||
|
|
"delta_A = 1e-10*np.random.normal(size=(n,n))\n",
|
||
|
|
"Ap = A + delta_A"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 66,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"100"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 66,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# according to matrix_rank the rank of Ap is almost as high as it can be\n",
|
||
|
|
"la.matrix_rank(Ap)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 67,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"absolute error: 1.0158658489730158e-08\n",
|
||
|
|
"relative error: 1.2119290778015715e-10\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# what error do we make if we replace Ap by a rank-1 matrix A?\n",
|
||
|
|
"print(\"absolute error:\", np.linalg.norm(A - Ap))\n",
|
||
|
|
"print(\"relative error:\", np.linalg.norm(A - Ap)/np.linalg.norm(A))"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"- Thus, rank-1 matrix $A$ provides a very good approximation of a rank-100 matrix $Ap$\n",
|
||
|
|
"- It is a very synthetic example. Our goal today is to understand how this may work in a general case"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Geometric considerations "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"- Consider an arbitrary $m\\times n$ matrix $A$. It maps arbitrary $n$-dimensional vector $v$ to $m$-dimensional vector $u$: $u=Av$.\n",
|
||
|
|
"- Lets inspect the geometric meaning of this fact\n",
|
||
|
|
""
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Singular value decomposition\n",
|
||
|
|
"\n",
|
||
|
|
"Any $m\\times n$ matrix $A$ can be written as a product of three matrices: \n",
|
||
|
|
"\n",
|
||
|
|
"$$ A = U \\Sigma V^*, $$\n",
|
||
|
|
"\n",
|
||
|
|
"where, for $K = \\min(m, n)$:\n",
|
||
|
|
"- $U$ is an $m \\times K$ matrix with orthonormal columns, \n",
|
||
|
|
"- $V$ is an $n \\times K$ matrix with orthonormal columns,\n",
|
||
|
|
"- $\\Sigma$ is a diagonal matrix with non-negative elements on the diagonal \n",
|
||
|
|
"\n",
|
||
|
|
"Reduced SVD ($m\\geq n$)\n",
|
||
|
|
""
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Sometimes its convenient to append the missing columns, so that both matrices ($U$ and $V$) become unitary."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Full SVD ($m\\geq n$)\n",
|
||
|
|
""
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"- We can write down $A$ element-wise: $$a_{ij} = \\displaystyle{\\sum_{\\alpha=1}^r} \\sigma_\\alpha u_{i\\alpha} v_{j\\alpha}^{*}$$\n",
|
||
|
|
"- One usually enumerates diagonal elements of $\\Sigma$ in a descending order $\\sigma_1 \\geq \\ldots, \\geq \\sigma_K$.\n",
|
||
|
|
"- If $\\text{rank}(A) = r$, then: $\\sigma_{r+1} = ... = \\sigma_K = 0.$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Columns of matrices $U$ and $V$ are known as left and right singular vectors.\n",
|
||
|
|
"They form orthonormal sets in the space of image ($u_i$) and pre-image ($v_i$) of the matrix $A$. Their geometrical meaning:\n",
|
||
|
|
"\n",
|
||
|
|
"$$\\text{ker}(A) = \\mathrm{span}\\{v_{r+1},\\dots,v_K\\}$$\n",
|
||
|
|
"\n",
|
||
|
|
"$$\\text{im}(A) = \\mathrm{span}\\{u_{1},\\dots,u_r\\}$$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Eckart-Young theorem\n",
|
||
|
|
"\n",
|
||
|
|
"The best low-rank approximation can be computed by SVD.\n",
|
||
|
|
"\n",
|
||
|
|
"Let $r < \\text{rank}(A)$, $B$ is a matrix of rank $r$ (to be optimized). Then\n",
|
||
|
|
"\n",
|
||
|
|
"$$ \\min_{\\text{rank}(B)=r} \\|A - B\\|_2 = \\sigma_{r+1}. $$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Low-rank approximation\n",
|
||
|
|
"\n",
|
||
|
|
"The best rank-$r$ approximation to $A$ is obtained by setting $\\sigma_{r+1}= 0, \\ldots, \\sigma_K = 0$:\n",
|
||
|
|
"\n",
|
||
|
|
"$$A_r = U_r \\Sigma_r V_r^*.$$\n",
|
||
|
|
"\n",
|
||
|
|
"The error \n",
|
||
|
|
"\n",
|
||
|
|
"$$ \\min_{\\text{rank}(B)=r} \\Vert A - B \\Vert_2 = \\Vert A - A_r \\Vert_2 = \\sigma_{r+1}$$\n",
|
||
|
|
"\n",
|
||
|
|
"Thus the faster the singular values decay, the better is accuracy of 'compression' of the matrix into low-rank form."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"The same holds for Frobenius norm:\n",
|
||
|
|
"$$\\min_{\\text{rank}(B)=r}\\Vert A - B \\Vert_F = \\Vert A - A_r \\Vert_F = \\sqrt{\\sigma_{r+1}^2 + \\dots + \\sigma_{K}^2}$$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Computing SVD\n",
|
||
|
|
"\n",
|
||
|
|
"- Algorithms for the computation of the SVD are tricky. Some approaches will be discussed later.\n",
|
||
|
|
"\n",
|
||
|
|
"- We are ready to use NumPy and SciPy implementations\n",
|
||
|
|
"\n",
|
||
|
|
"- Broadly speaking, one may think of (i) exact computation of complete SVD, (ii) exact computation of 'sparse' SVD, where we need only first $k \\ll n$ singular values and vectors and (iii) randomized SVD (approximate evaluation of the leading subset of singular values and vectors)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 68,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Rank of the matrix: 1\n",
|
||
|
|
"Rank of the matrix: 100\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"# Recall: Computing matrix rank\n",
|
||
|
|
"import numpy as np\n",
|
||
|
|
"print('Rank of the matrix:', np.linalg.matrix_rank(A))\n",
|
||
|
|
"print('Rank of the matrix:', np.linalg.matrix_rank(Ap))"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 69,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[ 0.62881492, -0.38142909],\n",
|
||
|
|
" [-0.38142909, 0.2313688 ]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 69,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"A[:2,:2]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 70,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[ 0.62881492, -0.38142909],\n",
|
||
|
|
" [-0.38142909, 0.2313688 ]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 70,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"Ap[:2,:2]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 74,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"u, s, vh = np.linalg.svd(A)\n",
|
||
|
|
"up, sp, vhp = np.linalg.svd(Ap)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"la.matrix_rank()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 75,
|
||
|
|
"metadata": {
|
||
|
|
"scrolled": true,
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([8.38222193e+01, 4.53158462e-14, 4.00097651e-14, 3.10833486e-14,\n",
|
||
|
|
" 2.59425086e-14, 2.13154110e-14, 2.07592143e-14, 1.90293475e-14,\n",
|
||
|
|
" 1.77971365e-14, 1.46656957e-14])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 75,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"s[:10]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 73,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([8.38222193e+01, 1.99954677e-09, 1.89619187e-09, 1.86088354e-09,\n",
|
||
|
|
" 1.84812585e-09, 1.81777344e-09, 1.75761017e-09, 1.75596734e-09,\n",
|
||
|
|
" 1.68481081e-09, 1.66847435e-09])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 73,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"sp[:10]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"SVD allows to properly define a numerical rank for a given tolerance $\\epsilon$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Singular values of a random Gaussian matrix\n",
|
||
|
|
"\n",
|
||
|
|
"What is the singular value decay of a random matrix?"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 23,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAENCAYAAADEw0mJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABBTElEQVR4nO3dd1iUV/YH8O+ZodcBBFTAMmDFiti7Eo1JTLVk07OxpWfXaOImv822bNaYnk1Bs+ldY2JM3GQBexfQaOyCvSBIEWwo3N8f80IQ6c7MO+X7eZ55hHdeZs6bDHcOd849V5RSICIiIiKiq2fQOwAiIiIiIlfB5JqIiIiIyEqYXBMRERERWQmTayIiIiIiK2FyTURERERkJUyuiYiIiIisxEPvAKylWbNmqk2bNnqHQUTUJBkZGXlKqXC947AnjttE5KzqGrNdJrlu06YN0tPT9Q6DiKhJROSg3jFYg4gkAchWSmXXdy7HbSJyVnWN2SwLISKiqyYiJhGZDSAdwFNakk1E5HaYXBMRkTVMAJCilCoEkAzgKX3DISLSh92TaxFJEpEMETHXcY5JRGZrt5S6ziUiIuurbayuY3zuBSBf+zobAMdtInJLdqu5FhETgAwAmQAS6jk9DcBIpVSh9tFiCoBY20ZIREQNGKtrG59DK07Q7rNDtEREjsduybU22PbS/lW1nacN1vnaR4tQSqWKiFlEEpRSmfaKl4jIHdU1Vtc1PsMyax2qnWeCZfaaiMjt2LUspGJArkcCrhyU+REjEZGd1DFW1zU+Z+C3mW4zLDPaRERuxxEXNIYBKKx2rBBVPnK0lvJyhcP5Z639sERErqrW8VkpNRdArIhMAZCklHrR2k9+/mIZck6ft/bDEhFZlSP2uTbVcjy/luNN9tQ3W7Fyby7WPj0SRgPrA4mI6mGq5Xg+ACilptryydN2nsSjX2Tilp7RmJAYjZ6tQuDl4YhzRETkzhwxuc4C0LuG44XVD2gzJFMAoFWrVo1+osHtwzE/4wi2HS1CjxhTo3+eiMjNNHh8rs3VjNvdooNx74A2+GLjIXyTeQS+nkYktglBr9Yh6NQiCP1jwxDk49moxyQisjZHTK6zAUysdsyslEqtfqL2MeRcAEhMTKx1kWRt4lsGAQD2nSxhck1EVL8Gj8+1uZpxOybUD8+Njcf0UR2wak8uNuzPx9qsPLyWurfynCiTLzo2D0T3GBM6twhCh+aBiDL5wsBPJ4nIThwmuRaRmQAWKKUWiMg8ETErpbK1VegNHrgbo1WoHwK8PbA++xTG9Yq2xVMQEbkMa43PIjIWwNi4uLgmxRHg7YExXVtgTNcWACy12JsPFSLzUAF2nyjGjuOnkbbrZOX5fl5GtIsMRPuIAHRoHoj2kYHo0DwQEYHeYMtAIrI2uybXWhun8dq3ySKSrJRaoH0/FZZZkWwAIwHMFpFNAMKUUuOvfLSr52k04LquzfHj1uP4+01d4OtltMXTEBE5lXrG6qsen5VSiwEsTkxMnGyNeH08jegfG4b+sWGVx4rPX8SenBLsySnG7hPF2JNTjGW7czE/40jlOcG+nugQGYh2kVWS7shAhPh7WSMsInJTolSjqykcUmJiokpPT2/0z63LOoXfzVuP12/vgZt6RNkgMiKi+olIhlIqUe847Kmp4/bVOFVy4bekO6cYe7Xk+/T5S5XnhAd6IzbcH9EhfmgV6oe4iADEhgegTTM/eHtwEoaI6h6zHaYsRC9924YiyuSLhZlHmVwTEdnB1ZaFXI2wAG/0D/C+bJZbKYWc0xewO6cYe7RZ7uy8M1i1Nxc5py9Unmc0CFqF+iE23B/m8AC0beaPts38YQ73R3gAS0yIyMLtk2uDQXBLzyj8e9k+rNqbi8HtwvUOiYjIpVm7LORqiQiaB/ugebAPhra//D3gXGkZsnJLkJVbgn0nf/t35d48lF4qrzwvwNujMtmuSLgrvg5kBxMit+L2yTUATB5ixrebj+Kvi3fgx8cG8WM/IiICAPh6GdElKhhdooIvO15ernCs6Bz2551Bdu4Zy795Z7D5cAEWbz2GqhWX4YHeloS7Mum2zHq3CvVjn24iF8TkGpZFLX+/OR6//zAdb6btw5OjO+gdEhEROTCDQRAd4ofoEL8rPvE8f7EMh/LPVibd+/NKkJ17Bik7cnDqTGnleUaDICbEV5vptiTcrcP80DrUHy1NPvAwMvEmckZMrjUjOkbi2vjm+GzDQTw2sh1nE4iIbETPmmt78PE0on2kpftIdUVnLyI7r0RLui2z3tl5Z7Au+xTOX/ytzMTDIIiulnibta8jg1jfTeTImFxXMaF3NH7afgIvp+zGrDGd9A6HiMglOVrNtT0F+3miZ6sQ9GwVctnx8nKFE6fP4+CpsziUfwYHT53FwVNna0y8fT2NaKMl21VrvM3hAQj2ZX03kd6YXFcxomMkbkuIRvKKbHRqHoSbe7J7CBER2Z7BIGhp8kVLk+9lnUyA3xJvy0x3CfbnncX+vBJsP1aEn7afQFn5bwXezQK8YG4WULmgMjY8ALERAYgJ8WWZCZGdMLmu5vlbumDNvjzM+Xk3bujWgoMRERHpqmriPTCu2WX3lV4q1+q7S5BdmXxfWd/taRS0DvOvbCPYPjIA7SIsG+hwET+RdTG5rsbH04i/3BiPaZ9mIHllNh4e7po1gUREenH1mmt78vIwIC4iAHERAVfcV3T2IrLySpB10pJ4Z520tBFM23kSl7TZbqNBEGXyReswy4Y5bZv5IzYiAO0iAtAy2BcGA2u7iRqLyXUNru3SHNd3a4HXUvdgcLtm6BZt0jskIiKX4c411/YU7OeJhFYhSKhW332xrBwH8s5gd04xdh0vxoFTZ3Ao/yx+2HocRecuVp7n52VEu4gAdGwehI4tAtEuIhCxEf5oHuTDBZVEdWByXYt/3NQFmw8W4M73NuCHRwehdZi/3iERERFdNU+jAe0iA9EuMhA3dLv8vvwzpdinzXDvPWnZrTJlZw6+Sj9ceY6/lxHm8ADEhvtbHiciAO0iA9Eq1A9GznQTMbmuTYi/F+bdm4jr31iNRz7fjO8fGci/1ImIyKWF+nuhT9tQ9GkbWnlMKYXc4gvYl2vp112xS+XG/fn4bsuxyvO8PAyIDvFFmzB/dG4RhE4tgtCheQBah/nDk+uXyI0wua5DfMtgPDQsFm8vz8J/Vu/HpMFmvUMiIiKyKxFBRJAPIoJ8MCD28gWVJRcuYd/JEuzJKca+kyU4UmDZPGfFntzKLiaeRoG5WQDiIgPQPiIQXaOD0D3ahLAAbz0uh8jmmFzX44mk9kg/WIB/LtmJgXHN0KlFkN4hERE5NS5odB0B3h7oEWNCjxjTZcfPXyyrTLr35JRgb04xth0pwpJtxyu3hm8W4IUOzQPRITIIHZsHokNzy8Y7vl7sXkLOTZRS9Z/lBBITE1V6erpNHvtUyQUMmr0MY7o2x8vju7M8hIisTkQylFKJesdhT7Yct8kxnS29hF8OF2H7sSLsOlGM3Scsdd0XLlk2yREBWof6oX1koJZwB6FD80C0CfNja1xyKHWN2Zy5boCwAG9MSIzGR+sOItjXE3++oTMTbCIiokby8/JA/9iwyzbKKStXOJR/FrtPnK5MuHefKEbqzhxU7I/j5WFA+8gAdI0yIb5lEDq3tMx2+3kxjSHHw1dlAz03Nh4igg/WHEDJ+UuYM7673iERERE5PaNBKrdxv7ZLi8rj5y+WYW9OCXbnFGP3idPYcfw0ftx6DF9sPAQAMAhgDg9Al5ZBiG8ZjPioIMS3CEawH7eAJ30xuW4gg0Hw5xs641xpGb5KP4yOLYLwwKC2eodFRETkknw8jegaHYyu0cGVx5RSOFp4DtuPncb2Y6ex41gR1mdf3rUkJtQXPWJC0D06GN2iLTPd/t5Md8h++GprBINB8KfrO+FY0Tn8/YcdMIf7Y3iHCL3DIiIicgsigugQP0SH+GF0fPPK43klF7SEuwi/Hi1C+oF8LP7lmPYzQFx4ALpGBaNdZCC6xwQjoVUIfDy5cJJsgwsam+DCpTJc9/oq5Jy+gIUPDUD7yEC7PC8RuS53WtB
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 864x288 with 2 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"n = 1000\n",
|
||
|
|
"a = np.random.randn(n, n)\n",
|
||
|
|
"s = np.linalg.svd(a, compute_uv=False)\n",
|
||
|
|
"fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
|
||
|
|
"ax[0].plot(s/s[0])\n",
|
||
|
|
"ax[1].plot(s/s[0])\n",
|
||
|
|
"ax[1].set_yscale('log')\n",
|
||
|
|
"ax[0].set_ylabel(r\"$\\sigma_i / \\sigma_0$\", fontsize=16)\n",
|
||
|
|
"ax[0].set_xlabel(r\"$i$\", fontsize=16)\n",
|
||
|
|
"ax[1].set_xlabel(r\"$i$\", fontsize=16);"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Applications of SVD"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Dense matrix compression\n",
|
||
|
|
"\n",
|
||
|
|
"Dense matrices require $N^2$ elements to be stored. A rank-$r$ approximation can reduces this number of $\\mathcal{O}(Nr)$. Let us consider image compression as an example."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 76,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"from skimage import data\n",
|
||
|
|
"from skimage.color import rgb2gray\n",
|
||
|
|
"from numpy.linalg import svd\n",
|
||
|
|
"from skimage import img_as_float"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 77,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"cat = rgb2gray(img_as_float(data.chelsea()))"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 80,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"(300, 451)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 80,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"cat.shape"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 78,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[0.48523098, 0.48523098, 0.47738784, 0.47738784, 0.47738784],\n",
|
||
|
|
" [0.49699569, 0.49307412, 0.48523098, 0.48130941, 0.48130941],\n",
|
||
|
|
" [0.50849255, 0.50457098, 0.49475569, 0.49083412, 0.49307412],\n",
|
||
|
|
" [0.52054 , 0.51269686, 0.50485373, 0.50485373, 0.50064941],\n",
|
||
|
|
" [0.52894863, 0.52110549, 0.51661843, 0.51269686, 0.50849255]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 78,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"cat[:5, :5]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 81,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.image.AxesImage at 0x2826347c400>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 81,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAAD/CAYAAAAHSua4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz9x65t2baeB37dDTPNcttF7Ig4lteIl5lKXIkJZCFLpB5BUDFr5CMIqmRVUBazJuoJEuAbkCpnApIIELTXnDguzHbLTTNcd1lofYw119pxzmXy5oHOAaMDa++55pxrmjH6+Hvrf/vb31TOme/H9+P78f34fvxhDv2/9wf4fnw/vh/fj+/Hf/z4HsS/H9+P78f34w94fA/i34/vx/fj+/EHPL4H8e/H9+P78f34Ax7fg/j34/vx/fh+/AGP70H8+/H9+H58P/6Ah/1dvbBS6gL478qvfw7845zzl7+r9/t+fD++H9+P/xTH7wzEgf8Z+Ac55zul1D8E/hnw09/h+30/vh/fj+/Hf3Ljd0KnFNC+yTnfAeSc/znwE6XUn/8u3u/78f34fnw//lMdvytO/M+Bp9TJl8BPfkfv9/34fnw/vh//SY7fFZ3yDLh7ct8dcPX0iUqpfwT8IwDduP+i/vw53+UEoJTcqZ7cln9Akeeby+OnI6PI+eH/lBU5K5QCoxNaZXKGmPTyWM4K8vweGVL5ndP7gQQqlbuXDwEo+RsVTx7XcjsboE4Yk3A6YlVaPmvImozCqohRmZQVPhl5SyXfNWVFygqtMlplYtLErNEqYXUiZcUUDTlqlM5YE+W7KnkMwEfz6HhplckolHzJ8hnK8xWEqElegwJlEzloOQ6ah+8ZynfPkCpQLpGjkudEsEeFSqBSJrSKZOUwlSOM9pBceU2T5c4kf0Mux81k1CSfKdsMCoyNxKhRo0YFSGs5tuloUQnSSn7XKuN0IpZjl7LCR0OaTDm3YAYwPqMSxEoR1nn5jrYOOB1JWRGzJkaJg/Kk0R7MBGaIZF0+n1KkSpEccnwiZAvJynFSQe5b5qkpP+V7keW4zp8tG0DLY9oklJJ5H5MmZ5kfwHINKSW3tS7nNKvlfGcgJ8XTodTH15Aq10Auc0epvMxF9Rt+f7iOZUJYndBkEg/vuRzHPF93j95VXu/ks8zvn8vtef5neHzdLs8/ebn5ej75WEpnmYBZkaNCBYWeQPuMjkD6DbYkjw728s/J9316399u7LtvP+ScX3zXY78rEL/4DfffPL0j5/xPgH8CsPnjT/Lf/X/+3wqQQiqACgK0SmWMzjgTUQjoGJ2WyTE/7+lI+fFE7XzFcXKM3qFUZl1PtM6jVabzjvtjy3CsyKMR4HZJzkssrxMUptOoqEhVRnmFGcrF6WRCxHUEl1G9prrTuL0ia5gu5HPaTtG/DrQvOj692PEn5+/Y2gGnIn103Pg1mszajtz5loOvsTpR6YjTkTFaroc1zkRWdqIxng/DhtuhxajMpho5qwZ+vb9g1zVs2xGjE5dNj4+GD92K+/2K2Fm2z48YlaldQKnMYahJSRGC4fn5AR8No7eMk2Xc1ajRUL/sGN632IMhPPOowWDvNdVOLQtZ/1kkWzk+KijMoHj5vyZsn0iV4sPfs/htprpXTGeZ+k6RDIwvElmBfj4SDxZzEIBNTca96JmOFebakT8d0F83hKvAi9d3vP/mAntrCZvIFz99zzfX59i/bpkuE5c/ukWXBfsn59ckFO/7DWOwfPvhHN7VpE3E3liaD4qrvwhUd57j65rbP9Vkk5meJS5/cEvjAikrOSZeLqHu3ZrmW8uzfxvZ/OJI1gp/URNrzf5zw3ihMBO4Q0YFmM4VsQHby8Ll9pnqmBjPNN0rxXSRiW1CJUV1o6n2UN9mpjNF90kmrDPuk451O2JNYt81TIOlbj3rZuI4VOX6UsSocC7K9ZTUAuopKWIwssgqBMRMxtWBqgpYnbBloQBZxEdvSUlL8GEi1qSHa7Bcj1plWfyRQCMmjTORT1f3DNExRctZ1VPrSK0DY7LsQ83d2DIExxgNIRpimq/9vFz/lOtegVyzZI6+4jhWxALgqfxdjHq5nZOWx+JDAKZtoqoCWmeGwZHeN6y+1Wx/lWiuA2ZMqJAkyHi6KmZAqweQL4GLypmsFOp09UhZnnsyVCwLqVrWELk/Qzbyusv9WvHP/9//919+BGxl/K5A/GfA3/+O++9+2x9plWlsICRNzIqYHlb++TjM0ZOZo4knAP0QYT88Pq/8uqzqtQ3LY/1YcaTC6ERjA7WJWBvRNhMnUFFBMMvimktEQ5IVe35AZSCCVgJEKipyncibiA8KPSp0lEg8bCNgcLeGoa3ZNQ27VYMmU2vPxo7UOtClCkPCqYeLBGBtJqxKdLZiZafl+13VHZWO3AwrcpYo+sdnN7x1W2LS7MeKD2nNeT0QYrl4J00IBuMCMWmMToRg8JNF68Sb9+cS3QQNU2HfgmI8VnJsygWRTSbVmfECzKhQGfSgqG4NqZbZWN3Lc4dLg98osoXqTmF7CC3ECobPPGqSaJ83NdqA6Uo0rOTCtO8d4YWHybD5VjH+eOA4VKgqEi7BrD2//vIF7nLA/3hgvRl5vup4f1xjNFxVHb/qLvnm9pyxd+j3FXGdUHUk1QYdFNpnVMpUu8jZzxX9C8X4aeLY1/hoyFnRHesFAKkTqcokowibiu5VxbSRHUhYKewA0xbGK1h9rah2mTjJXEkO7JhZ/6rDXdaEpiI2imwkgs9GgL79EKl3Cu0105nm6Br0i4y1ccEYrTOVDdyHhhiM7EaywnvBnxm8U1JonYkzBumy0+QBr6xJVDYsQDzvXAFC0ISgsTYtwK3Kz4xXSmVC1ISkWdcTX+cLbo4rrEncVS0rN7FxI1YlrEpc1D3BjQzR0QfHGCxTNPhgGIPsGq1OrCrPRdMDcPQVY7CEZYGav4N8vyVY1pGUFEpp+YxGAFwhi1McDe6oqO4zrhPgzlqB1cvmew6sVc6QKGS0EtCeQblcIjmdHEg5GKiUBZSVAvsA0stOVMkOVbbM/+Fx/O8KxL8E/psn9/2kJDh/4xAQ9+SsmJLBR7OsrmoBc4UG8skWKyT9UTQ+g31Cdp8Jec2ZLolZEZImeENKCl8ZKhOXBULrRATZtmfITkBERUW2WX6fHqJOKIBPBgMqKLLX6LUnbRSxV+RRSaTeJLwGd2vgaDluKna+YW1HSI42ezZmxGnZZ/epYh9qUlZ00aFVYmMnKiOLkVAeMnue1UcqE+hCxd3UsnEjr9o9h1CzHyvuDq0sVi7Qu0hqDOPgaCpPko8u30VljE2kr1uZUG3C3hviJqFiCUBclu+TAJeIa4RaujPkcgxcB7lXhBZ0gMOnsiDGWiat9mDGjDsqxotCIShw9xp7VIxXiWzBTAp/6dG/amWuVBH9TcN4CVUV2H9YCwglRbquUcDVWUfrPFYnfnl9SYqa18/ueTdu+Mu3L5jerTAHoUHSi0D2Gjso3CEv59V1gdBphmDQbaCuwkJdaJ0JhYYxbcCfG8YLTTYVd3+kcYeHyNt1meEZTJcRPVqaa7V8/2yhuo+Y+55KK1bvDWFtBETIhFXGbwQomvceM1qmrSEbS+ca7NpL0FKCxJg0YbLkzpJMFlqqCTSNx+iED0aCGRc4lGhcKcAmlM6kpJkmuwB4zoqY5JrRZacbgi67NY3P86wBUyL3OZBKSS9/c3NccTg2KJ3oK0dlayq7orGBxvqyoww0JrBxI10QgD74atkNWpOobaAPjv1YCwAXAJ8j75wECeb3PUVKpRPGZCoXqGxkCgbvDepgqXZqOV8C4EIpLhFkibrV6UWfC3DbQhMW4M5WcEOdROq50JIzLzqzMI/A+mSByGoOEH87nP9OQDzn/E+VUv+TUuonOecviyrltwI4CIhv3ChAGt3C7Z4C7+lYqJYTQNdPbqcymXSZ3L13DN4So2YaHfFoUU3EN4Zo9TJRXRUIlSFrs/CxKigB7/n9T/lcDWaEqEFHwRM1KZI18voXGrPXZJtxqwm9yYy6QXeGsXdc9ysuqg5nR8ZkMSoRs8aoxNqMjNbis2bA8WH
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"plt.imshow(cat)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 83,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"U, s, Vc = svd(cat, full_matrices=False)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 88,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"def compress_and_show(k):\n",
|
||
|
|
" os = cat.shape\n",
|
||
|
|
" cat_compressed = U[:,:k] @ np.diag(s[:k]) @ Vc[:k,:]\n",
|
||
|
|
" compression_ratio = 100.0* (k*(os[0] + os[1])+k)/(os[0]*os[1])\n",
|
||
|
|
" plt.title(\"compression ratio={:.2f}\".format(compression_ratio)+\"%\")\n",
|
||
|
|
" plt.imshow(cat_compressed)\n",
|
||
|
|
" plt.axis('off')\n",
|
||
|
|
" plt.tight_layout()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 90,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEQCAYAAABLMTQcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz92ZItW3aeB36z82Y10e69T5uJzATAriiKBlKXdQe8AVj1BALeQLAyK7O6lFFvAPAJKOhSV0XoAUolZFlRIiUykQ3ytLuJbnXezK4uxnRfK+Lsc04iQRCJrD3NwiJiNb58uc855hj/+Mc/VM6Zd+PdeDfejXfj12vov+0TeDfejXfj3Xg3/tOPd8b93Xg33o1349dwvDPu78a78W68G7+G451xfzfejXfj3fg1HO+M+7vxbrwb78av4Xhn3N+Nv/WhlPpzpdQf/22fxy8z/i6f+7vx6z3eGfd341dh/LfAn/5tn8QvOf6znbtS6veVUv/mb+r5v+p7lFIXSqk/VkrdlZ9/+Vc59rvxNzvs3/YJvBvvRs75f/jbPodvG0qp3wX+OOf8m6eP/+c492I0/wC4/Zt4/pc85gXw58AfA9/POd+Xx96NX5HxznN/N96NX/GRc/6jnPMl8FbP+K/7/C/5nn8J/FnO+b/LOd+X99z/osd/N/7mxzvj/is4lFL/jVLqx0qpXDDd3z157l+W5+6eYr1KqT8t751C5R8rpX63/EzH+9Mn7/njcszT8Pr3nxzzD8rPj5+cy+nn/MGTc7x7ev7f8Pi/eRrS/zLf8xe8tm/9PgV++PNybj9+eg2AfwP8oDyfJy/1r3ruv0bjD4B/U67nXbl2P/jbPql34zjeGfdfsVGMwf8V+BfAJfBHwEV57k+B3wF+D/g+cPUED71APKo/Lc//sPz9h8A/Kz+/f2qIgR8gC3V6z38P/OnJQr0o7/8j4I9yzn92ci4/KO/5PeBfKqV+pxjL30dCdVXed/t1j3/NNfhlvucvakTf+n2AK+C/Luf2h+Ua/A5AzvlfIPfjJzlnVX7uf8lz/zs/TubG/w25jt9H7uXf1bzJr+V4h7n/Co3iDf4B8Js555+Uhydj+juIcbw8MSz/onhNv3tipH54YoD/uLznj8t7fqiU+iHwCDcG/uTk/X9YDPFkAKEY8elzy+I+PZd7pdQfIZvS/xu4OgnVp3P5/bc9/pZr8Mt+z7+KAX30fcr5/MnJ33+mlPoJ8LvIxvELjV/k3P8Knvyfn57Tr9iYjPsfTfNUKfUvgDul1O/knH/ha/Zu/M2Nd8b9V2v8LnB/YthPxz9HPMf7J4//L4iX+Gcn/0/j9i2P/YQSCXzD+DOOCxgEWz393N8pv3+qlHp6Lv8t4qnncpw/LknHP/uax5+Ov873/EXH0+8DQIlofg/57r8MxPCt555z/sNf4ri/auMr86okVO+R6/bOuP8KjHewzN+dcfELvu7+6QP/CRJdb9tsfphzvnzy83s55/vCKPnDci5/qpT6b77u8bcc9+IXPKf7v/K3OI6vfB+l1J8j0Msf55z/Gb+cgbr4a5zT36XxtvnwbvyKjXee+6/W+CFwoZT6wVu89z9DcO2LJ8b6n/OL482/6PjdbznmD4Hfecu5zKNACn9SvOE/BP67b3r8ZPzn/J7ADDP9TsHb/zrjW8/91wGWKV76T4D/C/AnMF/DC9557b8y451x/xUaOeefKKX+BPFq/wUS/v4u8F/lnP9IKfVnwP908ty/QmCAvy7X+g+UUv8a8cj+FRJaf61heXKef1j+//3yvp88ef/vAafPP3r8Lcf+4d/g9/y6cQsCy+Sc/6Sc6+8A//rkNT9B2DIXHOGXR+f/i5z7ryIsc7K5/VWu779ENrL/Bbk2fwz8D18DKb4bfwvjHSzzKzbK4v8zJEF4h3i3/7o8N2HOfw78FLgtEMJfd/wZwnz4KWKA/9m3QTnlPH8I/LlSajrPP0MW+u+VY90h3tx//Q2Pv+3Yf1Pf8+u+yz0SQfxx+S7T59+fvOaHyPf9KcdE83+Wcy+0zYwY0ImO+eP/VM8jDsS/+qt8Zokq/gj4n8r3/ElhFb0bvyJDvevE9P/fo9D0fphz/lqD9W68G+/G373xznN/N96Nd+Pd+DUc74z7u/FuvBvvxq/heGfc34134914N34NxzvM/d14N96Nd+PXcLzz3N+Nd+PdeDd+Dcc74/5uvBvvxrvxazi+sYjpn/6P//cckyYDp8V7SmUUYHRCqYzRGaMTAE4ntBKox5mIUfK4PXlcjggJRcry45NhCJaYFd3oGIPF6MTlomPhRgAOvuKhaxi8xY+WGDQkBTqjdCaNBgYNGVRUkAEDZNC9wgzyHWKdyQayhuwSKijsTuN2Cj1CNpAqUBGShf57I+urPdfLA+8ttrTGE7MiZc3GN4SkWbmBy+pAypoH3/AwtCiV5+/sk+HgHU4nVtXAEC33XcsYDM9Wey7rA3fDgk9eX+IPDtNEri52NDbQWk9tA11wGJX49P4CpTLORAAWlScDMWkqE0lZoVVmVQ3sxhqjE/eHlvtXa9TB4N4/kKIm7J3czyrC1mF3mvqNIjmwHfg18F9uGDpHGgy6iaSd49n/y+AOGTskth9ZDu9nYpsxvbxXRWjeKLr3M8ll0iKCzajeoAZFtdFkBcOLADbTfOpILjO+F0BlVtcHYtT0XyxZfmLQA/j/8wZnI8O/vWDxBdz9c8/Fsx0haX77+jWHUD261l88nHH4dIUeFanOrH5mOPtZpHoIqAw3/6jm4Z+OKJvIe8vlxw989/yekDU+Gn76+ooYDDkq3Kc19Y1i/Ulk/Rc7stNkZ0iVZjyz3P+WJVlobjK2z3TPNcNlJjmo7hXNbUZ78AuFypnhQuHXGX8VoY7QG+yDQUeZk8ll4lVA2YQ2meWyx+iENYnNviF4K2uwfN8YNNpkWaMqU9de1lfSGJMwOpGBYXAEb2R9aMi5rA2TcVVA60xKskacjfP1BEhZYXSishGrj2u6sV5+Gz8/No2QxHdsTOBZvcOpyCY05XiafajYhZq9r9iPFTFpYjramZwVzkaciSg5VXw0xKTIWX5iVtQ2UtuAVpk+WLEPwZCSIiVFLjaMrMgJctTkqMR2ZMR+VAlbrkHwlrh1uDvD4gvF8otI/RAxh4AOCdJXoWyVIRs1P6dOX5IzTBpMT2Hwr3tuuoz6yf8ARslrk/z9//yf/x9vraz+RuO+qDxjNIQonzAZemC+8WL0c/nM48mpk+eVyrMR58n7U1bEpNFkahs4eDd/z5g0fZBTrE2Y3xu8IXpNHozcIJfAJjkNl+SxpNCDEgNugKzkYij5WwXAZrLJsqCqTLIKLXOU2GSyArdTmFvLoW5wJnJZH3hRb6m1TKaryvFmWAHgVCKV71rbgFWRxsjrxmSABZqMVYlVs6e1nttuwX6sSFmxdCP/4MOXfLE9k0UcNV5rXNYsVaTzC8Zo6HtH8prVeYcChmCpbGAMRiZ2NPPG4ss1HINBuQRK43uLrQN6Z+QeXCTUNIHKLYq1XIM4OPKmApXJncEeNPU2YfqMShkUxEXGdAo9KsIisXypSZUYKTKYVSB52XTNqBguE9klqque8b4mVZnxPQ9eo0bD2cc9r+9XqKgIDfS/6fl7l/f8xRfPqUfF/uPM4qKjdgEVDJWJaDWQsmLrG1JW9L1DBUVcJFSUeTCuNPU9qJTRIWNvHCpAWImTcj+0hDLH28aT8YyjJTQVsVGyCDWomEi1JVaaWMkFyxpCq1AJ7F7mUliCSmLUq02m2mZ0zGSjiI3C7DUpyLzUHqqtwh5kfnYHh19l4jISGj07T8Yk/ADaZppWnJ5xtGgtxj0lhVIQoy7GzRCUnh8jQ05i5CibgbJpXrNNFTD6uI5D1IzBkJLG2igOhI1YEyFp+uBQKjPE4zqd1nrOioSCqmPrG/axzHM70hrPRdXRGk9jamoTOPiKPlhi0rP9DFGTy8ZidBKHUavZhhgljqWPhv3oCNGQis3RWm5ZJKGyEiulVZnomqyz2FWdsS7iXJRNI4EaNG6jqB4y7pDQPoFSZKOLw/jESKdM1jJHVM6kYrBVzoA
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"compress_and_show(200)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 91,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAD7CAYAAACG0TnRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAyFUlEQVR4nO3de3xcdZ3/8ddnJmmatLknvaVNyiS9t1zStLbcCrSIoFVUCiwq7iq0q67Xn4r4W3/rLrpYXdbVXXVbXH+rLl6g7g+3AkJSLuUuaUBKr5BAofR+mRRoS9Pk+/tjzqRpSNKkmcs5k/fzYZzMOZM5nwyTk3e+/Xy/x5xziIiIiIgIhNJdgIiIiIiIXygci4iIiIh4FI5FRERERDwKxyIiIiIiHoVjERERERGPwrGIiIiIiCcr3QV0VVZW5iZOnJjuMkREBmzdunX7nHPl6a4jlXTOFpEg6+287atwPHHiRBobG9NdhojIgJnZtnTXkGo6Z4tIkPV23lZbhYiIiIiIR+FYRERERMSjcCwiIiIi4lE4FhERERHxKByLiMiAmNliM1vZ2tqa7lJERBJO4VhERAbEObfaObe0sLAw3aWIiCScwrGIiIiIiCfw4fj6n/2JL/zm2XSXISIi/XDbA1u4+J8eTncZIiK98tVFQE5H9PAxQpbuKkREpD+OHGtn96Gj6S5DRKRXgR85NsC5dFchIiL9kZ0Voq29I91liIj0Kvjh2AxlYxGRYMgOGW3tDqdRDRHxqQwIx+gkKyISEFnh2K+d9g6dt0XEnwIfjkNmaqsQEUkAM4uY2TozO2hm9ck4RlY4NknkuMKxiPhU4MOxAR1KxyIiibDIOTfbOVcMYGZLE32A7FDs1476jkXEr4Ifjk0T8kREemJmi7yR4Ei37UVmttz7qI/vd86t7PKwJuBAomvK9kaO29p14hYRfwr8Um6xCXk6yYqIxJlZEbCOWMCt7eEha4CFzrmomS0C6oHqbo+JOOduSnRt8Z7j4xo5FhGfCn44BtS6JiJyghd6Z3u3J50hvTB8wDkX9R7b4PUa1zrnmrzHLAduTEZtnSPHOnGLiE9lRFuFBo5FRE4WD789qAVaum1rASIAZvZV4NY+vn5QskIaORYRfwt8OA6prUJEZCBKgWi3bVGgxMxWADcD68ys2czWJfrg2VmakCci/hb8tgpTW4WIyAAU9bL9gHNuGbAsmQfPDmlCnoj4W2aMHGu5ChGR/moGSnrYHu3PF5vZUjNrNLPGvXv3DvjgJybk6bwtIv4U+HAMGjkWERmAzv7iLiLOuYb+fLFzbqVzrs45V1deXj7gg2d1TshTW4WI+FPgw3FsKTcREekP59wqIBJf29jMaoF+BeM4M1tsZitbW1sHfPzskEaORcTfAh+OQ4auAiIi0o13AZAV3t0VZnZVl90LgeXeyhTXOOeWDOS5nXOrnXNLCwsLB1zXiYuAaORYRPwp+BPyUFuFiEh3XptEAz1MsPPWMx5QIO7KzBYDi2tqagb8tfGeY4VjEfGrwI8c6wp5IiKplYiRY7VViIhf9Wvk2Lui0nJgiXOupdv2roq8frb45Utv9rbXAsu6fm2ihExdFSIiQdF5ERBNyBMRn+ozHHsBdx3QRCzgdreih22rvNs1wELv8qWLgHqg+vRL7bVKtVWIiKTQYNoq4iPHxzRyLCI+1WdbhXf50Nl9TNZocM5Vd/2AzhHlA/HLj3q9bxFvVnRCxUaOdZIVEUmVwbVV6PLRIuJvp+w5jgfcAaoltpZmVz2trTloprYKEZHAyFLPsYj43GAn5EXMrN7MnJk1d+lBLuWdV1uK0vNVmQbF0IQ8EZGgiI8c6yIgIuJXgw3HTcQm6RmxJYPqvT7l+Ed3BwZ5vHcIhTRyLCKSSoO5CEhWSCPHIuJvgwrHzrmbuvQVx9fSrAOa6XmUONp9g5ktNbNGM2vcu3fvgGswjA6lYxGRlBlUz3GW1jkWEX9L9DrHUaCRnvuLI97EvJM451Y65+qcc3Xl5eUDP6KhpgoRkYCIXz66TSPHIuJTpx2OzWx5t/uLgFudc1FvreOImUW8fbXE2i4SLmRKxyIiQXFiQp5GjkXEn055ERAv9MaXclthZiu88PtbM6snNkocBeqdc9/t8qULgeVm9gxQ2sdycIMSu3y00rGISKoM6vLRXs9xmxaoFxGfOmU49lohGoBl3bY3AZf28XVNnAjVSRPSwLGISEo551YDq+vq6m4c6NeaGVkhU8+xiPhWonuOU85ME/JERIIkOxxSW4WI+FbwwzFayk1EJEiywqYJeSLiW8EPx2YKxyIiAZIdDnFcFwEREZ/KgHAMTulYRCRlBnMREIhNytNFQETEr4IfjtGEPBGRVBrMRUAgNnJ8TD3HIuJTgQ/HIbVViIgESnZYI8ci4l+BD8dmWudYRCRIstRzLCI+lgHh2NRWISISILF1jnXmFhF/yoBwrAl5IiJBkh0O6SIgIuJbwQ/HaJ1jEZFUGuxqFeo5FhE/C3w4DqmtQkQkpQa7WkWWRo5FxMcCH441IU9EJFiyw8bxDp23RcSfgh+OUVuFiEiQZIU0ciwi/hX8cGymCXkiIgGSHdZqFSLiXxkQjjVyLCISJNnhEMc1ciwiPhX4cKwJeSIiwRK7CIjO3CLiT4EPx4Ym5ImIBMnInCwOHj6mljgR8aXgh2O1VYiIBMq0sflED7exs/VouksREXmHwIfjWFuF0rGISKoM9iIgM8bF1kd+4fXT+3oRkWQKfDjGQK1rIiKpM9iLgEwbm48ZbNhxKMGViYgMXuDDsWFo4FhEJDjyhmURKRuhcCwivhT4cBwy1FYhIhIwU8bk07L3zXSXISLyDoEPx6a2ChGRwBlfnMf26BE6dAIXEZ8JfDgO6Qp5IiKBM744l2PHO9j35tvpLkVE5CSBD8exdY7TXYWISGYws0Vmts7MIsk8zvjiXAC2R48k8zAiIgMW+HCMWborEBHJCGZ2FbACSGowBqgoygNg+0GFYxHxl8CH45CXjdVaISIyOM65Vc65aqAx2ceqiI8cHzyc7EOJiAxI4MOxEUvHaq0QETlZby0SZlZkZsu9j/pkt1D0ZGROFsV52Ro5FhHfyUp3AYNlJ40cq8VCRMTMioB1QBNQ28ND1gALnXNRM1sE1APVqaswprJ0BE3bDtLe4QiHdP4WEX8I/MhxZ1tFessQEfEN51wUmO2cW9J9nxeGD3iPwTnXAETMrKcQnVSfOG8im3e9wV2Nr6X60CIivQp8ODaLt1UoHouIxMXDbw9qgZZu21pIwSS87t5/1jjOHF/Ifz7xSqoPLSLSq8CH4zhlYxGRfikFot22RYES6Fyxog5Y5rVnJI2ZcdXs8Wze9QabdupS0iLiD4EPxyFv5FjhWESkX4q8j+4OQOeKFcXOuZv6GH1OmPfOGktWyLj72deTfSgRkX4JfDjunJCnrmMRkf5oxhsl7ibany82s6Vm1mhmjXv37h10MaUjc7hwcjm/f26HLiUtIr4Q+HB8Yp3j9NYhIhIQPfUXR7yJeafknFvpnKtzztWVl5cnpKArz6lg16GjPPXy/oQ8n4jIYAQ+HJ9Y51jpWETkVJxzq4itThEB8Fap6FcwTpZLp41mWFaIhzbvSWcZIiJAJq1znN4yRER8xVuyLb6U2wozW+EFY4CFwHIzewYo7WnJt1M892JgcU1NTUJqzR0WZtrYAta/3pqQ5xMRGYwMCMeakCci0p3XJtEALOthXxMngvPpPPdqYHVdXd2Np1/hyWZVFPD7Z2N9xyFdEERE0igD2ipinNKxiEhKmNliM1vZ2pq4kd5ZFYW88fZxth04nLDnFBE5HcEPx5qQJyKSUs651c65pYWFhQl7zpkVsedq2nYwYc8pInI6Ah+OQ7pCnohI4E0ZnU9VaR4r1jZzvL0j3eWIyBAW+HCsCXkiIqmVjLaKrHCIr71nKlt3v8mX7/ozbQrIIpImGRCONSFPRCSVktFWAfCemWP43CU13P3cDho27k7oc4uI9Ffww7F3qwl5IiLBZmZ8duEkhmeHePrlA+kuR0SGqOCHY7VViIhkjOxwiHMmFPPMKwrHIpI
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 864x288 with 2 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"fig, ax = plt.subplots(1, 2, figsize=(12, 4))\n",
|
||
|
|
"ax[0].plot(s)\n",
|
||
|
|
"ax[1].plot(s)\n",
|
||
|
|
"ax[1].set_yscale('log')"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 92,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEQCAYAAABLMTQcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAEAAElEQVR4nOz965IkWXbfi/3Wvrh7RGRmZVX3dM8ApIwEz5FJX0FKL6DBG4CmJxDwBqTxCY6RbzA8T8ADPYDMAOoBdIj5IjPJJBMBggQwl+6uyluEX/Zl6cPa7hGZVV3dA4DEcKy2WVVmRoR7uG/fe+21/+u//ktUlU/tU/vUPrVP7Terub/vC/jUPrVP7VP71P7u2yfj/ql9ap/ap/Yb2D4Z90/tU/vUPrXfwPbJuH9qn9qn9qn9BrZPxv1T+9Q+tU/tN7B9Mu6f2q9dE5E/FZGf/H1fx9+k/fd87Z/ab1b7ZNw/tV/H9j8Bf/T3fRF/w/Zf/dpF5FZE/khE3onIfxSRf33x3u+LiH7Lvz/4lvN9r2NE5F+073t3+Z2f2q9nk08890/tU/vVm4j8GPiJqv6Tv4fv/iPgj4H/BfhnwE+An6rqP2/v/86LQ26Bfw/8U1X9s28550ePEZF/AfyfgX8OvG3v/TtV/Td/F/f0qf3dt0/G/VP71P4G7e/LuDcj/BNV/b0X1/LHwGtVvfvAMX8M/JGq/ttf4Xu2Y0TkFngH/J6q/snFd/6Rqr7+29zPp/Zfr32CZf47aBfbYW2Y7o8v3vvXF1vln7w47o/asT+52ML/uP1bz/dHL475STvnesw7Efn9F+f8g/bvP764lsvvudzO/+v2+rPr/8jrf/xy2/83uc/v2bcfvJ8GVfxpu7b/+LIPMGP6Oxfwxe3f5Np/1dY87z988fJ/aD//2Qfu7/eBN7+iYX95zO+07/6Ti+v4E+BWRH73V7j8T+2/ZVPVT/9+jf9hW+4/BX4X2yr/GPj99t5mZNp7fwT88cWxfwxoO2Z9/137edvOqcAfvDjm3cUxP2mf+Z2L9/8U+I/rdby4ltt2Pe/a+X/cPnvbPvfjj71+8R3/+gPn/lXu8z9+z/79tvv5g4vr+XE7/+9evP/7H/qOX/Xa/47GyHp9tx947x3w41/xfM+OWc//gc/pZZ99+vfr9S/wqf3atuYN/gHwT/SMla7b4t/FDMzlVvyfN+/wx3r2sn6q5630T9oxP2nH/FREfgq8hBb+7cXxf9i82T8E/mV77XeAf7x+b4MKLq/lTkT+JYbR/q+YF3gHZ+/vwjt89voH+uBvep9//MFO/XB7dj/tev7txe9/IiJ/hhm5n37fk36fa/8VPPk/1W/3vv8l8G/0BSTT+vjtt/Xtt1zzh475D+2931XVn7bfV4z+zfc996f237Z9Mu6/3u3HwJ1+OAj2z4A/ezmhsYn4e7RFgPOWHSwQ9vK1P8M8yo+1P6Ftzde/X3zvujX/cxF5eS3/E/BWRLSd5yeq+n9tv3/o9Zftb3Of37e9vB8AGrT0e9i9vww4fp/2ndeuqi8hll+pNYjop6r6Lz/w9h8CH+rTj7X3jlHVOxH5E+Bfi8g/xwz6uih9MED7qf39t0+Y+3+/7fZ7fu7u5QsfMmS/YvvQhP6pqr5+8e/3VPVOLej4h+1a/khE/sW3vf6B895+z2u6+5Xv4tzeux8R+VOMGfITVf2n/Aoe+0W7/Vtc03e2Ztj/7EOGve361kDr9z3fx4755+3nn2OGfY0rfDLuv6btk+f+691+igWtfucD3vvqSd2+MNYrNe7vsv34O875U+B3P3AtW2uQwr9t3vAfAv/mY69ftP+W9wlskMPvqqp854c/3r7z2v+msMwa1P0IVLPuNP7Dt7z/Kx3Trv+SofMvsAX9k3H/NW2fjPuvcVPjF/9bzKtd+cU/Bv4Pqvov21b531+89z9jntyvuhV/2f5ARP4d5pX9z9ik/1a2xYvr/MP29++34/7sxfG/B1y+/+z1D5z7p/8V7/Pb2lswWEaNCvj7GPT07y4+82cYW+aWM/zy7Pq/z7X/TWCZZtj/V+B/WVk67Vx3Fx/7nQ+8th6/Ll4v++9jx/zuxT39PvCvgP/Tr3rtn9p/u/YJlvk1b23y/wlnFssf0oyMGtf5TzC2x59jgbB/+nfwtX+CTd4/xyb8P/0uKKdd50+BPxWR9Tr/BDOCv9fO9Q6DKv4vH3n9Q+f+r3Wf33Yvd9gO4iftXtbvv7v4zE+x+/1zzoHm/+rXfhG8/tdYv23/XtA/f4dvh6p+jC0yL9vHjvkdbPFW7Nn+0zW4+qn9erZPSUyf2rMmlrzybQG6T+1T+9T+O2mfPPdP7VP71D6138D2ybh/ap/ap/ap/Qa2T8b9U/vUPrVP7TewfcLcP7VP7VP71H4D2yfP/VP71D61T+03sH0y7p/ap/apfWq/ge2jSUz/x//bv9JSHaUKqkJVYQVxBHBifzlXcQIiihNFRPGiRF/sS1zFS6XzBYdSscS/qkKujqLOflZHqo4le3LxOFGuh5l9XHCijDnyNPec5khaAjl5qIL4inilZodmB0XsXwW8XaPMDj85qKBBt5uovaJO8UdPfBLCCNVDGUCdohHSZ5n9mxOvDyM/PDww+MxS/XYPAJ0rBGf3e8odU4nU1mdOlFQ9pxTpfCG4ypgic7bu38UEwOPU83gcqNnhQuWwn9l1idth5PPhyFQCT6nnZw83jFOk7zMiyr5LW99HVzmlyBAyb4YTY45kdbw77bh/2KP3He52wYdCmgOaHeIregzEB0+8F9SDFEjXCr9zJM8BrYLvKuUx8tmferoHRR0cf+Q4/UipQ0WK4GYhnASpkHdK2Sv1KoMDefKggksgVcjXBVkc3b2Q90q5LshQuH195DR1LD87sPu5s2fxv3uiFof78x1+EqZ/MtMfFgD+91/+HCfKUgNjjjzOPV+/u6Y8RnvOCof/FLj6q0p/V0Hg3f828PSPCxorOLj94pHPr47sQuJx6fn53Q3LZM9HftHTf+O4/i+V6/80IQrLq0iNwnzjOX0p1AhhhHBS8l4Yv1RUoH8n+AlcArUhw/Q55KtK2VUbnyo2PmfrN3VQvlhwsbLbz7zej0w5UKpwmnqWKSAOYpdxrpKznVgVRGA/LARf2zxylOIIoVKKIydPLefEW3HgQyFGG7u1Cs7Z3HWuom2OAkRfGGLGi1JUiK4SfcFLJbhKkMoQ0jbuc/U4qVzFmds4ApDVU1TsWZXIlCPH3DGmSK6OqkKpQm3f6Vwl+oqIoiqU6qgKl8nD3un2flXIxdvnqmz/tJ2vVgEFrYJWaX2guKD4YH2QF099jHTvPMNXwvCN0h0rfqpIUft3iWYryCW8ffl3bQ9m7e+XKPj60F78vd2eiLnglc0V18vPC/zf//2/+mAm9UeNu3eVqoI0w03r/A+h9CKKdxUvasbcVaIruGZ0glScVLtf9dsAAHDYMfaeEH3dHuRSPN4FM5DFs2RPSp6SHZoc1PW+KlrswVkPClIFVFBRJAtSWue0TlHfDL0otVNqhLrI+Sk4M3KyONISKHtpBrTQu7zde1JHlMrOp83A2/05KmILmsrWN9dx4jpOnHJHqp4gdTPOzlXmFFGFrg22KUfulh1zCTzOPdMcqcUjkvAXo2VKgcfit4H/9XjYJmEqnjp7JAtl9qiCPgXro+uKqCC59Us490/NHp1sES2TJ9x7ukclnipp76jRusqPzhZD3wangxpBoz7rT1TRIJShIPuMaiAfHOWqQBF09uTqUBW0q6RroQzKVZ9I2TPvzGi6WLnazYxLNCcgzoylMuaIiJoBq0BU+9keq18q6gQ/gz86ai/Uq0xVYUwRhxmJEAq1E0r21F6pXRs7wW0TTypIVWpnC2KNZpjdAt29UDobP7UHnBl+t0DpBfViI9+DW2zB8yfBL3aOWTvKoDwdzHDvukT0Fe8rziviKkPXDGl1dDGTsqeUCwOpzSmr0uaMpya32RNV60fnbO46UfA2l4Ovbc45arUFQtv5gq+bbVgdlNXRi7nDodv7wVVyzTzmgVQ9SzWn7XL+R1eoQXDF7IJ6tmtfbYsTm6eyOpS
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"compress_and_show(50)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Applications of SVD "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Principal Component Analysis (PCA)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Consider a dataset of points."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
""
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"The two 'properties' of each point in the sample ($x$ and $y$ on this plot) are correlated. Let us think of some new 'property', given by a linear combination $w_1 x + w_2 y$. PCA looks for properties that show as much variation across the data points as possible. Equivalently, PCA looks for properties that allow to reconstruct the original characteristics as well as possible."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Consider projection of the data points onto a certain (single) direction. Where should one project to maintain the largest variety of the datapoint? Related question: whats the best linear fit?"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
""
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"### Formally:\n",
|
||
|
|
"The first component should be defined in order to maximize variance. The projections of datapoints $\\mathbf {a}_{{(i)}}$ onto the vector ${\\mathbf {w}}$ are given by $\\mathbf {a}_{{(i)}}\\cdot {\\mathbf {w}}$. Suppose we've already standardized the data, then we need to solve the following optimization problem:\n",
|
||
|
|
"$$\n",
|
||
|
|
"{\\underset {\\Vert {\\mathbf {w}}\\Vert =1}{\\operatorname {\\arg \\,max}}}\\,\\left\\{\\sum _{i}\\left({\\mathbf {a}}_{{(i)}}\\cdot {\\mathbf {w}}\\right)^{2}\\right\\}\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"or, combining vectors $\\mathbf {a}_{{(i)}}$ into the rows of the matrix $\\mathbf {A}$, we may write:\n",
|
||
|
|
"\n",
|
||
|
|
"$$\n",
|
||
|
|
"{\\underset {\\Vert \\mathbf {w} \\Vert =1}{\\operatorname {\\arg \\,max} }}\\,\\{\\Vert \\mathbf {Aw} \\Vert ^{2}\\}={\\underset {\\Vert \\mathbf {w} \\Vert =1}{\\operatorname {\\arg \\,max} }}\\,\\left\\{\\mathbf {w} ^{\\top}\\mathbf {A^{\\top}} \\mathbf {Aw} \\right\\}\n",
|
||
|
|
"$$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"The key ingridient is the matrix $A^T A$, proportional to the empirical covariance matrix. It is known that for positive semidefinite matrix $A^\\top A$ the solution is delivered by the eigenvector of $A^\\top A$, corresponding to the largest eigenvalue. We can repeat this procedure and construct another direction, explaining as much variance in orthogonal directions. So, we can conclude, that the following linear operation ($n$ stays for the number of observation points and $d$ for the dimensions of observation vectors):\n",
|
||
|
|
"\n",
|
||
|
|
"$$\n",
|
||
|
|
"\\underset{n \\times d}{A} \\longrightarrow \\underset{n \\times k}{A^\\prime} = \\underset{n \\times d}{A} \\cdot \\underset{d \\times k}{W} \n",
|
||
|
|
"$$\n",
|
||
|
|
"\n",
|
||
|
|
"describes the projection of data onto the $k$ principal components, where $W$ contains first (by the size of eigenvalues) $k$ eigenvectors of $A^\\top A$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Formally, we may start with SVD decomposition of $A$:\n",
|
||
|
|
"\n",
|
||
|
|
"$$A = U \\Sigma W^\\top$$\n",
|
||
|
|
"\n",
|
||
|
|
"Then, consider matrix $A^\\top A$:\n",
|
||
|
|
"\n",
|
||
|
|
"$$A^\\top A=W \\Sigma^2 W^\\top$$\n",
|
||
|
|
"\n",
|
||
|
|
"Which corresponds to the eigendecomposition of matrix $A^\\top A$, where $W$ stands for the matrix of eigenvectors of $A^\\top A$, while $\\Sigma^2$ contains eigenvalues of $A^\\top A$. In other words, columns of $W$ are right singular vectors of $A$.\n",
|
||
|
|
"\n",
|
||
|
|
"In the end:\n",
|
||
|
|
"$$\n",
|
||
|
|
"A \\longrightarrow A \\cdot W = U \\Sigma W^\\top W = U \\Sigma\n",
|
||
|
|
"$$\n",
|
||
|
|
"\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"# Example of dimensionality reduction via PCA: Wine dataset"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 93,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import matplotlib.pyplot as plt\n",
|
||
|
|
"from sklearn.datasets import load_wine\n",
|
||
|
|
"from sklearn.preprocessing import StandardScaler"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 94,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"(178, 13) \n",
|
||
|
|
"\n",
|
||
|
|
"[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n",
|
||
|
|
" 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
|
||
|
|
" 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1\n",
|
||
|
|
" 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2\n",
|
||
|
|
" 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"wine = load_wine()\n",
|
||
|
|
"\n",
|
||
|
|
"A = wine['data']\n",
|
||
|
|
"print(wine['data'].shape, \"\\n\")\n",
|
||
|
|
"\n",
|
||
|
|
"labels = wine['target']\n",
|
||
|
|
"print(labels)\n",
|
||
|
|
"\n",
|
||
|
|
"classes = [0, 1, 2]\n",
|
||
|
|
"colors = ['red', 'green', 'blue']"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 95,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"# data standardization\n",
|
||
|
|
"A_std = StandardScaler().fit_transform(A)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 97,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"array([[ 1.51861254, -0.5622498 , 0.23205254],\n",
|
||
|
|
" [ 0.24628963, -0.49941338, -0.82799632],\n",
|
||
|
|
" [ 0.19687903, 0.02123125, 1.10933436],\n",
|
||
|
|
" [ 1.69154964, -0.34681064, 0.4879264 ],\n",
|
||
|
|
" [ 0.29570023, 0.22769377, 1.84040254],\n",
|
||
|
|
" [ 1.48155459, -0.51736664, 0.30515936],\n",
|
||
|
|
" [ 1.71625494, -0.4186237 , 0.30515936],\n",
|
||
|
|
" [ 1.3086175 , -0.16727801, 0.89001391],\n",
|
||
|
|
" [ 2.25977152, -0.62508622, -0.7183361 ],\n",
|
||
|
|
" [ 1.0615645 , -0.88540853, -0.352802 ]])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 97,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"A_std[:10,:3]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"How can we visually inspect the dataset? Recall that we have 178 datapoints in 13-dimensional space. What if we reduce dimension in some random way? Lets pick two random features."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 98,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"[0 3]\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"rank = 2\n",
|
||
|
|
"ix = np.random.choice(A.shape[1], 2, replace=False)\n",
|
||
|
|
"print(ix)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 99,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"def inspect_projections(projections):\n",
|
||
|
|
" for label, colour in zip(classes, colors):\n",
|
||
|
|
" plt.scatter(projections[labels == label, 0],\n",
|
||
|
|
" projections[labels == label, 1],\n",
|
||
|
|
" label = label,\n",
|
||
|
|
" c = colour)\n",
|
||
|
|
" plt.legend(loc='best')\n",
|
||
|
|
" plt.grid()\n",
|
||
|
|
" plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 100,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"rnd_projections = A[:, ix]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 102,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD7CAYAAACRxdTpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAApi0lEQVR4nO2dbWxc15nf/4cUVZmSQUpULDhSObOjZIECVQKT2i/dFJBAqi1UBEW7YrLecWFnUY+WAhy7QGDBIZq1sh0oNIzKMQp7SX0oXHAW2qUCFLDXxUIURMBCgcKikFpogN2GMoeRYyA2KdKmaFkyefrh3jucubzv99zX+f+AC5Lnvj334cwzZ87zJqSUIIQQki86khaAEEKIemjcCSEkh9C4E0JIDqFxJ4SQHELjTgghOYTGnRBCcsiOuG+4f/9+WSwWfZ1z79497N69OxqB2gDqLxzUXziov3AY+pubm/tUSvk1r+fFbtyLxSJu3Ljh65zZ2VkcO3YsGoHaAOovHNRfOKi/cBj6E0LU/ZzHZRlCCMkhNO6EEJJDaNwJISSHuBp3IcQpIcScEGJeCDEthOht2tcrhBjXtytCiFKk0hJCCPGEo0NVN9anAYwA6AUwDWBcHwOAqwCGpJQrQohhAFcAHI5MWkIIIZ5wm7n3AhiRUt6WUt4EMAFgHwDoxnxZSrkCAFLKGQAlIcRAdOISEg+1WzUUXyui41wHbv3uFmq3akmLRIgvHGfuukFv5jCA8/rvAwBum/bfBlACYD6PkMxQu1VD5e0K1h+uAwAebDxA5e0KAKB8pJykaIR4xpNDVQhREUJMQzPcK/pwX9PvBivQZ/aEZJWxq2MNw26w/nAdY1fHEpKIEP8IL806hBADUsqbQohxABUAvwdt7R1SytNNx80BOC+lvGw6v6KfhwMHDgxeunTJl5Bra2vYs2ePr3PIFtSfP+Y+nmv5+9A/OoQ7X94BAAw+PpiESJmGr79wGPo7fvz4nJTyqNfzPBn3xsFapMxdaA7WEoA/kFKONO2fA3BWX3+35OjRo5IZqvFC/fmj+FoR9dWtZMBXf/9V/OgffoRCTwELLywkJ1hG4esvHE0Zqr6Mu+OyjNk5ajhPoa2pG+vrzZScDDshWaA6VEV3V3fLWHdXN6pD1YQkIsQ/bmvuw82x60KIUwBe0aNnLkOLjinp+wYA0LCTzFM+UsbkdydR6ClAQGBn505MfneSzlSSKdwKh80AmBZC3AAwBy308WzT/iEA40KI9wH0NS/REJJlykfKDWM+OzuLY0eOJSsQIT7xEgpp60HS99OgE0JIymBtGUIIySE07oQQkkNo3AkhJIfQuBNCSA6hcSeEkBxC404IITmExl0xzaVii68VA5eKVXUdQkh74pbERHxgLhVbX60HKhWr6jqEkPaFM3eFqCoVy5KzhJCw0LgrZHF10dd41NchhLQvNO4K6e/p9zUe9XUIIe0LjbtCVJWKZclZQkhYaNwVYi4VW+gpBCoVq+o6hJD2hdEyimkuFZuG6xBC2hPO3AkhJIfQuBNCSA6hcSeEkBxC404IITmExp0QQnIIjTshhOQQGndCCMkhNO4kNmo1oFgEOjq0nzVWMSYkMpjERGKhVgMqFWBdL3ZZr2t/A0CZuVqEKIczdxILY2Nbht1gfV0bJ4Soh8adxMKiTbViu3FCSDho3Eks9NtUK7YbJ4SEg8adxEK1CnS3VjFGd7c2TghRD407iZzarRrGPili/V+V0bn3DiAkCgVgcpLOVEKigtEyJFJamn1/q46Nb/2V1niE9ekJiRTO3EmksNk3IclA404ihc2+CUkGGncSKWz2TUgy0LjHSO1WDcXXiug414Hia0XUbuU//57NvglJBhr3mDAci/XVOiQk6qt1VN6u5N7As9k3IcnAaJmYcHIs5t3Qsdk3IfHDmXtM0LFICIkTGveYoGOREBInNO4xQcdiPmGNepJWaNxjgo7F/GHUqK/XASm3atTTwJM0QIdqjNCxmC+catSzZg5JGlfjLoQYB6D3zMEMgGellCv6vmHT4b1SystKJSQkpbBGPUkzjsZdN94lAIP6z2kAFwGM6IdMWJxG407agv5+bSnGapyQpHFbcy9Bm6nfllLOADgLoHm2PiOlPNy8RSYpyTWZcEyahKyevB5pjfo06iStWdZp1FXSOM7cpZSTpqFlADeiE4e0I5lonm0hZPmtfwk8/XcYe/c7WFzUZuzVqhqZ06iTlvLNQCPLGkCivqQ06ioN+I2W+T5al2JKQogrQggphJi3WIMnxJVMNM+2EbL87lNYWAA2N4GFBXXGJI06SWv55jTqKg0IKaW3A4UYAPCSlHKkaWwcwHkp5YoQYgKa43Wv4XBtOq6i78OBAwcGL1265EvItbU17Nmzx9c5ZIu0629uzn7f4GB8ctixtraGPX//9/YHRCBkGnUy97G9UIOP2wsV9esvjbpSiaG/48ePz0kpj3o+UUrpugHoBTDt4TgJYNjpmMHBQemXa9eu+T6HbJF2/RUKUmqR4q1boZC0ZBrXrl2LXcg06qRwoSDxMrZthQvOQkX9+kujrlRi6A/ADenBXhub12WZi7Jpxu7ACrgm3/b4dbplonl2zEKGvV0UDsYos6zDyJuJ108SuFl/aOGPRkjkgL6VAIybjhsG8KLb9Thzj5849Tf1wZTsrna3zOy6q91y6oMp5/OmtJmWENrPKefDY6Whv5iFDHq7qSkpu7tbZ7Hd3WrEnfpgShYuFKR4WcjChYLr/1VK99efCnnT/PoJS9CZu+Oae9M6upmz0BKaxgHchjZjvyK1cElHjh49Km/c8De5n52dxbFjx3ydQ7aIU3/F14qor24P/i70FLDwwkIsMqgma6+/YtE6/r5Q0Jy+ceOmv7TJmzYM/QkhfK25u4VCngZw2uGQE15vRNoDljZOnqxlzmZN3qzAwmFEKSxtnDx2GbJpzZzNmrxZgca9DYkyyzBtpY0NR50QwI4d2s88ZDDWbtWw/6kfQvQuQIhN7P/6WuOZsuZgzJq8WYFVIduMqLMMjWuMXR3D4uoi+nv6UR2qJpLBaM5c3NjQfmY9g7F2q4Yf/MUMHv6P/wo83A0AWPp4D/70P3wFYEfjmcbGoDxzNgqyJm9W8JzEpAo6VOOnWX95dHjaYeeoM/DqsEvb66/4WhH1l2eB1eK2fWl0QqZNf1kjqEOVyzJtRjs5PN0ccll12C2uLgKr1gvSWX0moh4a9zajnRyebg65rDrs+nv6gR5rK57VZyLqoXFvM05+86Sv8Sxj5agzyLLDrjpURde/OAd03TPtkVhby76zmKiBxr3NePf/vetrPMuUy8DkpLYODQCdndrPQkEbz6rDrnykjP/2n4bR972XgEc+gVbSCQAElpbYx5Vo0Li3Ge205g5oBnxhQUtq/+or7afK0rxJUT5SxqdTr6Pw2NcAiJZ9LHdLABr3tqOd1tzbAWZ3Ejto3NuMtCUZkXAwu5PYQePeZpSPlDH53UkUegoQECj0FDD53clE26QFJa39PFVhLoN75sz2srjM7iR2MEO1DSkfKWfSmDeT1n6eqrDqC/rmm1v7jSzbyUltY3YnMcOZO2khK7PhtPbzVIVVX1AzhuPUcBqr7uNKsg1n7qRBlmbDeY/68eoQpeOU2MGZO2mQpdlw3qN+vDpE6TgldtC4tyF2Sy9RzYajWOqxivoBgLUHa5EtJcW5ZOWUXWsghLb2rqKEce3MdRR33EGH2ERxxx3UzlwPd0GSOFyWaTOcll76e/otK0aGmQ1HtdRjnPv8/3weS18sNcaXvliKZCkp7iUrqzK4J08C776rGXQhtIQsIHwJ49qZ66i8+QTWoZUPrm8cQuXNvQCuo/zGd8I/DEkEztzbDKellyhi4KNc6ikfKWPPzj3bxqNYSkpiycrsKH3jDe1nobBl2BuyhMhKHZssNgx743rYjbHJYrALklRA495mOC29RBEDH7XjMy7HapocuKqzUhc3vu5rnGQDLsu0GW5LL6pj4KNY6onz+nHfx5Ms/dZNSII6V/s7f4v6xiHLcWD7OMkGnLnnGMMBOPfxHPa/sh/7X9mP+mo
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"inspect_projections(rnd_projections)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"## Lets use PCA for (unsupervised) projection "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 103,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"u, s, wh = np.linalg.svd(A_std)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 115,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXUAAAD7CAYAAACVMATUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAcx0lEQVR4nO3df3BU1f3/8ecBFH8bCKGChmICxJJQZYlja4ApBbXORwUtlFYKYoWgFRg7aiID2B9qJWA7LWA1SRETtAyC1t+2Eq0i+KOGgOWX8iNfEIe2gmFRMGIg5/tH7m6XZbMbwiZ39+7rMXMn2XtPNu/jLq8czz33rrHWIiIi3tDB7QJERCR+FOoiIh6iUBcR8RCFuoiIhyjURUQ8RKEuIuIhndwuoFu3brZ3795ulyEikjTWrl27z1qbEemY66Heu3dvqqur3S5DRCRpGGN2NXdM0y8iIh6iUBcR8ZAWhboxZoQxZq0xJitKmzRjTImzrYzWVkRE2kbUOXVjTBqwFqgBfDGe6zVguLXWb4wZAawEsuNRpIiItEzUUHcCepDztdk7fzkhXmet9Ts/V2WMyTLG+Ky1NfEtWUREmhNz+iUQ1DH4gNqwfbWApmBERNpRvE6UpgP+sH1+oGucnl9ERFogXuvU05rZXxen5xcRSShHjhzB7/dTV1dHXV0d+/fvD35/+eWXM2jQIFfqileo7wAujbDfH6mxMaYQKATo1atXnEoQETkx1lrq6+sjBnP44/Bjn3/+ebPP++CDDyZ9qNcCY8P2ZVlrqyI1ttaWAWUA+fn5+uglEYkray379+9n9+7dx22ffPIJn376aTCoDx8+3OzzdOrUia5du9KlSxe6du1Kjx49yM3NDT4ObOGP09LS2q+z4TWfzA8bY4qAFdbaFcaYcmNMlrW21hjjAyIGuojIyfr888/55JNPIoZ2YPvyyy+P+ZmOHTty/vnnk5mZyUUXXRQ1lAOPzzrrLIwxLvWydWKGurNccYzzsNQYU2qtXeE8nkLTKL0WGA6UGGPeB9KttWOOfzYRkegOHz7Mrl27ogZ2+NSHMYbzzjuPzMxM8vLyuPrqq8nMzDxmO++88+jYsaNLvWo/xu0Pns7Pz7e6oZdIarHWsmfPHj766CM++ugjtm7dGvx+586dNDY2HtM+IyPjuJAO3Xr27Mkpp5ziUm/anzFmrbU2P9Ix1+/SKCLe9cUXXwQDOzS4t27dyqFDh4LtzjjjDPr168ell17KuHHj6NOnTzCwL7jgAk477TQXe5FcFOoiclKOHDnCzp07jwntQHDv2bMn2M4YQ+/evcnJyWHIkCHk5OSQk5NDv379OP/88+nQQfcXjAeFuoi0yKFDh9iyZQubNm1iy5YtwfDevn07DQ0NwXZdunQhJyeHK6644pjg7tOnj0bc7UChLiLHqK+vD4Z36LZz504C5+BOOeUU+vTpQ05ODtddd10wuHNycujWrZvLPUhtCnWRFPXVV1/x0UcfHRfeO3bsOCa8A3PdEydOJDc3l7y8PLKzs+nUSfGRiPSqiHjc4cOH2bp163HhvX379uAqk06dOtG3b18uueQSxo0bR25uLrm5ufTt2zelVpV4gUJdxCOstezevZt169axbt06Nm7cyKZNm9i2bRtHjx4Fmi7A6dOnD3l5eYwdOzYY3v369ePUU091uQcSDwp1kSTU2NjI9u3bqampYd26dcGvn332GdC00iQ7O5vc3FxuuOEG8vLyyM3NJScnh86dO7tcvbQlhbpIgmtoaGDz5s3HhPf69es5ePAg0DTvPWDAAK6//noGDhyIz+fj29/+NmeccYbLlYsbFOoiCaS+vp5//etfwQCvqalh48aNwZtOnXnmmVx88cVMnDgxGOD9+/fX1IkEKdRFXHLgwAHWr19/zBTKhx9+GJz/7tKlCz6fj2nTpuHz+Rg4cCB9+/ZNifuXSOsp1EXagbWWrVu3snr1alavXs2aNWvYtm1b8HiPHj3w+Xxcf/31+Hw+fD4fvXr1Sro7BIr7FOoibeDrr7+mpqbmmBDft28fAOnp6RQUFASnUAYOHMh5553ncsXiFQp1kTjw+/288847wQB/7733+OqrrwDo06cP11xzDQUFBQwePJicnByNwKXNKNRFWuHjjz8+ZhS+YcMGrLV07NgRn8/HbbfdRkFBAQUFBRqFS7tSqIvEcPToUTZu3HhMiO/evRuAs846i8svv5zRo0dTUFDAZZddxplnnulyxZLKFOoiYay1VFdX8/e//501a9bw9ttvBz9pp2fPngwZMoTBgwdTUFDAgAEDdA8USSh6N4o4PvnkE5544gkqKyvZsmULAHl5edx4443BEP/mN7+p+XBJaAp1SWmHDh3imWeeobKyktdeew1rLYMHD6a8vJwbbriBrl27ul2iyAlRqEvKaWxs5M0336SiooIVK1Zw6NAhLrzwQu69917Gjx9Pdna22yWKtJpCXVLG1q1bqaysZMmSJXz88cecc845/OQnP2HChAkUFBTo49TEExTq4ml1dXUsW7aMyspK3n33XTp06MCVV15JSUkJI0eO5PTTT3e7RJG4UqiL5zQ0NPC3v/2NyspKnn/+eb7++mvy8vKYN28e48aNo0ePHm6XKNJmFOriCdZa1q9fT0VFBX/5y1/Yu3cvGRkZ/PznP2fChAlccsklWrUiKUGhLknt3//+N08++SQVFRVs3LiRU089leuuu46bbrqJq666Sh/FJilHoS5J58iRI7zwwguUlZXx6quv0tjYyHe+8x0eeeQRxo4dS5cuXdwuUcQ1CnVJGnV1dSxatIiHH36YXbt2kZmZyYwZM5gwYQL9+vVzuzyRhKBQl4S3ceNGFixYwJIlS6ivr2fYsGH84Q9/4Nprr9UHRoiEUahLQjp69Cgvvvgi8+fP5/XXX+e0005j/PjxTJs2jQEDBrhdnkjCUqhLQtm/fz+PPfYYCxcuZOfOnWRmZjJnzhwmTZpEenq62+WJJDyFuiSEzZs3s2DBAiorK/nyyy8ZOnQoDz30ECNHjtRdEEVOgP61iGuOHj3Kyy+/zPz586mqqqJz586MGzeOadOmcckll7hdnkhSihnqxpg0YIbz0AdMsdbWRmjnA6YAfqddibW2Km6Vimf4/X4WL17MwoULqa2t5fzzz+e3v/0tkydPplu3bm6XJ5LUWjJSfw0Ybq31G2NGACuBSLexW26tzQYwxmQBawEtGJagLVu2sHDhQioqKjh06BCDBw9mzpw5jBo1ShcJicRJ1FB3QrzOWusHsNZWGWOyjDE+a21NSDsfELzxtLW21hhT10Y1SxJpbGzklVdeYf78+bz66quceuqp3HjjjUybNg2fz+d2eSKeE2uk7gPCp1pqgSwgGOrW2hpjTJoxZjkwGSgEiuNZqCSX+vp6ysrKWLhwIdu3b6dnz57cf//9TJ48me7du7tdnohnxbqBdDpNc+Sh/ISMykNcAYwG9jvfaz49BVlrWbZsGRdddBF33HEHGRkZLF26lJ07dzJz5kwFukgbizVST2tmf6SplTHOBlBO05x6xI+QMcYU0jSap1evXjGLlORQXV3NHXfcwZo1a7j44ot5/PHHGTZsmNtliaSUWCP1HUQelftDHxhjRgNYa1dYa1cAFwJZgf3hrLVl1tp8a21+RkbGiVctCWXPnj1MnDiRSy+9lG3btlFeXs7atWsV6CIuiBXqgfnzUFkRlipeStPIHADnxGoNx0/diIfU19fzwAMP0K9fP5YuXUpxcTHbtm1j0qRJuieLiEuihroz6s5yligGVrkEA90YU+Qce5+mefTA/jSaVs1oXt2DQufNZ82axVVXXcXmzZuZM2cO55xzjtvliaS0lqxTHw6UGGPeB9KttWNCjk0Baq21K5yljqU0Tdmk87/5dfGQ8HnziooKvve977ldlog4Yoa6sx49YkAHLjZyvp8bx7okwezZs4eZM2fy+OOP0717d8rLy7n55ps1zSKSYHTvF4mqvr6e3//+9zz44IM0NDRQVFTEzJkzNc0ikqAU6hKRtZbly5dTVFTErl27uOGGG5g7dy7Z2RFXqYpIgoi1+kVSUHV1NUO
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"cum_s = np.cumsum(s**2)\n",
|
||
|
|
"plt.plot(cum_s/cum_s[-1], color = 'black')\n",
|
||
|
|
"plt.bar(x=np.arange(len(s)), height=s**2/cum_s[-1]);"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Lets project onto the principal components: $A \\longrightarrow U \\Sigma$"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 110,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"projections = u[:,:rank] @ np.diag(s[:rank])"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 111,
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "slide"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD7CAYAAAB+B7/XAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAlVklEQVR4nO2dX2wVR5r2n7KNhQ2sjY/348JZH6+Z0V4xfzDfPazJxfBptLsRnrk4YhGsYiZIK3KxChr5YjY7siJHo8yglZLgSEEkPhIbs9kZjYZcBAR3ewORZrmINBOQzYJWM8EHMxDM2sH1XfQ55vicru7q7up/1c9PaoH7dHfVqVP19FtvvVUlpJQghBBiDx1pZ4AQQohZKOyEEGIZFHZCCLEMCjshhFgGhZ0QQiyDwk4IIZbRlXYGAGBwcFCOjIzEns5XX32Fbdu2xZ5OHmHZqGHZuMNyUZNE2dy4ceO+lPLP3T7LhLCPjIzg+vXrsadz7do17N+/P/Z08gjLRg3Lxh2Wi5okykYIsaj6jK4YQgixDAo7IYRYBoWdEEIsg8JOCCGWQWEnhBDLoLATkiOqVWBkBOjocP6tVtPOEckimQh3JIT4U60Ck5PAkyfO34uLzt8AUKmkly+SPWixE5ITpqaei3qDJ0+c84Q0Q2EnJCfcuRPsPCkuFHZCcsLwcLDzpLhQ2AnJCdPTQG/v5nO9vc55QpqhsBOSEyoVYHYWKJcBIZx/Z2c5cEraYVQMITmiUqGQE39osRNCiGVQ2AkhxDIo7IQQYhkUdkIIsQwKOyGEWAaFnRBCLIPCTgghlkFhJ4QQy6CwE0KIZVDYCSHEMijshBBiGRR2QgixDAp7A24mSUgmMNkUi9qsubojwM0kCckIJptikZs1LXaAm0kSkhFMNsUiN2sKO8DNJAnJCCabYpGbNYUd4GaShGQEk02xyM2awg5wM0lCEqBaBQYHnW39hHD+3zqYabIpFrlZU9gBbiZJSMxUq8Dx48DS0vNzS0vAsWObxd1kUyxys6awN6hUgIUFYH3d+bf51y9qzBQhhpiaAlZX28+vrbUPZno1xaCYfFaeYLijH0WOmSLEEF4DlkUYzEwaWux+FDlmihBDeA1YFmEwM2ko7H4UOWaKWE9SXsbpaaC7u/38li3FGMxMGgq7H0WOmSJW0/AyLi4CUj73MpoW92r1uY+9o0lxSiXg3Dl6NOOAwu6HLTFTHAAmLSThZWx+eQDOICbgiPqZMxT1uKCw+2FDzFRSphnJFUl4Gd1eHoAT6sgqGB8Udh3yHjPFAWDiQhJeRq+XBKtgfFDYiwAHgIkLSXgZ/V4Si4vpeAZt90xS2IsAB4CJC0l4Gd1eHq0k7RksgmeSwp4kaZkJtgwAE+PE7WVsvDxKJe/rVG6ZOJpMETyT+RX2vPWl0jQTbBgAJrnArVlWKsD9+8DcnFP1VLR6BuNqMoXwTEopjR4AZgA8qB/zAPr97hkbG5OBmJuTsrdXSuf3do7eXue8B1evXg2WjknK5c35bRzlcnp5aiLVssk4LBt3WstFt1mqmkKjOTSuV11XKjmfCbH5el2SaIpJ1BkA16VCU41a7EKIgwBGAYwBmABwEMB7JtMAkM++VCHMBFJkdJull9+92SpXNY2lpWhWfBE8k6ZdMaMAXpZS3pZSXgZwGo64myWPIskBTGI5us2y2TPoRuNloNs0gtp0RfBMGhV2KeWslHK56VQNwHWTaQDIp0gWwUwghSZIs2wM2grhfs+dO3oRNc3XByHvU1P8EI6rJqaHCzEP4N+klBddPpsEMAkAu3btGrtw4YL+g2s1pw/WmJ8MOKM15TIwMKC87fHjx9i+fbt+Oqap1YB795xFM7q7gaEhz/wmSeplk2FYNu60lkuYZnnzpvs67d3dwJ497U3m6683P79BVxfw7W9H/EIGSaLOHDhw4IaUcp/rhyrne9QDwF4A8zrXBh48ldIZMQk4glLIQTDNcipk2WjCsnHHrVyCNsugcRClknpANUukPXgay0YbQoh+AD+WUk7E8XwATt/Jtv6TabhJCEmYoM2yce3UlONOGR52XDCqZ9Rqwc4Xlbji2N+LVdSLjm4Mf5ajh/I2D4HERhB/d1LDa3mvnsaFve5XPyuEGBVC7K0fo6bTKQRutSvIrI2sRg8VYU43CYyOmCYRg2BF9VT5aMIcAM4CkC7Ha173hfKxhyBxX2mIcYBN97o5H1VORrfZFQFmYiRaNhmfrNUKfezumCyXIL72KM1KBxPVM20fu+lwxxNSSuFyvGkynVwQ9bWvcqMsLblf32yFN0yfxcX2eLIshFhmtSdBUiOI1zDuUEUbqmd+14rJOlH920FrUcPJ2LpljZTPxd1vJsbJk07cmBDOvydPBstD0LzqnifWkyUxtaF6Uth1CTqaErWmqmpRqeTtZHR7oUjpiHqzedP6fX73O+Cdd4Bnz5zPnz1z/o5D3DlZi7SQJTG1oXpS2HUI41aJWlNVtevMGe/50DovFLfv8+iR+32zs3r5DUIR5nSTQPiJqSqOII7IFSuqp8r5nuSR+cHTMKMpIVegbHtG0FEinby6XHP1Zz9zvw/Qz6+lFHXw1K/6mS4XVXpuTWnLFim7u6M1rzixavDUWsK4VUy89sOMEun0I4M4Ljs79a8l1pBGyJ+qurt5F9fW2pciyMoUjSxAYdfBz62i6hNGGb5vPLMxkCmEXn9T54USxHHZmKlKCkWW5rYFsUPyFLkSJxR2Hbys4DhMm9bIlsaApu6z/V4obt+nowMYH39uoXd2Aq+8Arz9dvjvQXJLHqJUol5rMxR2Hbys4DhMG7dnmno24P59ymXg8mVn+TwpnX8p6oUl61EqW7Y4qz02k7fIlTihsOuisoLjMG387m18HiUsoPX7ZGT5YJINshTy52aHnDsHvP9+ziNXYoTCHgQ3IY3DtPG7d3jYkgUtSFbJWsifm10VdQZq3hf68oLCrotKSA8dMm/aeG0d03h2lka3iJVkeZehqKLsZxflXfQp7LqohPTSJfOmTeumkI0BzeZnNwZWW1G5cfxmeNy8mb/aSwqJic6ql11kQ2eYwq6Lly89TtNGCOCFF4C5uefPrlbVm0W6+crdaurx48CxY8/Pra7mr/aSzBKnxWuis+rVnG3oDFPYdYnqSw9S0/1Mhqkp57wbf/pT+7PdaurqqjPLo5m81V6SSU6eBI4cic/iNRGv4NWcsxTqGRYKuy5ufm8hHB+7H0H7dn4mg1cNW1trF2fO8CAJUa0C777bbneYtBlMxCt4Rf1kKdQzLBR2XSoV4OjRzS4QKYHz5/1NkaB9Oz+Twa+Gtd7PGR4kIbw6k6ZsBhOhmF5RP1kK9QwLhT0Ily6FM0WC9u38TIZDh9Q+drf7vaJsmslb7SWZw0u8TdkMpkIxVUNjWQv1DAOFPQhhnW9B+3Z+SxicP682i9zEuVFTSyV1Hru68ld7SeZQVWkhzNoMcYdiZjnUUwcKexDCOt+C9u2CLmHQwM+0WFlR57GjI3+1l2QO1VDUj37E6pUkFPYghHW+henbBV3CAPA2LbxeCED7GqhpkfeZIQXHrap/+CGXHUoaCnsQojjfTPXtvPq6XiLo5y5qXVEpDWyYGUJy7cawxa6gsAcl7Vo7Pe0+cCql9yCul7uotxcYGoqet6icOpX/mSEkt6jsipMn8yf2FPa8UamEiydTRcaUSk6vI+3VHatVYGnJ/TPG1pMEUEUlv/tu/jqRFPa4CdO3U93TOK/Cyyp3cyPNzQH372ejrxy2t0GIIVT2Q5yTreKCwh4nuj7jZiEfHNy8hktzf7B5V6VWdAdxs+r89OttEBKAMPZUEPsh651ICnuc6Mw4bRX/pSX3NVxmZ8OHOSY9ImSyVZVK2XoBkdgwVU3DjsGrQjXdyHwnUkqZ+jE2NiaT4OrVq4mks4EQUjp1a/MhxPNrymX3a3SP5me5MTcnZW9v+z2Ak/bcnJTSo2zm5pzrhNh0feD0Xnkl+H29vf7pJUDi9SYnmCwXkz+/qkmVy3r5aK7ur7yil6/W+/79368Gz3hAAFyXCk1NXdSlzcKuU8NU4t96dHaGq61+L45
|
||
|
|
"text/plain": [
|
||
|
|
"<Figure size 432x288 with 1 Axes>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {
|
||
|
|
"needs_background": "light"
|
||
|
|
},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"inspect_projections(projections)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {
|
||
|
|
"slideshow": {
|
||
|
|
"slide_type": "fragment"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"source": [
|
||
|
|
"Its clear that projections obtained by SVD capture the data variety much better than just two random columns. The resulting projections can be fed to your favourite classification algorithm. In this case, we reduced (in an unsupervised manner) the data in 13-dimensional space to the data in 2-dimensional space!"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"metadata": {
|
||
|
|
"celltoolbar": "Slideshow",
|
||
|
|
"kernelspec": {
|
||
|
|
"display_name": "Python 3 (ipykernel)",
|
||
|
|
"language": "python",
|
||
|
|
"name": "python3"
|
||
|
|
},
|
||
|
|
"language_info": {
|
||
|
|
"codemirror_mode": {
|
||
|
|
"name": "ipython",
|
||
|
|
"version": 3
|
||
|
|
},
|
||
|
|
"file_extension": ".py",
|
||
|
|
"mimetype": "text/x-python",
|
||
|
|
"name": "python",
|
||
|
|
"nbconvert_exporter": "python",
|
||
|
|
"pygments_lexer": "ipython3",
|
||
|
|
"version": "3.9.12"
|
||
|
|
}
|
||
|
|
},
|
||
|
|
"nbformat": 4,
|
||
|
|
"nbformat_minor": 4
|
||
|
|
}
|