update on lin algebra
This commit is contained in:
mhjensen
+328
-363
File diff suppressed because it is too large
Load Diff
@@ -3,9 +3,9 @@
|
||||
|
||||
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
|
||||
<meta name="generator" content="DocOnce: https://github.com/hplgit/doconce/" />
|
||||
<meta name="description" content="Data analysis and Machine Learning Lectures: Linear Algebra methods">
|
||||
<meta name="description" content="Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays">
|
||||
|
||||
<title>Data analysis and Machine Learning Lectures: Linear Algebra methods</title>
|
||||
<title>Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays</title>
|
||||
|
||||
|
||||
|
||||
@@ -132,7 +132,7 @@ MathJax.Hub.Config({
|
||||
|
||||
|
||||
|
||||
<center><h1 style="text-align: center;">Data analysis and Machine Learning Lectures: Linear Algebra methods </h1></center> <!-- document title -->
|
||||
<center><h1 style="text-align: center;">Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays</h1></center> <!-- document title -->
|
||||
|
||||
<p>
|
||||
<!-- author(s): Morten Hjorth-Jensen -->
|
||||
@@ -148,7 +148,7 @@ MathJax.Hub.Config({
|
||||
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
|
||||
<br>
|
||||
<p> <br>
|
||||
<center><h4>May 22, 2018</h4></center> <!-- date -->
|
||||
<center><h4>May 24, 2018</h4></center> <!-- date -->
|
||||
<br>
|
||||
<p>
|
||||
|
||||
@@ -159,12 +159,19 @@ MathJax.Hub.Config({
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec0">To do </h2>
|
||||
<h2 id="___sec0">Introduction </h2>
|
||||
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.
|
||||
|
||||
<ul>
|
||||
<p><li> add material on python handling of matrices and vectors</li>
|
||||
<p><li> keep c++ material?</li>
|
||||
</ul>
|
||||
<p>
|
||||
As discussed in the introductory notes, these series of lectures focuses both on using
|
||||
central Python packages like <b>tensorflow</b> and <b>scikit-learn</b> 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.
|
||||
</section>
|
||||
|
||||
|
||||
@@ -172,10 +179,10 @@ MathJax.Hub.Config({
|
||||
<h2 id="___sec1">Important Matrix and vector handling packages </h2>
|
||||
|
||||
<p>
|
||||
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 <b>Numerical Recipes</b>. 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.
|
||||
|
||||
<ul>
|
||||
|
||||
@@ -187,7 +194,22 @@ developed in the 1970s, namely EISPACK and LINPACK.
|
||||
</ul>
|
||||
<p>
|
||||
|
||||
<b>Add python material on linear algebra and array handling, text on numpy etc</b>
|
||||
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,
|
||||
<b>Numpy</b> 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, <b>Armadillo</b> 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.
|
||||
</section>
|
||||
|
||||
|
||||
@@ -311,7 +333,217 @@ For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec7">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
|
||||
<h2 id="___sec7">Numpy and arrays </h2>
|
||||
<a href="http://www.numpy.org/" target="_blank">Numpy</a> provides an easy way to handle arrays in Python. The standard way to import this library is as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
x = np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.array([<span style="color: #B452CD">1</span>, <span style="color: #B452CD">2</span>, <span style="color: #B452CD">3</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4</span>, <span style="color: #B452CD">7</span>, <span style="color: #B452CD">8</span>]))
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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 <b>log</b> function
|
||||
from Python's <b>math</b> module. The looping is done explicitely by the
|
||||
<b>np.log</b> function. The alternative, and slower way to compute the
|
||||
logarithms of a vector would be to write
|
||||
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">math</span> <span style="color: #8B008B; font-weight: bold">import</span> log
|
||||
x = np.array([<span style="color: #B452CD">4</span>, <span style="color: #B452CD">7</span>, <span style="color: #B452CD">8</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">0</span>, <span style="color: #658b00">len</span>(x)):
|
||||
x[i] = log(x[i])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
We note that our code is much longer already and we need to import the <b>log</b> function from the <b>math</b> module.
|
||||
The attentive reader will also notice that the output is \( [1, 1, 2] \). Python interprets automacally our numbers as integers (like the <b>automatic</b> keyword in C++). To change this we could define our array elements to be double precision numbers as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4</span>, <span style="color: #B452CD">7</span>, <span style="color: #B452CD">8</span>], dtype = np.float64))
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
To check the number of bytes (remember that one byte contains eight bits for double precision variables), you can use simple use the <b>itemsize</b> functionality (the array \( x \) is actually an object which inherits the functionalities defined in Numpy) as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x.itemsize)
|
||||
</pre></div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec8">Matrices in Python </h2>
|
||||
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)
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
A = np.log(np.array([ [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>], [<span style="color: #B452CD">3.0</span>, <span style="color: #B452CD">10.0</span>, <span style="color: #B452CD">11.0</span>], [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">5.0</span>, <span style="color: #B452CD">7.0</span>] ]))
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
If we use the <b>shape</b> 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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
A = np.log(np.array([ [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>], [<span style="color: #B452CD">3.0</span>, <span style="color: #B452CD">10.0</span>, <span style="color: #B452CD">11.0</span>], [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">5.0</span>, <span style="color: #B452CD">7.0</span>] ]))
|
||||
<span style="color: #228B22"># print the first column, row-major order and elements start with 0</span>
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A[:,<span style="color: #B452CD">0</span>])
|
||||
</pre></div>
|
||||
<p>
|
||||
We can continue this was by printing out other columns or rows. The example here prints out the second column
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
A = np.log(np.array([ [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>], [<span style="color: #B452CD">3.0</span>, <span style="color: #B452CD">10.0</span>, <span style="color: #B452CD">11.0</span>], [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">5.0</span>, <span style="color: #B452CD">7.0</span>] ]))
|
||||
<span style="color: #228B22"># print the first column, row-major order and elements start with 0</span>
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A[<span style="color: #B452CD">1</span>,:])
|
||||
</pre></div>
|
||||
<p>
|
||||
Numpy contains many other functionalities that allow us to slice, subdivide etc etc arrays. We strongly recommend that you look up the <a href="http://www.numpy.org/" target="_blank">Numpy website for more details</a>. Useful functions when defining a matrix are the <b>np.zeros</b> function which declares a matrix of a given dimension and sets all elements to zero
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
<span style="color: #228B22"># define a matrix of dimension 10 x 10 and set all elements to zero</span>
|
||||
A = np.zeros( (n, n) )
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
or initializing all elements to
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
<span style="color: #228B22"># define a matrix of dimension 10 x 10 and set all elements to one</span>
|
||||
A = np.ones( (n, n) )
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
or as unitarily distributed random numbers (see the material on random number generators in the statistics part)
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
<span style="color: #228B22"># define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]</span>
|
||||
A = np.random.rand(n, n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p> <br>
|
||||
$$
|
||||
\hat{\Sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\
|
||||
\sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\
|
||||
\sigma_{zx} & \sigma_{zy} & \sigma_{zz}
|
||||
\end{bmatrix},
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
where for example
|
||||
<p> <br>
|
||||
$$
|
||||
\sigma_{xy} =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}).
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
The Numpy function <b>np.cov</b> 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 <b>np.vstack</b> function which takes each vector of dimension \( 1\times n \) and produces a $ 3\times n$ matrix \( \hat{W} \)
|
||||
<p> <br>
|
||||
$$
|
||||
\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},
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
<p>
|
||||
which in turn is converted into into the \( 3 times 3 \) covariance matrix
|
||||
\( \hat{\Sigma} \) via the Numpy function <b>np.cov()</b>. 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 <b>np.mean(x)</b>. We can also extract the eigenvalues of the
|
||||
covariance matrix through the <b>np.linalg.eig()</b> function.
|
||||
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #228B22"># Importing various packages</span>
|
||||
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
|
||||
n = <span style="color: #B452CD">100</span>
|
||||
x = np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(np.mean(x))
|
||||
y = <span style="color: #B452CD">4</span>+<span style="color: #B452CD">3</span>*x+np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(np.mean(y))
|
||||
z = x**<span style="color: #B452CD">3</span>+np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(np.mean(z))
|
||||
W = np.vstack((x, y, z))
|
||||
Sigma = np.cov(W)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(Sigma)
|
||||
Eigvals, Eigvecs = np.linalg.eig(Sigma)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(Eigvals)
|
||||
</pre></div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec9">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -337,7 +569,7 @@ Note the way the matrix is organized, row-major order.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec8">Matrix Handling in C/C++ </h2>
|
||||
<h2 id="___sec10">Matrix Handling in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -367,7 +599,7 @@ In C/C++ this would be coded like
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec9">Matrix Handling in C/C++ </h2>
|
||||
<h2 id="___sec11">Matrix Handling in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -398,7 +630,7 @@ In C/C++ this would be coded like
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec10">Dynamic memory allocation in C/C++ </h2>
|
||||
<h2 id="___sec12">Dynamic memory allocation in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
At least three possibilities in this course
|
||||
@@ -415,7 +647,7 @@ At least three possibilities in this course
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec11">Matrix Handling in C/C++, Dynamic Allocation </h2>
|
||||
<h2 id="___sec13">Matrix Handling in C/C++, Dynamic Allocation </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -445,7 +677,7 @@ Always free space when you don't need an array anymore.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec12">Armadillo, recommended!! </h2>
|
||||
<h2 id="___sec14">Armadillo, recommended!! </h2>
|
||||
|
||||
<ul>
|
||||
|
||||
@@ -463,7 +695,7 @@ Always free space when you don't need an array anymore.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec13">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec15">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -487,7 +719,7 @@ Always free space when you don't need an array anymore.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec14">Armadillo, how to compile and install </h2>
|
||||
<h2 id="___sec16">Armadillo, how to compile and install </h2>
|
||||
|
||||
<p>
|
||||
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
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec15">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec17">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -550,7 +782,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec16">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec18">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -582,7 +814,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec17">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec19">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -616,7 +848,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec18">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec20">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -645,7 +877,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec19">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec21">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -678,7 +910,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec20">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec22">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -709,7 +941,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec21">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec23">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -742,7 +974,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec22">Gaussian Elimination </h2>
|
||||
<h2 id="___sec24">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
We start with the linear set of equations
|
||||
@@ -781,7 +1013,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec23">Gaussian Elimination </h2>
|
||||
<h2 id="___sec25">Gaussian Elimination </h2>
|
||||
or
|
||||
|
||||
<p> <br>
|
||||
@@ -798,7 +1030,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec24">Gaussian Elimination </h2>
|
||||
<h2 id="___sec26">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
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.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec25">Gaussian Elimination </h2>
|
||||
<h2 id="___sec27">Gaussian Elimination </h2>
|
||||
This process can be expressed mathematically as
|
||||
|
||||
<p> <br>
|
||||
@@ -844,7 +1076,7 @@ the result from the $j$th equation. We assume obviously that \( a_{11}\ne 0 \) a
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec26">Gaussian Elimination </h2>
|
||||
<h2 id="___sec28">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
Our actual \( 4\times 4 \) example reads after the first operation
|
||||
@@ -888,7 +1120,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec27">Gaussian Elimination </h2>
|
||||
<h2 id="___sec29">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
The new coefficients are
|
||||
@@ -930,7 +1162,7 @@ We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \(
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec28">Gaussian Elimination </h2>
|
||||
<h2 id="___sec30">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
This step is called forward substitution.
|
||||
@@ -972,7 +1204,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec29">Linear Algebra Methods </h2>
|
||||
<h2 id="___sec31">Linear Algebra Methods </h2>
|
||||
|
||||
<ul>
|
||||
|
||||
@@ -992,7 +1224,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec30">LU Decomposition </h2>
|
||||
<h2 id="___sec32">LU Decomposition </h2>
|
||||
|
||||
<p>
|
||||
The LU decomposition method means that we can rewrite
|
||||
@@ -1025,7 +1257,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec31">LU Decomposition </h2>
|
||||
<h2 id="___sec33">LU Decomposition </h2>
|
||||
|
||||
<p>
|
||||
LU decomposition forms the backbone of other algorithms in linear algebra, such as the
|
||||
@@ -1058,7 +1290,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec32">LU Decomposition, why? </h2>
|
||||
<h2 id="___sec34">LU Decomposition, why? </h2>
|
||||
|
||||
<p>
|
||||
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
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec33">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec35">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
With the LU decomposition it is rather
|
||||
@@ -1110,7 +1342,7 @@ $$ \mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}. $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec34">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec36">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
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}, $$
|
||||
<p> <br>
|
||||
|
||||
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} \).
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec35">LU Decomposition, why? </h2>
|
||||
<h2 id="___sec37">LU Decomposition, why? </h2>
|
||||
|
||||
<p>
|
||||
For our four-dimentional example this takes the form
|
||||
@@ -1188,7 +1420,7 @@ needed to solve the set of \( n \) linear equations.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec36">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec38">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
The algorithm goes as follows
|
||||
@@ -1205,7 +1437,7 @@ The algorithm goes as follows
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec37">LU Decomposition, the inverse of a matrix </h2>
|
||||
<h2 id="___sec39">LU Decomposition, the inverse of a matrix </h2>
|
||||
|
||||
<p>
|
||||
If the inverse exists then
|
||||
@@ -1227,7 +1459,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec38">LU Decomposition, the inverse of a matrix </h2>
|
||||
<h2 id="___sec40">LU Decomposition, the inverse of a matrix </h2>
|
||||
|
||||
<p>
|
||||
If we assume that the first column (that is column 1) of the inverse matrix
|
||||
@@ -1265,7 +1497,7 @@ $$
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec39">LU Decomposition, the inverse </h2>
|
||||
<h2 id="___sec41">LU Decomposition, the inverse </h2>
|
||||
|
||||
<p>
|
||||
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.
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec40"><a href="https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp" target="_blank">Using Armadillo to perform an LU decomposition</a> </h2>
|
||||
<h2 id="___sec42"><a href="https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp" target="_blank">Using Armadillo to perform an LU decomposition</a> </h2>
|
||||
<p>
|
||||
|
||||
<!-- code=c++ (!bc cppcod) typeset with pygments style "perldoc" -->
|
||||
@@ -1326,260 +1558,6 @@ and continue till we have solved all \( n \) sets of linear equations.
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec41">Iterative methods, Chapter 6 </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<ul>
|
||||
<p><li> Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.</li>
|
||||
<p><li> 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.</li>
|
||||
<p><li> Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.</li>
|
||||
</ul>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec42">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
It is a simple method for solving
|
||||
<p> <br>
|
||||
$$
|
||||
\mathbf{A}\mathbf{x}=\mathbf{b},
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
where \( \mathbf{A} \) is a matrix and \( \mathbf{x} \) and \( \mathbf{b} \) are vectors. The vector \( \mathbf{x} \) is
|
||||
the unknown.
|
||||
|
||||
<p>
|
||||
It is an iterative scheme where we start with a guess for the unknown, and
|
||||
after \( k+1 \) iterations we have
|
||||
|
||||
<p> <br>
|
||||
$$
|
||||
\mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}),
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
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.
|
||||
|
||||
<p>
|
||||
If the matrix \( \mathbf{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution.
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec43">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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
|
||||
<p> <br>
|
||||
$$
|
||||
\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}
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
which after \( k+1 \) iterations reads
|
||||
<p> <br>
|
||||
$$
|
||||
\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}
|
||||
$$
|
||||
<p> <br>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec44">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
We can generalize the above equations to
|
||||
<p> <br>
|
||||
$$
|
||||
x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii}
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
or in an even more compact form as
|
||||
<p> <br>
|
||||
$$ \mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}),
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
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.
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec45">Iterative methods, Gauss-Seidel's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
Our \( 4\times 4 \) matrix problem
|
||||
<p> <br>
|
||||
$$
|
||||
\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}
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
can be rewritten as
|
||||
<p> <br>
|
||||
$$
|
||||
\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}
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence
|
||||
behavior and leads to the Gauss-Seidel method!
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec46">Iterative methods, Gauss-Seidel's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
We can generalize
|
||||
<p> <br>
|
||||
$$
|
||||
\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}
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
to the following form
|
||||
<p> <br>
|
||||
$$
|
||||
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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
The procedure is generally continued until the changes made by an iteration are below some tolerance.
|
||||
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec47">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
Given a square system of n linear equations with unknown \( \mathbf x \):
|
||||
<p> <br>
|
||||
$$
|
||||
\mathbf{A}\mathbf x = \mathbf b
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
where
|
||||
<p> <br>
|
||||
$$
|
||||
\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}.
|
||||
$$
|
||||
<p> <br>
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec48">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U:
|
||||
<p> <br>
|
||||
$$
|
||||
\mathbf{A} =\mathbf{D} + \mathbf{L} + \mathbf{U},
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
where
|
||||
<p> <br>
|
||||
$$
|
||||
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}.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
The system of linear equations may be rewritten as:
|
||||
<p> <br>
|
||||
$$
|
||||
(D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
for a constant \( \omega > 1 \).
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
<section>
|
||||
<h2 id="___sec49">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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:
|
||||
<p> <br>
|
||||
$$
|
||||
\mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big).
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
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:
|
||||
<p> <br>
|
||||
$$
|
||||
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.
|
||||
$$
|
||||
<p> <br>
|
||||
|
||||
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.
|
||||
</div>
|
||||
</section>
|
||||
|
||||
|
||||
|
||||
</div> <!-- class="slides" -->
|
||||
</div> <!-- class="reveal" -->
|
||||
|
||||
@@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source
|
||||
<head>
|
||||
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
|
||||
<meta name="generator" content="DocOnce: https://github.com/hplgit/doconce/" />
|
||||
<meta name="description" content="Data analysis and Machine Learning Lectures: Linear Algebra methods">
|
||||
<meta name="description" content="Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays">
|
||||
|
||||
<title>Data analysis and Machine Learning Lectures: Linear Algebra methods</title>
|
||||
<title>Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays</title>
|
||||
|
||||
|
||||
<link href="https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_styles/style_solarized_box/css/solarized_light_code.css" rel="stylesheet" type="text/css" title="light"/>
|
||||
@@ -60,7 +60,7 @@ div { text-align: justify; text-justify: inter-word; }
|
||||
|
||||
<!-- tocinfo
|
||||
{'highest level': 2,
|
||||
'sections': [('To do', 2, None, '___sec0'),
|
||||
'sections': [('Introduction', 2, None, '___sec0'),
|
||||
('Important Matrix and vector handling packages',
|
||||
2,
|
||||
None,
|
||||
@@ -70,74 +70,58 @@ div { text-align: justify; text-justify: inter-word; }
|
||||
('Basic Matrix Features', 2, None, '___sec4'),
|
||||
('Some famous Matrices', 2, None, '___sec5'),
|
||||
('Basic Matrix Features', 2, None, '___sec6'),
|
||||
('Numpy and arrays', 2, None, '___sec7'),
|
||||
('Matrices in Python', 2, None, '___sec8'),
|
||||
('Matrix Handling in C/C++, Static and Dynamical allocation',
|
||||
2,
|
||||
None,
|
||||
'___sec7'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec8'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec9'),
|
||||
('Dynamic memory allocation in C/C++', 2, None, '___sec10'),
|
||||
'___sec9'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec10'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec11'),
|
||||
('Dynamic memory allocation in C/C++', 2, None, '___sec12'),
|
||||
('Matrix Handling in C/C++, Dynamic Allocation',
|
||||
2,
|
||||
None,
|
||||
'___sec11'),
|
||||
('Armadillo, recommended!!', 2, None, '___sec12'),
|
||||
('Armadillo, simple examples', 2, None, '___sec13'),
|
||||
('Armadillo, how to compile and install', 2, None, '___sec14'),
|
||||
'___sec13'),
|
||||
('Armadillo, recommended!!', 2, None, '___sec14'),
|
||||
('Armadillo, simple examples', 2, None, '___sec15'),
|
||||
('Armadillo, simple examples', 2, None, '___sec16'),
|
||||
('Armadillo, how to compile and install', 2, None, '___sec16'),
|
||||
('Armadillo, simple examples', 2, None, '___sec17'),
|
||||
('Armadillo, simple examples', 2, None, '___sec18'),
|
||||
('Armadillo, simple examples', 2, None, '___sec19'),
|
||||
('Armadillo, simple examples', 2, None, '___sec20'),
|
||||
('Armadillo, simple examples', 2, None, '___sec21'),
|
||||
('Gaussian Elimination', 2, None, '___sec22'),
|
||||
('Gaussian Elimination', 2, None, '___sec23'),
|
||||
('Armadillo, simple examples', 2, None, '___sec22'),
|
||||
('Armadillo, simple examples', 2, None, '___sec23'),
|
||||
('Gaussian Elimination', 2, None, '___sec24'),
|
||||
('Gaussian Elimination', 2, None, '___sec25'),
|
||||
('Gaussian Elimination', 2, None, '___sec26'),
|
||||
('Gaussian Elimination', 2, None, '___sec27'),
|
||||
('Gaussian Elimination', 2, None, '___sec28'),
|
||||
('Linear Algebra Methods', 2, None, '___sec29'),
|
||||
('LU Decomposition', 2, None, '___sec30'),
|
||||
('LU Decomposition', 2, None, '___sec31'),
|
||||
('LU Decomposition, why?', 2, None, '___sec32'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec33'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec34'),
|
||||
('LU Decomposition, why?', 2, None, '___sec35'),
|
||||
('Gaussian Elimination', 2, None, '___sec29'),
|
||||
('Gaussian Elimination', 2, None, '___sec30'),
|
||||
('Linear Algebra Methods', 2, None, '___sec31'),
|
||||
('LU Decomposition', 2, None, '___sec32'),
|
||||
('LU Decomposition', 2, None, '___sec33'),
|
||||
('LU Decomposition, why?', 2, None, '___sec34'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec35'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec36'),
|
||||
('LU Decomposition, why?', 2, None, '___sec37'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec38'),
|
||||
('LU Decomposition, the inverse of a matrix',
|
||||
2,
|
||||
None,
|
||||
'___sec37'),
|
||||
'___sec39'),
|
||||
('LU Decomposition, the inverse of a matrix',
|
||||
2,
|
||||
None,
|
||||
'___sec38'),
|
||||
('LU Decomposition, the inverse', 2, None, '___sec39'),
|
||||
'___sec40'),
|
||||
('LU Decomposition, the inverse', 2, None, '___sec41'),
|
||||
('"Using Armadillo to perform an LU '
|
||||
'decomposition":"https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp"',
|
||||
2,
|
||||
None,
|
||||
'___sec40'),
|
||||
('Iterative methods, Chapter 6', 2, None, '___sec41'),
|
||||
("Iterative methods, Jacobi's method", 2, None, '___sec42'),
|
||||
("Iterative methods, Jacobi's method", 2, None, '___sec43'),
|
||||
("Iterative methods, Jacobi's method", 2, None, '___sec44'),
|
||||
("Iterative methods, Gauss-Seidel's method", 2, None, '___sec45'),
|
||||
("Iterative methods, Gauss-Seidel's method", 2, None, '___sec46'),
|
||||
('Iterative methods, Successive over-relaxation',
|
||||
2,
|
||||
None,
|
||||
'___sec47'),
|
||||
('Iterative methods, Successive over-relaxation',
|
||||
2,
|
||||
None,
|
||||
'___sec48'),
|
||||
('Iterative methods, Successive over-relaxation',
|
||||
2,
|
||||
None,
|
||||
'___sec49')]}
|
||||
'___sec42')]}
|
||||
end of tocinfo -->
|
||||
|
||||
<body>
|
||||
@@ -163,7 +147,7 @@ MathJax.Hub.Config({
|
||||
|
||||
|
||||
|
||||
<center><h1>Data analysis and Machine Learning Lectures: Linear Algebra methods </h1></center> <!-- document title -->
|
||||
<center><h1>Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays</h1></center> <!-- document title -->
|
||||
|
||||
<p>
|
||||
<!-- author(s): Morten Hjorth-Jensen -->
|
||||
@@ -179,27 +163,35 @@ MathJax.Hub.Config({
|
||||
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
|
||||
<br>
|
||||
<p>
|
||||
<center><h4>May 22, 2018</h4></center> <!-- date -->
|
||||
<center><h4>May 24, 2018</h4></center> <!-- date -->
|
||||
<br>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec0">To do </h2>
|
||||
<h2 id="___sec0">Introduction </h2>
|
||||
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.
|
||||
|
||||
<ul>
|
||||
<li> add material on python handling of matrices and vectors</li>
|
||||
<li> keep c++ material?</li>
|
||||
</ul>
|
||||
<p>
|
||||
As discussed in the introductory notes, these series of lectures focuses both on using
|
||||
central Python packages like <b>tensorflow</b> and <b>scikit-learn</b> 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.
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
<p>
|
||||
<!-- !split -->
|
||||
|
||||
<h2 id="___sec1">Important Matrix and vector handling packages </h2>
|
||||
|
||||
<p>
|
||||
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 <b>Numerical Recipes</b>. 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.
|
||||
|
||||
<ul>
|
||||
<li> LINPACK: package for linear equations and least square problems.</li>
|
||||
@@ -207,10 +199,25 @@ developed in the 1970s, namely EISPACK and LINPACK.
|
||||
<li> 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 <a href="http://www.netlib.org" target="_blank"><tt>http://www.netlib.org</tt></a>.</li>
|
||||
</ul>
|
||||
|
||||
<b>Add python material on linear algebra and array handling, text on numpy etc</b>
|
||||
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,
|
||||
<b>Numpy</b> 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, <b>Armadillo</b> 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.
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
<!-- !split -->
|
||||
|
||||
<h2 id="___sec2">Basic Matrix Features </h2>
|
||||
|
||||
@@ -239,6 +246,8 @@ $$
|
||||
<h2 id="___sec3">Basic Matrix Features </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec7">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
|
||||
<h2 id="___sec7">Numpy and arrays </h2>
|
||||
<a href="http://www.numpy.org/" target="_blank">Numpy</a> provides an easy way to handle arrays in Python. The standard way to import this library is as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
x = np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.array([<span style="color: #B452CD">1</span>, <span style="color: #B452CD">2</span>, <span style="color: #B452CD">3</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4</span>, <span style="color: #B452CD">7</span>, <span style="color: #B452CD">8</span>]))
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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 <b>log</b> function
|
||||
from Python's <b>math</b> module. The looping is done explicitely by the
|
||||
<b>np.log</b> function. The alternative, and slower way to compute the
|
||||
logarithms of a vector would be to write
|
||||
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">math</span> <span style="color: #8B008B; font-weight: bold">import</span> log
|
||||
x = np.array([<span style="color: #B452CD">4</span>, <span style="color: #B452CD">7</span>, <span style="color: #B452CD">8</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">0</span>, <span style="color: #658b00">len</span>(x)):
|
||||
x[i] = log(x[i])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
We note that our code is much longer already and we need to import the <b>log</b> function from the <b>math</b> module.
|
||||
The attentive reader will also notice that the output is \( [1, 1, 2] \). Python interprets automacally our numbers as integers (like the <b>automatic</b> keyword in C++). To change this we could define our array elements to be double precision numbers as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4</span>, <span style="color: #B452CD">7</span>, <span style="color: #B452CD">8</span>], dtype = np.float64))
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
To check the number of bytes (remember that one byte contains eight bits for double precision variables), you can use simple use the <b>itemsize</b> functionality (the array \( x \) is actually an object which inherits the functionalities defined in Numpy) as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
x = np.log(np.array([<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>])
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(x.itemsize)
|
||||
</pre></div>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec8">Matrices in Python </h2>
|
||||
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)
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
A = np.log(np.array([ [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>], [<span style="color: #B452CD">3.0</span>, <span style="color: #B452CD">10.0</span>, <span style="color: #B452CD">11.0</span>], [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">5.0</span>, <span style="color: #B452CD">7.0</span>] ]))
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
If we use the <b>shape</b> 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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
A = np.log(np.array([ [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>], [<span style="color: #B452CD">3.0</span>, <span style="color: #B452CD">10.0</span>, <span style="color: #B452CD">11.0</span>], [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">5.0</span>, <span style="color: #B452CD">7.0</span>] ]))
|
||||
<span style="color: #228B22"># print the first column, row-major order and elements start with 0</span>
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A[:,<span style="color: #B452CD">0</span>])
|
||||
</pre></div>
|
||||
<p>
|
||||
We can continue this was by printing out other columns or rows. The example here prints out the second column
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
A = np.log(np.array([ [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">7.0</span>, <span style="color: #B452CD">8.0</span>], [<span style="color: #B452CD">3.0</span>, <span style="color: #B452CD">10.0</span>, <span style="color: #B452CD">11.0</span>], [<span style="color: #B452CD">4.0</span>, <span style="color: #B452CD">5.0</span>, <span style="color: #B452CD">7.0</span>] ]))
|
||||
<span style="color: #228B22"># print the first column, row-major order and elements start with 0</span>
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A[<span style="color: #B452CD">1</span>,:])
|
||||
</pre></div>
|
||||
<p>
|
||||
Numpy contains many other functionalities that allow us to slice, subdivide etc etc arrays. We strongly recommend that you look up the <a href="http://www.numpy.org/" target="_blank">Numpy website for more details</a>. Useful functions when defining a matrix are the <b>np.zeros</b> function which declares a matrix of a given dimension and sets all elements to zero
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
<span style="color: #228B22"># define a matrix of dimension 10 x 10 and set all elements to zero</span>
|
||||
A = np.zeros( (n, n) )
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
or initializing all elements to
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
<span style="color: #228B22"># define a matrix of dimension 10 x 10 and set all elements to one</span>
|
||||
A = np.ones( (n, n) )
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
or as unitarily distributed random numbers (see the material on random number generators in the statistics part)
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
n = <span style="color: #B452CD">10</span>
|
||||
<span style="color: #228B22"># define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]</span>
|
||||
A = np.random.rand(n, n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
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 <b>np.cov</b> 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 <b>np.vstack</b> 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},
|
||||
$$
|
||||
|
||||
<p>
|
||||
which in turn is converted into into the \( 3 times 3 \) covariance matrix
|
||||
\( \hat{\Sigma} \) via the Numpy function <b>np.cov()</b>. 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 <b>np.mean(x)</b>. We can also extract the eigenvalues of the
|
||||
covariance matrix through the <b>np.linalg.eig()</b> function.
|
||||
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
|
||||
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #228B22"># Importing various packages</span>
|
||||
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
|
||||
|
||||
n = <span style="color: #B452CD">100</span>
|
||||
x = np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(np.mean(x))
|
||||
y = <span style="color: #B452CD">4</span>+<span style="color: #B452CD">3</span>*x+np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(np.mean(y))
|
||||
z = x**<span style="color: #B452CD">3</span>+np.random.normal(size=n)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(np.mean(z))
|
||||
W = np.vstack((x, y, z))
|
||||
Sigma = np.cov(W)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(Sigma)
|
||||
Eigvals, Eigvecs = np.linalg.eig(Sigma)
|
||||
<span style="color: #8B008B; font-weight: bold">print</span>(Eigvals)
|
||||
</pre></div>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec9">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -343,7 +554,7 @@ Note the way the matrix is organized, row-major order.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec8">Matrix Handling in C/C++ </h2>
|
||||
<h2 id="___sec10">Matrix Handling in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -372,7 +583,7 @@ In C/C++ this would be coded like
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec9">Matrix Handling in C/C++ </h2>
|
||||
<h2 id="___sec11">Matrix Handling in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -402,7 +613,7 @@ In C/C++ this would be coded like
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec10">Dynamic memory allocation in C/C++ </h2>
|
||||
<h2 id="___sec12">Dynamic memory allocation in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
At least three possibilities in this course
|
||||
@@ -415,7 +626,7 @@ At least three possibilities in this course
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec11">Matrix Handling in C/C++, Dynamic Allocation </h2>
|
||||
<h2 id="___sec13">Matrix Handling in C/C++, Dynamic Allocation </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -447,7 +658,7 @@ Always free space when you don't need an array anymore.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec12">Armadillo, recommended!! </h2>
|
||||
<h2 id="___sec14">Armadillo, recommended!! </h2>
|
||||
|
||||
<ul>
|
||||
<li> 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.</li>
|
||||
@@ -459,7 +670,7 @@ Always free space when you don't need an array anymore.
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec13">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec15">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -482,7 +693,7 @@ Always free space when you don't need an array anymore.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec14">Armadillo, how to compile and install </h2>
|
||||
<h2 id="___sec16">Armadillo, how to compile and install </h2>
|
||||
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec15">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec17">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -543,7 +754,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec16">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec18">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -574,7 +785,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec17">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec19">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -607,7 +818,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec18">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec20">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -635,7 +846,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec19">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec21">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -667,7 +878,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec20">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec22">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -697,7 +908,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec21">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec23">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -729,7 +940,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec22">Gaussian Elimination </h2>
|
||||
<h2 id="___sec24">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
We start with the linear set of equations
|
||||
@@ -764,7 +975,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec23">Gaussian Elimination </h2>
|
||||
<h2 id="___sec25">Gaussian Elimination </h2>
|
||||
or
|
||||
|
||||
$$
|
||||
@@ -779,7 +990,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec24">Gaussian Elimination </h2>
|
||||
<h2 id="___sec26">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
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.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec25">Gaussian Elimination </h2>
|
||||
<h2 id="___sec27">Gaussian Elimination </h2>
|
||||
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
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec26">Gaussian Elimination </h2>
|
||||
<h2 id="___sec28">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
Our actual \( 4\times 4 \) example reads after the first operation
|
||||
@@ -861,7 +1072,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec27">Gaussian Elimination </h2>
|
||||
<h2 id="___sec29">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
The new coefficients are
|
||||
@@ -897,7 +1108,7 @@ We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \(
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec28">Gaussian Elimination </h2>
|
||||
<h2 id="___sec30">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
This step is called forward substitution.
|
||||
@@ -935,7 +1146,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec29">Linear Algebra Methods </h2>
|
||||
<h2 id="___sec31">Linear Algebra Methods </h2>
|
||||
|
||||
<ul>
|
||||
<li> Gaussian elimination, \( O(2/3n^3) \) flops, general matrix</li>
|
||||
@@ -948,7 +1159,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec30">LU Decomposition </h2>
|
||||
<h2 id="___sec32">LU Decomposition </h2>
|
||||
|
||||
<p>
|
||||
The LU decomposition method means that we can rewrite
|
||||
@@ -979,7 +1190,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec31">LU Decomposition </h2>
|
||||
<h2 id="___sec33">LU Decomposition </h2>
|
||||
|
||||
<p>
|
||||
LU decomposition forms the backbone of other algorithms in linear algebra, such as the
|
||||
@@ -1008,7 +1219,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec32">LU Decomposition, why? </h2>
|
||||
<h2 id="___sec34">LU Decomposition, why? </h2>
|
||||
|
||||
<p>
|
||||
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
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec33">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec35">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
With the LU decomposition it is rather
|
||||
@@ -1050,7 +1261,7 @@ $$ \mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}. $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec34">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec36">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
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} \).
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec35">LU Decomposition, why? </h2>
|
||||
<h2 id="___sec37">LU Decomposition, why? </h2>
|
||||
|
||||
<p>
|
||||
For our four-dimentional example this takes the form
|
||||
@@ -1116,7 +1327,7 @@ needed to solve the set of \( n \) linear equations.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec36">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec38">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
The algorithm goes as follows
|
||||
@@ -1129,7 +1340,7 @@ The algorithm goes as follows
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec37">LU Decomposition, the inverse of a matrix </h2>
|
||||
<h2 id="___sec39">LU Decomposition, the inverse of a matrix </h2>
|
||||
|
||||
<p>
|
||||
If the inverse exists then
|
||||
@@ -1147,7 +1358,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec38">LU Decomposition, the inverse of a matrix </h2>
|
||||
<h2 id="___sec40">LU Decomposition, the inverse of a matrix </h2>
|
||||
|
||||
<p>
|
||||
If we assume that the first column (that is column 1) of the inverse matrix
|
||||
@@ -1181,7 +1392,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec39">LU Decomposition, the inverse </h2>
|
||||
<h2 id="___sec41">LU Decomposition, the inverse </h2>
|
||||
|
||||
<p>
|
||||
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.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec40"><a href="https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp" target="_blank">Using Armadillo to perform an LU decomposition</a> </h2>
|
||||
<h2 id="___sec42"><a href="https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp" target="_blank">Using Armadillo to perform an LU decomposition</a> </h2>
|
||||
<p>
|
||||
|
||||
<!-- code=c++ (!bc cppcod) typeset with pygments style "perldoc" -->
|
||||
@@ -1237,236 +1448,6 @@ and continue till we have solved all \( n \) sets of linear equations.
|
||||
<span style="color: #8B008B; font-weight: bold">return</span> <span style="color: #B452CD">0</span>;
|
||||
}
|
||||
</pre></div>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec41">Iterative methods, Chapter 6 </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
|
||||
<ul>
|
||||
<li> Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.</li>
|
||||
<li> 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.</li>
|
||||
<li> Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.</li>
|
||||
</ul>
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec42">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
If the matrix \( \mathbf{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec43">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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}
|
||||
$$
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec44">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec45">Iterative methods, Gauss-Seidel's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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!
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec46">Iterative methods, Gauss-Seidel's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec47">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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}.
|
||||
$$
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec48">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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 \).
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec49">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
|
||||
<!-- ------------------- end of main content --------------- -->
|
||||
|
||||
+302
-321
@@ -6,9 +6,9 @@ Automatically generated HTML file from DocOnce source
|
||||
<head>
|
||||
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
|
||||
<meta name="generator" content="DocOnce: https://github.com/hplgit/doconce/" />
|
||||
<meta name="description" content="Data analysis and Machine Learning Lectures: Linear Algebra methods">
|
||||
<meta name="description" content="Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays">
|
||||
|
||||
<title>Data analysis and Machine Learning Lectures: Linear Algebra methods</title>
|
||||
<title>Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays</title>
|
||||
|
||||
|
||||
<style type="text/css">
|
||||
@@ -65,7 +65,7 @@ div { text-align: justify; text-justify: inter-word; }
|
||||
|
||||
<!-- tocinfo
|
||||
{'highest level': 2,
|
||||
'sections': [('To do', 2, None, '___sec0'),
|
||||
'sections': [('Introduction', 2, None, '___sec0'),
|
||||
('Important Matrix and vector handling packages',
|
||||
2,
|
||||
None,
|
||||
@@ -75,74 +75,58 @@ div { text-align: justify; text-justify: inter-word; }
|
||||
('Basic Matrix Features', 2, None, '___sec4'),
|
||||
('Some famous Matrices', 2, None, '___sec5'),
|
||||
('Basic Matrix Features', 2, None, '___sec6'),
|
||||
('Numpy and arrays', 2, None, '___sec7'),
|
||||
('Matrices in Python', 2, None, '___sec8'),
|
||||
('Matrix Handling in C/C++, Static and Dynamical allocation',
|
||||
2,
|
||||
None,
|
||||
'___sec7'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec8'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec9'),
|
||||
('Dynamic memory allocation in C/C++', 2, None, '___sec10'),
|
||||
'___sec9'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec10'),
|
||||
('Matrix Handling in C/C++', 2, None, '___sec11'),
|
||||
('Dynamic memory allocation in C/C++', 2, None, '___sec12'),
|
||||
('Matrix Handling in C/C++, Dynamic Allocation',
|
||||
2,
|
||||
None,
|
||||
'___sec11'),
|
||||
('Armadillo, recommended!!', 2, None, '___sec12'),
|
||||
('Armadillo, simple examples', 2, None, '___sec13'),
|
||||
('Armadillo, how to compile and install', 2, None, '___sec14'),
|
||||
'___sec13'),
|
||||
('Armadillo, recommended!!', 2, None, '___sec14'),
|
||||
('Armadillo, simple examples', 2, None, '___sec15'),
|
||||
('Armadillo, simple examples', 2, None, '___sec16'),
|
||||
('Armadillo, how to compile and install', 2, None, '___sec16'),
|
||||
('Armadillo, simple examples', 2, None, '___sec17'),
|
||||
('Armadillo, simple examples', 2, None, '___sec18'),
|
||||
('Armadillo, simple examples', 2, None, '___sec19'),
|
||||
('Armadillo, simple examples', 2, None, '___sec20'),
|
||||
('Armadillo, simple examples', 2, None, '___sec21'),
|
||||
('Gaussian Elimination', 2, None, '___sec22'),
|
||||
('Gaussian Elimination', 2, None, '___sec23'),
|
||||
('Armadillo, simple examples', 2, None, '___sec22'),
|
||||
('Armadillo, simple examples', 2, None, '___sec23'),
|
||||
('Gaussian Elimination', 2, None, '___sec24'),
|
||||
('Gaussian Elimination', 2, None, '___sec25'),
|
||||
('Gaussian Elimination', 2, None, '___sec26'),
|
||||
('Gaussian Elimination', 2, None, '___sec27'),
|
||||
('Gaussian Elimination', 2, None, '___sec28'),
|
||||
('Linear Algebra Methods', 2, None, '___sec29'),
|
||||
('LU Decomposition', 2, None, '___sec30'),
|
||||
('LU Decomposition', 2, None, '___sec31'),
|
||||
('LU Decomposition, why?', 2, None, '___sec32'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec33'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec34'),
|
||||
('LU Decomposition, why?', 2, None, '___sec35'),
|
||||
('Gaussian Elimination', 2, None, '___sec29'),
|
||||
('Gaussian Elimination', 2, None, '___sec30'),
|
||||
('Linear Algebra Methods', 2, None, '___sec31'),
|
||||
('LU Decomposition', 2, None, '___sec32'),
|
||||
('LU Decomposition', 2, None, '___sec33'),
|
||||
('LU Decomposition, why?', 2, None, '___sec34'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec35'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec36'),
|
||||
('LU Decomposition, why?', 2, None, '___sec37'),
|
||||
('LU Decomposition, linear equations', 2, None, '___sec38'),
|
||||
('LU Decomposition, the inverse of a matrix',
|
||||
2,
|
||||
None,
|
||||
'___sec37'),
|
||||
'___sec39'),
|
||||
('LU Decomposition, the inverse of a matrix',
|
||||
2,
|
||||
None,
|
||||
'___sec38'),
|
||||
('LU Decomposition, the inverse', 2, None, '___sec39'),
|
||||
'___sec40'),
|
||||
('LU Decomposition, the inverse', 2, None, '___sec41'),
|
||||
('"Using Armadillo to perform an LU '
|
||||
'decomposition":"https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp"',
|
||||
2,
|
||||
None,
|
||||
'___sec40'),
|
||||
('Iterative methods, Chapter 6', 2, None, '___sec41'),
|
||||
("Iterative methods, Jacobi's method", 2, None, '___sec42'),
|
||||
("Iterative methods, Jacobi's method", 2, None, '___sec43'),
|
||||
("Iterative methods, Jacobi's method", 2, None, '___sec44'),
|
||||
("Iterative methods, Gauss-Seidel's method", 2, None, '___sec45'),
|
||||
("Iterative methods, Gauss-Seidel's method", 2, None, '___sec46'),
|
||||
('Iterative methods, Successive over-relaxation',
|
||||
2,
|
||||
None,
|
||||
'___sec47'),
|
||||
('Iterative methods, Successive over-relaxation',
|
||||
2,
|
||||
None,
|
||||
'___sec48'),
|
||||
('Iterative methods, Successive over-relaxation',
|
||||
2,
|
||||
None,
|
||||
'___sec49')]}
|
||||
'___sec42')]}
|
||||
end of tocinfo -->
|
||||
|
||||
<body>
|
||||
@@ -168,7 +152,7 @@ MathJax.Hub.Config({
|
||||
|
||||
|
||||
|
||||
<center><h1>Data analysis and Machine Learning Lectures: Linear Algebra methods </h1></center> <!-- document title -->
|
||||
<center><h1>Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays</h1></center> <!-- document title -->
|
||||
|
||||
<p>
|
||||
<!-- author(s): Morten Hjorth-Jensen -->
|
||||
@@ -184,27 +168,35 @@ MathJax.Hub.Config({
|
||||
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
|
||||
<br>
|
||||
<p>
|
||||
<center><h4>May 22, 2018</h4></center> <!-- date -->
|
||||
<center><h4>May 24, 2018</h4></center> <!-- date -->
|
||||
<br>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec0">To do </h2>
|
||||
<h2 id="___sec0">Introduction </h2>
|
||||
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.
|
||||
|
||||
<ul>
|
||||
<li> add material on python handling of matrices and vectors</li>
|
||||
<li> keep c++ material?</li>
|
||||
</ul>
|
||||
<p>
|
||||
As discussed in the introductory notes, these series of lectures focuses both on using
|
||||
central Python packages like <b>tensorflow</b> and <b>scikit-learn</b> 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.
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
<p>
|
||||
<!-- !split -->
|
||||
|
||||
<h2 id="___sec1">Important Matrix and vector handling packages </h2>
|
||||
|
||||
<p>
|
||||
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 <b>Numerical Recipes</b>. 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.
|
||||
|
||||
<ul>
|
||||
<li> LINPACK: package for linear equations and least square problems.</li>
|
||||
@@ -212,10 +204,25 @@ developed in the 1970s, namely EISPACK and LINPACK.
|
||||
<li> 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 <a href="http://www.netlib.org" target="_blank"><tt>http://www.netlib.org</tt></a>.</li>
|
||||
</ul>
|
||||
|
||||
<b>Add python material on linear algebra and array handling, text on numpy etc</b>
|
||||
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,
|
||||
<b>Numpy</b> 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, <b>Armadillo</b> 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.
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
<!-- !split -->
|
||||
|
||||
<h2 id="___sec2">Basic Matrix Features </h2>
|
||||
|
||||
@@ -244,6 +251,8 @@ $$
|
||||
<h2 id="___sec3">Basic Matrix Features </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
|
||||
<p>
|
||||
The inverse of a matrix is defined by
|
||||
|
||||
@@ -321,7 +330,209 @@ For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec7">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
|
||||
<h2 id="___sec7">Numpy and arrays </h2>
|
||||
<a href="http://www.numpy.org/" target="_blank">Numpy</a> provides an easy way to handle arrays in Python. The standard way to import this library is as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
n <span style="color: #666666">=</span> <span style="color: #666666">10</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(size<span style="color: #666666">=</span>n)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>array([<span style="color: #666666">1</span>, <span style="color: #666666">2</span>, <span style="color: #666666">3</span>])
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([<span style="color: #666666">4</span>, <span style="color: #666666">7</span>, <span style="color: #666666">8</span>]))
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
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 <b>log</b> function
|
||||
from Python's <b>math</b> module. The looping is done explicitely by the
|
||||
<b>np.log</b> function. The alternative, and slower way to compute the
|
||||
logarithms of a vector would be to write
|
||||
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">math</span> <span style="color: #008000; font-weight: bold">import</span> log
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>array([<span style="color: #666666">4</span>, <span style="color: #666666">7</span>, <span style="color: #666666">8</span>])
|
||||
<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">0</span>, <span style="color: #008000">len</span>(x)):
|
||||
x[i] <span style="color: #666666">=</span> log(x[i])
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
We note that our code is much longer already and we need to import the <b>log</b> function from the <b>math</b> module.
|
||||
The attentive reader will also notice that the output is \( [1, 1, 2] \). Python interprets automacally our numbers as integers (like the <b>automatic</b> keyword in C++). To change this we could define our array elements to be double precision numbers as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([<span style="color: #666666">4</span>, <span style="color: #666666">7</span>, <span style="color: #666666">8</span>], dtype <span style="color: #666666">=</span> np<span style="color: #666666">.</span>float64))
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([<span style="color: #666666">4.0</span>, <span style="color: #666666">7.0</span>, <span style="color: #666666">8.0</span>])
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x)
|
||||
</pre></div>
|
||||
<p>
|
||||
To check the number of bytes (remember that one byte contains eight bits for double precision variables), you can use simple use the <b>itemsize</b> functionality (the array \( x \) is actually an object which inherits the functionalities defined in Numpy) as
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([<span style="color: #666666">4.0</span>, <span style="color: #666666">7.0</span>, <span style="color: #666666">8.0</span>])
|
||||
<span style="color: #008000; font-weight: bold">print</span>(x<span style="color: #666666">.</span>itemsize)
|
||||
</pre></div>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec8">Matrices in Python </h2>
|
||||
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)
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
A <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([ [<span style="color: #666666">4.0</span>, <span style="color: #666666">7.0</span>, <span style="color: #666666">8.0</span>], [<span style="color: #666666">3.0</span>, <span style="color: #666666">10.0</span>, <span style="color: #666666">11.0</span>], [<span style="color: #666666">4.0</span>, <span style="color: #666666">5.0</span>, <span style="color: #666666">7.0</span>] ]))
|
||||
<span style="color: #008000; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
If we use the <b>shape</b> 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
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
A <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([ [<span style="color: #666666">4.0</span>, <span style="color: #666666">7.0</span>, <span style="color: #666666">8.0</span>], [<span style="color: #666666">3.0</span>, <span style="color: #666666">10.0</span>, <span style="color: #666666">11.0</span>], [<span style="color: #666666">4.0</span>, <span style="color: #666666">5.0</span>, <span style="color: #666666">7.0</span>] ]))
|
||||
<span style="color: #408080; font-style: italic"># print the first column, row-major order and elements start with 0</span>
|
||||
<span style="color: #008000; font-weight: bold">print</span>(A[:,<span style="color: #666666">0</span>])
|
||||
</pre></div>
|
||||
<p>
|
||||
We can continue this was by printing out other columns or rows. The example here prints out the second column
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
A <span style="color: #666666">=</span> np<span style="color: #666666">.</span>log(np<span style="color: #666666">.</span>array([ [<span style="color: #666666">4.0</span>, <span style="color: #666666">7.0</span>, <span style="color: #666666">8.0</span>], [<span style="color: #666666">3.0</span>, <span style="color: #666666">10.0</span>, <span style="color: #666666">11.0</span>], [<span style="color: #666666">4.0</span>, <span style="color: #666666">5.0</span>, <span style="color: #666666">7.0</span>] ]))
|
||||
<span style="color: #408080; font-style: italic"># print the first column, row-major order and elements start with 0</span>
|
||||
<span style="color: #008000; font-weight: bold">print</span>(A[<span style="color: #666666">1</span>,:])
|
||||
</pre></div>
|
||||
<p>
|
||||
Numpy contains many other functionalities that allow us to slice, subdivide etc etc arrays. We strongly recommend that you look up the <a href="http://www.numpy.org/" target="_blank">Numpy website for more details</a>. Useful functions when defining a matrix are the <b>np.zeros</b> function which declares a matrix of a given dimension and sets all elements to zero
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
n <span style="color: #666666">=</span> <span style="color: #666666">10</span>
|
||||
<span style="color: #408080; font-style: italic"># define a matrix of dimension 10 x 10 and set all elements to zero</span>
|
||||
A <span style="color: #666666">=</span> np<span style="color: #666666">.</span>zeros( (n, n) )
|
||||
<span style="color: #008000; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
or initializing all elements to
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
n <span style="color: #666666">=</span> <span style="color: #666666">10</span>
|
||||
<span style="color: #408080; font-style: italic"># define a matrix of dimension 10 x 10 and set all elements to one</span>
|
||||
A <span style="color: #666666">=</span> np<span style="color: #666666">.</span>ones( (n, n) )
|
||||
<span style="color: #008000; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
or as unitarily distributed random numbers (see the material on random number generators in the statistics part)
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
n <span style="color: #666666">=</span> <span style="color: #666666">10</span>
|
||||
<span style="color: #408080; font-style: italic"># define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]</span>
|
||||
A <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(n, n)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(A)
|
||||
</pre></div>
|
||||
<p>
|
||||
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 <b>np.cov</b> 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 <b>np.vstack</b> 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},
|
||||
$$
|
||||
|
||||
<p>
|
||||
which in turn is converted into into the \( 3 times 3 \) covariance matrix
|
||||
\( \hat{\Sigma} \) via the Numpy function <b>np.cov()</b>. 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 <b>np.mean(x)</b>. We can also extract the eigenvalues of the
|
||||
covariance matrix through the <b>np.linalg.eig()</b> function.
|
||||
|
||||
<p>
|
||||
|
||||
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
|
||||
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #408080; font-style: italic"># Importing various packages</span>
|
||||
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
|
||||
|
||||
n <span style="color: #666666">=</span> <span style="color: #666666">100</span>
|
||||
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(size<span style="color: #666666">=</span>n)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(np<span style="color: #666666">.</span>mean(x))
|
||||
y <span style="color: #666666">=</span> <span style="color: #666666">4+3*</span>x<span style="color: #666666">+</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(size<span style="color: #666666">=</span>n)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(np<span style="color: #666666">.</span>mean(y))
|
||||
z <span style="color: #666666">=</span> x<span style="color: #666666">**3+</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>normal(size<span style="color: #666666">=</span>n)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(np<span style="color: #666666">.</span>mean(z))
|
||||
W <span style="color: #666666">=</span> np<span style="color: #666666">.</span>vstack((x, y, z))
|
||||
Sigma <span style="color: #666666">=</span> np<span style="color: #666666">.</span>cov(W)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(Sigma)
|
||||
Eigvals, Eigvecs <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>eig(Sigma)
|
||||
<span style="color: #008000; font-weight: bold">print</span>(Eigvals)
|
||||
</pre></div>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec9">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -348,7 +559,7 @@ Note the way the matrix is organized, row-major order.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec8">Matrix Handling in C/C++ </h2>
|
||||
<h2 id="___sec10">Matrix Handling in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -377,7 +588,7 @@ In C/C++ this would be coded like
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec9">Matrix Handling in C/C++ </h2>
|
||||
<h2 id="___sec11">Matrix Handling in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -407,7 +618,7 @@ In C/C++ this would be coded like
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec10">Dynamic memory allocation in C/C++ </h2>
|
||||
<h2 id="___sec12">Dynamic memory allocation in C/C++ </h2>
|
||||
|
||||
<p>
|
||||
At least three possibilities in this course
|
||||
@@ -420,7 +631,7 @@ At least three possibilities in this course
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec11">Matrix Handling in C/C++, Dynamic Allocation </h2>
|
||||
<h2 id="___sec13">Matrix Handling in C/C++, Dynamic Allocation </h2>
|
||||
|
||||
<p>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
@@ -452,7 +663,7 @@ Always free space when you don't need an array anymore.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec12">Armadillo, recommended!! </h2>
|
||||
<h2 id="___sec14">Armadillo, recommended!! </h2>
|
||||
|
||||
<ul>
|
||||
<li> 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.</li>
|
||||
@@ -464,7 +675,7 @@ Always free space when you don't need an array anymore.
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec13">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec15">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -487,7 +698,7 @@ Always free space when you don't need an array anymore.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec14">Armadillo, how to compile and install </h2>
|
||||
<h2 id="___sec16">Armadillo, how to compile and install </h2>
|
||||
|
||||
<p>
|
||||
For people using Ubuntu, Debian, Linux Mint, simply go to the synaptic package manager and install
|
||||
@@ -515,7 +726,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec15">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec17">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -548,7 +759,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec16">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec18">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -579,7 +790,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec17">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec19">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -612,7 +823,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec18">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec20">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -640,7 +851,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec19">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec21">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -672,7 +883,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec20">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec22">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -702,7 +913,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec21">Armadillo, simple examples </h2>
|
||||
<h2 id="___sec23">Armadillo, simple examples </h2>
|
||||
|
||||
<p>
|
||||
|
||||
@@ -734,7 +945,7 @@ For OS X users you may have to declare the paths to the include files and the li
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec22">Gaussian Elimination </h2>
|
||||
<h2 id="___sec24">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
We start with the linear set of equations
|
||||
@@ -769,7 +980,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec23">Gaussian Elimination </h2>
|
||||
<h2 id="___sec25">Gaussian Elimination </h2>
|
||||
or
|
||||
|
||||
$$
|
||||
@@ -784,7 +995,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec24">Gaussian Elimination </h2>
|
||||
<h2 id="___sec26">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown \( x_1 \)
|
||||
@@ -808,7 +1019,7 @@ what is called a backward substitution.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec25">Gaussian Elimination </h2>
|
||||
<h2 id="___sec27">Gaussian Elimination </h2>
|
||||
This process can be expressed mathematically as
|
||||
|
||||
$$
|
||||
@@ -826,7 +1037,7 @@ the result from the $j$th equation. We assume obviously that \( a_{11}\ne 0 \) a
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec26">Gaussian Elimination </h2>
|
||||
<h2 id="___sec28">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
Our actual \( 4\times 4 \) example reads after the first operation
|
||||
@@ -866,7 +1077,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec27">Gaussian Elimination </h2>
|
||||
<h2 id="___sec29">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
The new coefficients are
|
||||
@@ -902,7 +1113,7 @@ We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \(
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec28">Gaussian Elimination </h2>
|
||||
<h2 id="___sec30">Gaussian Elimination </h2>
|
||||
|
||||
<p>
|
||||
This step is called forward substitution.
|
||||
@@ -940,7 +1151,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec29">Linear Algebra Methods </h2>
|
||||
<h2 id="___sec31">Linear Algebra Methods </h2>
|
||||
|
||||
<ul>
|
||||
<li> Gaussian elimination, \( O(2/3n^3) \) flops, general matrix</li>
|
||||
@@ -953,7 +1164,7 @@ adding \( 10^7+1 \). With single precision this results in \( 10^7 \).
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec30">LU Decomposition </h2>
|
||||
<h2 id="___sec32">LU Decomposition </h2>
|
||||
|
||||
<p>
|
||||
The LU decomposition method means that we can rewrite
|
||||
@@ -984,7 +1195,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec31">LU Decomposition </h2>
|
||||
<h2 id="___sec33">LU Decomposition </h2>
|
||||
|
||||
<p>
|
||||
LU decomposition forms the backbone of other algorithms in linear algebra, such as the
|
||||
@@ -1013,7 +1224,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec32">LU Decomposition, why? </h2>
|
||||
<h2 id="___sec34">LU Decomposition, why? </h2>
|
||||
|
||||
<p>
|
||||
There are at least three main advantages with LU decomposition compared with standard Gaussian elimination:
|
||||
@@ -1026,7 +1237,7 @@ There are at least three main advantages with LU decomposition compared with sta
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec33">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec35">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
With the LU decomposition it is rather
|
||||
@@ -1055,7 +1266,7 @@ $$ \mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}. $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec34">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec36">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
The previous equation can be calculated in two steps
|
||||
@@ -1068,7 +1279,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
|
||||
|
||||
@@ -1088,7 +1299,7 @@ through \( \mathbf{Ux}=\mathbf{y} \).
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec35">LU Decomposition, why? </h2>
|
||||
<h2 id="___sec37">LU Decomposition, why? </h2>
|
||||
|
||||
<p>
|
||||
For our four-dimentional example this takes the form
|
||||
@@ -1121,7 +1332,7 @@ needed to solve the set of \( n \) linear equations.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec36">LU Decomposition, linear equations </h2>
|
||||
<h2 id="___sec38">LU Decomposition, linear equations </h2>
|
||||
|
||||
<p>
|
||||
The algorithm goes as follows
|
||||
@@ -1134,7 +1345,7 @@ The algorithm goes as follows
|
||||
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec37">LU Decomposition, the inverse of a matrix </h2>
|
||||
<h2 id="___sec39">LU Decomposition, the inverse of a matrix </h2>
|
||||
|
||||
<p>
|
||||
If the inverse exists then
|
||||
@@ -1152,7 +1363,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec38">LU Decomposition, the inverse of a matrix </h2>
|
||||
<h2 id="___sec40">LU Decomposition, the inverse of a matrix </h2>
|
||||
|
||||
<p>
|
||||
If we assume that the first column (that is column 1) of the inverse matrix
|
||||
@@ -1186,7 +1397,7 @@ $$
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec39">LU Decomposition, the inverse </h2>
|
||||
<h2 id="___sec41">LU Decomposition, the inverse </h2>
|
||||
|
||||
<p>
|
||||
In a similar way we can compute the unknow entries of the second column,
|
||||
@@ -1210,7 +1421,7 @@ and continue till we have solved all \( n \) sets of linear equations.
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec40"><a href="https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp" target="_blank">Using Armadillo to perform an LU decomposition</a> </h2>
|
||||
<h2 id="___sec42"><a href="https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp" target="_blank">Using Armadillo to perform an LU decomposition</a> </h2>
|
||||
<p>
|
||||
|
||||
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
|
||||
@@ -1242,236 +1453,6 @@ and continue till we have solved all \( n \) sets of linear equations.
|
||||
<span style="color: #008000; font-weight: bold">return</span> <span style="color: #666666">0</span>;
|
||||
}
|
||||
</pre></div>
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec41">Iterative methods, Chapter 6 </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
|
||||
<ul>
|
||||
<li> Direct solvers such as Gauss elimination and LU decomposition discussed in connection with project 1.</li>
|
||||
<li> 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.</li>
|
||||
<li> Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.</li>
|
||||
</ul>
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec42">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
If the matrix \( \mathbf{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec43">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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}
|
||||
$$
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec44">Iterative methods, Jacobi's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec45">Iterative methods, Gauss-Seidel's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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!
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec46">Iterative methods, Gauss-Seidel's method </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec47">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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}.
|
||||
$$
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec48">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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 \).
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
|
||||
|
||||
<h2 id="___sec49">Iterative methods, Successive over-relaxation </h2>
|
||||
<div class="alert alert-block alert-block alert-text-normal">
|
||||
<b></b>
|
||||
<p>
|
||||
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.
|
||||
</div>
|
||||
|
||||
|
||||
<p>
|
||||
|
||||
<!-- ------------------- end of main content --------------- -->
|
||||
|
||||
+545
-484
File diff suppressed because it is too large
Load Diff
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
+200
-203
@@ -1,30 +1,53 @@
|
||||
TITLE: Data analysis and Machine Learning Lectures: Linear Algebra methods
|
||||
TITLE: Data analysis and Machine Learning Lectures: Linear Algebra and Handling of Arrays
|
||||
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
|
||||
DATE: today
|
||||
|
||||
|
||||
!split
|
||||
===== 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.
|
||||
|
||||
|
||||
!split
|
||||
!split
|
||||
===== 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.
|
||||
* LAPACK:package for solving symmetric, unsymmetric and generalized eigenvalue problems. From LAPACK's website URL: "http://www.netlib.org" it is possible to download for free all source codes from this library. Both C/C++ and Fortran versions are available.
|
||||
* 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 URL: "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.
|
||||
|
||||
!split
|
||||
!split
|
||||
===== Basic Matrix Features =====
|
||||
|
||||
!bblock Matrix properties reminder
|
||||
@@ -48,6 +71,7 @@ _Add python material on linear algebra and array handling, text on numpy etc_
|
||||
!split
|
||||
===== Basic Matrix Features =====
|
||||
!bblock
|
||||
|
||||
The inverse of a matrix is defined by
|
||||
|
||||
!bt
|
||||
@@ -101,6 +125,171 @@ For an $N\times N$ matrix $\mathbf{A}$ the following properties are all equival
|
||||
* $0$ is not eigenvalue of $\mathbf{A}$.
|
||||
!eblock
|
||||
|
||||
!split
|
||||
===== Numpy and arrays =====
|
||||
"Numpy":"http://www.numpy.org/" provides an easy way to handle arrays in Python. The standard way to import this library is as
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
n = 10
|
||||
x = np.random.normal(size=n)
|
||||
print(x)
|
||||
!ec
|
||||
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
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
x = np.array([1, 2, 3])
|
||||
print(x)
|
||||
!ec
|
||||
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
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
x = np.log(np.array([4, 7, 8]))
|
||||
print(x)
|
||||
!ec
|
||||
|
||||
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
|
||||
|
||||
!bc pycod
|
||||
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)
|
||||
!ec
|
||||
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
|
||||
!bc pycod
|
||||
import numpy as np
|
||||
x = np.log(np.array([4, 7, 8], dtype = np.float64))
|
||||
print(x)
|
||||
!ec
|
||||
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])
|
||||
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])
|
||||
print(x.itemsize)
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== 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)
|
||||
!bc pycod
|
||||
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)
|
||||
!ec
|
||||
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
|
||||
!bc pycod
|
||||
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])
|
||||
!ec
|
||||
We can continue this was by printing out other columns or rows. The example here prints out the second column
|
||||
!bc pycod
|
||||
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,:])
|
||||
!ec
|
||||
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
|
||||
!bc pycod
|
||||
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)
|
||||
!ec
|
||||
or initializing all elements to
|
||||
!bc pycod
|
||||
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)
|
||||
!ec
|
||||
or as unitarily distributed random numbers (see the material on random number generators in the statistics part)
|
||||
!bc pycod
|
||||
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)
|
||||
!ec
|
||||
|
||||
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
|
||||
!bt
|
||||
\[
|
||||
\hat{\Sigma} = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\
|
||||
\sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\
|
||||
\sigma_{zx} & \sigma_{zy} & \sigma_{zz}
|
||||
\end{bmatrix},
|
||||
\]
|
||||
!et
|
||||
where for example
|
||||
!bt
|
||||
\[
|
||||
\sigma_{xy} =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}).
|
||||
\]
|
||||
!et
|
||||
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}$
|
||||
!bt
|
||||
\[
|
||||
\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},
|
||||
\]
|
||||
!et
|
||||
|
||||
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.
|
||||
|
||||
!bc pycod
|
||||
# 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)
|
||||
!ec
|
||||
|
||||
|
||||
!split
|
||||
===== Matrix Handling in C/C++, Static and Dynamical allocation =====
|
||||
@@ -744,7 +933,7 @@ to rewrite our system of linear equations as
|
||||
!bt
|
||||
\[ \mathbf{LUx}=\mathbf{w}, \]
|
||||
!et
|
||||
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
|
||||
|
||||
@@ -910,198 +1099,6 @@ int main()
|
||||
}
|
||||
!ec
|
||||
|
||||
!split
|
||||
===== Iterative methods, Chapter 6 =====
|
||||
!bblock
|
||||
* 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.
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Iterative methods, Jacobi's method =====
|
||||
!bblock
|
||||
It is a simple method for solving
|
||||
!bt
|
||||
\[
|
||||
\mathbf{A}\mathbf{x}=\mathbf{b},
|
||||
\]
|
||||
!et
|
||||
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
|
||||
!bt
|
||||
\[
|
||||
\mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}),
|
||||
\]
|
||||
!et
|
||||
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.
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Iterative methods, Jacobi's method =====
|
||||
!bblock
|
||||
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
|
||||
!bt
|
||||
\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}
|
||||
!et
|
||||
which after $k+1$ iterations reads
|
||||
!bt
|
||||
\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}
|
||||
!et
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Iterative methods, Jacobi's method =====
|
||||
!bblock
|
||||
We can generalize the above equations to
|
||||
!bt
|
||||
\[
|
||||
x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii}
|
||||
\]
|
||||
!et
|
||||
or in an even more compact form as
|
||||
!bt
|
||||
\[ \mathbf{x}^{(k+1)}= \mathbf{D}^{-1}(\mathbf{b}-(\mathbf{L}+\mathbf{U})\mathbf{x}^{(k)}),
|
||||
\]
|
||||
!et
|
||||
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.
|
||||
!eblock
|
||||
|
||||
!split
|
||||
===== Iterative methods, Gauss-Seidel's method =====
|
||||
!bblock
|
||||
Our $4\times 4$ matrix problem
|
||||
!bt
|
||||
\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}
|
||||
!et
|
||||
can be rewritten as
|
||||
!bt
|
||||
\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}
|
||||
!et
|
||||
which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence
|
||||
behavior and leads to the Gauss-Seidel method!
|
||||
!eblock
|
||||
|
||||
!split
|
||||
===== Iterative methods, Gauss-Seidel's method =====
|
||||
!bblock
|
||||
We can generalize
|
||||
!bt
|
||||
\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}
|
||||
!et
|
||||
to the following form
|
||||
!bt
|
||||
\[
|
||||
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.
|
||||
\]
|
||||
!et
|
||||
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.
|
||||
!eblock
|
||||
|
||||
!split
|
||||
===== Iterative methods, Successive over-relaxation =====
|
||||
!bblock
|
||||
Given a square system of n linear equations with unknown $\mathbf x$:
|
||||
!bt
|
||||
\[
|
||||
\mathbf{A}\mathbf x = \mathbf b
|
||||
\]
|
||||
!et
|
||||
where
|
||||
!bt
|
||||
\[
|
||||
\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}.
|
||||
\]
|
||||
!et
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Iterative methods, Successive over-relaxation =====
|
||||
!bblock
|
||||
Then A can be decomposed into a diagonal component D, and strictly lower and upper triangular components L and U:
|
||||
!bt
|
||||
\[
|
||||
\mathbf{A} =\mathbf{D} + \mathbf{L} + \mathbf{U},
|
||||
\]
|
||||
!et
|
||||
where
|
||||
!bt
|
||||
\[
|
||||
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}.
|
||||
\]
|
||||
!et
|
||||
The system of linear equations may be rewritten as:
|
||||
!bt
|
||||
\[
|
||||
(D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}
|
||||
\]
|
||||
!et
|
||||
for a constant $\omega > 1$.
|
||||
!eblock
|
||||
|
||||
|
||||
!split
|
||||
===== Iterative methods, Successive over-relaxation =====
|
||||
!bblock
|
||||
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:
|
||||
!bt
|
||||
\[
|
||||
\mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big).
|
||||
\]
|
||||
!et
|
||||
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:
|
||||
!bt
|
||||
\[
|
||||
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.
|
||||
\]
|
||||
!et
|
||||
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.
|
||||
!eblock
|
||||
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user