From 5d2bca6a62ff6207f3455efe8b37e055299aaff2 Mon Sep 17 00:00:00 2001 From: Morten Hjorth-Jensen Date: Tue, 23 Aug 2022 21:21:31 +0200 Subject: [PATCH] update --- doc/BookChapters/linalg.do.txt | 14 +- doc/LectureNotes/linalg.ipynb | 714 ++++++++++++++++++++++++++------- 2 files changed, 573 insertions(+), 155 deletions(-) diff --git a/doc/BookChapters/linalg.do.txt b/doc/BookChapters/linalg.do.txt index c4e890255..71215be66 100644 --- a/doc/BookChapters/linalg.do.txt +++ b/doc/BookChapters/linalg.do.txt @@ -153,13 +153,13 @@ print(x) or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is !bc pycod import numpy as np -x = np.log(np.array([4.0, 7.0, 8.0]) +x = np.log(np.array([4.0, 7.0, 8.0])) print(x) !ec 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 !bc pycod import numpy as np -x = np.log(np.array([4.0, 7.0, 8.0]) +x = np.log(np.array([4.0, 7.0, 8.0])) print(x.itemsize) !ec @@ -305,7 +305,10 @@ print(f"This is the mean for each column:{np.mean(a, axis=0, keepdims=True,dtype print(f"This is the mean value for each row:{np.mean(a, axis=1, keepdims=True,dtype=np.float64)}") # We print then the mean value for each row by setting keepdims=False print(f"This is the mean value for each row with keepdims false:{np.mean(a, axis=1, keepdims=False,dtype=np.float64)}") +!ec +Another useful function is the _ravel_ function, which returns a flattened array as shown in the example here. +!bc pycod # Ravel return a contiguous flattened array. print(f"Flatten the matrix:{np.ravel(a)}") # It is the same as reshaping the matrix into a one-dimensional array @@ -726,6 +729,11 @@ and continue till we have solved all $n$ sets of linear equations. - +The calculation of the inverse here assumes that it actually +exists. In many machine learning applications there may be strong +linear dependencies among the various columns and/or rows. In our +discussions of linear regression we will dive into the mathematics of +the singular value decomposition, an algorithm which will allow us to calculate the so-called pseudo-inverse. +These details will be presented in our linear regression chapter. diff --git a/doc/LectureNotes/linalg.ipynb b/doc/LectureNotes/linalg.ipynb index e915954aa..056e50939 100644 --- a/doc/LectureNotes/linalg.ipynb +++ b/doc/LectureNotes/linalg.ipynb @@ -2,10 +2,32 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "id": "43efcac0", + "metadata": { + "editable": true + }, + "source": [ + "" + ] + }, + { + "cell_type": "markdown", + "id": "074ac7c2", + "metadata": { + "editable": true + }, + "source": [ + "# Linear Algebra, Handling of Arrays and more Python Features" + ] + }, + { + "cell_type": "markdown", + "id": "0be7d59e", + "metadata": { + "editable": true + }, "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", @@ -18,9 +40,16 @@ "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", + "attention to Python, C++ and Fortran." + ] + }, + { + "cell_type": "markdown", + "id": "5e77e6c1", + "metadata": { + "editable": true + }, + "source": [ "## 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", @@ -49,9 +78,16 @@ "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", + " matrices and vectors." + ] + }, + { + "cell_type": "markdown", + "id": "2eed76da", + "metadata": { + "editable": true + }, + "source": [ "## Basic Matrix Features\n", "\n", "Matrix properties reminder" @@ -59,7 +95,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "109ff37c", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{A} =\n", @@ -79,14 +118,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2c87c75d", + "metadata": { + "editable": true + }, "source": [ "The inverse of a matrix is defined by" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "aab763d1", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{A}^{-1} \\cdot \\mathbf{A} = I\n", @@ -95,9 +140,12 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "f382e0f8", + "metadata": { + "editable": true + }, "source": [ - "\n", + "
\n", "\n", "\n", "\n", @@ -108,7 +156,16 @@ "\n", "\n", "\n", - "
Relations Name matrix elements
$A = A^{\\dagger}$ hermitian $a_{ij} = a_{ji}^{ * }$
$A = \\left (A^{\\dagger} \\right )^{-1}$ unitary $\\sum_k a_{ik} a_{jk}^{ * } = \\sum_k a_{ki}^{ * } a_{kj} = \\delta_{ij}$
\n", + "" + ] + }, + { + "cell_type": "markdown", + "id": "5d394b4e", + "metadata": { + "editable": true + }, + "source": [ "### Some famous Matrices\n", "\n", " * Diagonal if $a_{ij}=0$ for $i\\ne j$\n", @@ -141,15 +198,24 @@ "\n", " * $\\mathbf{A}$ is a product of elementary matrices.\n", "\n", - " * $0$ is not eigenvalue of $\\mathbf{A}$.\n", - "\n", + " * $0$ is not eigenvalue of $\\mathbf{A}$." + ] + }, + { + "cell_type": "markdown", + "id": "0a85f2b2", + "metadata": { + "editable": true + }, + "source": [ "## 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, + "execution_count": 1, + "id": "7043c92e", "metadata": { "collapsed": false, "editable": true @@ -164,7 +230,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2b5abf29", + "metadata": { + "editable": true + }, "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" @@ -172,7 +241,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, + "id": "4193cb75", "metadata": { "collapsed": false, "editable": true @@ -186,7 +256,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "b112c576", + "metadata": { + "editable": true + }, "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" @@ -194,7 +267,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, + "id": "18167a1f", "metadata": { "collapsed": false, "editable": true @@ -208,7 +282,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "ee482a57", + "metadata": { + "editable": true + }, "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", @@ -221,7 +298,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, + "id": "d8305632", "metadata": { "collapsed": false, "editable": true @@ -238,7 +316,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "a50fea8e", + "metadata": { + "editable": true + }, "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" @@ -246,7 +327,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, + "id": "e9263103", "metadata": { "collapsed": false, "editable": true @@ -260,14 +342,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "d68631ad", + "metadata": { + "editable": true + }, "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, + "execution_count": 6, + "id": "074cfbda", "metadata": { "collapsed": false, "editable": true @@ -275,20 +361,24 @@ "outputs": [], "source": [ "import numpy as np\n", - "x = np.log(np.array([4.0, 7.0, 8.0])\n", + "x = np.log(np.array([4.0, 7.0, 8.0]))\n", "print(x)" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "dab5cce4", + "metadata": { + "editable": true + }, "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, + "execution_count": 7, + "id": "9d335488", "metadata": { "collapsed": false, "editable": true @@ -296,13 +386,16 @@ "outputs": [], "source": [ "import numpy as np\n", - "x = np.log(np.array([4.0, 7.0, 8.0])\n", + "x = np.log(np.array([4.0, 7.0, 8.0]))\n", "print(x.itemsize)" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "64c9f3e5", + "metadata": { + "editable": true + }, "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)" @@ -310,7 +403,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, + "id": "a3b32334", "metadata": { "collapsed": false, "editable": true @@ -324,14 +418,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4ff9ce6b", + "metadata": { + "editable": true + }, "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, + "execution_count": 9, + "id": "d19f85e2", "metadata": { "collapsed": false, "editable": true @@ -346,14 +444,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4e7e8796", + "metadata": { + "editable": true + }, "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, + "execution_count": 10, + "id": "4733c8b7", "metadata": { "collapsed": false, "editable": true @@ -368,14 +470,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "df9f205b", + "metadata": { + "editable": true + }, "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, + "execution_count": 11, + "id": "19757d00", "metadata": { "collapsed": false, "editable": true @@ -391,14 +497,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "f1911274", + "metadata": { + "editable": true + }, "source": [ "or initializing all elements to" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, + "id": "4f737773", "metadata": { "collapsed": false, "editable": true @@ -414,14 +524,18 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "cd241572", + "metadata": { + "editable": true + }, "source": [ "or as unitarily distributed random numbers (see the material on random number generators in the statistics part)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, + "id": "54645585", "metadata": { "collapsed": false, "editable": true @@ -437,7 +551,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "3e993bf6", + "metadata": { + "editable": true + }, "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", @@ -446,7 +563,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4fe66190", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\hat{\\Sigma} = \\begin{bmatrix} \\sigma_{xx} & \\sigma_{xy} & \\sigma_{xz} \\\\\n", @@ -458,14 +578,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "cdd238d5", + "metadata": { + "editable": true + }, "source": [ "where for example" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "41663b63", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\sigma_{xy} =\\frac{1}{n} \\sum_{i=0}^{n-1}(x_i- \\overline{x})(y_i- \\overline{y}).\n", @@ -474,7 +600,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "e7bd73e9", + "metadata": { + "editable": true + }, "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}$" @@ -482,7 +611,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2c4d0878", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\hat{W} = \\begin{bmatrix} x_0 & y_0 & z_0 \\\\\n", @@ -497,7 +629,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "32eb00e7", + "metadata": { + "editable": true + }, "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", @@ -510,7 +645,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, + "id": "6b91d50d", "metadata": { "collapsed": false, "editable": true @@ -536,7 +672,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, + "id": "6a0aa964", "metadata": { "collapsed": false, "editable": true @@ -560,7 +697,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "ed00bf63", + "metadata": { + "editable": true + }, "source": [ "## Other Matrix and Vector Operations\n", "\n", @@ -569,7 +709,8 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, + "id": "a12e8986", "metadata": { "collapsed": false, "editable": true @@ -593,8 +734,29 @@ "# Try setting keepdims=False\n", "print(f\"This is the mean value for each row:{np.mean(a, axis=1, keepdims=True,dtype=np.float64)}\")\n", "# We print then the mean value for each row by setting keepdims=False\n", - "print(f\"This is the mean value for each row with keepdims false:{np.mean(a, axis=1, keepdims=False,dtype=np.float64)}\")\n", - "\n", + "print(f\"This is the mean value for each row with keepdims false:{np.mean(a, axis=1, keepdims=False,dtype=np.float64)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "596d78a6", + "metadata": { + "editable": true + }, + "source": [ + "Another useful function is the **ravel** function, which returns a flattened array as shown in the example here." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "fe01c225", + "metadata": { + "collapsed": false, + "editable": true + }, + "outputs": [], + "source": [ "# Ravel return a contiguous flattened array.\n", "print(f\"Flatten the matrix:{np.ravel(a)}\")\n", "# It is the same as reshaping the matrix into a one-dimensional array\n", @@ -612,7 +774,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "545df59a", + "metadata": { + "editable": true + }, "source": [ "## Gaussian Elimination\n", "\n", @@ -621,7 +786,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "796d7554", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{A}\\mathbf{x} = \\mathbf{w}.\n", @@ -630,7 +798,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "35a9c235", + "metadata": { + "editable": true + }, "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" @@ -638,7 +809,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "432d4f8e", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\begin{bmatrix}\n", @@ -663,14 +837,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "c06fa57b", + "metadata": { + "editable": true + }, "source": [ "or" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "4e8d9bd5", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \\nonumber\n", @@ -679,7 +859,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "112c2488", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \\nonumber\n", @@ -688,7 +871,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4774bcd7", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \\nonumber\n", @@ -697,7 +883,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "5bf16627", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \\nonumber\n", @@ -706,7 +895,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4da1a1c6", + "metadata": { + "editable": true + }, "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", @@ -716,7 +908,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "041684ac", + "metadata": { + "editable": true + }, "source": [ "$$\n", "b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=y_1 \\nonumber\n", @@ -725,7 +920,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "ba6214cb", + "metadata": { + "editable": true + }, "source": [ "$$\n", "b_{22}x_2 + b_{23}x_3 + b_{24}x_4=y_2 \\nonumber\n", @@ -734,7 +932,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2a1d0555", + "metadata": { + "editable": true + }, "source": [ "$$\n", "b_{33}x_3 + b_{34}x_4=y_3 \\nonumber\n", @@ -743,7 +944,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "51570380", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -756,18 +960,23 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "e51b3aa8", + "metadata": { + "editable": true + }, "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": {}, + "id": "2c78b600", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -782,20 +991,25 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "1cd1015f", + "metadata": { + "editable": true + }, "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": {}, + "id": "34093ce7", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\begin{bmatrix}\n", @@ -820,14 +1034,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "f269239d", + "metadata": { + "editable": true + }, "source": [ "or" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "d4477133", + "metadata": { + "editable": true + }, "source": [ "$$\n", "b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=y_1 \\nonumber\n", @@ -836,7 +1056,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "25bdfc8c", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a^{(2)}_{22}x_2 + a^{(2)}_{23}x_3 + a^{(2)}_{24}x_4=w^{(2)}_2 \\nonumber\n", @@ -845,7 +1068,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "34b8408c", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a^{(2)}_{32}x_2 + a^{(2)}_{33}x_3 + a^{(2)}_{34}x_4=w^{(2)}_3 \\nonumber\n", @@ -854,7 +1080,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "b8185f8d", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a^{(2)}_{42}x_2 + a^{(2)}_{43}x_3 + a^{(2)}_{44}x_4=w^{(2)}_4, \\nonumber\n", @@ -863,7 +1092,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4a218e4e", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -877,14 +1109,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "ba424207", + "metadata": { + "editable": true + }, "source": [ "The new coefficients are" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "10f17fb1", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -899,14 +1137,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "e22879fc", + "metadata": { + "editable": true + }, "source": [ "where each $a_{1k}^{(1)}$ is equal to the original $a_{1k}$ element. The other coefficients are" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "3b765484", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -921,14 +1165,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "9383783c", + "metadata": { + "editable": true + }, "source": [ "with a new right-hand side given by" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "c0a3994f", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -943,13 +1193,14 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "b5fab15e", + "metadata": { + "editable": true + }, "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" @@ -957,7 +1208,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "4d2d01dd", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -972,7 +1226,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8c1e5ff9", + "metadata": { + "editable": true + }, "source": [ "with $m=1,\\dots,n-1$ and a\n", "right-hand side given by" @@ -980,7 +1237,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "7165e693", + "metadata": { + "editable": true + }, "source": [ "\n", "
\n", @@ -995,7 +1255,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2c6646d7", + "metadata": { + "editable": true + }, "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", @@ -1007,9 +1270,6 @@ "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", @@ -1029,7 +1289,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "cc3c0cf7", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\begin{bmatrix}\n", @@ -1055,7 +1318,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "7256c8d3", + "metadata": { + "editable": true + }, "source": [ "LU decomposition forms the backbone of other algorithms in linear algebra, such as the\n", "solution of linear equations given by" @@ -1063,7 +1329,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "634741e9", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \\nonumber\n", @@ -1072,7 +1341,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8a16abbb", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \\nonumber\n", @@ -1081,7 +1353,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "442b2dd1", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \\nonumber\n", @@ -1090,7 +1365,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2be3e530", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \\nonumber\n", @@ -1099,7 +1377,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8f743387", + "metadata": { + "editable": true + }, "source": [ "The above set of equations is conveniently solved by using LU decomposition as an intermediate step.\n", "\n", @@ -1110,7 +1391,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "0fe06953", + "metadata": { + "editable": true + }, "source": [ "$$\n", "det\\{\\mathbf{A}\\}=det\\{\\mathbf{LU}\\}= det\\{\\mathbf{L}\\}det\\{\\mathbf{U}\\}=u_{11}u_{22}\\dots u_{nn}.\n", @@ -1119,7 +1403,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8232172f", + "metadata": { + "editable": true + }, "source": [ "There are at least three main advantages with LU decomposition compared with standard Gaussian elimination:\n", "\n", @@ -1135,7 +1422,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "87c2e89a", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \\nonumber\n", @@ -1144,7 +1434,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "50fc7a76", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \\nonumber\n", @@ -1153,7 +1446,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "e0379f33", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \\nonumber\n", @@ -1162,7 +1458,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "11e67b7e", + "metadata": { + "editable": true + }, "source": [ "$$\n", "a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \\nonumber\n", @@ -1171,14 +1470,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8cda9070", + "metadata": { + "editable": true + }, "source": [ "This can be written in matrix form as" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "908f05f7", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{Ax}=\\mathbf{w}.\n", @@ -1187,7 +1492,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "53053832", + "metadata": { + "editable": true + }, "source": [ "where $\\mathbf{A}$ and $\\mathbf{w}$ are known and we have to solve for\n", "$\\mathbf{x}$. Using the LU dcomposition we write" @@ -1195,7 +1503,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "f3314fa8", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{A} \\mathbf{x} \\equiv \\mathbf{L} \\mathbf{U} \\mathbf{x} =\\mathbf{w}.\n", @@ -1204,14 +1515,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2043fafc", + "metadata": { + "editable": true + }, "source": [ "The previous equation can be calculated in two steps" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "6fff8406", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{L} \\mathbf{y} = \\mathbf{w};\\qquad \\mathbf{Ux}=\\mathbf{y}.\n", @@ -1220,7 +1537,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "9a60901d", + "metadata": { + "editable": true + }, "source": [ "To show that this is correct we use to the LU decomposition\n", "to rewrite our system of linear equations as" @@ -1228,7 +1548,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "daacf54d", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{LUx}=\\mathbf{w},\n", @@ -1237,7 +1560,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "da702290", + "metadata": { + "editable": true + }, "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", @@ -1246,7 +1572,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "c2a04b00", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{Ux}=\\mathbf{L^{-1}w}=\\mathbf{y},\n", @@ -1255,14 +1584,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "d1214f7f", + "metadata": { + "editable": true + }, "source": [ "which yields the intermediate step" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "e8572956", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{L^{-1}w}=\\mathbf{y}\n", @@ -1271,18 +1606,23 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "6920b812", + "metadata": { + "editable": true + }, "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": {}, + "id": "379ca640", + "metadata": { + "editable": true + }, "source": [ "$$\n", "y_1=w_1 \\nonumber\n", @@ -1291,7 +1631,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8c753134", + "metadata": { + "editable": true + }, "source": [ "$$\n", "l_{21}y_1 + y_2=w_2\\nonumber\n", @@ -1300,7 +1643,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "12e32410", + "metadata": { + "editable": true + }, "source": [ "$$\n", "l_{31}y_1 + l_{32}y_2 + y_3 =w_3\\nonumber\n", @@ -1309,7 +1655,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "01e9e4ae", + "metadata": { + "editable": true + }, "source": [ "$$\n", "l_{41}y_1 + l_{42}y_2 + l_{43}y_3 + y_4=w_4. \\nonumber\n", @@ -1318,14 +1667,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "a27d4632", + "metadata": { + "editable": true + }, "source": [ "and" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "4dae7cfb", + "metadata": { + "editable": true + }, "source": [ "$$\n", "u_{11}x_1 +u_{12}x_2 +u_{13}x_3 + u_{14}x_4=y_1 \\nonumber\n", @@ -1334,7 +1689,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "5241a055", + "metadata": { + "editable": true + }, "source": [ "$$\n", "u_{22}x_2 + u_{23}x_3 + u_{24}x_4=y_2\\nonumber\n", @@ -1343,7 +1701,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "f590f30b", + "metadata": { + "editable": true + }, "source": [ "$$\n", "u_{33}x_3 + u_{34}x_4=y_3\\nonumber\n", @@ -1352,7 +1713,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "dbddd25a", + "metadata": { + "editable": true + }, "source": [ "$$\n", "u_{44}x_4=y_4 \\nonumber\n", @@ -1361,21 +1725,30 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8d5dc72b", + "metadata": { + "editable": true + }, "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", + " * 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." + ] + }, + { + "cell_type": "markdown", + "id": "fcf615b7", + "metadata": { + "editable": true + }, + "source": [ "### LU Decomposition, the inverse of a matrix\n", "\n", "If the inverse exists then" @@ -1383,7 +1756,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "3e2be427", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{A}^{-1}\\mathbf{A}=\\mathbf{I},\n", @@ -1392,14 +1768,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "519b78bb", + "metadata": { + "editable": true + }, "source": [ "the identity matrix. With an LU decomposed matrix we can rewrite the last equation as" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "a730fd5f", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{LU}\\mathbf{A}^{-1}=\\mathbf{I}.\n", @@ -1408,7 +1790,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "2cc4e115", + "metadata": { + "editable": true + }, "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" @@ -1416,7 +1801,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "03c4a0ae", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{A}_1^{-1}= \\begin{bmatrix}\n", @@ -1430,14 +1818,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "93ef84cc", + "metadata": { + "editable": true + }, "source": [ "then we have a linear set of equations" ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "04c367a7", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{LU}\\begin{bmatrix}\n", @@ -1456,14 +1850,20 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "8f903ad8", + "metadata": { + "editable": true + }, "source": [ "In a similar way we can compute the unknow entries of the second column," ] }, { "cell_type": "markdown", - "metadata": {}, + "id": "c343638b", + "metadata": { + "editable": true + }, "source": [ "$$\n", "\\mathbf{LU}\\begin{bmatrix}\n", @@ -1482,13 +1882,23 @@ }, { "cell_type": "markdown", - "metadata": {}, + "id": "1c013532", + "metadata": { + "editable": true + }, "source": [ - "and continue till we have solved all $n$ sets of linear equations." + "and continue till we have solved all $n$ sets of linear equations.\n", + "\n", + "The calculation of the inverse here assumes that it actually\n", + "exists. In many machine learning applications there may be strong\n", + "linear dependencies among the various columns and/or rows. In our\n", + "discussions of linear regression we will dive into the mathematics of\n", + "the singular value decomposition, an algorithm which will allow us to calculate the so-called pseudo-inverse.\n", + "These details will be presented in our linear regression chapter." ] } ], "metadata": {}, "nbformat": 4, - "nbformat_minor": 4 + "nbformat_minor": 5 }