Commit aa90ae26 authored by Éamonn Murray's avatar Éamonn Murray

Add note on matrix inversion to lesson 2

parent 1f49160d
......@@ -604,7 +604,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.3"
"version": "3.5.3"
}
},
"nbformat": 4,
......
......@@ -692,6 +692,38 @@
"brand = np.random.rand(333334)\n",
"%timeit scipy.linalg.solve_banded((1, 1), Mrand, brand)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### An Application of LU decomposition: Matrix Inversion\n",
"\n",
"One of those most common applications of of these approaches is to matrix inversion. In general, to invert a matrix $M$, we want to find the matrix $M^{-1}$ such that $MM^{-1}=I$, where $I$ is the identity matrix.\n",
"\n",
"Say we have a two by two matrix, we could write this as\n",
"\n",
"$$\n",
"\\begin{pmatrix} m_{11} & m_{12} \\\\ m_{21} & m_{22} \\end{pmatrix}\n",
"\\begin{pmatrix} m_{11}' & m_{12}' \\\\ m_{21}' & m_{22}' \\end{pmatrix} =\n",
"\\begin{pmatrix} 1 & 0 \\\\ 0 & 1 \\end{pmatrix},\n",
"$$\n",
"\n",
"where we want to solve for $m_{ij}'$. We can break this down and solve for one column of $M^{-1}$ at a time, giving us 2 linear systems to solve here, but $n$ for an $n\\times n$ matrix.\n",
"E.g the first would be:\n",
"\n",
"$$\n",
"\\begin{pmatrix} m_{11} & m_{12} \\\\ m_{21} & m_{22} \\end{pmatrix}\n",
"\\begin{pmatrix} m_{11}' \\\\ m_{21}' \\end{pmatrix} =\n",
"\\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}.\n",
"$$\n",
"\n",
"So we have a linear system, of the type $\\mathbf M \\mathbf x = \\mathbf b$, where we want to solve for a set of different $\\mathbf b$ but the same $\\mathbf M$. We saw here that LU decomposition gives a good general way to solve this set of systems and find the elements of the inverse.\n",
"\n",
"This is exactly how the the [scipy.linalg.inv](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.inv.html) function operates (via an interface to LAPACK routines that actually do the work). Note that special types of matrices can have other approaches that are more optimal for computing an inverse, but there is no numpy or scipy function that calls these LAPACK routines implementing these methods currently. New features are added to SciPy all the time, so if you do know your matrix may have a better algorithm for inversion and it's large enough that it may be a time consuming part of your code I suggest double checking the SciPy documentation, or considering building an interface to a more optimal method.\n",
"\n",
"Note while you may be tempted to solve a linear system $\\mathbf M \\mathbf x = \\mathbf b$ by calculating $\\mathbf M^{-1}$ and evaluating $\\mathbf x = \\mathbf M^{-1}\\mathbf b$, it should be clear that this is a very sub-optimal approach as you will end up solving $N$ linear systems to invert the matrix, and accruing all the numerical errors associated with that, rather than solving the single linear system you had to begin with."
]
}
],
"metadata": {
......@@ -710,7 +742,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.3"
"version": "3.5.3"
}
},
"nbformat": 4,
......
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