From bfc6e385237f827ef5766226db80c06fd90111ec Mon Sep 17 00:00:00 2001 From: mhjensen Date: Thu, 24 May 2018 17:41:49 -0400 Subject: [PATCH] update on lin algebra --- doc/pub/Linalg/html/Linalg-bs.html | 691 ++++++----- doc/pub/Linalg/html/Linalg-reveal.html | 582 +++++----- doc/pub/Linalg/html/Linalg-solarized.html | 623 +++++----- doc/pub/Linalg/html/Linalg.html | 623 +++++----- doc/pub/Linalg/ipynb/Linalg.ipynb | 1029 +++++++++-------- doc/pub/Linalg/ipynb/ipynb-Linalg-src.tar.gz | Bin 207 -> 208 bytes .../Linalg/pdf/Linalg-beamer-handouts2x3.pdf | Bin 375052 -> 365727 bytes doc/pub/Linalg/pdf/Linalg-beamer.pdf | Bin 359199 -> 357095 bytes doc/pub/Linalg/pdf/Linalg-minted.pdf | Bin 384441 -> 387781 bytes doc/src/Linalg/Linalg.do.txt | 403 ++++--- 10 files changed, 1957 insertions(+), 1994 deletions(-) diff --git a/doc/pub/Linalg/html/Linalg-bs.html b/doc/pub/Linalg/html/Linalg-bs.html index f908a1b08..d16de27a2 100644 --- a/doc/pub/Linalg/html/Linalg-bs.html +++ b/doc/pub/Linalg/html/Linalg-bs.html @@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source - + -Data analysis and Machine Learning Lectures: Linear Algebra methods +Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays @@ -40,7 +40,7 @@ Automatically generated HTML file from DocOnce source @@ -147,7 +131,7 @@ MathJax.Hub.Config({ - Data analysis and Machine Learning Lectures: Linear Algebra methods + Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays
@@ -155,56 +139,49 @@ MathJax.Hub.Config({
  • Contents
    • @@ -222,7 +199,7 @@ MathJax.Hub.Config({
    -

    Data analysis and Machine Learning Lectures: Linear Algebra methods

    +

    Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays

    @@ -238,7 +215,7 @@ MathJax.Hub.Config({

    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

    -

    May 22, 2018

    +

    May 24, 2018


    @@ -246,22 +223,30 @@ MathJax.Hub.Config({ -

    To do

    +

    Introduction

    +The aim of this set of lectures is to review some central linear algebra algorithms that we will need in our +data analysis part and in the construction of Machine Learning algorithms (ML). +This will allow us to introduce some central programming features of high-level languages like Python and +compiled languages like C++ and/or Fortran. -
      -
    • add material on python handling of matrices and vectors
      • -
    • keep c++ material?
      • -
    +

    +As discussed in the introductory notes, these series of lectures focuses both on using +central Python packages like tensorflow and scikit-learn as well +as writing your own codes for some central ML algorithms. The +latter can be written in a language of your choice, be it Python, Julia, R, +Rust, C++, Fortran etc. In order to avoid confusion however, in these lectures we will limit our +attention to Python, C++ and Fortran. - +

    +

    Important Matrix and vector handling packages

    -The Numerical Recipes codes have been rewritten in Fortran 90/95 and -C/C++ by us. The original source codes are taken from the widely used +There are several central software packages for linear algebra and eigenvalue problems. Several of the more +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 software package LAPACK, which follows two other popular packages -developed in the 1970s, namely EISPACK and LINPACK. +developed in the 1970s, namely EISPACK and LINPACK. We describe them shortly here.

    • LINPACK: package for linear equations and least square problems.
      • @@ -269,10 +254,25 @@ developed in the 1970s, namely EISPACK and LINPACK.
    • 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.
    -Add python material on linear algebra and array handling, text on numpy etc +When dealing with matrices and vectors a central issue is memory +handling and allocation. If our code is written in Python the way we +declare these objects and the way they are handled, interpreted and +used by say a linear algebra library, requires codes that interface +our Python program with such libraries. For Python programmers, +Numpy is by now the standard Python package for numerical arrays in +Python as well as the source of functions which act on these +arrays. These functions span from eigenvalue solvers to functions that +compute the mean value, variance or the covariance matrix. If you are +not familiar with how arrays are handled in say Python or compiled +languages like C++ and Fortran, the sections in this chapter may be +useful. For C++ programmer, Armadillo is widely used library for +linear algebra and eigenvalue problems. In addition it offers a +convenient way to handle and organize arrays. We discuss this library +as well. Before we proceed we believe it may be convenient to repeat some basic features of + matrices and vectors.

    - +

    Basic Matrix Features

    @@ -303,6 +303,8 @@ $$

    + +

    The inverse of a matrix is defined by $$ @@ -387,7 +389,209 @@ For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all

    -

    Matrix Handling in C/C++, Static and Dynamical allocation

    +

    Numpy and arrays

    +Numpy provides an easy way to handle arrays in Python. The standard way to import this library is as +

    + + +

    import numpy as np
+n = 10
+x = np.random.normal(size=n)
+print(x)
+
    +

    +Here we have defined a vector \( x \) with \( n=10 \) elements with its values given by the Normal distribution \( N(0,1) \). +Another alternative is to declare a vector as follows +

    + + +

    import numpy as np
+x = np.array([1, 2, 3])
+print(x)
+
    +

    +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++ +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4, 7, 8]))
+print(x)
+
    +

    +Here we have used Numpy's unary function \( np.log \). This function is +highly tuned to compute array elements since the code is vectorized +and does not require looping. We normaly recommend that you use the +Numpy intrinsic functions instead of the corresponding log function +from Python's math module. The looping is done explicitely by the +np.log function. The alternative, and slower way to compute the +logarithms of a vector would be to write + +

    + + +

    import numpy as np
+from math import log
+x = np.array([4, 7, 8])
+for i in range(0, len(x)):
+    x[i] = log(x[i])
+print(x)
+
    +

    +We note that our code is much longer already and we need to import the log function from the math module. +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4, 7, 8], dtype = np.float64))
+print(x)
+
    +

    +or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is +

    + + +

    import numpy as np
+x = np.log(np.array([4.0, 7.0, 8.0])
+print(x)
+
    +

    +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4.0, 7.0, 8.0])
+print(x.itemsize)
+
    +

    + + +

    Matrices in Python

    +Having defined vectors, we are now ready to try out matrices. We can define a \( 3 \times 3 \) real matrix \( \hat{A} \) +as (recall that we user lowercase letters for vectors and uppercase letters for matrices) +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+print(A)
+
    +

    +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 +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+# print the first column, row-major order and elements start with 0
+print(A[:,0]) 
+
    +

    +We can continue this was by printing out other columns or rows. The example here prints out the second column +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+# print the first column, row-major order and elements start with 0
+print(A[1,:]) 
+
    +

    +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. 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 +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to zero
+A = np.zeros( (n, n) )
+print(A) 
+
    +

    +or initializing all elements to +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to one
+A = np.ones( (n, n) )
+print(A) 
+
    +

    +or as unitarily distributed random numbers (see the material on random number generators in the statistics part) +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]
+A = np.random.rand(n, n)
+print(A) 
+
    +

    +As we will see throughout these lectures, there are several extremely useful functionalities in Numpy. +As an example, consider the discussion of the covariance matrix. Suppose we have defined three vectors +\( \hat{x}, \hat{y}, \hat{z} \) with \( n \) elements each. The covariance matrix is defined as +$$ +\hat{\Sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ + \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ + \sigma_{zx} & \sigma_{zy} & \sigma_{zz} + \end{bmatrix}, +$$ + +where for example +$$ +\sigma_{xy} =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}). +$$ + +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. +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} \) +$$ +\hat{W} = \begin{bmatrix} x_0 & y_0 & z_0 \\ + x_1 & y_1 & z_1 \\ + x_2 & y_2 & z_2 \\ + \dots & \dots & \dots \\ + x_{n-2} & y_{n-2} & z_{n-2} \\ + x_{n-1} & y_{n-1} & z_{n-1} + \end{bmatrix}, +$$ + +

    +which in turn is converted into into the \( 3 times 3 \) covariance matrix +\( \hat{\Sigma} \) via the Numpy function np.cov(). In our review of +statistical functions and quantities we will discuss more about the +meaning of the covariance matrix. Here we note that we can calculate +the mean value of each set of samples \( \hat{x} \) etc using the Numpy +function np.mean(x). We can also extract the eigenvalues of the +covariance matrix through the np.linalg.eig() function. + +

    + + +

    # Importing various packages
+import numpy as np
+
+n = 100
+x = np.random.normal(size=n)
+print(np.mean(x))
+y = 4+3*x+np.random.normal(size=n)
+print(np.mean(y))
+z = x**3+np.random.normal(size=n)
+print(np.mean(z))
+W = np.vstack((x, y, z))
+Sigma = np.cov(W)
+print(Sigma)
+Eigvals, Eigvecs = np.linalg.eig(Sigma)
+print(Eigvals)
+
    +

    + + +

    Matrix Handling in C/C++, Static and Dynamical allocation

    @@ -415,7 +619,7 @@ Note the way the matrix is organized, row-major order.

    -

    Matrix Handling in C/C++

    +

    Matrix Handling in C/C++

    @@ -445,7 +649,7 @@ In C/C++ this would be coded like

    -

    Matrix Handling in C/C++

    +

    Matrix Handling in C/C++

    @@ -476,7 +680,7 @@ In C/C++ this would be coded like

    -

    Dynamic memory allocation in C/C++

    +

    Dynamic memory allocation in C/C++

    At least three possibilities in this course @@ -489,7 +693,7 @@ At least three possibilities in this course -

    Matrix Handling in C/C++, Dynamic Allocation

    +

    Matrix Handling in C/C++, Dynamic Allocation

    @@ -522,7 +726,7 @@ Always free space when you don't need an array anymore.

    -

    Armadillo, recommended!!

    +

    Armadillo, recommended!!

    • Armadillo is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use. The syntax is deliberately similar to Matlab.
      • @@ -534,7 +738,7 @@ Always free space when you don't need an array anymore. -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -557,7 +761,7 @@ Always free space when you don't need an array anymore.

      -

      Armadillo, how to compile and install

      +

      Armadillo, how to compile and install

      For people using Ubuntu, Debian, Linux Mint, simply go to the synaptic package manager and install @@ -585,7 +789,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -618,7 +822,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -649,7 +853,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -682,7 +886,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -710,7 +914,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -742,7 +946,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -772,7 +976,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -804,7 +1008,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      We start with the linear set of equations @@ -839,7 +1043,7 @@ $$

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      or $$ @@ -854,7 +1058,7 @@ $$

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown \( x_1 \) @@ -878,7 +1082,7 @@ what is called a backward substitution.

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      This process can be expressed mathematically as $$ @@ -896,7 +1100,7 @@ the result from the $j$th equation. We assume obviously that \( a_{11}\ne 0 \) a

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      Our actual \( 4\times 4 \) example reads after the first operation @@ -936,7 +1140,7 @@ $$

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      The new coefficients are @@ -972,7 +1176,7 @@ We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \(

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      This step is called forward substitution. @@ -1010,7 +1214,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).

      -

      Linear Algebra Methods

      +

      Linear Algebra Methods

      • Gaussian elimination, \( O(2/3n^3) \) flops, general matrix
        • @@ -1023,7 +1227,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \). -

        LU Decomposition

        +

        LU Decomposition

        The LU decomposition method means that we can rewrite @@ -1054,7 +1258,7 @@ $$

        -

        LU Decomposition

        +

        LU Decomposition

        LU decomposition forms the backbone of other algorithms in linear algebra, such as the @@ -1083,7 +1287,7 @@ $$

        -

        LU Decomposition, why?

        +

        LU Decomposition, why?

        There are at least three main advantages with LU decomposition compared with standard Gaussian elimination: @@ -1096,7 +1300,7 @@ There are at least three main advantages with LU decomposition compared with sta -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        With the LU decomposition it is rather @@ -1125,7 +1329,7 @@ $$ \mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}. $$

        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        The previous equation can be calculated in two steps @@ -1138,7 +1342,7 @@ to rewrite our system of linear equations as $$ \mathbf{LUx}=\mathbf{w}, $$ -and since the determinat of \( \mathbf{L} \) is equal to 1 (by construction +and since the determinant of \( \mathbf{L} \) is equal to 1 (by construction since the diagonals of \( \mathbf{L} \) equal 1) we can use the inverse of \( \mathbf{L} \) to obtain @@ -1158,7 +1362,7 @@ through \( \mathbf{Ux}=\mathbf{y} \).

        -

        LU Decomposition, why?

        +

        LU Decomposition, why?

        For our four-dimentional example this takes the form @@ -1191,7 +1395,7 @@ needed to solve the set of \( n \) linear equations.

        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        The algorithm goes as follows @@ -1204,7 +1408,7 @@ The algorithm goes as follows -

        LU Decomposition, the inverse of a matrix

        +

        LU Decomposition, the inverse of a matrix

        If the inverse exists then @@ -1222,7 +1426,7 @@ $$

        -

        LU Decomposition, the inverse of a matrix

        +

        LU Decomposition, the inverse of a matrix

        If we assume that the first column (that is column 1) of the inverse matrix @@ -1256,7 +1460,7 @@ $$

        -

        LU Decomposition, the inverse

        +

        LU Decomposition, the inverse

        In a similar way we can compute the unknow entries of the second column, @@ -1280,7 +1484,7 @@ and continue till we have solved all \( n \) sets of linear equations.

        -

        Using Armadillo to perform an LU decomposition

        +

        Using Armadillo to perform an LU decomposition

        @@ -1312,245 +1516,6 @@ and continue till we have solved all \( n \) sets of linear equations. return 0; }

    -

    - - -

    Iterative methods, Chapter 6

    -
    -
    -

    - -

      -
    • Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.
      • -
    • Iterative solvers such as Basic iterative solvers, Jacobi, Gauss-Seidel, Successive over-relaxation. These methods are easy to parallelize, as we will se later. Much used in solutions of partial differential equations.
      • -
    • Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.
      • -
    -
    -
    - - -

    - - -

    Iterative methods, Jacobi's method

    -
    -
    -

    -It is a simple method for solving -$$ -\mathbf{A}\mathbf{x}=\mathbf{b}, -$$ - -where \( \mathbf{A} \) is a matrix and \( \mathbf{x} \) and \( \mathbf{b} \) are vectors. The vector \( \mathbf{x} \) is -the unknown. - -

    -It is an iterative scheme where we start with a guess for the unknown, and -after \( k+1 \) iterations we have -$$ -\mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}), -$$ - -with \( \mathbf{A}=\mathbf{D}+\mathbf{U}+\mathbf{L} \) and -\( \mathbf{D} \) being a diagonal matrix, \( \mathbf{U} \) an upper triangular matrix and \( \mathbf{L} \) a lower triangular -matrix. - -

    -If the matrix \( \mathbf{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution. -

    -
    - - -

    - - -

    Iterative methods, Jacobi's method

    -
    -
    -

    -We can demonstrate Jacobi's method by this \( 4\times 4 \) matrix problem. We assume a guess -for the vector elements \( x_i^{(0)} \), a guess which represents our first iteration. The new -values are obtained by substitution -$$ -\begin{align} - x_1^{(1)} =&(b_1-a_{12}x_2^{(0)} -a_{13}x_3^{(0)} - a_{14}x_4^{(0)})/a_{11} \nonumber \\ - x_2^{(1)} =&(b_2-a_{21}x_1^{(0)} - a_{23}x_3^{(0)} - a_{24}x_4^{(0)})/a_{22} \nonumber \\ - x_3^{(1)} =&(b_3- a_{31}x_1^{(0)} -a_{32}x_2^{(0)} -a_{34}x_4^{(0)})/a_{33} \nonumber \\ - x_4^{(1)}=&(b_4-a_{41}x_1^{(0)} -a_{42}x_2^{(0)} - a_{43}x_3^{(0)})/a_{44}, \nonumber -\end{align} -$$ - -which after \( k+1 \) iterations reads -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber -\end{align} -$$ -

    -
    - - -

    - - -

    Iterative methods, Jacobi's method

    -
    -
    -

    -We can generalize the above equations to -$$ - x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii} -$$ - -or in an even more compact form as -$$ \mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}), -$$ - -with \( \mathbf{A}=\mathbf{D}+\mathbf{U}+\mathbf{L} \) and -\( \mathbf{D} \) being a diagonal matrix, \( \mathbf{U} \) an upper triangular matrix and \( \mathbf{L} \) a lower triangular -matrix. -

    -
    - - -

    - - -

    Iterative methods, Gauss-Seidel's method

    -
    -
    -

    -Our \( 4\times 4 \) matrix problem -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber -\end{align} -$$ - -can be rewritten as -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber -\end{align} -$$ - -which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence -behavior and leads to the Gauss-Seidel method! -

    -
    - - -

    - - -

    Iterative methods, Gauss-Seidel's method

    -
    -
    -

    -We can generalize -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber -\end{align} -$$ - -to the following form -$$ - x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j > i}a_{ij}x^{(k)}_j - \sum_{j < i}a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. -$$ - -The procedure is generally continued until the changes made by an iteration are below some tolerance. - -

    -The convergence properties of the Jacobi method and the -Gauss-Seidel method are dependent on the matrix \( \mathbf{A} \). These methods converge when -the matrix is symmetric positive-definite, or is strictly or irreducibly diagonally dominant. -Both methods sometimes converge even if these conditions are not satisfied. -

    -
    - - -

    - - -

    Iterative methods, Successive over-relaxation

    -
    -
    -

    -Given a square system of n linear equations with unknown \( \mathbf x \): -$$ - \mathbf{A}\mathbf x = \mathbf b -$$ - -where -$$ - \mathbf{A}=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}. -$$ -

    -
    - - -

    - - -

    Iterative methods, Successive over-relaxation

    -
    -
    -

    -Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U: -$$ - \mathbf{A} =\mathbf{D} + \mathbf{L} + \mathbf{U}, -$$ - -where -$$ - D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \quad L = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix}, \quad U = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & 0 \end{bmatrix}. -$$ - -The system of linear equations may be rewritten as: -$$ - (D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x} -$$ - -for a constant \( \omega > 1 \). -

    -
    - - -

    - - -

    Iterative methods, Successive over-relaxation

    -
    -
    -

    -The method of successive over-relaxation is an iterative technique that solves the left hand side of this expression for \( x \), using previous value for \( x \) on the right hand side. Analytically, this may be written as: -$$ - \mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big). -$$ - -However, by taking advantage of the triangular form of \( (D+\omega L) \), the elements of \( x^{(k+1)} \) can be computed sequentially using forward substitution: -$$ - x^{(k+1)}_i = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j > i} a_{ij}x^{(k)}_j - \sum_{j < i} a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. -$$ - -The choice of relaxation factor is not necessarily easy, and depends upon the properties of the coefficient matrix. For symmetric, positive-definite matrices it can be proven that \( 0 < \omega < 2 \) will lead to convergence, but we are generally interested in faster convergence rather than just convergence. -

    -
    - -

    diff --git a/doc/pub/Linalg/html/Linalg-reveal.html b/doc/pub/Linalg/html/Linalg-reveal.html index c20710a77..1e36e4cc0 100644 --- a/doc/pub/Linalg/html/Linalg-reveal.html +++ b/doc/pub/Linalg/html/Linalg-reveal.html @@ -3,9 +3,9 @@ - + -Data analysis and Machine Learning Lectures: Linear Algebra methods +Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays @@ -132,7 +132,7 @@ MathJax.Hub.Config({ -

    Data analysis and Machine Learning Lectures: Linear Algebra methods

    +

    Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays

    @@ -148,7 +148,7 @@ MathJax.Hub.Config({

    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

     
    -

    May 22, 2018

    +

    May 24, 2018


    @@ -159,12 +159,19 @@ MathJax.Hub.Config({

    -

    To do

    +

    Introduction

    +The aim of this set of lectures is to review some central linear algebra algorithms that we will need in our +data analysis part and in the construction of Machine Learning algorithms (ML). +This will allow us to introduce some central programming features of high-level languages like Python and +compiled languages like C++ and/or Fortran. -
      -

    • add material on python handling of matrices and vectors
      • -

    • keep c++ material?
      • -
    +

    +As discussed in the introductory notes, these series of lectures focuses both on using +central Python packages like tensorflow and scikit-learn as well +as writing your own codes for some central ML algorithms. The +latter can be written in a language of your choice, be it Python, Julia, R, +Rust, C++, Fortran etc. In order to avoid confusion however, in these lectures we will limit our +attention to Python, C++ and Fortran.

    @@ -172,10 +179,10 @@ MathJax.Hub.Config({

    Important Matrix and vector handling packages

    -The Numerical Recipes codes have been rewritten in Fortran 90/95 and -C/C++ by us. The original source codes are taken from the widely used +There are several central software packages for linear algebra and eigenvalue problems. Several of the more +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 software package LAPACK, which follows two other popular packages -developed in the 1970s, namely EISPACK and LINPACK. +developed in the 1970s, namely EISPACK and LINPACK. We describe them shortly here.

      @@ -187,7 +194,22 @@ developed in the 1970s, namely EISPACK and LINPACK.

    -Add python material on linear algebra and array handling, text on numpy etc +When dealing with matrices and vectors a central issue is memory +handling and allocation. If our code is written in Python the way we +declare these objects and the way they are handled, interpreted and +used by say a linear algebra library, requires codes that interface +our Python program with such libraries. For Python programmers, +Numpy is by now the standard Python package for numerical arrays in +Python as well as the source of functions which act on these +arrays. These functions span from eigenvalue solvers to functions that +compute the mean value, variance or the covariance matrix. If you are +not familiar with how arrays are handled in say Python or compiled +languages like C++ and Fortran, the sections in this chapter may be +useful. For C++ programmer, Armadillo is widely used library for +linear algebra and eigenvalue problems. In addition it offers a +convenient way to handle and organize arrays. We discuss this library +as well. Before we proceed we believe it may be convenient to repeat some basic features of + matrices and vectors. @@ -311,7 +333,217 @@ For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all

    -

    Matrix Handling in C/C++, Static and Dynamical allocation

    +

    Numpy and arrays

    +Numpy provides an easy way to handle arrays in Python. The standard way to import this library is as +

    + + +

    import numpy as np
+n = 10
+x = np.random.normal(size=n)
+print(x)
+
    +

    +Here we have defined a vector \( x \) with \( n=10 \) elements with its values given by the Normal distribution \( N(0,1) \). +Another alternative is to declare a vector as follows +

    + + +

    import numpy as np
+x = np.array([1, 2, 3])
+print(x)
+
    +

    +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++ +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4, 7, 8]))
+print(x)
+
    +

    +Here we have used Numpy's unary function \( np.log \). This function is +highly tuned to compute array elements since the code is vectorized +and does not require looping. We normaly recommend that you use the +Numpy intrinsic functions instead of the corresponding log function +from Python's math module. The looping is done explicitely by the +np.log function. The alternative, and slower way to compute the +logarithms of a vector would be to write + +

    + + +

    import numpy as np
+from math import log
+x = np.array([4, 7, 8])
+for i in range(0, len(x)):
+    x[i] = log(x[i])
+print(x)
+
    +

    +We note that our code is much longer already and we need to import the log function from the math module. +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4, 7, 8], dtype = np.float64))
+print(x)
+
    +

    +or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is +

    + + +

    import numpy as np
+x = np.log(np.array([4.0, 7.0, 8.0])
+print(x)
+
    +

    +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4.0, 7.0, 8.0])
+print(x.itemsize)
+
    +
    + + +
    +

    Matrices in Python

    +Having defined vectors, we are now ready to try out matrices. We can define a \( 3 \times 3 \) real matrix \( \hat{A} \) +as (recall that we user lowercase letters for vectors and uppercase letters for matrices) +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+print(A)
+
    +

    +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 +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+# print the first column, row-major order and elements start with 0
+print(A[:,0]) 
+
    +

    +We can continue this was by printing out other columns or rows. The example here prints out the second column +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+# print the first column, row-major order and elements start with 0
+print(A[1,:]) 
+
    +

    +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. 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 +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to zero
+A = np.zeros( (n, n) )
+print(A) 
+
    +

    +or initializing all elements to +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to one
+A = np.ones( (n, n) )
+print(A) 
+
    +

    +or as unitarily distributed random numbers (see the material on random number generators in the statistics part) +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]
+A = np.random.rand(n, n)
+print(A) 
+
    +

    +As we will see throughout these lectures, there are several extremely useful functionalities in Numpy. +As an example, consider the discussion of the covariance matrix. Suppose we have defined three vectors +\( \hat{x}, \hat{y}, \hat{z} \) with \( n \) elements each. The covariance matrix is defined as +

     
    +$$ +\hat{\Sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ + \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ + \sigma_{zx} & \sigma_{zy} & \sigma_{zz} + \end{bmatrix}, +$$ +

     
    + +where for example +

     
    +$$ +\sigma_{xy} =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}). +$$ +

     
    + +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. +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} \) +

     
    +$$ +\hat{W} = \begin{bmatrix} x_0 & y_0 & z_0 \\ + x_1 & y_1 & z_1 \\ + x_2 & y_2 & z_2 \\ + \dots & \dots & \dots \\ + x_{n-2} & y_{n-2} & z_{n-2} \\ + x_{n-1} & y_{n-1} & z_{n-1} + \end{bmatrix}, +$$ +

     
    + +

    +which in turn is converted into into the \( 3 times 3 \) covariance matrix +\( \hat{\Sigma} \) via the Numpy function np.cov(). In our review of +statistical functions and quantities we will discuss more about the +meaning of the covariance matrix. Here we note that we can calculate +the mean value of each set of samples \( \hat{x} \) etc using the Numpy +function np.mean(x). We can also extract the eigenvalues of the +covariance matrix through the np.linalg.eig() function. + +

    + + +

    # Importing various packages
+import numpy as np
+
+n = 100
+x = np.random.normal(size=n)
+print(np.mean(x))
+y = 4+3*x+np.random.normal(size=n)
+print(np.mean(y))
+z = x**3+np.random.normal(size=n)
+print(np.mean(z))
+W = np.vstack((x, y, z))
+Sigma = np.cov(W)
+print(Sigma)
+Eigvals, Eigvecs = np.linalg.eig(Sigma)
+print(Eigvals)
+
    +
    + + +
    +

    Matrix Handling in C/C++, Static and Dynamical allocation

    @@ -337,7 +569,7 @@ Note the way the matrix is organized, row-major order.
    -

    Matrix Handling in C/C++

    +

    Matrix Handling in C/C++

    @@ -367,7 +599,7 @@ In C/C++ this would be coded like
    -

    Matrix Handling in C/C++

    +

    Matrix Handling in C/C++

    @@ -398,7 +630,7 @@ In C/C++ this would be coded like
    -

    Dynamic memory allocation in C/C++

    +

    Dynamic memory allocation in C/C++

    At least three possibilities in this course @@ -415,7 +647,7 @@ At least three possibilities in this course

    -

    Matrix Handling in C/C++, Dynamic Allocation

    +

    Matrix Handling in C/C++, Dynamic Allocation

    @@ -445,7 +677,7 @@ Always free space when you don't need an array anymore.
    -

    Armadillo, recommended!!

    +

    Armadillo, recommended!!

      @@ -463,7 +695,7 @@ Always free space when you don't need an array anymore.
      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -487,7 +719,7 @@ Always free space when you don't need an array anymore.

      -

      Armadillo, how to compile and install

      +

      Armadillo, how to compile and install

      For people using Ubuntu, Debian, Linux Mint, simply go to the synaptic package manager and install @@ -516,7 +748,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -550,7 +782,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -582,7 +814,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -616,7 +848,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -645,7 +877,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -678,7 +910,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -709,7 +941,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -742,7 +974,7 @@ For OS X users you may have to declare the paths to the include files and the li

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      We start with the linear set of equations @@ -781,7 +1013,7 @@ $$

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      or

       
      @@ -798,7 +1030,7 @@ $$

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown \( x_1 \) @@ -824,7 +1056,7 @@ what is called a backward substitution.

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      This process can be expressed mathematically as

       
      @@ -844,7 +1076,7 @@ the result from the $j$th equation. We assume obviously that \( a_{11}\ne 0 \) a

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      Our actual \( 4\times 4 \) example reads after the first operation @@ -888,7 +1120,7 @@ $$

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      The new coefficients are @@ -930,7 +1162,7 @@ We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \(

      -

      Gaussian Elimination

      +

      Gaussian Elimination

      This step is called forward substitution. @@ -972,7 +1204,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).

      -

      Linear Algebra Methods

      +

      Linear Algebra Methods

        @@ -992,7 +1224,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
        -

        LU Decomposition

        +

        LU Decomposition

        The LU decomposition method means that we can rewrite @@ -1025,7 +1257,7 @@ $$

        -

        LU Decomposition

        +

        LU Decomposition

        LU decomposition forms the backbone of other algorithms in linear algebra, such as the @@ -1058,7 +1290,7 @@ $$

        -

        LU Decomposition, why?

        +

        LU Decomposition, why?

        There are at least three main advantages with LU decomposition compared with standard Gaussian elimination: @@ -1075,7 +1307,7 @@ There are at least three main advantages with LU decomposition compared with sta

        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        With the LU decomposition it is rather @@ -1110,7 +1342,7 @@ $$ \mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}. $$

        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        The previous equation can be calculated in two steps @@ -1127,7 +1359,7 @@ to rewrite our system of linear equations as $$ \mathbf{LUx}=\mathbf{w}, $$

         
        -and since the determinat of \( \mathbf{L} \) is equal to 1 (by construction +and since the determinant of \( \mathbf{L} \) is equal to 1 (by construction since the diagonals of \( \mathbf{L} \) equal 1) we can use the inverse of \( \mathbf{L} \) to obtain @@ -1151,7 +1383,7 @@ through \( \mathbf{Ux}=\mathbf{y} \).

        -

        LU Decomposition, why?

        +

        LU Decomposition, why?

        For our four-dimentional example this takes the form @@ -1188,7 +1420,7 @@ needed to solve the set of \( n \) linear equations.

        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        The algorithm goes as follows @@ -1205,7 +1437,7 @@ The algorithm goes as follows

        -

        LU Decomposition, the inverse of a matrix

        +

        LU Decomposition, the inverse of a matrix

        If the inverse exists then @@ -1227,7 +1459,7 @@ $$

        -

        LU Decomposition, the inverse of a matrix

        +

        LU Decomposition, the inverse of a matrix

        If we assume that the first column (that is column 1) of the inverse matrix @@ -1265,7 +1497,7 @@ $$

        -

        LU Decomposition, the inverse

        +

        LU Decomposition, the inverse

        In a similar way we can compute the unknow entries of the second column, @@ -1291,7 +1523,7 @@ and continue till we have solved all \( n \) sets of linear equations.

        -

        Using Armadillo to perform an LU decomposition

        +

        Using Armadillo to perform an LU decomposition

        @@ -1326,260 +1558,6 @@ and continue till we have solved all \( n \) sets of linear equations.

        -
        -

        Iterative methods, Chapter 6

        -
        - -
          -

        • Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.
          • -

        • Iterative solvers such as Basic iterative solvers, Jacobi, Gauss-Seidel, Successive over-relaxation. These methods are easy to parallelize, as we will se later. Much used in solutions of partial differential equations.
          • -

        • Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.
          • -
        -
        -
        - - -
        -

        Iterative methods, Jacobi's method

        -
        - -

        -It is a simple method for solving -

         
        -$$ -\mathbf{A}\mathbf{x}=\mathbf{b}, -$$ -

         
        - -where \( \mathbf{A} \) is a matrix and \( \mathbf{x} \) and \( \mathbf{b} \) are vectors. The vector \( \mathbf{x} \) is -the unknown. - -

        -It is an iterative scheme where we start with a guess for the unknown, and -after \( k+1 \) iterations we have - -

         
        -$$ -\mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}), -$$ -

         
        - -with \( \mathbf{A}=\mathbf{D}+\mathbf{U}+\mathbf{L} \) and -\( \mathbf{D} \) being a diagonal matrix, \( \mathbf{U} \) an upper triangular matrix and \( \mathbf{L} \) a lower triangular -matrix. - -

        -If the matrix \( \mathbf{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution. -

        -
        - - -
        -

        Iterative methods, Jacobi's method

        -
        - -

        -We can demonstrate Jacobi's method by this \( 4\times 4 \) matrix problem. We assume a guess -for the vector elements \( x_i^{(0)} \), a guess which represents our first iteration. The new -values are obtained by substitution -

         
        -$$ -\begin{align} - x_1^{(1)} =&(b_1-a_{12}x_2^{(0)} -a_{13}x_3^{(0)} - a_{14}x_4^{(0)})/a_{11} \nonumber \\ - x_2^{(1)} =&(b_2-a_{21}x_1^{(0)} - a_{23}x_3^{(0)} - a_{24}x_4^{(0)})/a_{22} \nonumber \\ - x_3^{(1)} =&(b_3- a_{31}x_1^{(0)} -a_{32}x_2^{(0)} -a_{34}x_4^{(0)})/a_{33} \nonumber \\ - x_4^{(1)}=&(b_4-a_{41}x_1^{(0)} -a_{42}x_2^{(0)} - a_{43}x_3^{(0)})/a_{44}, \nonumber -\end{align} -$$ -

         
        - -which after \( k+1 \) iterations reads -

         
        -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber -\end{align} -$$ -

         
        -

        -
        - - -
        -

        Iterative methods, Jacobi's method

        -
        - -

        -We can generalize the above equations to -

         
        -$$ - x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii} -$$ -

         
        - -or in an even more compact form as -

         
        -$$ \mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}), -$$ -

         
        - -with \( \mathbf{A}=\mathbf{D}+\mathbf{U}+\mathbf{L} \) and -\( \mathbf{D} \) being a diagonal matrix, \( \mathbf{U} \) an upper triangular matrix and \( \mathbf{L} \) a lower triangular -matrix. -

        -
        - - -
        -

        Iterative methods, Gauss-Seidel's method

        -
        - -

        -Our \( 4\times 4 \) matrix problem -

         
        -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber -\end{align} -$$ -

         
        - -can be rewritten as -

         
        -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber -\end{align} -$$ -

         
        - -which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence -behavior and leads to the Gauss-Seidel method! -

        -
        - - -
        -

        Iterative methods, Gauss-Seidel's method

        -
        - -

        -We can generalize -

         
        -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber -\end{align} -$$ -

         
        - -to the following form -

         
        -$$ - x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j > i}a_{ij}x^{(k)}_j - \sum_{j < i}a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. -$$ -

         
        - -The procedure is generally continued until the changes made by an iteration are below some tolerance. - -

        -The convergence properties of the Jacobi method and the -Gauss-Seidel method are dependent on the matrix \( \mathbf{A} \). These methods converge when -the matrix is symmetric positive-definite, or is strictly or irreducibly diagonally dominant. -Both methods sometimes converge even if these conditions are not satisfied. -

        -
        - - -
        -

        Iterative methods, Successive over-relaxation

        -
        - -

        -Given a square system of n linear equations with unknown \( \mathbf x \): -

         
        -$$ - \mathbf{A}\mathbf x = \mathbf b -$$ -

         
        - -where -

         
        -$$ - \mathbf{A}=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}. -$$ -

         
        -

        -
        - - -
        -

        Iterative methods, Successive over-relaxation

        -
        - -

        -Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U: -

         
        -$$ - \mathbf{A} =\mathbf{D} + \mathbf{L} + \mathbf{U}, -$$ -

         
        - -where -

         
        -$$ - D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \quad L = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix}, \quad U = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & 0 \end{bmatrix}. -$$ -

         
        - -The system of linear equations may be rewritten as: -

         
        -$$ - (D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x} -$$ -

         
        - -for a constant \( \omega > 1 \). -

        -
        - - -
        -

        Iterative methods, Successive over-relaxation

        -
        - -

        -The method of successive over-relaxation is an iterative technique that solves the left hand side of this expression for \( x \), using previous value for \( x \) on the right hand side. Analytically, this may be written as: -

         
        -$$ - \mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big). -$$ -

         
        - -However, by taking advantage of the triangular form of \( (D+\omega L) \), the elements of \( x^{(k+1)} \) can be computed sequentially using forward substitution: -

         
        -$$ - x^{(k+1)}_i = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j > i} a_{ij}x^{(k)}_j - \sum_{j < i} a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. -$$ -

         
        - -The choice of relaxation factor is not necessarily easy, and depends upon the properties of the coefficient matrix. For symmetric, positive-definite matrices it can be proven that \( 0 < \omega < 2 \) will lead to convergence, but we are generally interested in faster convergence rather than just convergence. -

        -
        - -
    diff --git a/doc/pub/Linalg/html/Linalg-solarized.html b/doc/pub/Linalg/html/Linalg-solarized.html index a6e7673b5..de7a34d25 100644 --- a/doc/pub/Linalg/html/Linalg-solarized.html +++ b/doc/pub/Linalg/html/Linalg-solarized.html @@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source - + -Data analysis and Machine Learning Lectures: Linear Algebra methods +Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays @@ -60,7 +60,7 @@ div { text-align: justify; text-justify: inter-word; } @@ -163,7 +147,7 @@ MathJax.Hub.Config({ -

    Data analysis and Machine Learning Lectures: Linear Algebra methods

    +

    Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays

    @@ -179,27 +163,35 @@ MathJax.Hub.Config({

    [2] Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

    -

    May 22, 2018

    +

    May 24, 2018












    -

    To do

    +

    Introduction

    +The aim of this set of lectures is to review some central linear algebra algorithms that we will need in our +data analysis part and in the construction of Machine Learning algorithms (ML). +This will allow us to introduce some central programming features of high-level languages like Python and +compiled languages like C++ and/or Fortran. -
      -
    • add material on python handling of matrices and vectors
      • -
    • keep c++ material?
      • -
    +

    +As discussed in the introductory notes, these series of lectures focuses both on using +central Python packages like tensorflow and scikit-learn as well +as writing your own codes for some central ML algorithms. The +latter can be written in a language of your choice, be it Python, Julia, R, +Rust, C++, Fortran etc. In order to avoid confusion however, in these lectures we will limit our +attention to Python, C++ and Fortran. -









    +

    +

    Important Matrix and vector handling packages

    -The Numerical Recipes codes have been rewritten in Fortran 90/95 and -C/C++ by us. The original source codes are taken from the widely used +There are several central software packages for linear algebra and eigenvalue problems. Several of the more +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 software package LAPACK, which follows two other popular packages -developed in the 1970s, namely EISPACK and LINPACK. +developed in the 1970s, namely EISPACK and LINPACK. We describe them shortly here.

    • LINPACK: package for linear equations and least square problems.
      • @@ -207,10 +199,25 @@ developed in the 1970s, namely EISPACK and LINPACK.
    • 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.
    -Add python material on linear algebra and array handling, text on numpy etc +When dealing with matrices and vectors a central issue is memory +handling and allocation. If our code is written in Python the way we +declare these objects and the way they are handled, interpreted and +used by say a linear algebra library, requires codes that interface +our Python program with such libraries. For Python programmers, +Numpy is by now the standard Python package for numerical arrays in +Python as well as the source of functions which act on these +arrays. These functions span from eigenvalue solvers to functions that +compute the mean value, variance or the covariance matrix. If you are +not familiar with how arrays are handled in say Python or compiled +languages like C++ and Fortran, the sections in this chapter may be +useful. For C++ programmer, Armadillo is widely used library for +linear algebra and eigenvalue problems. In addition it offers a +convenient way to handle and organize arrays. We discuss this library +as well. Before we proceed we believe it may be convenient to repeat some basic features of + matrices and vectors.

    -









    +

    Basic Matrix Features

    @@ -239,6 +246,8 @@ $$

    Basic Matrix Features

    +

    +

    The inverse of a matrix is defined by @@ -316,7 +325,209 @@ For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all











    -

    Matrix Handling in C/C++, Static and Dynamical allocation

    +

    Numpy and arrays

    +Numpy provides an easy way to handle arrays in Python. The standard way to import this library is as +

    + + +

    import numpy as np
+n = 10
+x = np.random.normal(size=n)
+print(x)
+
    +

    +Here we have defined a vector \( x \) with \( n=10 \) elements with its values given by the Normal distribution \( N(0,1) \). +Another alternative is to declare a vector as follows +

    + + +

    import numpy as np
+x = np.array([1, 2, 3])
+print(x)
+
    +

    +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++ +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4, 7, 8]))
+print(x)
+
    +

    +Here we have used Numpy's unary function \( np.log \). This function is +highly tuned to compute array elements since the code is vectorized +and does not require looping. We normaly recommend that you use the +Numpy intrinsic functions instead of the corresponding log function +from Python's math module. The looping is done explicitely by the +np.log function. The alternative, and slower way to compute the +logarithms of a vector would be to write + +

    + + +

    import numpy as np
+from math import log
+x = np.array([4, 7, 8])
+for i in range(0, len(x)):
+    x[i] = log(x[i])
+print(x)
+
    +

    +We note that our code is much longer already and we need to import the log function from the math module. +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4, 7, 8], dtype = np.float64))
+print(x)
+
    +

    +or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is +

    + + +

    import numpy as np
+x = np.log(np.array([4.0, 7.0, 8.0])
+print(x)
+
    +

    +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 +

    + + +

    import numpy as np
+x = np.log(np.array([4.0, 7.0, 8.0])
+print(x.itemsize)
+
    +

    +









    + +

    Matrices in Python

    +Having defined vectors, we are now ready to try out matrices. We can define a \( 3 \times 3 \) real matrix \( \hat{A} \) +as (recall that we user lowercase letters for vectors and uppercase letters for matrices) +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+print(A)
+
    +

    +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 +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+# print the first column, row-major order and elements start with 0
+print(A[:,0]) 
+
    +

    +We can continue this was by printing out other columns or rows. The example here prints out the second column +

    + + +

    import numpy as np
+A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
+# print the first column, row-major order and elements start with 0
+print(A[1,:]) 
+
    +

    +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. 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 +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to zero
+A = np.zeros( (n, n) )
+print(A) 
+
    +

    +or initializing all elements to +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to one
+A = np.ones( (n, n) )
+print(A) 
+
    +

    +or as unitarily distributed random numbers (see the material on random number generators in the statistics part) +

    + + +

    import numpy as np
+n = 10
+# define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]
+A = np.random.rand(n, n)
+print(A) 
+
    +

    +As we will see throughout these lectures, there are several extremely useful functionalities in Numpy. +As an example, consider the discussion of the covariance matrix. Suppose we have defined three vectors +\( \hat{x}, \hat{y}, \hat{z} \) with \( n \) elements each. The covariance matrix is defined as +$$ +\hat{\Sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ + \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ + \sigma_{zx} & \sigma_{zy} & \sigma_{zz} + \end{bmatrix}, +$$ + +where for example +$$ +\sigma_{xy} =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}). +$$ + +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. +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} \) +$$ +\hat{W} = \begin{bmatrix} x_0 & y_0 & z_0 \\ + x_1 & y_1 & z_1 \\ + x_2 & y_2 & z_2 \\ + \dots & \dots & \dots \\ + x_{n-2} & y_{n-2} & z_{n-2} \\ + x_{n-1} & y_{n-1} & z_{n-1} + \end{bmatrix}, +$$ + +

    +which in turn is converted into into the \( 3 times 3 \) covariance matrix +\( \hat{\Sigma} \) via the Numpy function np.cov(). In our review of +statistical functions and quantities we will discuss more about the +meaning of the covariance matrix. Here we note that we can calculate +the mean value of each set of samples \( \hat{x} \) etc using the Numpy +function np.mean(x). We can also extract the eigenvalues of the +covariance matrix through the np.linalg.eig() function. + +

    + + +

    # Importing various packages
+import numpy as np
+
+n = 100
+x = np.random.normal(size=n)
+print(np.mean(x))
+y = 4+3*x+np.random.normal(size=n)
+print(np.mean(y))
+z = x**3+np.random.normal(size=n)
+print(np.mean(z))
+W = np.vstack((x, y, z))
+Sigma = np.cov(W)
+print(Sigma)
+Eigvals, Eigvecs = np.linalg.eig(Sigma)
+print(Eigvals)
+
    +

    +









    + +

    Matrix Handling in C/C++, Static and Dynamical allocation

    @@ -343,7 +554,7 @@ Note the way the matrix is organized, row-major order.











    -

    Matrix Handling in C/C++

    +

    Matrix Handling in C/C++

    @@ -372,7 +583,7 @@ In C/C++ this would be coded like











    -

    Matrix Handling in C/C++

    +

    Matrix Handling in C/C++

    @@ -402,7 +613,7 @@ In C/C++ this would be coded like











    -

    Dynamic memory allocation in C/C++

    +

    Dynamic memory allocation in C/C++

    At least three possibilities in this course @@ -415,7 +626,7 @@ At least three possibilities in this course









    -

    Matrix Handling in C/C++, Dynamic Allocation

    +

    Matrix Handling in C/C++, Dynamic Allocation

    @@ -447,7 +658,7 @@ Always free space when you don't need an array anymore.











    -

    Armadillo, recommended!!

    +

    Armadillo, recommended!!

    • Armadillo is a C++ linear algebra library (matrix maths) aiming towards a good balance between speed and ease of use. The syntax is deliberately similar to Matlab.
      • @@ -459,7 +670,7 @@ Always free space when you don't need an array anymore.









      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -482,7 +693,7 @@ Always free space when you don't need an array anymore.











      -

      Armadillo, how to compile and install

      +

      Armadillo, how to compile and install

      For people using Ubuntu, Debian, Linux Mint, simply go to the synaptic package manager and install @@ -510,7 +721,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -543,7 +754,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -574,7 +785,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -607,7 +818,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -635,7 +846,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -667,7 +878,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -697,7 +908,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Armadillo, simple examples

      +

      Armadillo, simple examples

      @@ -729,7 +940,7 @@ For OS X users you may have to declare the paths to the include files and the li











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      We start with the linear set of equations @@ -764,7 +975,7 @@ $$











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      or $$ @@ -779,7 +990,7 @@ $$











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown \( x_1 \) @@ -803,7 +1014,7 @@ what is called a backward substitution.











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      This process can be expressed mathematically as $$ @@ -821,7 +1032,7 @@ the result from the $j$th equation. We assume obviously that \( a_{11}\ne 0 \) a











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      Our actual \( 4\times 4 \) example reads after the first operation @@ -861,7 +1072,7 @@ $$











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      The new coefficients are @@ -897,7 +1108,7 @@ We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \(











      -

      Gaussian Elimination

      +

      Gaussian Elimination

      This step is called forward substitution. @@ -935,7 +1146,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).











      -

      Linear Algebra Methods

      +

      Linear Algebra Methods

      • Gaussian elimination, \( O(2/3n^3) \) flops, general matrix
        • @@ -948,7 +1159,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).









        -

        LU Decomposition

        +

        LU Decomposition

        The LU decomposition method means that we can rewrite @@ -979,7 +1190,7 @@ $$











        -

        LU Decomposition

        +

        LU Decomposition

        LU decomposition forms the backbone of other algorithms in linear algebra, such as the @@ -1008,7 +1219,7 @@ $$











        -

        LU Decomposition, why?

        +

        LU Decomposition, why?

        There are at least three main advantages with LU decomposition compared with standard Gaussian elimination: @@ -1021,7 +1232,7 @@ There are at least three main advantages with LU decomposition compared with sta









        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        With the LU decomposition it is rather @@ -1050,7 +1261,7 @@ $$ \mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}. $$











        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        The previous equation can be calculated in two steps @@ -1063,7 +1274,7 @@ to rewrite our system of linear equations as $$ \mathbf{LUx}=\mathbf{w}, $$ -and since the determinat of \( \mathbf{L} \) is equal to 1 (by construction +and since the determinant of \( \mathbf{L} \) is equal to 1 (by construction since the diagonals of \( \mathbf{L} \) equal 1) we can use the inverse of \( \mathbf{L} \) to obtain @@ -1083,7 +1294,7 @@ through \( \mathbf{Ux}=\mathbf{y} \).











        -

        LU Decomposition, why?

        +

        LU Decomposition, why?

        For our four-dimentional example this takes the form @@ -1116,7 +1327,7 @@ needed to solve the set of \( n \) linear equations.











        -

        LU Decomposition, linear equations

        +

        LU Decomposition, linear equations

        The algorithm goes as follows @@ -1129,7 +1340,7 @@ The algorithm goes as follows









        -

        LU Decomposition, the inverse of a matrix

        +

        LU Decomposition, the inverse of a matrix

        If the inverse exists then @@ -1147,7 +1358,7 @@ $$











        -

        LU Decomposition, the inverse of a matrix

        +

        LU Decomposition, the inverse of a matrix

        If we assume that the first column (that is column 1) of the inverse matrix @@ -1181,7 +1392,7 @@ $$











        -

        LU Decomposition, the inverse

        +

        LU Decomposition, the inverse

        In a similar way we can compute the unknow entries of the second column, @@ -1205,7 +1416,7 @@ and continue till we have solved all \( n \) sets of linear equations.











        -

        Using Armadillo to perform an LU decomposition

        +

        Using Armadillo to perform an LU decomposition

        @@ -1237,236 +1448,6 @@ and continue till we have solved all \( n \) sets of linear equations. return 0; }

    -

    -









    - -

    Iterative methods, Chapter 6

    -
    - -

    - -

      -
    • Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.
      • -
    • Iterative solvers such as Basic iterative solvers, Jacobi, Gauss-Seidel, Successive over-relaxation. These methods are easy to parallelize, as we will se later. Much used in solutions of partial differential equations.
      • -
    • Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.
      • -
    -
    - - -

    -









    - -

    Iterative methods, Jacobi's method

    -
    - -

    -It is a simple method for solving -$$ -\mathbf{A}\mathbf{x}=\mathbf{b}, -$$ - -where \( \mathbf{A} \) is a matrix and \( \mathbf{x} \) and \( \mathbf{b} \) are vectors. The vector \( \mathbf{x} \) is -the unknown. - -

    -It is an iterative scheme where we start with a guess for the unknown, and -after \( k+1 \) iterations we have -$$ -\mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}), -$$ - -with \( \mathbf{A}=\mathbf{D}+\mathbf{U}+\mathbf{L} \) and -\( \mathbf{D} \) being a diagonal matrix, \( \mathbf{U} \) an upper triangular matrix and \( \mathbf{L} \) a lower triangular -matrix. - -

    -If the matrix \( \mathbf{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution. -

    - - -

    -









    - -

    Iterative methods, Jacobi's method

    -
    - -

    -We can demonstrate Jacobi's method by this \( 4\times 4 \) matrix problem. We assume a guess -for the vector elements \( x_i^{(0)} \), a guess which represents our first iteration. The new -values are obtained by substitution -$$ -\begin{align} - x_1^{(1)} =&(b_1-a_{12}x_2^{(0)} -a_{13}x_3^{(0)} - a_{14}x_4^{(0)})/a_{11} \nonumber \\ - x_2^{(1)} =&(b_2-a_{21}x_1^{(0)} - a_{23}x_3^{(0)} - a_{24}x_4^{(0)})/a_{22} \nonumber \\ - x_3^{(1)} =&(b_3- a_{31}x_1^{(0)} -a_{32}x_2^{(0)} -a_{34}x_4^{(0)})/a_{33} \nonumber \\ - x_4^{(1)}=&(b_4-a_{41}x_1^{(0)} -a_{42}x_2^{(0)} - a_{43}x_3^{(0)})/a_{44}, \nonumber -\end{align} -$$ - -which after \( k+1 \) iterations reads -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber -\end{align} -$$ -

    - - -

    -









    - -

    Iterative methods, Jacobi's method

    -
    - -

    -We can generalize the above equations to -$$ - x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii} -$$ - -or in an even more compact form as -$$ \mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}), -$$ - -with \( \mathbf{A}=\mathbf{D}+\mathbf{U}+\mathbf{L} \) and -\( \mathbf{D} \) being a diagonal matrix, \( \mathbf{U} \) an upper triangular matrix and \( \mathbf{L} \) a lower triangular -matrix. -

    - - -

    -









    - -

    Iterative methods, Gauss-Seidel's method

    -
    - -

    -Our \( 4\times 4 \) matrix problem -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber -\end{align} -$$ - -can be rewritten as -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber -\end{align} -$$ - -which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence -behavior and leads to the Gauss-Seidel method! -

    - - -

    -









    - -

    Iterative methods, Gauss-Seidel's method

    -
    - -

    -We can generalize -$$ -\begin{align} - x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ - x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ - x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ - x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber -\end{align} -$$ - -to the following form -$$ - x^{(k+1)}_i = \frac{1}{a_{ii}} \left(b_i - \sum_{j > i}a_{ij}x^{(k)}_j - \sum_{j < i}a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. -$$ - -The procedure is generally continued until the changes made by an iteration are below some tolerance. - -

    -The convergence properties of the Jacobi method and the -Gauss-Seidel method are dependent on the matrix \( \mathbf{A} \). These methods converge when -the matrix is symmetric positive-definite, or is strictly or irreducibly diagonally dominant. -Both methods sometimes converge even if these conditions are not satisfied. -

    - - -

    -









    - -

    Iterative methods, Successive over-relaxation

    -
    - -

    -Given a square system of n linear equations with unknown \( \mathbf x \): -$$ - \mathbf{A}\mathbf x = \mathbf b -$$ - -where -$$ - \mathbf{A}=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}. -$$ -

    - - -

    -









    - -

    Iterative methods, Successive over-relaxation

    -
    - -

    -Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U: -$$ - \mathbf{A} =\mathbf{D} + \mathbf{L} + \mathbf{U}, -$$ - -where -$$ - D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \quad L = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix}, \quad U = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & 0 \end{bmatrix}. -$$ - -The system of linear equations may be rewritten as: -$$ - (D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x} -$$ - -for a constant \( \omega > 1 \). -

    - - -

    -









    - -

    Iterative methods, Successive over-relaxation

    -
    - -

    -The method of successive over-relaxation is an iterative technique that solves the left hand side of this expression for \( x \), using previous value for \( x \) on the right hand side. Analytically, this may be written as: -$$ - \mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big). -$$ - -However, by taking advantage of the triangular form of \( (D+\omega L) \), the elements of \( x^{(k+1)} \) can be computed sequentially using forward substitution: -$$ - x^{(k+1)}_i = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j > i} a_{ij}x^{(k)}_j - \sum_{j < i} a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. -$$ - -The choice of relaxation factor is not necessarily easy, and depends upon the properties of the coefficient matrix. For symmetric, positive-definite matrices it can be proven that \( 0 < \omega < 2 \) will lead to convergence, but we are generally interested in faster convergence rather than just convergence. -

    - -

    diff --git a/doc/pub/Linalg/html/Linalg.html b/doc/pub/Linalg/html/Linalg.html index c28ebdded..80ec79dd2 100644 --- a/doc/pub/Linalg/html/Linalg.html +++ b/doc/pub/Linalg/html/Linalg.html @@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source - + -Data analysis and Machine Learning Lectures: Linear Algebra methods +Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays