Morten Hjorth-Jensen
1443 lines
42 KiB
Plaintext
1443 lines
42 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Linear Algebra, Handling of Arrays and more Python Features\n",
|
|
"\n",
|
|
"## Introduction\n",
|
|
"\n",
|
|
"The aim of this set of lectures is to review some central linear algebra algorithms that we will need in our \n",
|
|
"data analysis part and in the construction of Machine Learning algorithms (ML). \n",
|
|
"This will allow us to introduce some central programming features of high-level languages like Python and \n",
|
|
"compiled languages like C++ and/or Fortran. \n",
|
|
"\n",
|
|
"As discussed in the introductory notes, these series of lectures focuses both on using\n",
|
|
"central Python packages like **tensorflow** and **scikit-learn** as well\n",
|
|
"as writing your own codes for some central ML algorithms. The\n",
|
|
"latter can be written in a language of your choice, be it Python, Julia, R,\n",
|
|
"Rust, C++, Fortran etc. In order to avoid confusion however, in these lectures we will limit our\n",
|
|
"attention to Python, C++ and Fortran. \n",
|
|
"\n",
|
|
"\n",
|
|
"## Important Matrix and vector handling packages\n",
|
|
"\n",
|
|
"There are several central software packages for linear algebra and eigenvalue problems. Several of the more\n",
|
|
"popular ones have been wrapped into ofter software packages like those from the widely used text **Numerical Recipes**. The original source codes in many of the available packages are often taken from the widely used\n",
|
|
"software package LAPACK, which follows two other popular packages\n",
|
|
"developed in the 1970s, namely EISPACK and LINPACK. We describe them shortly here.\n",
|
|
"\n",
|
|
" * LINPACK: package for linear equations and least square problems.\n",
|
|
"\n",
|
|
" * LAPACK:package for solving symmetric, unsymmetric and generalized eigenvalue problems. From LAPACK's website <http://www.netlib.org> it is possible to download for free all source codes from this library. Both C/C++ and Fortran versions are available.\n",
|
|
"\n",
|
|
" * BLAS (I, II and III): (Basic Linear Algebra Subprograms) are routines that provide standard building blocks for performing basic vector and matrix operations. Blas I is vector operations, II vector-matrix operations and III matrix-matrix operations. Highly parallelized and efficient codes, all available for download from <http://www.netlib.org>.\n",
|
|
"\n",
|
|
"When dealing with matrices and vectors a central issue is memory\n",
|
|
"handling and allocation. If our code is written in Python the way we\n",
|
|
"declare these objects and the way they are handled, interpreted and\n",
|
|
"used by say a linear algebra library, requires codes that interface\n",
|
|
"our Python program with such libraries. For Python programmers,\n",
|
|
"**Numpy** is by now the standard Python package for numerical arrays in\n",
|
|
"Python as well as the source of functions which act on these\n",
|
|
"arrays. These functions span from eigenvalue solvers to functions that\n",
|
|
"compute the mean value, variance or the covariance matrix. If you are\n",
|
|
"not familiar with how arrays are handled in say Python or compiled\n",
|
|
"languages like C++ and Fortran, the sections in this chapter may be\n",
|
|
"useful. For C++ programmer, **Armadillo** is widely used library for\n",
|
|
"linear algebra and eigenvalue problems. In addition it offers a\n",
|
|
"convenient way to handle and organize arrays. We discuss this library\n",
|
|
"as well. Before we proceed we believe it may be convenient to repeat some basic features of \n",
|
|
" matrices and vectors.\n",
|
|
"\n",
|
|
"\n",
|
|
"## Basic Matrix Features\n",
|
|
"\n",
|
|
"Matrix properties reminder"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{A} =\n",
|
|
" \\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\\\\n",
|
|
" a_{21} & a_{22} & a_{23} & a_{24} \\\\\n",
|
|
" a_{31} & a_{32} & a_{33} & a_{34} \\\\\n",
|
|
" a_{41} & a_{42} & a_{43} & a_{44}\n",
|
|
" \\end{bmatrix}\\qquad\n",
|
|
"\\mathbf{I} =\n",
|
|
" \\begin{bmatrix} 1 & 0 & 0 & 0 \\\\\n",
|
|
" 0 & 1 & 0 & 0 \\\\\n",
|
|
" 0 & 0 & 1 & 0 \\\\\n",
|
|
" 0 & 0 & 0 & 1\n",
|
|
" \\end{bmatrix}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The inverse of a matrix is defined by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{A}^{-1} \\cdot \\mathbf{A} = I\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<table border=\"1\">\n",
|
|
"<thead>\n",
|
|
"<tr><th align=\"center\"> Relations </th> <th align=\"center\"> Name </th> <th align=\"center\"> matrix elements </th> </tr>\n",
|
|
"</thead>\n",
|
|
"<tbody>\n",
|
|
"<tr><td align=\"center\"> $A = A^{T}$ </td> <td align=\"center\"> symmetric </td> <td align=\"center\"> $a_{ij} = a_{ji}$ </td> </tr>\n",
|
|
"<tr><td align=\"center\"> $A = \\left (A^{T} \\right )^{-1}$ </td> <td align=\"center\"> real orthogonal </td> <td align=\"center\"> $\\sum_k a_{ik} a_{jk} = \\sum_k a_{ki} a_{kj} = \\delta_{ij}$ </td> </tr>\n",
|
|
"<tr><td align=\"center\"> $A = A^{ * }$ </td> <td align=\"center\"> real matrix </td> <td align=\"center\"> $a_{ij} = a_{ij}^{ * }$ </td> </tr>\n",
|
|
"<tr><td align=\"center\"> $A = A^{\\dagger}$ </td> <td align=\"center\"> hermitian </td> <td align=\"center\"> $a_{ij} = a_{ji}^{ * }$ </td> </tr>\n",
|
|
"<tr><td align=\"center\"> $A = \\left (A^{\\dagger} \\right )^{-1}$ </td> <td align=\"center\"> unitary </td> <td align=\"center\"> $\\sum_k a_{ik} a_{jk}^{ * } = \\sum_k a_{ki}^{ * } a_{kj} = \\delta_{ij}$ </td> </tr>\n",
|
|
"</tbody>\n",
|
|
"</table>\n",
|
|
"### Some famous Matrices\n",
|
|
"\n",
|
|
" * Diagonal if $a_{ij}=0$ for $i\\ne j$\n",
|
|
"\n",
|
|
" * Upper triangular if $a_{ij}=0$ for $i > j$\n",
|
|
"\n",
|
|
" * Lower triangular if $a_{ij}=0$ for $i < j$\n",
|
|
"\n",
|
|
" * Upper Hessenberg if $a_{ij}=0$ for $i > j+1$\n",
|
|
"\n",
|
|
" * Lower Hessenberg if $a_{ij}=0$ for $i < j+1$\n",
|
|
"\n",
|
|
" * Tridiagonal if $a_{ij}=0$ for $|i -j| > 1$\n",
|
|
"\n",
|
|
" * Lower banded with bandwidth $p$: $a_{ij}=0$ for $i > j+p$\n",
|
|
"\n",
|
|
" * Upper banded with bandwidth $p$: $a_{ij}=0$ for $i < j+p$\n",
|
|
"\n",
|
|
" * Banded, block upper triangular, block lower triangular....\n",
|
|
"\n",
|
|
"Some Equivalent Statements. For an $N\\times N$ matrix $\\mathbf{A}$ the following properties are all equivalent\n",
|
|
"\n",
|
|
" * If the inverse of $\\mathbf{A}$ exists, $\\mathbf{A}$ is nonsingular.\n",
|
|
"\n",
|
|
" * The equation $\\mathbf{Ax}=0$ implies $\\mathbf{x}=0$.\n",
|
|
"\n",
|
|
" * The rows of $\\mathbf{A}$ form a basis of $R^N$.\n",
|
|
"\n",
|
|
" * The columns of $\\mathbf{A}$ form a basis of $R^N$.\n",
|
|
"\n",
|
|
" * $\\mathbf{A}$ is a product of elementary matrices.\n",
|
|
"\n",
|
|
" * $0$ is not eigenvalue of $\\mathbf{A}$.\n",
|
|
"\n",
|
|
"## Numpy and arrays\n",
|
|
"[Numpy](http://www.numpy.org/) provides an easy way to handle arrays in Python. The standard way to import this library is as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"n = 10\n",
|
|
"x = np.random.normal(size=n)\n",
|
|
"print(x)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we have defined a vector $x$ with $n=10$ elements with its values given by the Normal distribution $N(0,1)$.\n",
|
|
"Another alternative is to declare a vector as follows"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"x = np.array([1, 2, 3])\n",
|
|
"print(x)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we have defined a vector with three elements, with $x_0=1$, $x_1=2$ and $x_2=3$. Note that both Python and C++\n",
|
|
"start numbering array elements from $0$ and on. This means that a vector with $n$ elements has a sequence of entities $x_0, x_1, x_2, \\dots, x_{n-1}$. We could also let (recommended) Numpy to compute the logarithms of a specific array as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"x = np.log(np.array([4, 7, 8]))\n",
|
|
"print(x)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Here we have used Numpy's unary function $np.log$. This function is\n",
|
|
"highly tuned to compute array elements since the code is vectorized\n",
|
|
"and does not require looping. We normaly recommend that you use the\n",
|
|
"Numpy intrinsic functions instead of the corresponding **log** function\n",
|
|
"from Python's **math** module. The looping is done explicitely by the\n",
|
|
"**np.log** function. The alternative, and slower way to compute the\n",
|
|
"logarithms of a vector would be to write"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"from math import log\n",
|
|
"x = np.array([4, 7, 8])\n",
|
|
"for i in range(0, len(x)):\n",
|
|
" x[i] = log(x[i])\n",
|
|
"print(x)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We note that our code is much longer already and we need to import the **log** function from the **math** module. \n",
|
|
"The attentive reader will also notice that the output is $[1, 1, 2]$. Python interprets automacally our numbers as integers (like the **automatic** keyword in C++). To change this we could define our array elements to be double precision numbers as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"x = np.log(np.array([4, 7, 8], dtype = np.float64))\n",
|
|
"print(x)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"x = np.log(np.array([4.0, 7.0, 8.0])\n",
|
|
"print(x)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"To check the number of bytes (remember that one byte contains eight bits for double precision variables), you can use simple use the **itemsize** functionality (the array $x$ is actually an object which inherits the functionalities defined in Numpy) as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"x = np.log(np.array([4.0, 7.0, 8.0])\n",
|
|
"print(x.itemsize)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Having defined vectors, we are now ready to try out matrices. We can define a $3 \\times 3 $ real matrix $\\hat{A}$\n",
|
|
"as (recall that we user lowercase letters for vectors and uppercase letters for matrices)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))\n",
|
|
"print(A)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"If we use the **shape** function we would get $(3, 3)$ as output, that is verifying that our matrix is a $3\\times 3$ matrix. We can slice the matrix and print for example the first column (Python organized matrix elements in a row-major order, see below) as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))\n",
|
|
"# print the first column, row-major order and elements start with 0\n",
|
|
"print(A[:,0])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can continue this was by printing out other columns or rows. The example here prints out the second column"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))\n",
|
|
"# print the first column, row-major order and elements start with 0\n",
|
|
"print(A[1,:])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Numpy contains many other functionalities that allow us to slice, subdivide etc etc arrays. We strongly recommend that you look up the [Numpy website for more details](http://www.numpy.org/). Useful functions when defining a matrix are the **np.zeros** function which declares a matrix of a given dimension and sets all elements to zero"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"n = 10\n",
|
|
"# define a matrix of dimension 10 x 10 and set all elements to zero\n",
|
|
"A = np.zeros( (n, n) )\n",
|
|
"print(A)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"or initializing all elements to"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"n = 10\n",
|
|
"# define a matrix of dimension 10 x 10 and set all elements to one\n",
|
|
"A = np.ones( (n, n) )\n",
|
|
"print(A)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"or as unitarily distributed random numbers (see the material on random number generators in the statistics part)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"n = 10\n",
|
|
"# define a matrix of dimension 10 x 10 and set all elements to random numbers with x \\in [0, 1]\n",
|
|
"A = np.random.rand(n, n)\n",
|
|
"print(A)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"As we will see throughout these lectures, there are several extremely useful functionalities in Numpy.\n",
|
|
"As an example, consider the discussion of the covariance matrix. Suppose we have defined three vectors\n",
|
|
"$\\hat{x}, \\hat{y}, \\hat{z}$ with $n$ elements each. The covariance matrix is defined as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\hat{\\Sigma} = \\begin{bmatrix} \\sigma_{xx} & \\sigma_{xy} & \\sigma_{xz} \\\\\n",
|
|
" \\sigma_{yx} & \\sigma_{yy} & \\sigma_{yz} \\\\\n",
|
|
" \\sigma_{zx} & \\sigma_{zy} & \\sigma_{zz} \n",
|
|
" \\end{bmatrix},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where for example"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\sigma_{xy} =\\frac{1}{n} \\sum_{i=0}^{n-1}(x_i- \\overline{x})(y_i- \\overline{y}).\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The Numpy function **np.cov** calculates the covariance elements using the factor $1/(n-1)$ instead of $1/n$ since it assumes we do not have the exact mean values. For a more in-depth discussion of the covariance and covariance matrix and its meaning, we refer you to the lectures on statistics. \n",
|
|
"The following simple function uses the **np.vstack** function which takes each vector of dimension $1\\times n$ and produces a $ 3\\times n$ matrix $\\hat{W}$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\hat{W} = \\begin{bmatrix} x_0 & y_0 & z_0 \\\\\n",
|
|
" x_1 & y_1 & z_1 \\\\\n",
|
|
" x_2 & y_2 & z_2 \\\\\n",
|
|
" \\dots & \\dots & \\dots \\\\\n",
|
|
" x_{n-2} & y_{n-2} & z_{n-2} \\\\\n",
|
|
" x_{n-1} & y_{n-1} & z_{n-1}\n",
|
|
" \\end{bmatrix},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"which in turn is converted into into the $3 times 3$ covariance matrix\n",
|
|
"$\\hat{\\Sigma}$ via the Numpy function **np.cov()**. In our review of\n",
|
|
"statistical functions and quantities we will discuss more about the\n",
|
|
"meaning of the covariance matrix. Here we note that we can calculate\n",
|
|
"the mean value of each set of samples $\\hat{x}$ etc using the Numpy\n",
|
|
"function **np.mean(x)**. We can also extract the eigenvalues of the\n",
|
|
"covariance matrix through the **np.linalg.eig()** function."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Importing various packages\n",
|
|
"import numpy as np\n",
|
|
"\n",
|
|
"n = 100\n",
|
|
"x = np.random.normal(size=n)\n",
|
|
"print(np.mean(x))\n",
|
|
"y = 4+3*x+np.random.normal(size=n)\n",
|
|
"print(np.mean(y))\n",
|
|
"z = x**3+np.random.normal(size=n)\n",
|
|
"print(np.mean(z))\n",
|
|
"W = np.vstack((x, y, z))\n",
|
|
"Sigma = np.cov(W)\n",
|
|
"print(Sigma)\n",
|
|
"Eigvals, Eigvecs = np.linalg.eig(Sigma)\n",
|
|
"print(Eigvals)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"collapsed": false,
|
|
"editable": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"%matplotlib inline\n",
|
|
"\n",
|
|
"import numpy as np\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from scipy import sparse\n",
|
|
"eye = np.eye(4)\n",
|
|
"print(eye)\n",
|
|
"sparse_mtx = sparse.csr_matrix(eye)\n",
|
|
"print(sparse_mtx)\n",
|
|
"x = np.linspace(-10,10,100)\n",
|
|
"y = np.sin(x)\n",
|
|
"plt.plot(x,y,marker='x')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Gaussian Elimination\n",
|
|
"\n",
|
|
"We start with the linear set of equations"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{A}\\mathbf{x} = \\mathbf{w}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We assume also that the matrix $\\mathbf{A}$ is non-singular and that the\n",
|
|
"matrix elements along the diagonal satisfy $a_{ii} \\ne 0$. Simple $4\\times 4 $ example"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{bmatrix}\n",
|
|
" a_{11}& a_{12} &a_{13}& a_{14}\\\\\n",
|
|
" a_{21}& a_{22} &a_{23}& a_{24}\\\\\n",
|
|
" a_{31}& a_{32} &a_{33}& a_{34}\\\\\n",
|
|
" a_{41}& a_{42} &a_{43}& a_{44}\\\\\n",
|
|
" \\end{bmatrix} \\begin{bmatrix}\n",
|
|
" x_1\\\\\n",
|
|
" x_2\\\\\n",
|
|
" x_3 \\\\\n",
|
|
" x_4 \\\\\n",
|
|
" \\end{bmatrix}\n",
|
|
" =\\begin{bmatrix}\n",
|
|
" w_1\\\\\n",
|
|
" w_2\\\\\n",
|
|
" w_3 \\\\\n",
|
|
" w_4\\\\\n",
|
|
" \\end{bmatrix}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"or"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown $x_1$\n",
|
|
"from the remaining $n-1$ equations. Then we use the new second equation to eliminate the second unknown\n",
|
|
"$x_2$ from the remaining $n-2$ equations. With $n-1$ such eliminations\n",
|
|
"we obtain a so-called upper triangular set of equations of the form"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=y_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"b_{22}x_2 + b_{23}x_3 + b_{24}x_4=y_2 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"b_{33}x_3 + b_{34}x_4=y_3 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"eq:gaussbacksub\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"b_{44}x_4=y_4. \\nonumber\n",
|
|
"\\label{eq:gaussbacksub} \\tag{1}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We can solve this system of equations recursively starting from $x_n$ (in our case $x_4$) and proceed with\n",
|
|
"what is called a backward substitution. \n",
|
|
"\n",
|
|
"\n",
|
|
"This process can be expressed mathematically as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto1\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" x_m = \\frac{1}{b_{mm}}\\left(y_m-\\sum_{k=m+1}^nb_{mk}x_k\\right)\\quad m=n-1,n-2,\\dots,1.\n",
|
|
"\\label{_auto1} \\tag{2}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"To arrive at such an upper triangular system of equations, we start by eliminating\n",
|
|
"the unknown $x_1$ for $j=2,n$. We achieve this by multiplying the first equation by $a_{j1}/a_{11}$ and then subtract\n",
|
|
"the result from the $j$th equation. We assume obviously that $a_{11}\\ne 0$ and that\n",
|
|
"$\\mathbf{A}$ is not singular.\n",
|
|
"\n",
|
|
"\n",
|
|
"Our actual $4\\times 4$ example reads after the first operation"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{bmatrix}\n",
|
|
" a_{11}& a_{12} &a_{13}& a_{14}\\\\\n",
|
|
" 0& (a_{22}-\\frac{a_{21}a_{12}}{a_{11}}) &(a_{23}-\\frac{a_{21}a_{13}}{a_{11}}) & (a_{24}-\\frac{a_{21}a_{14}}{a_{11}})\\\\\n",
|
|
"0& (a_{32}-\\frac{a_{31}a_{12}}{a_{11}})& (a_{33}-\\frac{a_{31}a_{13}}{a_{11}})& (a_{34}-\\frac{a_{31}a_{14}}{a_{11}})\\\\\n",
|
|
"0&(a_{42}-\\frac{a_{41}a_{12}}{a_{11}}) &(a_{43}-\\frac{a_{41}a_{13}}{a_{11}}) & (a_{44}-\\frac{a_{41}a_{14}}{a_{11}}) \\\\\n",
|
|
" \\end{bmatrix} \\begin{bmatrix}\n",
|
|
" x_1\\\\\n",
|
|
" x_2\\\\\n",
|
|
" x_3 \\\\\n",
|
|
" x_4 \\\\\n",
|
|
" \\end{bmatrix} \n",
|
|
" =\\begin{bmatrix}\n",
|
|
" y_1\\\\\n",
|
|
" w_2^{(2)}\\\\\n",
|
|
" w_3^{(2)} \\\\\n",
|
|
" w_4^{(2)}\\\\\n",
|
|
" \\end{bmatrix},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"or"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=y_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a^{(2)}_{22}x_2 + a^{(2)}_{23}x_3 + a^{(2)}_{24}x_4=w^{(2)}_2 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a^{(2)}_{32}x_2 + a^{(2)}_{33}x_3 + a^{(2)}_{34}x_4=w^{(2)}_3 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a^{(2)}_{42}x_2 + a^{(2)}_{43}x_3 + a^{(2)}_{44}x_4=w^{(2)}_4, \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto2\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation} \n",
|
|
"\\label{_auto2} \\tag{3}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The new coefficients are"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto3\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" b_{1k} = a_{1k}^{(1)} \\quad k=1,\\dots,n,\n",
|
|
"\\label{_auto3} \\tag{4}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where each $a_{1k}^{(1)}$ is equal to the original $a_{1k}$ element. The other coefficients are"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto4\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
"a_{jk}^{(2)} = a_{jk}^{(1)}-\\frac{a_{j1}^{(1)}a_{1k}^{(1)}}{a_{11}^{(1)}} \\quad j,k=2,\\dots,n,\n",
|
|
"\\label{_auto4} \\tag{5}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"with a new right-hand side given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto5\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
"y_{1}=w_1^{(1)}, \\quad w_j^{(2)} =w_j^{(1)}-\\frac{a_{j1}^{(1)}w_1^{(1)}}{a_{11}^{(1)}} \\quad j=2,\\dots,n.\n",
|
|
"\\label{_auto5} \\tag{6}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We have also set $w_1^{(1)}=w_1$, the original vector element.\n",
|
|
"We see that the system of unknowns $x_1,\\dots,x_n$ is transformed into an $(n-1)\\times (n-1)$ problem.\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"This step is called forward substitution.\n",
|
|
"Proceeding with these substitutions, we obtain the\n",
|
|
"general expressions for the new coefficients"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto6\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" a_{jk}^{(m+1)} = a_{jk}^{(m)}-\\frac{a_{jm}^{(m)}a_{mk}^{(m)}}{a_{mm}^{(m)}} \\quad j,k=m+1,\\dots,n,\n",
|
|
"\\label{_auto6} \\tag{7}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"with $m=1,\\dots,n-1$ and a\n",
|
|
"right-hand side given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"<!-- Equation labels as ordinary links -->\n",
|
|
"<div id=\"_auto7\"></div>\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\begin{equation}\n",
|
|
" w_j^{(m+1)} =w_j^{(m)}-\\frac{a_{jm}^{(m)}w_m^{(m)}}{a_{mm}^{(m)}}\\quad j=m+1,\\dots,n.\n",
|
|
"\\label{_auto7} \\tag{8}\n",
|
|
"\\end{equation}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"This set of $n-1$ elimations leads us to an equations which is solved by back substitution.\n",
|
|
"If the arithmetics is exact and the matrix $\\mathbf{A}$ is not singular, then the computed answer will be exact.\n",
|
|
"\n",
|
|
"Even though the matrix elements along the diagonal are not zero,\n",
|
|
"numerically small numbers may appear and subsequent divisions may lead to large numbers, which, if added\n",
|
|
"to a small number may yield losses of precision. Suppose for example that our first division in $(a_{22}-a_{21}a_{12}/a_{11})$\n",
|
|
"results in $-10^{-7}$ and that $a_{22}$ is one.\n",
|
|
"one. We are then\n",
|
|
"adding $10^7+1$. With single precision this results in $10^7$.\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
" * Gaussian elimination, $O(2/3n^3)$ flops, general matrix\n",
|
|
"\n",
|
|
" * LU decomposition, upper triangular and lower tridiagonal matrices, $O(2/3n^3)$ flops, general matrix. Get easily the inverse, determinant and can solve linear equations with back-substitution only, $O(n^2)$ flops\n",
|
|
"\n",
|
|
" * Cholesky decomposition. Real symmetric or hermitian positive definite matrix, $O(1/3n^3)$ flops.\n",
|
|
"\n",
|
|
" * Tridiagonal linear systems, important for differential equations. Normally positive definite and non-singular. $O(8n)$ flops for symmetric. Special case of banded matrices.\n",
|
|
"\n",
|
|
" * Singular value decomposition\n",
|
|
"\n",
|
|
" * the QR method will be discussed in chapter 7 in connection with eigenvalue systems. $O(4/3n^3)$ flops.\n",
|
|
"\n",
|
|
"The LU decomposition method means that we can rewrite\n",
|
|
"this matrix as the product of two matrices $\\mathbf{L}$ and $\\mathbf{U}$\n",
|
|
"where"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\begin{bmatrix}\n",
|
|
" a_{11} & a_{12} & a_{13} & a_{14} \\\\\n",
|
|
" a_{21} & a_{22} & a_{23} & a_{24} \\\\\n",
|
|
" a_{31} & a_{32} & a_{33} & a_{34} \\\\\n",
|
|
" a_{41} & a_{42} & a_{43} & a_{44}\n",
|
|
" \\end{bmatrix}\n",
|
|
" = \\begin{bmatrix}\n",
|
|
" 1 & 0 & 0 & 0 \\\\\n",
|
|
" l_{21} & 1 & 0 & 0 \\\\\n",
|
|
" l_{31} & l_{32} & 1 & 0 \\\\\n",
|
|
" l_{41} & l_{42} & l_{43} & 1\n",
|
|
" \\end{bmatrix}\n",
|
|
" \\begin{bmatrix}\n",
|
|
" u_{11} & u_{12} & u_{13} & u_{14} \\\\\n",
|
|
" 0 & u_{22} & u_{23} & u_{24} \\\\\n",
|
|
" 0 & 0 & u_{33} & u_{34} \\\\\n",
|
|
" 0 & 0 & 0 & u_{44}\n",
|
|
" \\end{bmatrix}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"LU decomposition forms the backbone of other algorithms in linear algebra, such as the\n",
|
|
"solution of linear equations given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The above set of equations is conveniently solved by using LU decomposition as an intermediate step.\n",
|
|
"\n",
|
|
"The matrix $\\mathbf{A}\\in \\mathbb{R}^{n\\times n}$ has an LU factorization if the determinant\n",
|
|
"is different from zero. If the LU factorization exists and $\\mathbf{A}$ is non-singular, then the LU factorization\n",
|
|
"is unique and the determinant is given by"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"det\\{\\mathbf{A}\\}=det\\{\\mathbf{LU}\\}= det\\{\\mathbf{L}\\}det\\{\\mathbf{U}\\}=u_{11}u_{22}\\dots u_{nn}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"There are at least three main advantages with LU decomposition compared with standard Gaussian elimination:\n",
|
|
"\n",
|
|
" * It is straightforward to compute the determinant of a matrix\n",
|
|
"\n",
|
|
" * If we have to solve sets of linear equations with the same matrix but with different vectors $\\mathbf{y}$, the number of FLOPS is of the order $n^3$.\n",
|
|
"\n",
|
|
" * The inverse is such an operation \n",
|
|
"\n",
|
|
"With the LU decomposition it is rather\n",
|
|
"simple to solve a system of linear equations"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"This can be written in matrix form as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{Ax}=\\mathbf{w}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"where $\\mathbf{A}$ and $\\mathbf{w}$ are known and we have to solve for\n",
|
|
"$\\mathbf{x}$. Using the LU dcomposition we write"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{A} \\mathbf{x} \\equiv \\mathbf{L} \\mathbf{U} \\mathbf{x} =\\mathbf{w}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The previous equation can be calculated in two steps"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{L} \\mathbf{y} = \\mathbf{w};\\qquad \\mathbf{Ux}=\\mathbf{y}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"To show that this is correct we use to the LU decomposition\n",
|
|
"to rewrite our system of linear equations as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{LUx}=\\mathbf{w},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and since the determinant of $\\mathbf{L}$ is equal to 1 (by construction\n",
|
|
"since the diagonals of $\\mathbf{L}$ equal 1) we can use the inverse of\n",
|
|
"$\\mathbf{L}$ to obtain"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{Ux}=\\mathbf{L^{-1}w}=\\mathbf{y},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"which yields the intermediate step"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{L^{-1}w}=\\mathbf{y}\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and as soon as we have $\\mathbf{y}$ we can obtain $\\mathbf{x}$\n",
|
|
"through $\\mathbf{Ux}=\\mathbf{y}$.\n",
|
|
"\n",
|
|
"\n",
|
|
"For our four-dimentional example this takes the form"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"y_1=w_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"l_{21}y_1 + y_2=w_2\\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"l_{31}y_1 + l_{32}y_2 + y_3 =w_3\\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"l_{41}y_1 + l_{42}y_2 + l_{43}y_3 + y_4=w_4. \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"u_{11}x_1 +u_{12}x_2 +u_{13}x_3 + u_{14}x_4=y_1 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"u_{22}x_2 + u_{23}x_3 + u_{24}x_4=y_2\\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"u_{33}x_3 + u_{34}x_4=y_3\\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"u_{44}x_4=y_4 \\nonumber\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"This example shows the basis for the algorithm\n",
|
|
"needed to solve the set of $n$ linear equations.\n",
|
|
"\n",
|
|
"\n",
|
|
"\n",
|
|
"The algorithm goes as follows\n",
|
|
"\n",
|
|
" * Set up the matrix $\\bf A$ and the vector $\\bf w$ with their correct dimensions. This determines the dimensionality of the unknown vector $\\bf x$.\n",
|
|
"\n",
|
|
" * Then LU decompose the matrix $\\bf A$ through a call to the function `ludcmp(double a, int n, int indx, double &d)`. This functions returns the LU decomposed matrix $\\bf A$, its determinant and the vector indx which keeps track of the number of interchanges of rows. If the determinant is zero, the solution is malconditioned.\n",
|
|
"\n",
|
|
" * Thereafter you call the function `lubksb(double a, int n, int indx, double w)` which uses the LU decomposed matrix $\\bf A$ and the vector $\\bf w$ and returns $\\bf x$ in the same place as $\\bf w$. Upon exit the original content in $\\bf w$ is destroyed. If you wish to keep this information, you should make a backup of it in your calling function.\n",
|
|
"\n",
|
|
"### LU Decomposition, the inverse of a matrix\n",
|
|
"\n",
|
|
"If the inverse exists then"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{A}^{-1}\\mathbf{A}=\\mathbf{I},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"the identity matrix. With an LU decomposed matrix we can rewrite the last equation as"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{LU}\\mathbf{A}^{-1}=\\mathbf{I}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"If we assume that the first column (that is column 1) of the inverse matrix\n",
|
|
"can be written as a vector with unknown entries"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{A}_1^{-1}= \\begin{bmatrix}\n",
|
|
" a_{11}^{-1} \\\\\n",
|
|
" a_{21}^{-1} \\\\\n",
|
|
" \\dots \\\\\n",
|
|
" a_{n1}^{-1} \\\\\n",
|
|
" \\end{bmatrix},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"then we have a linear set of equations"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{LU}\\begin{bmatrix}\n",
|
|
" a_{11}^{-1} \\\\\n",
|
|
" a_{21}^{-1} \\\\\n",
|
|
" \\dots \\\\\n",
|
|
" a_{n1}^{-1} \\\\\n",
|
|
" \\end{bmatrix} =\\begin{bmatrix}\n",
|
|
" 1 \\\\\n",
|
|
" 0 \\\\\n",
|
|
" \\dots \\\\\n",
|
|
" 0 \\\\\n",
|
|
" \\end{bmatrix}.\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"In a similar way we can compute the unknow entries of the second column,"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"$$\n",
|
|
"\\mathbf{LU}\\begin{bmatrix}\n",
|
|
" a_{12}^{-1} \\\\\n",
|
|
" a_{22}^{-1} \\\\\n",
|
|
" \\dots \\\\\n",
|
|
" a_{n2}^{-1} \\\\\n",
|
|
" \\end{bmatrix}=\\begin{bmatrix}\n",
|
|
" 0 \\\\\n",
|
|
" 1 \\\\\n",
|
|
" \\dots \\\\\n",
|
|
" 0 \\\\\n",
|
|
" \\end{bmatrix},\n",
|
|
"$$"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"and continue till we have solved all $n$ sets of linear equations."
|
|
]
|
|
}
|
|
],
|
|
"metadata": {},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|