Commit 3442f025 authored by Éamonn Murray's avatar Éamonn Murray

Simplify lu decompostion code, fix some text

parent aa90ae26
......@@ -3,9 +3,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
......@@ -242,7 +240,8 @@
"So combining the various operations we have\n",
"$$ \\mathbf{P_{23}N_0P_{12}M_0} = \\mathbf{U}.$$\n",
"\n",
"Permutation matrices have the useful property that they are their own inverse, so we can write\n",
"Permutation matrices have the useful property that their inverse is their transpose, and beyond that, permutations involving a single swap are their own inverse (e.g. a permutation that swaps row 1 and 2 and keeps row 3 the same, will be inverted by swapping rows 1 and 2 again).\n",
"So we can write\n",
"$$ \\mathbf{P_{23}N_0P_{23}^{-1}P_{23}P_{12}M_0} = (\\mathbf{P_{23}N_0P_{23})P_{23}P_{12}M_0} = \\mathbf{U}$$\n",
"\n",
"Pre-multiplying by a permutation vector exchanges rows, while post-multiplying exchanges columns, so\n",
......@@ -264,9 +263,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"def pivot_matrix(M):\n",
......@@ -294,34 +291,37 @@
" The function returns P, L and U.\n",
" \"\"\"\n",
"\n",
" # Initialize L as the identity and U with zeros\n",
" L = np.identity(len(M))\n",
" U = np.zeros((len(M), len(M)))\n",
" n = len(M)\n",
" # Initialize L as the identity\n",
" L = np.identity(n)\n",
" \n",
" # Create the pivot matrix P\n",
" PM, P = pivot_matrix(M)\n",
"\n",
" # Perform the LU decomposition on the permuted matrix\n",
" for j in range(len(M)):\n",
"\n",
" # LaTeX: u_{ij} = m_{ij} - \\sum_{k=1}^{i-1} u_{kj} l_{ik}\n",
" for i in range(j+1):\n",
" s1 = sum(U[k, j] * L[i, k] for k in range(i))\n",
" U[i, j] = PM[i, j] - s1\n",
"\n",
" # LaTeX: l_{ij} = \\frac{1}{u_{jj}} (m_{ij} - \\sum_{k=1}^{j-1} u_{kj} l_{ik} )\n",
" for i in range(j, len(M)):\n",
" s2 = sum(U[k, j] * L[i, k] for k in range(j))\n",
" L[i, j] = (PM[i, j] - s2) / U[j, j]\n",
"\n",
" # Initialize U as a copy of the permuted matrix, and\n",
" # ensure it's storing floats rather than integers as\n",
" # this can lead to some unexpected results otherwise.\n",
" U = 1.0 * PM.copy()\n",
" \n",
" for k in range(n-1):\n",
" for j in range(k+1, n):\n",
" L[j, k] = U[j, k] / U[k, k]\n",
" U[j, k:] = U[j, k:] - L[j, k] * U[k, k:]\n",
" return (P, L, U)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you're still somewhat confused by how the code here corresponds to what we did previously, you should be able to find more detailed descriptions in a textbook on numerical linear algebra. For example, the chapter at <http://www.math.iit.edu/~fass/477577_Chapter_7.pdf> covers this nicely."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
"scrolled": true
},
"outputs": [],
"source": [
......@@ -329,15 +329,14 @@
"M1 = np.array([[1, -1, 2], [-4, 4, 1], [-2, 4, 3]])\n",
"\n",
"P1, L1, U1 = lu_decomp(M1)\n",
"print(\"P =\", P1, \"\\n\\nL =\", L1, \"\\n\\nU =\", U1)"
"print(\"P =\", P1, \"\\n\\nL =\", L1, \"\\n\\nU =\", U1)\n",
"print(\"\\nP^{-1}LU = \\n\", np.dot(P1.transpose(), np.dot(L1, U1)))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"%timeit lu_decomp(M1)"
......@@ -361,9 +360,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"M1 = np.array([[1, -1, 2], [-4, 4, 1], [-2, 4, 3]])\n",
......@@ -374,9 +371,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"%timeit scipy.linalg.lu_factor(M1)"
......@@ -385,9 +380,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"M1 = np.array([[3, 5, 2], [2, -3, -3], [1, 1, 1]])\n",
......@@ -413,9 +406,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"print(np.linalg.solve(M1,b1))\n",
......@@ -428,9 +419,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"%timeit np.linalg.solve(M1, b1)\n",
......@@ -447,9 +436,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"dimrand = 2000\n",
......@@ -465,14 +452,13 @@
"source": [
"## Scaling\n",
"\n",
"If you examine our Python implementation above, you'll see we have a triple nested loop. This means we likely have $O(n^3)$ scaling. Let's use `%timeit -o` to generate output we can save and try to see this ourselves. This generates a \"TimeitResult\" object which has few methods available for examining the output. We'll save the average of each run (these are in seconds)."
"If you examine our Python implementation above, you'll see we have a triple nested loop (two explicit and one implicit). This means we likely have $O(n^3)$ scaling. Let's use `%timeit -o` to generate output we can save and try to see this ourselves. This generates a \"TimeitResult\" object which has few methods available for examining the output. We'll save the average of each run (these are in seconds)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true,
"scrolled": true
},
"outputs": [],
......@@ -489,9 +475,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Set a scaling so the final points match up\n",
......@@ -527,7 +511,7 @@
"0 & \\cdots & \\dots & 0 & m_{N, N-1} & m_{N, N}\n",
"\\end{pmatrix} $$\n",
"\n",
"For example, any problem with a 1D system where we approximate some piece of it as only interacting with its left and right neighbours can usually be cast as the solution of a tridiagonal matrix problem. (And by extension if we also included interactions with next-nearest neighbours it would be a banded matrix with $k_l = k_u = 2$.\n",
"For example, any problem with a 1D system where we approximate some piece of it as only interacting with its left and right neighbours can usually be cast as the solution of a tridiagonal matrix problem. (And by extension if we also included interactions with next-nearest neighbours it would be a banded matrix with $k_l = k_u = 2$.)\n",
"\n",
"For these types of matrices there are better approaches we can use than doing e.g. a full LU factorisation using Gaussian elimination. If the matrix is diagonal, we can obtain the solutions directly. The usual method used for tridiagonal matrices is known as the Thomas algorithm, which is composed of a forward sweep which converts the matrix to upper triangular in a single pass, and is followed by backward substitution sweep to produce the result.\n",
"\n",
......@@ -546,9 +530,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"help(scipy.linalg.solve_banded)"
......@@ -581,9 +563,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"m_packed = np.empty((3, 10)) # Initialize our 3x10 array\n",
......@@ -604,9 +584,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# We can compare this to the solution we get from the full matrix:\n",
......@@ -637,9 +615,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"bsolve_timing = []\n",
......@@ -655,9 +631,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"plt.plot(bsize, np.array(bsolve_timing), \"bo\")\n",
......@@ -678,9 +652,7 @@
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"metadata": {},
"outputs": [],
"source": [
"# Timing for large matrix and large tridiagonal matrix with similar numbers of non-zero elements\n",
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment