diff --git a/LU_decomposition/LU_Decomposition.ipynb b/LU_decomposition/LU_Decomposition.ipynb index e370abe..55439e8 100644 --- a/LU_decomposition/LU_Decomposition.ipynb +++ b/LU_decomposition/LU_Decomposition.ipynb @@ -7,21 +7,19 @@ "Copyright © 2015 Ondrej Martinsky, All rights reserved\n", "\n", "[www.quantandfinancial.com](http://www.quantandfinancial.com)\n", - "# LU Decomposition" + "# Solving Linear Equations using LU Decomposition" ] }, { "cell_type": "code", "execution_count": 1, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Using matplotlib backend: Qt4Agg\n", + "Using matplotlib backend: Qt5Agg\n", "Populating the interactive namespace from numpy and matplotlib\n" ] } @@ -31,24 +29,22 @@ "%matplotlib inline\n", "import matplotlib.pyplot\n", "from mpl_toolkits.mplot3d import Axes3D\n", - "from numpy import array, multiply, dot\n", - "from scipy.linalg import lu" + "import numpy\n", + "import scipy.linalg" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Gaussian Elimination\n", - "Code for calculating **X** from **M*X=Y** where *M* is either lower or upper diagonal matrix" + "## Gaussian Elimination\n", + "Code for calculating $\\bf{X}$ from $\\bf{M \\cdot X = Y}$ where $\\bf{M}$ is either lower or upper diagonal matrix" ] }, { "cell_type": "code", "execution_count": 2, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def solve_l(m, y): # solves x from m*x = y\n", @@ -83,15 +79,20 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Solving using L U decomposition" + "## Example" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Unknowns ($\\bf{X^{org}}$), we will later try to solve these" ] }, { "cell_type": "code", "execution_count": 3, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -102,17 +103,21 @@ } ], "source": [ - "# Unknowns\n", "x_org = array([2, 4, 1])\n", "print(x_org)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Coefficients ($\\bf{M}$)" + ] + }, { "cell_type": "code", "execution_count": 4, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -125,17 +130,21 @@ } ], "source": [ - "# Coefficients\n", "m = array([[2,-1,1],[3,3,9],[3,3,5]])\n", "print(m)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Results ($\\bf{Y}$)" + ] + }, { "cell_type": "code", "execution_count": 5, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -146,33 +155,42 @@ } ], "source": [ - "# Results\n", "y = dot(m,x_org)\n", "print(y)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Decompose matrix $\\bf{M}$" + ] + }, { "cell_type": "code", "execution_count": 6, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "# Note: matrix dot-product is not commutative, but is associative\n", - "p, l, u = lu(m, permute_l=False)\n", - "pl, u = lu(m, permute_l=True)\n", + "p, l, u = scipy.linalg.lu(m, permute_l=False)\n", + "pl, u = scipy.linalg.lu(m, permute_l=True)\n", "assert (dot(p,l)==pl).all()\n", "assert (dot(pl,u)==m).all()\n", "assert (pinv(p)==p).all()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Lower diagonal matrix ($\\bf{L}$), zero element above the principal diagonal" + ] + }, { "cell_type": "code", "execution_count": 7, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -185,15 +203,20 @@ } ], "source": [ - "print(l) # Lower diagonal matrix, zero element above the principal diagonal" + "print(l)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Upper diagnonal matrix ($\\bf{U}$), zero elements below the principal diagonal" ] }, { "cell_type": "code", "execution_count": 8, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -206,15 +229,20 @@ } ], "source": [ - "print(u) # Upper diagnonal matrix, zero elements below the principal diagonal" + "print(u)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Permutation matrix ($\\bf{P}$)" ] }, { "cell_type": "code", "execution_count": 9, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -227,45 +255,69 @@ } ], "source": [ - "print(p) # Permutation matrix for \"l\"" + "print(p)" ] }, { "cell_type": "code", "execution_count": 10, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "assert (l*u==multiply(l,u)).all() # memberwise multiplication\n", - "assert (m==dot(dot(p,l),u)).all() # matrix multiplication, M=LU" + "assert (m==dot(dot(p,l),u)).all() # matrix multiplication, M=LU\n", + "assert (pinv(p)==p).all() # Assert that permutation matrix is identical to it's inverse" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Solution" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\\bf{M \\cdot X = Y}$$\n", + "$$\\bf{P \\cdot L \\cdot U \\cdot X = Y}$$\n", + "$$\\bf{L \\cdot U \\cdot X = P^{-1} \\cdot Y}$$\n", + "set $$ \\bf{Z = U \\cdot X} $$\n", + "then\n", + "$$ \\bf{ L \\cdot Z = P \\cdot Y } $$\n", + "Solve $\\bf{Z}$ in $\\bf{ L \\cdot Z = P \\cdot Y }$ using Gaussian elimination" ] }, { "cell_type": "code", "execution_count": 11, - "metadata": { - "collapsed": true - }, - "outputs": [], + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[ 27. -17. -4.]\n" + ] + } + ], "source": [ - "assert (pinv(p)==p).all()\n", - "# P*L*U*X = Y\n", - "# L*U*X = pinv(P)*Y\n", - "# set Z=U*X\n", - "# L*Z = P*Y (solve Z)\n", "z = solve(l, dot(p,y))\n", - "# solve X from U*X=Z\n", - "x = solve(u, z)" + "print(z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Solve $\\bf{X}$ in $ \\bf{U \\cdot X = Z} $ using Gaussian elimination" ] }, { "cell_type": "code", "execution_count": 12, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [ { "name": "stdout", @@ -276,25 +328,17 @@ } ], "source": [ + "x = solve(u, z)\n", "assert allclose(x_org,x)\n", "print(x)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true - }, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python [Root]", + "display_name": "Python 3", "language": "python", - "name": "Python [Root]" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -306,9 +350,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.0" } }, "nbformat": 4, - "nbformat_minor": 0 + "nbformat_minor": 1 } diff --git a/autodiff_templated/README.md b/autodiff_templated/README.md index 955b777..7761b08 100644 --- a/autodiff_templated/README.md +++ b/autodiff_templated/README.md @@ -1,3 +1,3 @@ -### C++ Template-based approach to Automatic Differentiation +### Automatic differentiation via C++ operators overloading Please refer to the corresponding article on http://www.quantandfinancial.com/2017/02/automatic-differentiation-templated.html