Skip to content

Commit

Permalink
change title
Browse files Browse the repository at this point in the history 
LU decomposition
  • Loading branch information
@omartinsky
omartinsky committed Feb 14, 2018
1 parent 10abba0 commit 98dcd3a
Show file tree
Hide file tree
Showing 2 changed files with 118 additions and 74 deletions.
190 changes: 117 additions & 73 deletions LU_decomposition/LU_Decomposition.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
}
Expand All @@ -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",
Expand Down Expand Up @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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",
Expand All @@ -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": {
Expand All @@ -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
}
2 changes: 1 addition & 1 deletion autodiff_templated/README.md
View file
Original file line number Diff line number Diff line change
@@ -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

0 comments on commit 98dcd3a

Please sign in to comment.