lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
f459b195f4e979247e79fdbffe93ca5fee1f466d
FYS-STK4155/doc/pub/Linalg/html/Linalg.html
T
mhjensen 065331c18d Update some material
2018-01-27 09:55:46 -05:00

1479 lines
71 KiB
HTML
Raw Blame History

<!--
Automatically generated HTML file from DocOnce source
(https://github.com/hplgit/doconce/)
-->
<html>
<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">
<title>Data analysis and Machine Learning Lectures: Linear Algebra methods</title>
<style type="text/css">
/* bloodish style */
body {
font-family: Helvetica, Verdana, Arial, Sans-serif;
color: #404040;
background: #ffffff;
}
h1 { font-size: 1.8em; color: #8A0808; }
h2 { font-size: 1.6em; color: #8A0808; }
h3 { font-size: 1.4em; color: #8A0808; }
h4 { color: #8A0808; }
a { color: #8A0808; text-decoration:none; }
tt { font-family: "Courier New", Courier; }
/* pre style removed because it will interfer with pygments */
p { text-indent: 0px; }
hr { border: 0; width: 80%; border-bottom: 1px solid #aaa}
p.caption { width: 80%; font-style: normal; text-align: left; }
hr.figure { border: 0; width: 80%; border-bottom: 1px solid #aaa}
.alert-text-small { font-size: 80%; }
.alert-text-large { font-size: 130%; }
.alert-text-normal { font-size: 90%; }
.alert {
padding:8px 35px 8px 14px; margin-bottom:18px;
text-shadow:0 1px 0 rgba(255,255,255,0.5);
border:1px solid #bababa;
border-radius: 4px;
-webkit-border-radius: 4px;
-moz-border-radius: 4px;
color: #555;
background-color: #f8f8f8;
background-position: 10px 5px;
background-repeat: no-repeat;
background-size: 38px;
padding-left: 55px;
width: 75%;
}
.alert-block {padding-top:14px; padding-bottom:14px}
.alert-block > p, .alert-block > ul {margin-bottom:1em}
.alert li {margin-top: 1em}
.alert-block p+p {margin-top:5px}
.alert-notice { background-image: url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_notice.png); }
.alert-summary { background-image:url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_summary.png); }
.alert-warning { background-image: url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_warning.png); }
.alert-question {background-image:url(https://cdn.rawgit.com/hplgit/doconce/master/bundled/html_images/small_gray_question.png); }
div { text-align: justify; text-justify: inter-word; }
</style>
</head>
<!-- tocinfo
{'highest level': 2,
'sections': [('Important Matrix and vector handling packages',
2,
None,
'___sec0'),
('Basic Matrix Features', 2, None, '___sec1'),
('Basic Matrix Features', 2, None, '___sec2'),
('Basic Matrix Features', 2, None, '___sec3'),
('Some famous Matrices', 2, None, '___sec4'),
('Basic Matrix Features', 2, None, '___sec5'),
('Matrix Handling in C/C++, Static and Dynamical allocation',
2,
None,
'___sec6'),
('Matrix Handling in C/C++', 2, None, '___sec7'),
('Matrix Handling in C/C++', 2, None, '___sec8'),
('Dynamic memory allocation in C/C++', 2, None, '___sec9'),
('Matrix Handling in C/C++, Dynamic Allocation',
2,
None,
'___sec10'),
('Armadillo, recommended!!', 2, None, '___sec11'),
('Armadillo, simple examples', 2, None, '___sec12'),
('Armadillo, how to compile and install', 2, None, '___sec13'),
('Armadillo, simple examples', 2, None, '___sec14'),
('Armadillo, simple examples', 2, None, '___sec15'),
('Armadillo, simple examples', 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'),
('Gaussian Elimination', 2, None, '___sec21'),
('Gaussian Elimination', 2, None, '___sec22'),
('Gaussian Elimination', 2, None, '___sec23'),
('Gaussian Elimination', 2, None, '___sec24'),
('Gaussian Elimination', 2, None, '___sec25'),
('Gaussian Elimination', 2, None, '___sec26'),
('Gaussian Elimination', 2, None, '___sec27'),
('Linear Algebra Methods', 2, None, '___sec28'),
('LU Decomposition', 2, None, '___sec29'),
('LU Decomposition', 2, None, '___sec30'),
('LU Decomposition, why?', 2, None, '___sec31'),
('LU Decomposition, linear equations', 2, None, '___sec32'),
('LU Decomposition, linear equations', 2, None, '___sec33'),
('LU Decomposition, why?', 2, None, '___sec34'),
('LU Decomposition, linear equations', 2, None, '___sec35'),
('LU Decomposition, the inverse of a matrix',
2,
None,
'___sec36'),
('LU Decomposition, the inverse of a matrix',
2,
None,
'___sec37'),
('LU Decomposition, the inverse', 2, None, '___sec38'),
('"Using Armadillo to perform an LU '
'decomposition":"https://github.com/CompPhysics/ComputationalPhysicsMSU/blob/master/doc/Programs/CppQtCodesLectures/MatrixTest/main.cpp"',
2,
None,
'___sec39'),
('Iterative methods, Chapter 6', 2, None, '___sec40'),
("Iterative methods, Jacobi's method", 2, None, '___sec41'),
("Iterative methods, Jacobi's method", 2, None, '___sec42'),
("Iterative methods, Jacobi's method", 2, None, '___sec43'),
("Iterative methods, Gauss-Seidel's method", 2, None, '___sec44'),
("Iterative methods, Gauss-Seidel's method", 2, None, '___sec45'),
('Iterative methods, Successive over-relaxation',
2,
None,
'___sec46'),
('Iterative methods, Successive over-relaxation',
2,
None,
'___sec47'),
('Iterative methods, Successive over-relaxation',
2,
None,
'___sec48')]}
end of tocinfo -->
<body>
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
TeX: {
equationNumbers: { autoNumber: "AMS" },
extensions: ["AMSmath.js", "AMSsymbols.js", "autobold.js", "color.js"]
}
});
</script>
<script type="text/javascript" async
src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.1/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
<!-- ------------------- main content ---------------------- -->
<center><h1>Data analysis and Machine Learning Lectures: Linear Algebra methods </h1></center> <!-- document title -->
<p>
<!-- author(s): Morten Hjorth-Jensen -->
<center>
<b>Morten Hjorth-Jensen</b> [1, 2]
</center>
<p>
<!-- institution(s) -->
<center>[1] <b>Department of Physics, University of Oslo</b></center>
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Jan 27, 2018</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec0">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
software package LAPACK, which follows two other popular packages
developed in the 1970s, namely EISPACK and LINPACK.
<ul>
<li> LINPACK: package for linear equations and least square problems.</li>
<li> LAPACK:package for solving symmetric, unsymmetric and generalized eigenvalue problems. From LAPACK's website <a href="http://www.netlib.org" target="_blank"><tt>http://www.netlib.org</tt></a> it is possible to download for free all source codes from this library. Both C/C++ and Fortran versions are available.</li>
<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>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec1">Basic Matrix Features </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Matrix properties reminder.</b>
<p>
$$
\mathbf{A} =
\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\
a_{21} & a_{22} & a_{23} & a_{24} \\
a_{31} & a_{32} & a_{33} & a_{34} \\
a_{41} & a_{42} & a_{43} & a_{44}
\end{bmatrix}\qquad
\mathbf{I} =
\begin{bmatrix} 1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix}
$$
</div>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec2">Basic Matrix Features </h2>
<div class="alert alert-block alert-block alert-text-normal">
<b></b>
<p>
The inverse of a matrix is defined by
$$
\mathbf{A}^{-1} \cdot \mathbf{A} = I
$$
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec3">Basic Matrix Features </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Matrix Properties Reminder.</b>
<p>
<p>
<table border="1">
<thead>
<tr><th align="center"> Relations </th> <th align="center"> Name </th> <th align="center"> matrix elements </th> </tr>
</thead>
<tbody>
<tr><td align="center"> \( A = A^{T} \) </td> <td align="center"> symmetric </td> <td align="center"> \( a_{ij} = a_{ji} \) </td> </tr>
<tr><td align="center"> \( A = \left (A^{T} \right )^{-1} \) </td> <td align="center"> real orthogonal </td> <td align="center"> \( \sum_k a_{ik} a_{jk} = \sum_k a_{ki} a_{kj} = \delta_{ij} \) </td> </tr>
<tr><td align="center"> \( A = A^{ * } \) </td> <td align="center"> real matrix </td> <td align="center"> \( a_{ij} = a_{ij}^{ * } \) </td> </tr>
<tr><td align="center"> \( A = A^{\dagger} \) </td> <td align="center"> hermitian </td> <td align="center"> \( a_{ij} = a_{ji}^{ * } \) </td> </tr>
<tr><td align="center"> \( A = \left (A^{\dagger} \right )^{-1} \) </td> <td align="center"> unitary </td> <td align="center"> \( \sum_k a_{ik} a_{jk}^{ * } = \sum_k a_{ki}^{ * } a_{kj} = \delta_{ij} \) </td> </tr>
</tbody>
</table>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec4">Some famous Matrices </h2>
<ul>
<li> Diagonal if \( a_{ij}=0 \) for \( i\ne j \)</li>
<li> Upper triangular if \( a_{ij}=0 \) for \( i > j \)</li>
<li> Lower triangular if \( a_{ij}=0 \) for \( i < j \)</li>
<li> Upper Hessenberg if \( a_{ij}=0 \) for \( i > j+1 \)</li>
<li> Lower Hessenberg if \( a_{ij}=0 \) for \( i < j+1 \)</li>
<li> Tridiagonal if \( a_{ij}=0 \) for \( |i -j| > 1 \)</li>
<li> Lower banded with bandwidth \( p \): \( a_{ij}=0 \) for \( i > j+p \)</li>
<li> Upper banded with bandwidth \( p \): \( a_{ij}=0 \) for \( i < j+p \)</li>
<li> Banded, block upper triangular, block lower triangular....</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec5">Basic Matrix Features </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Some Equivalent Statements.</b>
<p>
For an \( N\times N \) matrix \( \mathbf{A} \) the following properties are all equivalent
<ul>
<li> If the inverse of \( \mathbf{A} \) exists, \( \mathbf{A} \) is nonsingular.</li>
<li> The equation \( \mathbf{Ax}=0 \) implies \( \mathbf{x}=0 \).</li>
<li> The rows of \( \mathbf{A} \) form a basis of \( R^N \).</li>
<li> The columns of \( \mathbf{A} \) form a basis of \( R^N \).</li>
<li> \( \mathbf{A} \) is a product of elementary matrices.</li>
<li> \( 0 \) is not eigenvalue of \( \mathbf{A} \).</li>
</ul>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec6">Matrix Handling in C/C++, Static and Dynamical allocation </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Static.</b>
<p>
We have an \( N\times N \) matrix A with \( N=100 \)
In C/C++ this would be defined as
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span> <span style="color: #B00040">int</span> N <span style="color: #666666">=</span> <span style="color: #666666">100</span>;
<span style="color: #B00040">double</span> A[<span style="color: #666666">100</span>][<span style="color: #666666">100</span>];
<span style="color: #408080; font-style: italic">// initialize all elements to zero</span>
<span style="color: #008000; font-weight: bold">for</span>(i<span style="color: #666666">=0</span> ; i <span style="color: #666666">&lt;</span> N ; i<span style="color: #666666">++</span>) {
<span style="color: #008000; font-weight: bold">for</span>(j<span style="color: #666666">=0</span> ; j <span style="color: #666666">&lt;</span> N ; j<span style="color: #666666">++</span>) {
A[i][j] <span style="color: #666666">=</span> <span style="color: #666666">0.0</span>;
</pre></div>
<p>
Note the way the matrix is organized, row-major order.
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec7">Matrix Handling in C/C++ </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Row Major Order, Addition.</b>
<p>
We have \( N\times N \) matrices A, B and C and we wish to
evaluate \( A=B+C \).
$$
\mathbf{A}= \mathbf{B}\pm\mathbf{C} \Longrightarrow a_{ij} = b_{ij}\pm c_{ij},
$$
In C/C++ this would be coded like
<p>
<!-- code=c++ (!bc cppcod) 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">for</span>(i<span style="color: #666666">=0</span> ; i <span style="color: #666666">&lt;</span> N ; i<span style="color: #666666">++</span>) {
<span style="color: #008000; font-weight: bold">for</span>(j<span style="color: #666666">=0</span> ; j <span style="color: #666666">&lt;</span> N ; j<span style="color: #666666">++</span>) {
a[i][j] <span style="color: #666666">=</span> b[i][j]<span style="color: #666666">+</span>c[i][j]
</pre></div>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec8">Matrix Handling in C/C++ </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Row Major Order, Multiplication.</b>
<p>
We have \( N\times N \) matrices A, B and C and we wish to
evaluate \( A=BC \).
$$
\mathbf{A}=\mathbf{BC} \Longrightarrow a_{ij} = \sum_{k=1}^{n} b_{ik}c_{kj},
$$
In C/C++ this would be coded like
<p>
<!-- code=c++ (!bc cppcod) 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">for</span>(i<span style="color: #666666">=0</span> ; i <span style="color: #666666">&lt;</span> N ; i<span style="color: #666666">++</span>) {
<span style="color: #008000; font-weight: bold">for</span>(j<span style="color: #666666">=0</span> ; j <span style="color: #666666">&lt;</span> N ; j<span style="color: #666666">++</span>) {
<span style="color: #008000; font-weight: bold">for</span>(k<span style="color: #666666">=0</span> ; k <span style="color: #666666">&lt;</span> N ; k<span style="color: #666666">++</span>) {
a[i][j]<span style="color: #666666">+=</span>b[i][k]<span style="color: #666666">*</span>c[k][j];
</pre></div>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec9">Dynamic memory allocation in C/C++ </h2>
<p>
At least three possibilities in this course
<ul>
<li> Do it yourself</li>
<li> Use the functions provided in the library package lib.cpp</li>
<li> Use Armadillo <a href="http://arma.sourceforgenet" target="_blank"><tt>http://arma.sourceforgenet</tt></a> (a C++ linear algebra library, discussion both here and at lab).</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec10">Matrix Handling in C/C++, Dynamic Allocation </h2>
<p>
<div class="alert alert-block alert-block alert-text-normal">
<b>Do it yourself.</b>
<p>
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #B00040">int</span> N;
<span style="color: #B00040">double</span> <span style="color: #666666">**</span> A;
A <span style="color: #666666">=</span> <span style="color: #008000; font-weight: bold">new</span> <span style="color: #B00040">double</span><span style="color: #666666">*</span>[N]
<span style="color: #008000; font-weight: bold">for</span> ( i <span style="color: #666666">=</span> <span style="color: #666666">0</span>; i <span style="color: #666666">&lt;</span> N; i<span style="color: #666666">++</span>)
A[i] <span style="color: #666666">=</span> <span style="color: #008000; font-weight: bold">new</span> <span style="color: #B00040">double</span>[N];
</pre></div>
<p>
Always free space when you don't need an array anymore.
<p>
<!-- code=c++ (!bc cppcod) 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">for</span> ( i <span style="color: #666666">=</span> <span style="color: #666666">0</span>; i <span style="color: #666666">&lt;</span> N; i<span style="color: #666666">++</span>)
<span style="color: #008000; font-weight: bold">delete</span>[] A[i];
<span style="color: #008000; font-weight: bold">delete</span>[] A;
</pre></div>
</div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec11">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>
<li> Integer, floating point and complex numbers are supported, as well as a subset of trigonometric and statistics functions. Various matrix decompositions are provided through optional integration with LAPACK, or one of its high performance drop-in replacements (such as the multi-threaded MKL or ACML libraries).</li>
<li> A delayed evaluation approach is employed (at compile-time) to combine several operations into one and reduce (or eliminate) the need for temporaries. This is accomplished through recursive templates and template meta-programming.</li>
<li> Useful for conversion of research code into production environments, or if C++ has been decided as the language of choice, due to speed and/or integration capabilities.</li>
<li> The library is open-source software, and is distributed under a license that is useful in both open-source and commercial/proprietary contexts.</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec12">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&lt;iostream&gt;</span><span style="color: #BC7A00"></span>
<span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&lt;armadillo&gt;</span><span style="color: #BC7A00"></span>
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> std;
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> arma;
<span style="color: #B00040">int</span> <span style="color: #0000FF">main</span>(<span style="color: #B00040">int</span> argc, <span style="color: #B00040">char</span><span style="color: #666666">**</span> argv)
{
mat A <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>,<span style="color: #666666">5</span>);
mat B <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>,<span style="color: #666666">5</span>);
cout <span style="color: #666666">&lt;&lt;</span> A<span style="color: #666666">*</span>B <span style="color: #666666">&lt;&lt;</span> endl;
<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="___sec13">Armadillo, how to compile and install </h2>
<p>
For people using Ubuntu, Debian, Linux Mint, simply go to the synaptic package manager and install
armadillo from there.
You may have to install Lapack as well.
For Mac and Windows users, follow the instructions from the webpage
<a href="http://arma.sourceforge.net" target="_blank"><tt>http://arma.sourceforge.net</tt></a>.
To compile, use for example (linux/ubuntu)
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>c<span style="color: #666666">++</span> <span style="color: #666666">-</span>O2 <span style="color: #666666">-</span>o program.x program.cpp <span style="color: #666666">-</span>larmadillo <span style="color: #666666">-</span>llapack <span style="color: #666666">-</span>lblas
</pre></div>
<p>
where the <code>-l</code> option indicates the library you wish to link to.
<p>
For OS X users you may have to declare the paths to the include files and the libraries as
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span>c<span style="color: #666666">++</span> <span style="color: #666666">-</span>O2 <span style="color: #666666">-</span>o program.x program.cpp <span style="color: #666666">-</span>L<span style="color: #666666">/</span>usr<span style="color: #666666">/</span>local<span style="color: #666666">/</span>lib <span style="color: #666666">-</span>I<span style="color: #666666">/</span>usr<span style="color: #666666">/</span>local<span style="color: #666666">/</span>include <span style="color: #666666">-</span>larmadillo <span style="color: #666666">-</span>llapack <span style="color: #666666">-</span>lblas
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec14">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&lt;iostream&gt;</span><span style="color: #BC7A00"></span>
<span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&quot;armadillo&quot;</span><span style="color: #BC7A00"></span>
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> arma;
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> std;
<span style="color: #B00040">int</span> <span style="color: #0000FF">main</span>(<span style="color: #B00040">int</span> argc, <span style="color: #B00040">char</span><span style="color: #666666">**</span> argv)
{
<span style="color: #408080; font-style: italic">// directly specify the matrix size (elements are uninitialised)</span>
mat A(<span style="color: #666666">2</span>,<span style="color: #666666">3</span>);
<span style="color: #408080; font-style: italic">// .n_rows = number of rows (read only)</span>
<span style="color: #408080; font-style: italic">// .n_cols = number of columns (read only)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;A.n_rows = &quot;</span> <span style="color: #666666">&lt;&lt;</span> A.n_rows <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;A.n_cols = &quot;</span> <span style="color: #666666">&lt;&lt;</span> A.n_cols <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// directly access an element (indexing starts at 0)</span>
A(<span style="color: #666666">1</span>,<span style="color: #666666">2</span>) <span style="color: #666666">=</span> <span style="color: #666666">456.0</span>;
A.print(<span style="color: #BA2121">&quot;A:&quot;</span>);
<span style="color: #408080; font-style: italic">// scalars are treated as a 1x1 matrix,</span>
<span style="color: #408080; font-style: italic">// hence the code below will set A to have a size of 1x1</span>
A <span style="color: #666666">=</span> <span style="color: #666666">5.0</span>;
A.print(<span style="color: #BA2121">&quot;A:&quot;</span>);
<span style="color: #408080; font-style: italic">// if you want a matrix with all elements set to a particular value</span>
<span style="color: #408080; font-style: italic">// the .fill() member function can be used</span>
A.set_size(<span style="color: #666666">3</span>,<span style="color: #666666">3</span>);
A.fill(<span style="color: #666666">5.0</span>); A.print(<span style="color: #BA2121">&quot;A:&quot;</span>);
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec15">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span> mat B;
<span style="color: #408080; font-style: italic">// endr indicates &quot;end of row&quot;</span>
B <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.555950</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.274690</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.540605</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.798938</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.108929</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.830123</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.891726</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.895283</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.948014</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.973234</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.216504</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.883152</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.023787</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.675382</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.231751</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.450332</span> <span style="color: #666666">&lt;&lt;</span> endr;
<span style="color: #408080; font-style: italic">// print to the cout stream</span>
<span style="color: #408080; font-style: italic">// with an optional string before the contents of the matrix</span>
B.print(<span style="color: #BA2121">&quot;B:&quot;</span>);
<span style="color: #408080; font-style: italic">// the &lt;&lt; operator can also be used to print the matrix</span>
<span style="color: #408080; font-style: italic">// to an arbitrary stream (cout in this case)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;B:&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl <span style="color: #666666">&lt;&lt;</span> B <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// save to disk</span>
B.save(<span style="color: #BA2121">&quot;B.txt&quot;</span>, raw_ascii);
<span style="color: #408080; font-style: italic">// load from disk</span>
mat C;
C.load(<span style="color: #BA2121">&quot;B.txt&quot;</span>);
C <span style="color: #666666">+=</span> <span style="color: #666666">2.0</span> <span style="color: #666666">*</span> B;
C.print(<span style="color: #BA2121">&quot;C:&quot;</span>);
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec16">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) 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">// submatrix types:</span>
<span style="color: #408080; font-style: italic">//</span>
<span style="color: #408080; font-style: italic">// .submat(first_row, first_column, last_row, last_column)</span>
<span style="color: #408080; font-style: italic">// .row(row_number)</span>
<span style="color: #408080; font-style: italic">// .col(column_number)</span>
<span style="color: #408080; font-style: italic">// .cols(first_column, last_column)</span>
<span style="color: #408080; font-style: italic">// .rows(first_row, last_row)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;C.submat(0,0,3,1) =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> C.submat(<span style="color: #666666">0</span>,<span style="color: #666666">0</span>,<span style="color: #666666">3</span>,<span style="color: #666666">1</span>) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// generate the identity matrix</span>
mat D <span style="color: #666666">=</span> eye<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">4</span>,<span style="color: #666666">4</span>);
D.submat(<span style="color: #666666">0</span>,<span style="color: #666666">0</span>,<span style="color: #666666">3</span>,<span style="color: #666666">1</span>) <span style="color: #666666">=</span> C.cols(<span style="color: #666666">1</span>,<span style="color: #666666">2</span>);
D.print(<span style="color: #BA2121">&quot;D:&quot;</span>);
<span style="color: #408080; font-style: italic">// transpose</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;trans(B) =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> trans(B) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// maximum from each column (traverse along rows)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;max(B) =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> max(B) <span style="color: #666666">&lt;&lt;</span> endl;
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec17">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) 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">// maximum from each row (traverse along columns)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;max(B,1) =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> max(B,<span style="color: #666666">1</span>) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// maximum value in B</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;max(max(B)) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> max(max(B)) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// sum of each column (traverse along rows)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;sum(B) =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> sum(B) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// sum of each row (traverse along columns)</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;sum(B,1) =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> sum(B,<span style="color: #666666">1</span>) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// sum of all elements</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;sum(sum(B)) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> sum(sum(B)) <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;accu(B) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> accu(B) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// trace = sum along diagonal</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;trace(B) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> trace(B) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// random matrix -- values are uniformly distributed in the [0,1] interval</span>
mat E <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">4</span>,<span style="color: #666666">4</span>);
E.print(<span style="color: #BA2121">&quot;E:&quot;</span>);
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec18">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) 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">// row vectors are treated like a matrix with one row</span>
rowvec r;
r <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.59499</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.88807</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.88532</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.19968</span>;
r.print(<span style="color: #BA2121">&quot;r:&quot;</span>);
<span style="color: #408080; font-style: italic">// column vectors are treated like a matrix with one column</span>
colvec q;
q <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.81114</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.06256</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.95989</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.73628</span>;
q.print(<span style="color: #BA2121">&quot;q:&quot;</span>);
<span style="color: #408080; font-style: italic">// dot or inner product</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;as_scalar(r*q) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> as_scalar(r<span style="color: #666666">*</span>q) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// outer product</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;q*r =&quot;</span> <span style="color: #666666">&lt;&lt;</span> endl;
cout <span style="color: #666666">&lt;&lt;</span> q<span style="color: #666666">*</span>r <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// sum of three matrices (no temporary matrices are created)</span>
mat F <span style="color: #666666">=</span> B <span style="color: #666666">+</span> C <span style="color: #666666">+</span> D;
F.print(<span style="color: #BA2121">&quot;F:&quot;</span>);
<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="___sec19">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&lt;iostream&gt;</span><span style="color: #BC7A00"></span>
<span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&quot;armadillo&quot;</span><span style="color: #BC7A00"></span>
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> arma;
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> std;
<span style="color: #B00040">int</span> <span style="color: #0000FF">main</span>(<span style="color: #B00040">int</span> argc, <span style="color: #B00040">char</span><span style="color: #666666">**</span> argv)
{
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;Armadillo version: &quot;</span> <span style="color: #666666">&lt;&lt;</span> arma_version<span style="color: #666666">::</span>as_string() <span style="color: #666666">&lt;&lt;</span> endl;
mat A;
A <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.165300</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.454037</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.995795</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.124098</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.047084</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.688782</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.036549</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.552848</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.937664</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.866401</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.348740</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.479388</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.506228</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.145673</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.491547</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.148678</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.682258</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.571154</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.874724</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.444632</span> <span style="color: #666666">&lt;&lt;</span> endr
<span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.245726</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.595218</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.409327</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.367827</span> <span style="color: #666666">&lt;&lt;</span> <span style="color: #666666">0.385736</span> <span style="color: #666666">&lt;&lt;</span> endr;
A.print(<span style="color: #BA2121">&quot;A =&quot;</span>);
<span style="color: #408080; font-style: italic">// determinant</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;det(A) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> det(A) <span style="color: #666666">&lt;&lt;</span> endl;
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec20">Armadillo, simple examples </h2>
<p>
<!-- code=c++ (!bc cppcod) 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">// inverse</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;inv(A) = &quot;</span> <span style="color: #666666">&lt;&lt;</span> endl <span style="color: #666666">&lt;&lt;</span> inv(A) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #B00040">double</span> k <span style="color: #666666">=</span> <span style="color: #666666">1.23</span>;
mat B <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>,<span style="color: #666666">5</span>);
mat C <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>,<span style="color: #666666">5</span>);
rowvec r <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>rowvec<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>);
colvec q <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>colvec<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>);
<span style="color: #408080; font-style: italic">// examples of some expressions</span>
<span style="color: #408080; font-style: italic">// for which optimised implementations exist</span>
<span style="color: #408080; font-style: italic">// optimised implementation of a trinary expression</span>
<span style="color: #408080; font-style: italic">// that results in a scalar</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;as_scalar( r*inv(diagmat(B))*q ) = &quot;</span>;
cout <span style="color: #666666">&lt;&lt;</span> as_scalar( r<span style="color: #666666">*</span>inv(diagmat(B))<span style="color: #666666">*</span>q ) <span style="color: #666666">&lt;&lt;</span> endl;
<span style="color: #408080; font-style: italic">// example of an expression which is optimised</span>
<span style="color: #408080; font-style: italic">// as a call to the dgemm() function in BLAS:</span>
cout <span style="color: #666666">&lt;&lt;</span> <span style="color: #BA2121">&quot;k*trans(B)*C = &quot;</span> <span style="color: #666666">&lt;&lt;</span> endl <span style="color: #666666">&lt;&lt;</span> k<span style="color: #666666">*</span>trans(B)<span style="color: #666666">*</span>C;
<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="___sec21">Gaussian Elimination </h2>
<p>
We start with the linear set of equations
$$
\mathbf{A}\mathbf{x} = \mathbf{w}.
$$
We assume also that the matrix \( \mathbf{A} \) is non-singular and that the
matrix elements along the diagonal satisfy \( a_{ii} \ne 0 \). Simple \( 4\times 4 \) example
$$
\begin{bmatrix}
a_{11}& a_{12} &a_{13}& a_{14}\\
a_{21}& a_{22} &a_{23}& a_{24}\\
a_{31}& a_{32} &a_{33}& a_{34}\\
a_{41}& a_{42} &a_{43}& a_{44}\\
\end{bmatrix} \begin{bmatrix}
x_1\\
x_2\\
x_3 \\
x_4 \\
\end{bmatrix}
=\begin{bmatrix}
w_1\\
w_2\\
w_3 \\
w_4\\
\end{bmatrix}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec22">Gaussian Elimination </h2>
or
$$
\begin{align}
a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=&w_1 \nonumber \\
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=&w_2 \nonumber \\
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=&w_3 \nonumber \\
a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=&w_4. \nonumber
\end{align}
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec23">Gaussian Elimination </h2>
<p>
The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown \( x_1 \)
from the remaining \( n-1 \) equations. Then we use the new second equation to eliminate the second unknown
\( x_2 \) from the remaining \( n-2 \) equations. With \( n-1 \) such eliminations
we obtain a so-called upper triangular set of equations of the form
$$
\begin{align}
b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=&y_1 \nonumber \\
b_{22}x_2 + b_{23}x_3 + b_{24}x_4=&y_2 \nonumber \\
b_{33}x_3 + b_{34}x_4=&y_3 \nonumber \\
b_{44}x_4=&y_4. \nonumber
\label{eq:gaussbacksub}
\end{align}
$$
We can solve this system of equations recursively starting from \( x_n \) (in our case \( x_4 \)) and proceed with
what is called a backward substitution.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec24">Gaussian Elimination </h2>
This process can be expressed mathematically as
$$
\begin{equation}
x_m = \frac{1}{b_{mm}}\left(y_m-\sum_{k=m+1}^nb_{mk}x_k\right)\quad m=n-1,n-2,\dots,1.
\label{_auto1}
\end{equation}
$$
To arrive at such an upper triangular system of equations, we start by eliminating
the unknown \( x_1 \) for \( j=2,n \). We achieve this by multiplying the first equation by \( a_{j1}/a_{11} \) and then subtract
the result from the $j$th equation. We assume obviously that \( a_{11}\ne 0 \) and that
\( \mathbf{A} \) is not singular.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec25">Gaussian Elimination </h2>
<p>
Our actual \( 4\times 4 \) example reads after the first operation
$$
\begin{bmatrix}
a_{11}& a_{12} &a_{13}& a_{14}\\
0& (a_{22}-\frac{a_{21}a_{12}}{a_{11}}) &(a_{23}-\frac{a_{21}a_{13}}{a_{11}}) & (a_{24}-\frac{a_{21}a_{14}}{a_{11}})\\
0& (a_{32}-\frac{a_{31}a_{12}}{a_{11}})& (a_{33}-\frac{a_{31}a_{13}}{a_{11}})& (a_{34}-\frac{a_{31}a_{14}}{a_{11}})\\
0&(a_{42}-\frac{a_{41}a_{12}}{a_{11}}) &(a_{43}-\frac{a_{41}a_{13}}{a_{11}}) & (a_{44}-\frac{a_{41}a_{14}}{a_{11}}) \\
\end{bmatrix} \begin{bmatrix}
x_1\\
x_2\\
x_3 \\
x_4 \\
\end{bmatrix}
=\begin{bmatrix}
y_1\\
w_2^{(2)}\\
w_3^{(2)} \\
w_4^{(2)}\\
\end{bmatrix},
$$
or
$$
\begin{align}
b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=&y_1 \nonumber \\
a^{(2)}_{22}x_2 + a^{(2)}_{23}x_3 + a^{(2)}_{24}x_4=&w^{(2)}_2 \nonumber \\
a^{(2)}_{32}x_2 + a^{(2)}_{33}x_3 + a^{(2)}_{34}x_4=&w^{(2)}_3 \nonumber \\
a^{(2)}_{42}x_2 + a^{(2)}_{43}x_3 + a^{(2)}_{44}x_4=&w^{(2)}_4, \nonumber \\
\label{_auto2}
\end{align}
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec26">Gaussian Elimination </h2>
<p>
The new coefficients are
$$
\begin{equation}
b_{1k} = a_{1k}^{(1)} \quad k=1,\dots,n,
\label{_auto3}
\end{equation}
$$
where each \( a_{1k}^{(1)} \) is equal to the original \( a_{1k} \) element. The other coefficients are
$$
\begin{equation}
a_{jk}^{(2)} = a_{jk}^{(1)}-\frac{a_{j1}^{(1)}a_{1k}^{(1)}}{a_{11}^{(1)}} \quad j,k=2,\dots,n,
\label{_auto4}
\end{equation}
$$
with a new right-hand side given by
$$
\begin{equation}
y_{1}=w_1^{(1)}, \quad w_j^{(2)} =w_j^{(1)}-\frac{a_{j1}^{(1)}w_1^{(1)}}{a_{11}^{(1)}} \quad j=2,\dots,n.
\label{_auto5}
\end{equation}
$$
We have also set \( w_1^{(1)}=w_1 \), the original vector element.
We see that the system of unknowns \( x_1,\dots,x_n \) is transformed into an \( (n-1)\times (n-1) \) problem.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec27">Gaussian Elimination </h2>
<p>
This step is called forward substitution.
Proceeding with these substitutions, we obtain the
general expressions for the new coefficients
$$
\begin{equation}
a_{jk}^{(m+1)} = a_{jk}^{(m)}-\frac{a_{jm}^{(m)}a_{mk}^{(m)}}{a_{mm}^{(m)}} \quad j,k=m+1,\dots,n,
\label{_auto6}
\end{equation}
$$
with \( m=1,\dots,n-1 \) and a
right-hand side given by
$$
\begin{equation}
w_j^{(m+1)} =w_j^{(m)}-\frac{a_{jm}^{(m)}w_m^{(m)}}{a_{mm}^{(m)}}\quad j=m+1,\dots,n.
\label{_auto7}
\end{equation}
$$
This set of \( n-1 \) elimations leads us to an equations which is solved by back substitution.
If the arithmetics is exact and the matrix \( \mathbf{A} \) is not singular, then the computed answer will be exact.
<p>
Even though the matrix elements along the diagonal are not zero,
numerically small numbers may appear and subsequent divisions may lead to large numbers, which, if added
to a small number may yield losses of precision. Suppose for example that our first division in \( (a_{22}-a_{21}a_{12}/a_{11}) \)
results in \( -10^{-7} \) and that \( a_{22} \) is one.
one. We are then
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="___sec28">Linear Algebra Methods </h2>
<ul>
<li> Gaussian elimination, \( O(2/3n^3) \) flops, general matrix</li>
<li> LU decomposition, upper triangular and lower tridiagonal matrices, \( O(2/3n^3) \) flops, general matrix. Get easily the inverse, determinant and can solve linear equations with back-substitution only, \( O(n^2) \) flops</li>
<li> Cholesky decomposition. Real symmetric or hermitian positive definite matrix, \( O(1/3n^3) \) flops.</li>
<li> Tridiagonal linear systems, important for differential equations. Normally positive definite and non-singular. \( O(8n) \) flops for symmetric. Special case of banded matrices.</li>
<li> Singular value decomposition</li>
<li> the QR method will be discussed in chapter 7 in connection with eigenvalue systems. \( O(4/3n^3) \) flops.</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec29">LU Decomposition </h2>
<p>
The LU decomposition method means that we can rewrite
this matrix as the product of two matrices \( \mathbf{L} \) and \( \mathbf{U} \)
where
$$
\begin{bmatrix}
a_{11} & a_{12} & a_{13} & a_{14} \\
a_{21} & a_{22} & a_{23} & a_{24} \\
a_{31} & a_{32} & a_{33} & a_{34} \\
a_{41} & a_{42} & a_{43} & a_{44}
\end{bmatrix}
= \begin{bmatrix}
1 & 0 & 0 & 0 \\
l_{21} & 1 & 0 & 0 \\
l_{31} & l_{32} & 1 & 0 \\
l_{41} & l_{42} & l_{43} & 1
\end{bmatrix}
\begin{bmatrix}
u_{11} & u_{12} & u_{13} & u_{14} \\
0 & u_{22} & u_{23} & u_{24} \\
0 & 0 & u_{33} & u_{34} \\
0 & 0 & 0 & u_{44}
\end{bmatrix}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec30">LU Decomposition </h2>
<p>
LU decomposition forms the backbone of other algorithms in linear algebra, such as the
solution of linear equations given by
$$
\begin{align}
a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=&w_1 \nonumber \\
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=&w_2 \nonumber \\
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=&w_3 \nonumber \\
a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=&w_4. \nonumber
\end{align}
$$
The above set of equations is conveniently solved by using LU decomposition as an intermediate step.
<p>
The matrix \( \mathbf{A}\in \mathbb{R}^{n\times n} \) has an LU factorization if the determinant
is different from zero. If the LU factorization exists and \( \mathbf{A} \) is non-singular, then the LU factorization
is unique and the determinant is given by
$$
det\{\mathbf{A}\}=det\{\mathbf{LU}\}= det\{\mathbf{L}\}det\{\mathbf{U}\}=u_{11}u_{22}\dots u_{nn}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec31">LU Decomposition, why? </h2>
<p>
There are at least three main advantages with LU decomposition compared with standard Gaussian elimination:
<ul>
<li> It is straightforward to compute the determinant of a matrix</li>
<li> If we have to solve sets of linear equations with the same matrix but with different vectors \( \mathbf{y} \), the number of FLOPS is of the order \( n^3 \).</li>
<li> The inverse is such an operation</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec32">LU Decomposition, linear equations </h2>
<p>
With the LU decomposition it is rather
simple to solve a system of linear equations
$$
\begin{align}
a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=&w_1 \nonumber \\
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=&w_2 \nonumber \\
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=&w_3 \nonumber \\
a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=&w_4. \nonumber
\end{align}
$$
<p>
This can be written in matrix form as
$$ \mathbf{Ax}=\mathbf{w}. $$
<p>
where \( \mathbf{A} \) and \( \mathbf{w} \) are known and we have to solve for
\( \mathbf{x} \). Using the LU dcomposition we write
$$ \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="___sec33">LU Decomposition, linear equations </h2>
<p>
The previous equation can be calculated in two steps
$$ \mathbf{L} \mathbf{y} = \mathbf{w};\qquad \mathbf{Ux}=\mathbf{y}. $$
<p>
To show that this is correct we use to the LU decomposition
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
since the diagonals of \( \mathbf{L} \) equal 1) we can use the inverse of
\( \mathbf{L} \) to obtain
$$
\mathbf{Ux}=\mathbf{L^{-1}w}=\mathbf{y},
$$
which yields the intermediate step
$$
\mathbf{L^{-1}w}=\mathbf{y}
$$
and as soon as we have \( \mathbf{y} \) we can obtain \( \mathbf{x} \)
through \( \mathbf{Ux}=\mathbf{y} \).
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec34">LU Decomposition, why? </h2>
<p>
For our four-dimentional example this takes the form
$$
\begin{align}
y_1=&w_1 \nonumber\\
l_{21}y_1 + y_2=&w_2\nonumber \\
l_{31}y_1 + l_{32}y_2 + y_3 =&w_3\nonumber \\
l_{41}y_1 + l_{42}y_2 + l_{43}y_3 + y_4=&w_4. \nonumber
\end{align}
$$
<p>
and
$$
\begin{align}
u_{11}x_1 +u_{12}x_2 +u_{13}x_3 + u_{14}x_4=&y_1 \nonumber\\
u_{22}x_2 + u_{23}x_3 + u_{24}x_4=&y_2\nonumber \\
u_{33}x_3 + u_{34}x_4=&y_3\nonumber \\
u_{44}x_4=&y_4 \nonumber
\end{align}
$$
<p>
This example shows the basis for the algorithm
needed to solve the set of \( n \) linear equations.
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec35">LU Decomposition, linear equations </h2>
<p>
The algorithm goes as follows
<ul>
<li> Set up the matrix \( \bf A \) and the vector \( \bf w \) with their correct dimensions. This determines the dimensionality of the unknown vector \( \bf x \).</li>
<li> Then LU decompose the matrix \( \bf A \) through a call to the function <code>ludcmp(double a, int n, int indx, double &d)</code>. This functions returns the LU decomposed matrix \( \bf A \), its determinant and the vector indx which keeps track of the number of interchanges of rows. If the determinant is zero, the solution is malconditioned.</li>
<li> Thereafter you call the function <code>lubksb(double a, int n, int indx, double w)</code> which uses the LU decomposed matrix \( \bf A \) and the vector \( \bf w \) and returns \( \bf x \) in the same place as \( \bf w \). Upon exit the original content in \( \bf w \) is destroyed. If you wish to keep this information, you should make a backup of it in your calling function.</li>
</ul>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec36">LU Decomposition, the inverse of a matrix </h2>
<p>
If the inverse exists then
$$
\mathbf{A}^{-1}\mathbf{A}=\mathbf{I},
$$
the identity matrix. With an LU decomposed matrix we can rewrite the last equation as
$$
\mathbf{LU}\mathbf{A}^{-1}=\mathbf{I}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec37">LU Decomposition, the inverse of a matrix </h2>
<p>
If we assume that the first column (that is column 1) of the inverse matrix
can be written as a vector with unknown entries
$$
\mathbf{A}_1^{-1}= \begin{bmatrix}
a_{11}^{-1} \\
a_{21}^{-1} \\
\dots \\
a_{n1}^{-1} \\
\end{bmatrix},
$$
then we have a linear set of equations
$$
\mathbf{LU}\begin{bmatrix}
a_{11}^{-1} \\
a_{21}^{-1} \\
\dots \\
a_{n1}^{-1} \\
\end{bmatrix} =\begin{bmatrix}
1 \\
0 \\
\dots \\
0 \\
\end{bmatrix}.
$$
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
<h2 id="___sec38">LU Decomposition, the inverse </h2>
<p>
In a similar way we can compute the unknow entries of the second column,
$$
\mathbf{LU}\begin{bmatrix}
a_{12}^{-1} \\
a_{22}^{-1} \\
\dots \\
a_{n2}^{-1} \\
\end{bmatrix}=\begin{bmatrix}
0 \\
1 \\
\dots \\
0 \\
\end{bmatrix},
$$
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="___sec39"><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" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&lt;iostream&gt;</span><span style="color: #BC7A00"></span>
<span style="color: #BC7A00">#include</span> <span style="color: #408080; font-style: italic">&quot;armadillo&quot;</span><span style="color: #BC7A00"></span>
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> arma;
<span style="color: #008000; font-weight: bold">using</span> <span style="color: #008000; font-weight: bold">namespace</span> std;
<span style="color: #B00040">int</span> <span style="color: #0000FF">main</span>()
{
mat A <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>mat<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>,<span style="color: #666666">5</span>);
vec b <span style="color: #666666">=</span> randu<span style="color: #666666">&lt;</span>vec<span style="color: #666666">&gt;</span>(<span style="color: #666666">5</span>);
A.print(<span style="color: #BA2121">&quot;A =&quot;</span>);
b.print(<span style="color: #BA2121">&quot;b=&quot;</span>);
<span style="color: #408080; font-style: italic">// solve Ax = b</span>
vec x <span style="color: #666666">=</span> solve(A,b);
<span style="color: #408080; font-style: italic">// print x</span>
x.print(<span style="color: #BA2121">&quot;x=&quot;</span>);
<span style="color: #408080; font-style: italic">// find LU decomp of A, if needed, P is the permutation matrix</span>
mat L, U;
lu(L,U,A);
<span style="color: #408080; font-style: italic">// print l</span>
L.print(<span style="color: #BA2121">&quot; L= &quot;</span>);
<span style="color: #408080; font-style: italic">// print U</span>
U.print(<span style="color: #BA2121">&quot; U= &quot;</span>);
<span style="color: #408080; font-style: italic">//Check that A = LU</span>
(A<span style="color: #666666">-</span>L<span style="color: #666666">*</span>U).print(<span style="color: #BA2121">&quot;Test of LU decomposition&quot;</span>);
<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="___sec40">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="___sec41">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="___sec42">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="___sec43">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="___sec44">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="___sec45">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="___sec46">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="___sec47">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="___sec48">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 --------------- -->
<center style="font-size:80%">
<!-- copyright --> &copy; 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license
</center>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink