{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Linear Algebra Examples\n", "====\n", "\n", "This just shows the machanics of linear algebra calculations with python. See Lecture 5 for motivation and understanding." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "import scipy.linalg as la\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "plt.style.use('ggplot')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Resources\n", "----\n", "\n", "- [Tutorial for `scipy.linalg`](http://docs.scipy.org/doc/scipy/reference/tutorial/linalg.html)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Exact solution of linear system of equations\n", "----" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{align}\n", "x + 2y &= 3 \\\\\n", "3x + 4y &= 17\n", "\\end{align}\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[1, 2],\n", " [3, 4]])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = np.array([[1,2],[3,4]])\n", "A" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 3, 17])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b = np.array([3,17])\n", "b" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 11., -4.])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x = la.solve(A, b)\n", "x" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(A @ x, b)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [], "source": [ "A1 = np.random.random((1000,1000))\n", "b1 = np.random.random(1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Using solve is faster and more stable numerically than using matrix inversion" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The slowest run took 5.24 times longer than the fastest. This could mean that an intermediate result is being cached \n", "1 loops, best of 3: 90.3 ms per loop\n" ] } ], "source": [ "%timeit la.solve(A1, b1)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1 loops, best of 3: 140 ms per loop\n" ] } ], "source": [ "%timeit la.inv(A1) @ b1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Under the hood (Optional)\n", "\n", "The `solve` function uses the `dgesv` fortran function to do the actual work. Here is an example of how to do this directly with the `lapack` function. There is rarely any reason to use `blas` or `lapack` functions directly becuase the `linalg` package provides more convenient functions that also perfrom error checking, but you can use Python to experiment with `lapack` or `blass` before using them in a language like C or Fortran.\n", "\n", "- [How to interpret lapack function names](http://www.netlib.org/lapack/lug/node24.html)\n", "- [Summary of BLAS functions](http://cvxopt.org/userguide/blas.html)\n", "- [Sumary of Lapack functions](http://cvxopt.org/userguide/lapack.html)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import scipy.linalg.lapack as lapack" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 11., -4.])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lu, piv, x, info = lapack.dgesv(A, b)\n", "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Basic information about a matrix\n", "----" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1.+0.j, 2.+3.j],\n", " [ 3.-2.j, 4.+0.j]])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "C = np.array([[1, 2+3j], [3-2j, 4]])\n", "C" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1.-0.j, 2.-3.j],\n", " [ 3.+2.j, 4.-0.j]])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "C.conjugate()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Trace" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def trace(M):\n", " return np.diag(M).sum()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(5+0j)" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "trace(C)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(trace(C), la.eigvals(C).sum())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Determinant" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "(-8-5j)" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.det(C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Rank" ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "2" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.linalg.matrix_rank(C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Norm" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "6.5574385243020004" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.norm(C, None) # Frobenius (default)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "6.3890280236012158" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.norm(C, 2) # largest sinular value" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "1.4765909770949921" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.norm(C, -2) # smallest singular value" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 6.38902802, 1.47659098])" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.svdvals(C)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Least-squares solution\n", "----" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 11., -4.])" ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "la.solve(A, b)" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 11., -4.])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, resid, rank, s = la.lstsq(A, b)\n", "x" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[1, 2],\n", " [2, 4]])" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A1 = np.array([[1,2],[2,4]])\n", "A1" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 3, 17])" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "b1 = np.array([3, 17])\n", "b1" ] }, { "cell_type": "code", "execution_count": 27, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "singular matrix\n" ] } ], "source": [ "try:\n", " la.solve(A1, b1)\n", "except la.LinAlgError as e:\n", " print(e)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 1.48, 2.96])" ] }, "execution_count": 28, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, resid, rank, s = la.lstsq(A1, b1)\n", "x" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "collapsed": true }, "outputs": [], "source": [ "A2 = np.random.random((10,3))\n", "b2 = np.random.random(10)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "expected square matrix\n" ] } ], "source": [ "try:\n", " la.solve(A2, b2)\n", "except ValueError as e:\n", " print(e)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0.4036226 , 0.38604513, 0.40359296])" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, resid, rank, s = la.lstsq(A2, b2)\n", "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Normal equations\n", "\n", "One way to solve least squares equations $X\\beta = y$ for $\\beta$ is by using the formula $\\beta = (X^TX)^{-1}X^Ty$ as you may have learnt in statistical theory classes (or can derive yourself with a bit of calculus). This is implemented below.\n", "\n", "Note: This is not how the `la.lstsq` function solves least square problems as it can be inefficent for large matrices." ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def least_squares(X, y):\n", " return la.solve(X.T @ X, X.T @ y)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0.4036226 , 0.38604513, 0.40359296])" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "least_squares(A2, b2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Matrix Decompositinos\n", "----" ] }, { "cell_type": "code", "execution_count": 35, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0.6],\n", " [ 0.6, 4. ]])" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "A = np.array([[1,0.6],[0.6,4]])\n", "A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### LU" ] }, { "cell_type": "code", "execution_count": 36, "metadata": { "collapsed": false }, "outputs": [], "source": [ "p, l, u = la.lu(A)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1., 0.],\n", " [ 0., 1.]])" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p" ] }, { "cell_type": "code", "execution_count": 38, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0. ],\n", " [ 0.6, 1. ]])" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "l" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0.6 ],\n", " [ 0. , 3.64]])" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u" ] }, { "cell_type": "code", "execution_count": 40, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(p@l@u, A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Choleskey" ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 1. , 0.6 ],\n", " [ 0. , 1.9078784]])" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "U = la.cholesky(A)\n", "U" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 42, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(U.T @ U, A)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# If workiing wiht complex matrices\n", "np.allclose(U.T.conj() @ U, A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### QR" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "collapsed": false }, "outputs": [], "source": [ "Q, R = la.qr(A)" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-0.85749293, -0.51449576],\n", " [-0.51449576, 0.85749293]])" ] }, "execution_count": 45, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Q" ] }, { "cell_type": "code", "execution_count": 46, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose((la.norm(Q[:,0]), la.norm(Q[:,1])), (1,1))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": 47, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(Q@R, A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Spectral" ] }, { "cell_type": "code", "execution_count": 48, "metadata": { "collapsed": true }, "outputs": [], "source": [ "u, v = la.eig(A)" ] }, { "cell_type": "code", "execution_count": 49, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 0.88445056+0.j, 4.11554944+0.j])" ] }, "execution_count": 49, "metadata": {}, "output_type": "execute_result" } ], "source": [ "u" ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[-0.98195639, -0.18910752],\n", " [ 0.18910752, -0.98195639]])" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "v" ] }, { "cell_type": "code", "execution_count": 51, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 51, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose((la.norm(v[:,0]), la.norm(v[:,1])), (1,1))" ] }, { "cell_type": "code", "execution_count": 52, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(v @ np.diag(u) @ v.T, A)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Inverting A" ] }, { "cell_type": "code", "execution_count": 53, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(v @ np.diag(1/u) @ v.T, la.inv(A))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Powers of A" ] }, { "cell_type": "code", "execution_count": 54, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(v @ np.diag(u**5) @ v.T, np.linalg.matrix_power(A, 5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### SVD" ] }, { "cell_type": "code", "execution_count": 55, "metadata": { "collapsed": true }, "outputs": [], "source": [ "U, s, V = la.svd(A)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.18910752, 0.98195639],\n", " [ 0.98195639, -0.18910752]])" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "U" ] }, { "cell_type": "code", "execution_count": 57, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose((la.norm(U[:,0]), la.norm(U[:,1])), (1,1))" ] }, { "cell_type": "code", "execution_count": 58, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 4.11554944, 0.88445056])" ] }, "execution_count": 58, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s" ] }, { "cell_type": "code", "execution_count": 59, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([[ 0.18910752, 0.98195639],\n", " [ 0.98195639, -0.18910752]])" ] }, "execution_count": 59, "metadata": {}, "output_type": "execute_result" } ], "source": [ "V" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose((la.norm(V[:,0]), la.norm(V[:,1])), (1,1))" ] }, { "cell_type": "code", "execution_count": 61, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.allclose(U @ np.diag(s) @ V, A)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.5.1" } }, "nbformat": 4, "nbformat_minor": 0 }