lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
20fa1096b61c60115b68365e0bfbfbc2a842df2f
FYS-STK4155/doc/LectureNotes/_build/html/linalg.html
T
Morten Hjorth-Jensen 20fa1096b6 udpate
2023-10-02 07:01:55 +02:00

1526 lines
89 KiB
HTML
Raw Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!DOCTYPE html>
<html>
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1.0" />
<title>2. Linear Algebra, Handling of Arrays and more Python Features &#8212; Applied Data Analysis and Machine Learning</title>
<link href="_static/css/theme.css" rel="stylesheet">
<link href="_static/css/index.ff1ffe594081f20da1ef19478df9384b.css" rel="stylesheet">
<link rel="stylesheet"
href="_static/vendor/fontawesome/5.13.0/css/all.min.css">
<link rel="preload" as="font" type="font/woff2" crossorigin
href="_static/vendor/fontawesome/5.13.0/webfonts/fa-solid-900.woff2">
<link rel="preload" as="font" type="font/woff2" crossorigin
href="_static/vendor/fontawesome/5.13.0/webfonts/fa-brands-400.woff2">
<link rel="stylesheet" type="text/css" href="_static/pygments.css" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-book-theme.css?digest=c3fdc42140077d1ad13ad2f1588a4309" />
<link rel="stylesheet" type="text/css" href="_static/togglebutton.css" />
<link rel="stylesheet" type="text/css" href="_static/copybutton.css" />
<link rel="stylesheet" type="text/css" href="_static/mystnb.css" />
<link rel="stylesheet" type="text/css" href="_static/sphinx-thebe.css" />
<link rel="stylesheet" type="text/css" href="_static/panels-main.c949a650a448cc0ae9fd3441c0e17fb0.css" />
<link rel="stylesheet" type="text/css" href="_static/panels-variables.06eb56fa6e07937060861dad626602ad.css" />
<link rel="preload" as="script" href="_static/js/index.be7d3bbb2ef33a8344ce.js">
<script data-url_root="./" id="documentation_options" src="_static/documentation_options.js"></script>
<script src="_static/jquery.js"></script>
<script src="_static/underscore.js"></script>
<script src="_static/doctools.js"></script>
<script src="_static/clipboard.min.js"></script>
<script src="_static/copybutton.js"></script>
<script>let toggleHintShow = 'Click to show';</script>
<script>let toggleHintHide = 'Click to hide';</script>
<script>let toggleOpenOnPrint = 'true';</script>
<script src="_static/togglebutton.js"></script>
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown, .tag_hide_input div.cell_input, .tag_hide-input div.cell_input, .tag_hide_output div.cell_output, .tag_hide-output div.cell_output, .tag_hide_cell.cell, .tag_hide-cell.cell';</script>
<script src="_static/sphinx-book-theme.d59cb220de22ca1c485ebbdc042f0030.js"></script>
<script>const THEBE_JS_URL = "https://unpkg.com/thebe@0.8.2/lib/index.js"
const thebe_selector = ".thebe,.cell"
const thebe_selector_input = "pre"
const thebe_selector_output = ".output, .cell_output"
</script>
<script async="async" src="_static/sphinx-thebe.js"></script>
<script>window.MathJax = {"options": {"processHtmlClass": "tex2jax_process|mathjax_process|math|output_area"}}</script>
<script defer="defer" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
<link rel="index" title="Index" href="genindex.html" />
<link rel="search" title="Search" href="search.html" />
<link rel="next" title="3. Linear Regression" href="chapter1.html" />
<link rel="prev" title="1. Elements of Probability Theory and Statistical Data Analysis" href="statistics.html" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<meta name="docsearch:language" content="None">
<!-- Google Analytics -->
</head>
<body data-spy="scroll" data-target="#bd-toc-nav" data-offset="80">
<div class="container-fluid" id="banner"></div>
<div class="container-xl">
<div class="row">
<div class="col-12 col-md-3 bd-sidebar site-navigation show" id="site-navigation">
<div class="navbar-brand-box">
<a class="navbar-brand text-wrap" href="index.html">
<!-- `logo` is deprecated in Sphinx 4.0, so remove this when we stop supporting 3 -->
<img src="_static/logo.png" class="logo" alt="logo">
<h1 class="site-logo" id="site-title">Applied Data Analysis and Machine Learning</h1>
</a>
</div><form class="bd-search d-flex align-items-center" action="search.html" method="get">
<i class="icon fas fa-search"></i>
<input type="search" class="form-control" name="q" id="search-input" placeholder="Search this book..." aria-label="Search this book..." autocomplete="off" >
</form><nav class="bd-links" id="bd-docs-nav" aria-label="Main">
<div class="bd-toc-item active">
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="intro.html">
Applied Data Analysis and Machine Learning
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
About the course
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="schedule.html">
Teaching schedule with links to material
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="teachers.html">
Teachers and Grading
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="textbooks.html">
Textbooks
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
Review of Statistics with Resampling Techniques and Linear Algebra
</span>
</p>
<ul class="current nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="statistics.html">
1. Elements of Probability Theory and Statistical Data Analysis
</a>
</li>
<li class="toctree-l1 current active">
<a class="current reference internal" href="#">
2. Linear Algebra, Handling of Arrays and more Python Features
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
From Regression to Support Vector Machines
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="chapter1.html">
3. Linear Regression
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter2.html">
4. Ridge and Lasso Regression
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter3.html">
5. Resampling Methods
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter4.html">
6. Logistic Regression
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapteroptimization.html">
7. Optimization, the central part of any Machine Learning algortithm
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter5.html">
8. Support Vector Machines, overarching aims
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
Decision Trees, Ensemble Methods and Boosting
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="chapter6.html">
9. Decision trees, overarching aims
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter7.html">
10. Ensemble Methods: From a Single Tree to Many Trees and Extreme Boosting, Meet the Jungle of Methods
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
Dimensionality Reduction
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="chapter8.html">
11. Basic ideas of the Principal Component Analysis (PCA)
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="clustering.html">
12. Clustering and Unsupervised Learning
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
Deep Learning Methods
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="chapter9.html">
13. Neural networks
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter10.html">
14. Building a Feed Forward Neural Network
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter11.html">
15. Solving Differential Equations with Deep Learning
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter12.html">
16. Convolutional Neural Networks
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="chapter13.html">
17. Recurrent neural networks: Overarching view
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
Weekly material, notes and exercises
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek34.html">
Exercises week 34
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week34.html">
Week 34: Introduction to the course, Logistics and Practicalities
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek35.html">
Exercises week 35
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week35.html">
Week 35: From Ordinary Linear Regression to Ridge and Lasso Regression
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek36.html">
Exercises week 36
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week36.html">
Week 36: Statistical interpretation of Linear Regression and Resampling techniques
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek37.html">
Exercises week 37
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week37.html">
Week 37: Statistical interpretations and Resampling Methods
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek38.html">
Exercises week 38
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week38.html">
Week 38: Logistic Regression and Optimization
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="exercisesweek39.html">
Exercises week 39
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week39.html">
Week 39: Optimization and Gradient Methods
</a>
</li>
<li class="toctree-l1">
<a class="reference internal" href="week40.html">
Week 40: Gradient descent methods (continued) and start Neural networks
</a>
</li>
</ul>
<p aria-level="2" class="caption" role="heading">
<span class="caption-text">
Projects
</span>
</p>
<ul class="nav bd-sidenav">
<li class="toctree-l1">
<a class="reference internal" href="project1.html">
Project 1 on Machine Learning, deadline October 9 (midnight), 2023
</a>
</li>
</ul>
</div>
</nav> <!-- To handle the deprecated key -->
<div class="navbar_extra_footer">
Powered by <a href="https://jupyterbook.org">Jupyter Book</a>
</div>
</div>
<main class="col py-md-3 pl-md-4 bd-content overflow-auto" role="main">
<div class="topbar container-xl fixed-top">
<div class="topbar-contents row">
<div class="col-12 col-md-3 bd-topbar-whitespace site-navigation show"></div>
<div class="col pl-md-4 topbar-main">
<button id="navbar-toggler" class="navbar-toggler ml-0" type="button" data-toggle="collapse"
data-toggle="tooltip" data-placement="bottom" data-target=".site-navigation" aria-controls="navbar-menu"
aria-expanded="true" aria-label="Toggle navigation" aria-controls="site-navigation"
title="Toggle navigation" data-toggle="tooltip" data-placement="left">
<i class="fas fa-bars"></i>
<i class="fas fa-arrow-left"></i>
<i class="fas fa-arrow-up"></i>
</button>
<div class="dropdown-buttons-trigger">
<button id="dropdown-buttons-trigger" class="btn btn-secondary topbarbtn" aria-label="Download this page"><i
class="fas fa-download"></i></button>
<div class="dropdown-buttons">
<!-- ipynb file if we had a myst markdown file -->
<!-- Download raw file -->
<a class="dropdown-buttons" href="_sources/linalg.ipynb"><button type="button"
class="btn btn-secondary topbarbtn" title="Download source file" data-toggle="tooltip"
data-placement="left">.ipynb</button></a>
<!-- Download PDF via print -->
<button type="button" id="download-print" class="btn btn-secondary topbarbtn" title="Print to PDF"
onclick="printPdf(this)" data-toggle="tooltip" data-placement="left">.pdf</button>
</div>
</div>
<!-- Source interaction buttons -->
<!-- Full screen (wrap in <a> to have style consistency -->
<a class="full-screen-button"><button type="button" class="btn btn-secondary topbarbtn" data-toggle="tooltip"
data-placement="bottom" onclick="toggleFullScreen()" aria-label="Fullscreen mode"
title="Fullscreen mode"><i
class="fas fa-expand"></i></button></a>
<!-- Launch buttons -->
</div>
<!-- Table of contents -->
<div class="d-none d-md-block col-md-2 bd-toc show noprint">
<div class="tocsection onthispage pt-5 pb-3">
<i class="fas fa-list"></i> Contents
</div>
<nav id="bd-toc-nav" aria-label="Page">
<ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#introduction">
2.1. Introduction
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#important-matrix-and-vector-handling-packages">
2.2. Important Matrix and vector handling packages
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#basic-matrix-features">
2.3. Basic Matrix Features
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#some-famous-matrices">
2.3.1. Some famous Matrices
</a>
</li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#numpy-and-arrays">
2.4. Numpy and arrays
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#other-matrix-and-vector-operations">
2.5. Other Matrix and Vector Operations
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#gaussian-elimination">
2.6. Gaussian Elimination
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#lu-decomposition-the-inverse-of-a-matrix">
2.6.1. LU Decomposition, the inverse of a matrix
</a>
</li>
</ul>
</li>
</ul>
</nav>
</div>
</div>
</div>
<div id="main-content" class="row">
<div class="col-12 col-md-9 pl-md-3 pr-md-0">
<!-- Table of contents that is only displayed when printing the page -->
<div id="jb-print-docs-body" class="onlyprint">
<h1>Linear Algebra, Handling of Arrays and more Python Features</h1>
<!-- Table of contents -->
<div id="print-main-content">
<div id="jb-print-toc">
<div>
<h2> Contents </h2>
</div>
<nav aria-label="Page">
<ul class="visible nav section-nav flex-column">
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#introduction">
2.1. Introduction
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#important-matrix-and-vector-handling-packages">
2.2. Important Matrix and vector handling packages
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#basic-matrix-features">
2.3. Basic Matrix Features
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#some-famous-matrices">
2.3.1. Some famous Matrices
</a>
</li>
</ul>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#numpy-and-arrays">
2.4. Numpy and arrays
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#other-matrix-and-vector-operations">
2.5. Other Matrix and Vector Operations
</a>
</li>
<li class="toc-h2 nav-item toc-entry">
<a class="reference internal nav-link" href="#gaussian-elimination">
2.6. Gaussian Elimination
</a>
<ul class="nav section-nav flex-column">
<li class="toc-h3 nav-item toc-entry">
<a class="reference internal nav-link" href="#lu-decomposition-the-inverse-of-a-matrix">
2.6.1. LU Decomposition, the inverse of a matrix
</a>
</li>
</ul>
</li>
</ul>
</nav>
</div>
</div>
</div>
<div>
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
doconce format html linalg.do.txt --><div class="tex2jax_ignore mathjax_ignore section" id="linear-algebra-handling-of-arrays-and-more-python-features">
<h1><span class="section-number">2. </span>Linear Algebra, Handling of Arrays and more Python Features<a class="headerlink" href="#linear-algebra-handling-of-arrays-and-more-python-features" title="Permalink to this headline"></a></h1>
<div class="section" id="introduction">
<h2><span class="section-number">2.1. </span>Introduction<a class="headerlink" href="#introduction" title="Permalink to this headline"></a></h2>
<p>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.</p>
<p>As discussed in the introductory notes, these series of lectures focuses both on using
central Python packages like <strong>tensorflow</strong> and <strong>scikit-learn</strong> 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.</p>
</div>
<div class="section" id="important-matrix-and-vector-handling-packages">
<h2><span class="section-number">2.2. </span>Important Matrix and vector handling packages<a class="headerlink" href="#important-matrix-and-vector-handling-packages" title="Permalink to this headline"></a></h2>
<p>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 <strong>Numerical Recipes</strong>. 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. We describe them shortly here.</p>
<ul class="simple">
<li><p>LINPACK: package for linear equations and least square problems.</p></li>
<li><p>LAPACK:package for solving symmetric, unsymmetric and generalized eigenvalue problems. From LAPACK’s website <a class="reference external" href="http://www.netlib.org">http://www.netlib.org</a> it is possible to download for free all source codes from this library. Both C/C++ and Fortran versions are available.</p></li>
<li><p>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 class="reference external" href="http://www.netlib.org">http://www.netlib.org</a>.</p></li>
</ul>
<p>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,
<strong>Numpy</strong> 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, <strong>Armadillo</strong> 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>
</div>
<div class="section" id="basic-matrix-features">
<h2><span class="section-number">2.3. </span>Basic Matrix Features<a class="headerlink" href="#basic-matrix-features" title="Permalink to this headline"></a></h2>
<p>Matrix properties reminder</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\mathbf{A} =
\begin{bmatrix} a_{11} &amp; a_{12} &amp; a_{13} &amp; a_{14} \\
a_{21} &amp; a_{22} &amp; a_{23} &amp; a_{24} \\
a_{31} &amp; a_{32} &amp; a_{33} &amp; a_{34} \\
a_{41} &amp; a_{42} &amp; a_{43} &amp; a_{44}
\end{bmatrix}\qquad
\mathbf{I} =
\begin{bmatrix} 1 &amp; 0 &amp; 0 &amp; 0 \\
0 &amp; 1 &amp; 0 &amp; 0 \\
0 &amp; 0 &amp; 1 &amp; 0 \\
0 &amp; 0 &amp; 0 &amp; 1
\end{bmatrix}
\end{split}\]</div>
<p>The inverse of a matrix is defined by</p>
<div class="math notranslate nohighlight">
\[
\mathbf{A}^{-1} \cdot \mathbf{A} = I
\]</div>
<table class="dotable" 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 class="section" id="some-famous-matrices">
<h3><span class="section-number">2.3.1. </span>Some famous Matrices<a class="headerlink" href="#some-famous-matrices" title="Permalink to this headline"></a></h3>
<ul class="simple">
<li><p>Diagonal if <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i\ne j\)</span></p></li>
<li><p>Upper triangular if <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i &gt; j\)</span></p></li>
<li><p>Lower triangular if <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i &lt; j\)</span></p></li>
<li><p>Upper Hessenberg if <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i &gt; j+1\)</span></p></li>
<li><p>Lower Hessenberg if <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i &lt; j+1\)</span></p></li>
<li><p>Tridiagonal if <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(|i -j| &gt; 1\)</span></p></li>
<li><p>Lower banded with bandwidth <span class="math notranslate nohighlight">\(p\)</span>: <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i &gt; j+p\)</span></p></li>
<li><p>Upper banded with bandwidth <span class="math notranslate nohighlight">\(p\)</span>: <span class="math notranslate nohighlight">\(a_{ij}=0\)</span> for <span class="math notranslate nohighlight">\(i &lt; j+p\)</span></p></li>
<li><p>Banded, block upper triangular, block lower triangular….</p></li>
</ul>
<p>Some Equivalent Statements. For an <span class="math notranslate nohighlight">\(N\times N\)</span> matrix <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> the following properties are all equivalent</p>
<ul class="simple">
<li><p>If the inverse of <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> exists, <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> is nonsingular.</p></li>
<li><p>The equation <span class="math notranslate nohighlight">\(\mathbf{Ax}=0\)</span> implies <span class="math notranslate nohighlight">\(\mathbf{x}=0\)</span>.</p></li>
<li><p>The rows of <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> form a basis of <span class="math notranslate nohighlight">\(R^N\)</span>.</p></li>
<li><p>The columns of <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> form a basis of <span class="math notranslate nohighlight">\(R^N\)</span>.</p></li>
<li><p><span class="math notranslate nohighlight">\(\mathbf{A}\)</span> is a product of elementary matrices.</p></li>
<li><p><span class="math notranslate nohighlight">\(0\)</span> is not eigenvalue of <span class="math notranslate nohighlight">\(\mathbf{A}\)</span>.</p></li>
</ul>
</div>
</div>
<div class="section" id="numpy-and-arrays">
<h2><span class="section-number">2.4. </span>Numpy and arrays<a class="headerlink" href="#numpy-and-arrays" title="Permalink to this headline"></a></h2>
<p><a class="reference external" href="http://www.numpy.org/">Numpy</a> provides an easy way to handle arrays in Python. The standard way to import this library is as</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">10</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[ 0.31579721 1.74724767 -0.89609007 2.48464841 -0.36051635 -2.41246325
0.28008933 0.45134965 -0.52204004 -0.48145226]
</pre></div>
</div>
</div>
</div>
<p>Here we have defined a vector <span class="math notranslate nohighlight">\(x\)</span> with <span class="math notranslate nohighlight">\(n=10\)</span> elements with its values given by the Normal distribution <span class="math notranslate nohighlight">\(N(0,1)\)</span>.
Another alternative is to declare a vector as follows</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1 2 3]
</pre></div>
</div>
</div>
</div>
<p>Here we have defined a vector with three elements, with <span class="math notranslate nohighlight">\(x_0=1\)</span>, <span class="math notranslate nohighlight">\(x_1=2\)</span> and <span class="math notranslate nohighlight">\(x_2=3\)</span>. Note that both Python and C++
start numbering array elements from <span class="math notranslate nohighlight">\(0\)</span> and on. This means that a vector with <span class="math notranslate nohighlight">\(n\)</span> elements has a sequence of entities <span class="math notranslate nohighlight">\(x_0, x_1, x_2, \dots, x_{n-1}\)</span>. We could also let (recommended) Numpy to compute the logarithms of a specific array as</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">4</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">]))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1.38629436 1.94591015 2.07944154]
</pre></div>
</div>
</div>
</div>
<p>Here we have used Numpy’s unary function <span class="math notranslate nohighlight">\(np.log\)</span>. 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 <strong>log</strong> function
from Python’s <strong>math</strong> module. The looping is done explicitely by the
<strong>np.log</strong> function. The alternative, and slower way to compute the
logarithms of a vector would be to write</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">math</span> <span class="kn">import</span> <span class="n">log</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">4</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">])</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="nb">len</span><span class="p">(</span><span class="n">x</span><span class="p">)):</span>
<span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">log</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1 1 2]
</pre></div>
</div>
</div>
</div>
<p>We note that our code is much longer already and we need to import the <strong>log</strong> function from the <strong>math</strong> module.
The attentive reader will also notice that the output is <span class="math notranslate nohighlight">\([1, 1, 2]\)</span>. Python interprets automacally our numbers as integers (like the <strong>automatic</strong> keyword in C++). To change this we could define our array elements to be double precision numbers as</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mi">4</span><span class="p">,</span> <span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">],</span> <span class="n">dtype</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1.38629436 1.94591015 2.07944154]
</pre></div>
</div>
</div>
</div>
<p>or simply write them as double precision numbers (Python uses 64 bits as default for floating point type variables), that is</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">,</span> <span class="mf">8.0</span><span class="p">]))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1.38629436 1.94591015 2.07944154]
</pre></div>
</div>
</div>
</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 <strong>itemsize</strong> functionality (the array <span class="math notranslate nohighlight">\(x\)</span> is actually an object which inherits the functionalities defined in Numpy) as</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">,</span> <span class="mf">8.0</span><span class="p">]))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">x</span><span class="o">.</span><span class="n">itemsize</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>8
</pre></div>
</div>
</div>
</div>
<p>Having defined vectors, we are now ready to try out matrices. We can define a <span class="math notranslate nohighlight">\(3 \times 3 \)</span> real matrix <span class="math notranslate nohighlight">\(\hat{A}\)</span>
as (recall that we user lowercase letters for vectors and uppercase letters for matrices)</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span> <span class="p">[</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">,</span> <span class="mf">8.0</span><span class="p">],</span> <span class="p">[</span><span class="mf">3.0</span><span class="p">,</span> <span class="mf">10.0</span><span class="p">,</span> <span class="mf">11.0</span><span class="p">],</span> <span class="p">[</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">5.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">]</span> <span class="p">]))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[1.38629436 1.94591015 2.07944154]
[1.09861229 2.30258509 2.39789527]
[1.38629436 1.60943791 1.94591015]]
</pre></div>
</div>
</div>
</div>
<p>If we use the <strong>shape</strong> function we would get <span class="math notranslate nohighlight">\((3, 3)\)</span> as output, that is verifying that our matrix is a <span class="math notranslate nohighlight">\(3\times 3\)</span> 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>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span> <span class="p">[</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">,</span> <span class="mf">8.0</span><span class="p">],</span> <span class="p">[</span><span class="mf">3.0</span><span class="p">,</span> <span class="mf">10.0</span><span class="p">,</span> <span class="mf">11.0</span><span class="p">],</span> <span class="p">[</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">5.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">]</span> <span class="p">]))</span>
<span class="c1"># print the first column, row-major order and elements start with 0</span>
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">[:,</span><span class="mi">0</span><span class="p">])</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1.38629436 1.09861229 1.38629436]
</pre></div>
</div>
</div>
</div>
<p>We can continue this was by printing out other columns or rows. The example here prints out the second column</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span> <span class="p">[</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">,</span> <span class="mf">8.0</span><span class="p">],</span> <span class="p">[</span><span class="mf">3.0</span><span class="p">,</span> <span class="mf">10.0</span><span class="p">,</span> <span class="mf">11.0</span><span class="p">],</span> <span class="p">[</span><span class="mf">4.0</span><span class="p">,</span> <span class="mf">5.0</span><span class="p">,</span> <span class="mf">7.0</span><span class="p">]</span> <span class="p">]))</span>
<span class="c1"># print the first column, row-major order and elements start with 0</span>
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">[</span><span class="mi">1</span><span class="p">,:])</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[1.09861229 2.30258509 2.39789527]
</pre></div>
</div>
</div>
</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 class="reference external" href="http://www.numpy.org/">Numpy website for more details</a>. Useful functions when defining a matrix are the <strong>np.zeros</strong> function which declares a matrix of a given dimension and sets all elements to zero</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">10</span>
<span class="c1"># define a matrix of dimension 10 x 10 and set all elements to zero</span>
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">(</span> <span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span> <span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
</pre></div>
</div>
</div>
</div>
<p>or initializing all elements to</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">10</span>
<span class="c1"># define a matrix of dimension 10 x 10 and set all elements to one</span>
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">ones</span><span class="p">(</span> <span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span> <span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
</pre></div>
</div>
</div>
</div>
<p>or as unitarily distributed random numbers (see the material on random number generators in the statistics part)</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">10</span>
<span class="c1"># define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]</span>
<span class="n">A</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">rand</span><span class="p">(</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">A</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[0.69873514 0.39457095 0.31927572 0.01704432 0.11837671 0.15047127
0.56465688 0.50274255 0.69818111 0.82296251]
[0.45308692 0.37112277 0.77794957 0.43043913 0.95327702 0.89793609
0.29401213 0.08245909 0.1520039 0.59511582]
[0.98004227 0.101058 0.57572321 0.61394448 0.963198 0.49616116
0.83488328 0.05432856 0.12814914 0.03856554]
[0.283078 0.19166136 0.29275129 0.39300201 0.59004971 0.29199381
0.40644745 0.9036573 0.44729805 0.34447052]
[0.57670824 0.6568551 0.84380376 0.86221134 0.28908491 0.25663096
0.62862896 0.1937079 0.15673992 0.44921888]
[0.25139357 0.11197884 0.26514544 0.11896755 0.13404683 0.21059098
0.86012593 0.67890723 0.97948913 0.30567713]
[0.51845286 0.76010633 0.71640333 0.75282841 0.47447472 0.78882958
0.84159521 0.3729492 0.80152684 0.04084872]
[0.38088413 0.63567272 0.82909728 0.13829298 0.26037366 0.92772833
0.74577867 0.42239354 0.20594513 0.72328506]
[0.38831624 0.44520102 0.17305512 0.0106014 0.94866246 0.84929103
0.81753152 0.06664867 0.96750421 0.50462474]
[0.32584888 0.58948347 0.22729927 0.01919702 0.74280244 0.73448544
0.33995567 0.61197218 0.5320148 0.34517495]]
</pre></div>
</div>
</div>
</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
<span class="math notranslate nohighlight">\(\hat{x}, \hat{y}, \hat{z}\)</span> with <span class="math notranslate nohighlight">\(n\)</span> elements each. The covariance matrix is defined as</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\hat{\Sigma} = \begin{bmatrix} \sigma_{xx} &amp; \sigma_{xy} &amp; \sigma_{xz} \\
\sigma_{yx} &amp; \sigma_{yy} &amp; \sigma_{yz} \\
\sigma_{zx} &amp; \sigma_{zy} &amp; \sigma_{zz}
\end{bmatrix},
\end{split}\]</div>
<p>where for example</p>
<div class="math notranslate nohighlight">
\[
\sigma_{xy} =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}).
\]</div>
<p>The Numpy function <strong>np.cov</strong> calculates the covariance elements using the factor <span class="math notranslate nohighlight">\(1/(n-1)\)</span> instead of <span class="math notranslate nohighlight">\(1/n\)</span> 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 <strong>np.vstack</strong> function which takes each vector of dimension <span class="math notranslate nohighlight">\(1\times n\)</span> and produces a <span class="math notranslate nohighlight">\( 3\times n\)</span> matrix <span class="math notranslate nohighlight">\(\hat{W}\)</span></p>
<div class="math notranslate nohighlight">
\[\begin{split}
\hat{W} = \begin{bmatrix} x_0 &amp; y_0 &amp; z_0 \\
x_1 &amp; y_1 &amp; z_1 \\
x_2 &amp; y_2 &amp; z_2 \\
\dots &amp; \dots &amp; \dots \\
x_{n-2} &amp; y_{n-2} &amp; z_{n-2} \\
x_{n-1} &amp; y_{n-1} &amp; z_{n-1}
\end{bmatrix},
\end{split}\]</div>
<p>which in turn is converted into into the <span class="math notranslate nohighlight">\(3 times 3\)</span> covariance matrix
<span class="math notranslate nohighlight">\(\hat{\Sigma}\)</span> via the Numpy function <strong>np.cov()</strong>. 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 <span class="math notranslate nohighlight">\(\hat{x}\)</span> etc using the Numpy
function <strong>np.mean(x)</strong>. We can also extract the eigenvalues of the
covariance matrix through the <strong>np.linalg.eig()</strong> function.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Importing various packages</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="n">n</span> <span class="o">=</span> <span class="mi">100</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
<span class="n">y</span> <span class="o">=</span> <span class="mi">4</span><span class="o">+</span><span class="mi">3</span><span class="o">*</span><span class="n">x</span><span class="o">+</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">y</span><span class="p">))</span>
<span class="n">z</span> <span class="o">=</span> <span class="n">x</span><span class="o">**</span><span class="mi">3</span><span class="o">+</span><span class="n">np</span><span class="o">.</span><span class="n">random</span><span class="o">.</span><span class="n">normal</span><span class="p">(</span><span class="n">size</span><span class="o">=</span><span class="n">n</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">z</span><span class="p">))</span>
<span class="n">W</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">vstack</span><span class="p">((</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">z</span><span class="p">))</span>
<span class="n">Sigma</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">cov</span><span class="p">(</span><span class="n">W</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">Sigma</span><span class="p">)</span>
<span class="n">Eigvals</span><span class="p">,</span> <span class="n">Eigvecs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linalg</span><span class="o">.</span><span class="n">eig</span><span class="p">(</span><span class="n">Sigma</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">Eigvals</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.1388976715362099
4.3703468543933255
0.08844723450419088
[[ 1.06638817 3.39483726 3.28607817]
[ 3.39483726 11.55955126 10.56366546]
[ 3.28607817 10.56366546 16.19343949]]
[25.58818643 0.05982961 3.17136288]
</pre></div>
</div>
</div>
</div>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="o">%</span><span class="k">matplotlib</span> inline
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">from</span> <span class="nn">scipy</span> <span class="kn">import</span> <span class="n">sparse</span>
<span class="n">eye</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">eye</span><span class="p">(</span><span class="mi">4</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">eye</span><span class="p">)</span>
<span class="n">sparse_mtx</span> <span class="o">=</span> <span class="n">sparse</span><span class="o">.</span><span class="n">csr_matrix</span><span class="p">(</span><span class="n">eye</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="n">sparse_mtx</span><span class="p">)</span>
<span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="o">-</span><span class="mi">10</span><span class="p">,</span><span class="mi">10</span><span class="p">,</span><span class="mi">100</span><span class="p">)</span>
<span class="n">y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sin</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">marker</span><span class="o">=</span><span class="s1">&#39;x&#39;</span><span class="p">)</span>
<span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[1. 0. 0. 0.]
[0. 1. 0. 0.]
[0. 0. 1. 0.]
[0. 0. 0. 1.]]
(0, 0) 1.0
(1, 1) 1.0
(2, 2) 1.0
(3, 3) 1.0
</pre></div>
</div>
<img alt="_images/linalg_44_1.png" src="_images/linalg_44_1.png" />
</div>
</div>
</div>
<div class="section" id="other-matrix-and-vector-operations">
<h2><span class="section-number">2.5. </span>Other Matrix and Vector Operations<a class="headerlink" href="#other-matrix-and-vector-operations" title="Permalink to this headline"></a></h2>
<p>The following examples show how to compute various quantities like the <strong>mean</strong> value of a matrix or a vector and how to use functions like <strong>reshape</strong> and <strong>ravel</strong>. These are all useful quantities when scaling the data and preparing the data for various machine learning algorithms and when calculating quantities like the mean squared error or the variance.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Simple code that tests various numpy functions</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="c1"># Simple test-matrix of dim 3 x 4</span>
<span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span> <span class="p">[</span><span class="mi">1</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">],</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span> <span class="mi">5</span><span class="p">,</span> <span class="mi">6</span><span class="p">],</span> <span class="p">[</span><span class="mi">7</span><span class="p">,</span> <span class="mi">8</span><span class="p">,</span> <span class="mi">9</span><span class="p">],[</span><span class="mi">10</span><span class="p">,</span> <span class="mi">11</span><span class="p">,</span> <span class="mi">12</span><span class="p">]],</span><span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">)</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;The test matrix:</span><span class="si">{</span><span class="n">a</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="c1"># This is the total mean summed over all elements, which here has to be 6.5</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;This is the total mean summed over all elements:</span><span class="si">{</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">)</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="c1"># This is the mean for each column, it returns an array with the mean values for each column. It returns a row-like vector</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;This is the mean for each column:</span><span class="si">{</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="w"> </span><span class="n">axis</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span><span class="w"> </span><span class="n">keepdims</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span><span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">)</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="c1"># This is the mean value for each row, it returns an array via the keepdims option which is a column-like vector if</span>
<span class="c1"># keepdims=True. Else it return a row-like vector</span>
<span class="c1"># Try setting keepdims=False</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;This is the mean value for each row:</span><span class="si">{</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="w"> </span><span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">keepdims</span><span class="o">=</span><span class="kc">True</span><span class="p">,</span><span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">)</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="c1"># We print then the mean value for each row by setting keepdims=False</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;This is the mean value for each row with keepdims false:</span><span class="si">{</span><span class="n">np</span><span class="o">.</span><span class="n">mean</span><span class="p">(</span><span class="n">a</span><span class="p">,</span><span class="w"> </span><span class="n">axis</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span><span class="w"> </span><span class="n">keepdims</span><span class="o">=</span><span class="kc">False</span><span class="p">,</span><span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">)</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>The test matrix:[[ 1. 2. 3.]
[ 4. 5. 6.]
[ 7. 8. 9.]
[10. 11. 12.]]
This is the total mean summed over all elements:6.5
This is the mean for each column:[[5.5 6.5 7.5]]
This is the mean value for each row:[[ 2.]
[ 5.]
[ 8.]
[11.]]
This is the mean value for each row with keepdims false:[ 2. 5. 8. 11.]
</pre></div>
</div>
</div>
</div>
<p>Another useful function is the <strong>ravel</strong> function, which returns a flattened array as shown in the example here.</p>
<div class="cell docutils container">
<div class="cell_input docutils container">
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="c1"># Ravel return a contiguous flattened array.</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;Flatten the matrix:</span><span class="si">{</span><span class="n">np</span><span class="o">.</span><span class="n">ravel</span><span class="p">(</span><span class="n">a</span><span class="p">)</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="c1"># It is the same as reshaping the matrix into a one-dimensional array</span>
<span class="nb">print</span><span class="p">(</span><span class="sa">f</span><span class="s2">&quot;Reshape the matrix to a one-dim array:</span><span class="si">{</span><span class="n">a</span><span class="o">.</span><span class="n">reshape</span><span class="p">(</span><span class="o">-</span><span class="mi">1</span><span class="p">)</span><span class="si">}</span><span class="s2">&quot;</span><span class="p">)</span>
<span class="c1"># ‘C’ means to index the elements in row-major, C-style order, with the last axis index changing fastest, back to the first axis index changing slowest.</span>
<span class="c1"># ‘F’ means to index the elements in column-major, Fortran-style order, with the first index changing fastest, and the last index changing slowest </span>
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">ravel</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="s1">&#39;F&#39;</span><span class="p">))</span>
<span class="c1"># When order is ‘A’, it will preserve the array’s ‘C’ or ‘F’ ordering</span>
<span class="c1"># ‘A’ means to read the elements in Fortran-like index order if a is Fortran contiguous in memory, C-like order otherwise.</span>
<span class="c1"># ‘K’ means to read the elements in the order they occur in memory, except for reversing the data when strides are negative. By default, ‘C’ index order is used.</span>
<span class="c1"># Transposing it</span>
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">ravel</span><span class="p">(</span><span class="n">a</span><span class="o">.</span><span class="n">T</span><span class="p">))</span>
<span class="nb">print</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">ravel</span><span class="p">(</span><span class="n">a</span><span class="o">.</span><span class="n">T</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="s1">&#39;A&#39;</span><span class="p">))</span>
</pre></div>
</div>
</div>
<div class="cell_output docutils container">
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>Flatten the matrix:[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.]
Reshape the matrix to a one-dim array:[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.]
[ 1. 4. 7. 10. 2. 5. 8. 11. 3. 6. 9. 12.]
[ 1. 4. 7. 10. 2. 5. 8. 11. 3. 6. 9. 12.]
[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12.]
</pre></div>
</div>
</div>
</div>
</div>
<div class="section" id="gaussian-elimination">
<h2><span class="section-number">2.6. </span>Gaussian Elimination<a class="headerlink" href="#gaussian-elimination" title="Permalink to this headline"></a></h2>
<p>We start with the linear set of equations</p>
<div class="math notranslate nohighlight">
\[
\mathbf{A}\mathbf{x} = \mathbf{w}.
\]</div>
<p>We assume also that the matrix <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> is non-singular and that the
matrix elements along the diagonal satisfy <span class="math notranslate nohighlight">\(a_{ii} \ne 0\)</span>. Simple <span class="math notranslate nohighlight">\(4\times 4 \)</span> example</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{bmatrix}
a_{11}&amp; a_{12} &amp;a_{13}&amp; a_{14}\\
a_{21}&amp; a_{22} &amp;a_{23}&amp; a_{24}\\
a_{31}&amp; a_{32} &amp;a_{33}&amp; a_{34}\\
a_{41}&amp; a_{42} &amp;a_{43}&amp; 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}.
\end{split}\]</div>
<p>or</p>
<div class="math notranslate nohighlight">
\[
a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \nonumber
\]</div>
<p>The basic idea of Gaussian elimination is to use the first equation to eliminate the first unknown <span class="math notranslate nohighlight">\(x_1\)</span>
from the remaining <span class="math notranslate nohighlight">\(n-1\)</span> equations. Then we use the new second equation to eliminate the second unknown
<span class="math notranslate nohighlight">\(x_2\)</span> from the remaining <span class="math notranslate nohighlight">\(n-2\)</span> equations. With <span class="math notranslate nohighlight">\(n-1\)</span> such eliminations
we obtain a so-called upper triangular set of equations of the form</p>
<div class="math notranslate nohighlight">
\[
b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=y_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
b_{22}x_2 + b_{23}x_3 + b_{24}x_4=y_2 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
b_{33}x_3 + b_{34}x_4=y_3 \nonumber
\]</div>
<!-- Equation labels as ordinary links -->
<div id="eq:gaussbacksub"></div>
<div class="math notranslate nohighlight">
\[
b_{44}x_4=y_4. \nonumber
\label{eq:gaussbacksub} \tag{1}
\]</div>
<p>We can solve this system of equations recursively starting from <span class="math notranslate nohighlight">\(x_n\)</span> (in our case <span class="math notranslate nohighlight">\(x_4\)</span>) and proceed with
what is called a backward substitution.</p>
<p>This process can be expressed mathematically as</p>
<!-- Equation labels as ordinary links -->
<div id="_auto1"></div>
<div class="math notranslate nohighlight">
\[
\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} \tag{2}
\end{equation}
\]</div>
<p>To arrive at such an upper triangular system of equations, we start by eliminating
the unknown <span class="math notranslate nohighlight">\(x_1\)</span> for <span class="math notranslate nohighlight">\(j=2,n\)</span>. We achieve this by multiplying the first equation by <span class="math notranslate nohighlight">\(a_{j1}/a_{11}\)</span> and then subtract
the result from the <span class="math notranslate nohighlight">\(j\)</span>th equation. We assume obviously that <span class="math notranslate nohighlight">\(a_{11}\ne 0\)</span> and that
<span class="math notranslate nohighlight">\(\mathbf{A}\)</span> is not singular.</p>
<p>Our actual <span class="math notranslate nohighlight">\(4\times 4\)</span> example reads after the first operation</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{bmatrix}
a_{11}&amp; a_{12} &amp;a_{13}&amp; a_{14}\\
0&amp; (a_{22}-\frac{a_{21}a_{12}}{a_{11}}) &amp;(a_{23}-\frac{a_{21}a_{13}}{a_{11}}) &amp; (a_{24}-\frac{a_{21}a_{14}}{a_{11}})\\
0&amp; (a_{32}-\frac{a_{31}a_{12}}{a_{11}})&amp; (a_{33}-\frac{a_{31}a_{13}}{a_{11}})&amp; (a_{34}-\frac{a_{31}a_{14}}{a_{11}})\\
0&amp;(a_{42}-\frac{a_{41}a_{12}}{a_{11}}) &amp;(a_{43}-\frac{a_{41}a_{13}}{a_{11}}) &amp; (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},
\end{split}\]</div>
<p>or</p>
<div class="math notranslate nohighlight">
\[
b_{11}x_1 +b_{12}x_2 +b_{13}x_3 + b_{14}x_4=y_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a^{(2)}_{22}x_2 + a^{(2)}_{23}x_3 + a^{(2)}_{24}x_4=w^{(2)}_2 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a^{(2)}_{32}x_2 + a^{(2)}_{33}x_3 + a^{(2)}_{34}x_4=w^{(2)}_3 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a^{(2)}_{42}x_2 + a^{(2)}_{43}x_3 + a^{(2)}_{44}x_4=w^{(2)}_4, \nonumber
\]</div>
<!-- Equation labels as ordinary links -->
<div id="_auto2"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
\label{_auto2} \tag{3}
\end{equation}
\]</div>
<p>The new coefficients are</p>
<!-- Equation labels as ordinary links -->
<div id="_auto3"></div>
<div class="math notranslate nohighlight">
\[
\begin{equation}
b_{1k} = a_{1k}^{(1)} \quad k=1,\dots,n,
\label{_auto3} \tag{4}
\end{equation}
\]</div>
<p>where each <span class="math notranslate nohighlight">\(a_{1k}^{(1)}\)</span> is equal to the original <span class="math notranslate nohighlight">\(a_{1k}\)</span> element. The other coefficients are</p>
<!-- Equation labels as ordinary links -->
<div id="_auto4"></div>
<div class="math notranslate nohighlight">
\[
\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} \tag{5}
\end{equation}
\]</div>
<p>with a new right-hand side given by</p>
<!-- Equation labels as ordinary links -->
<div id="_auto5"></div>
<div class="math notranslate nohighlight">
\[
\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} \tag{6}
\end{equation}
\]</div>
<p>We have also set <span class="math notranslate nohighlight">\(w_1^{(1)}=w_1\)</span>, the original vector element.
We see that the system of unknowns <span class="math notranslate nohighlight">\(x_1,\dots,x_n\)</span> is transformed into an <span class="math notranslate nohighlight">\((n-1)\times (n-1)\)</span> problem.</p>
<p>This step is called forward substitution.
Proceeding with these substitutions, we obtain the
general expressions for the new coefficients</p>
<!-- Equation labels as ordinary links -->
<div id="_auto6"></div>
<div class="math notranslate nohighlight">
\[
\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} \tag{7}
\end{equation}
\]</div>
<p>with <span class="math notranslate nohighlight">\(m=1,\dots,n-1\)</span> and a
right-hand side given by</p>
<!-- Equation labels as ordinary links -->
<div id="_auto7"></div>
<div class="math notranslate nohighlight">
\[
\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} \tag{8}
\end{equation}
\]</div>
<p>This set of <span class="math notranslate nohighlight">\(n-1\)</span> elimations leads us to an equations which is solved by back substitution.
If the arithmetics is exact and the matrix <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> is not singular, then the computed answer will be exact.</p>
<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 <span class="math notranslate nohighlight">\((a_{22}-a_{21}a_{12}/a_{11})\)</span>
results in <span class="math notranslate nohighlight">\(-10^{-7}\)</span> and that <span class="math notranslate nohighlight">\(a_{22}\)</span> is one.
one. We are then
adding <span class="math notranslate nohighlight">\(10^7+1\)</span>. With single precision this results in <span class="math notranslate nohighlight">\(10^7\)</span>.</p>
<ul class="simple">
<li><p>Gaussian elimination, <span class="math notranslate nohighlight">\(O(2/3n^3)\)</span> flops, general matrix</p></li>
<li><p>LU decomposition, upper triangular and lower tridiagonal matrices, <span class="math notranslate nohighlight">\(O(2/3n^3)\)</span> flops, general matrix. Get easily the inverse, determinant and can solve linear equations with back-substitution only, <span class="math notranslate nohighlight">\(O(n^2)\)</span> flops</p></li>
<li><p>Cholesky decomposition. Real symmetric or hermitian positive definite matrix, <span class="math notranslate nohighlight">\(O(1/3n^3)\)</span> flops.</p></li>
<li><p>Tridiagonal linear systems, important for differential equations. Normally positive definite and non-singular. <span class="math notranslate nohighlight">\(O(8n)\)</span> flops for symmetric. Special case of banded matrices.</p></li>
<li><p>Singular value decomposition</p></li>
<li><p>the QR method will be discussed in chapter 7 in connection with eigenvalue systems. <span class="math notranslate nohighlight">\(O(4/3n^3)\)</span> flops.</p></li>
</ul>
<p>The LU decomposition method means that we can rewrite
this matrix as the product of two matrices <span class="math notranslate nohighlight">\(\mathbf{L}\)</span> and <span class="math notranslate nohighlight">\(\mathbf{U}\)</span>
where</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\begin{bmatrix}
a_{11} &amp; a_{12} &amp; a_{13} &amp; a_{14} \\
a_{21} &amp; a_{22} &amp; a_{23} &amp; a_{24} \\
a_{31} &amp; a_{32} &amp; a_{33} &amp; a_{34} \\
a_{41} &amp; a_{42} &amp; a_{43} &amp; a_{44}
\end{bmatrix}
= \begin{bmatrix}
1 &amp; 0 &amp; 0 &amp; 0 \\
l_{21} &amp; 1 &amp; 0 &amp; 0 \\
l_{31} &amp; l_{32} &amp; 1 &amp; 0 \\
l_{41} &amp; l_{42} &amp; l_{43} &amp; 1
\end{bmatrix}
\begin{bmatrix}
u_{11} &amp; u_{12} &amp; u_{13} &amp; u_{14} \\
0 &amp; u_{22} &amp; u_{23} &amp; u_{24} \\
0 &amp; 0 &amp; u_{33} &amp; u_{34} \\
0 &amp; 0 &amp; 0 &amp; u_{44}
\end{bmatrix}.
\end{split}\]</div>
<p>LU decomposition forms the backbone of other algorithms in linear algebra, such as the
solution of linear equations given by</p>
<div class="math notranslate nohighlight">
\[
a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \nonumber
\]</div>
<p>The above set of equations is conveniently solved by using LU decomposition as an intermediate step.</p>
<p>The matrix <span class="math notranslate nohighlight">\(\mathbf{A}\in \mathbb{R}^{n\times n}\)</span> has an LU factorization if the determinant
is different from zero. If the LU factorization exists and <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> is non-singular, then the LU factorization
is unique and the determinant is given by</p>
<div class="math notranslate nohighlight">
\[
det\{\mathbf{A}\}=det\{\mathbf{LU}\}= det\{\mathbf{L}\}det\{\mathbf{U}\}=u_{11}u_{22}\dots u_{nn}.
\]</div>
<p>There are at least three main advantages with LU decomposition compared with standard Gaussian elimination:</p>
<ul class="simple">
<li><p>It is straightforward to compute the determinant of a matrix</p></li>
<li><p>If we have to solve sets of linear equations with the same matrix but with different vectors <span class="math notranslate nohighlight">\(\mathbf{y}\)</span>, the number of FLOPS is of the order <span class="math notranslate nohighlight">\(n^3\)</span>.</p></li>
<li><p>The inverse is such an operation</p></li>
</ul>
<p>With the LU decomposition it is rather
simple to solve a system of linear equations</p>
<div class="math notranslate nohighlight">
\[
a_{11}x_1 +a_{12}x_2 +a_{13}x_3 + a_{14}x_4=w_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + a_{24}x_4=w_2 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 + a_{34}x_4=w_3 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
a_{41}x_1 + a_{42}x_2 + a_{43}x_3 + a_{44}x_4=w_4. \nonumber
\]</div>
<p>This can be written in matrix form as</p>
<div class="math notranslate nohighlight">
\[
\mathbf{Ax}=\mathbf{w}.
\]</div>
<p>where <span class="math notranslate nohighlight">\(\mathbf{A}\)</span> and <span class="math notranslate nohighlight">\(\mathbf{w}\)</span> are known and we have to solve for
<span class="math notranslate nohighlight">\(\mathbf{x}\)</span>. Using the LU dcomposition we write</p>
<div class="math notranslate nohighlight">
\[
\mathbf{A} \mathbf{x} \equiv \mathbf{L} \mathbf{U} \mathbf{x} =\mathbf{w}.
\]</div>
<p>The previous equation can be calculated in two steps</p>
<div class="math notranslate nohighlight">
\[
\mathbf{L} \mathbf{y} = \mathbf{w};\qquad \mathbf{Ux}=\mathbf{y}.
\]</div>
<p>To show that this is correct we use to the LU decomposition
to rewrite our system of linear equations as</p>
<div class="math notranslate nohighlight">
\[
\mathbf{LUx}=\mathbf{w},
\]</div>
<p>and since the determinant of <span class="math notranslate nohighlight">\(\mathbf{L}\)</span> is equal to 1 (by construction
since the diagonals of <span class="math notranslate nohighlight">\(\mathbf{L}\)</span> equal 1) we can use the inverse of
<span class="math notranslate nohighlight">\(\mathbf{L}\)</span> to obtain</p>
<div class="math notranslate nohighlight">
\[
\mathbf{Ux}=\mathbf{L^{-1}w}=\mathbf{y},
\]</div>
<p>which yields the intermediate step</p>
<div class="math notranslate nohighlight">
\[
\mathbf{L^{-1}w}=\mathbf{y}
\]</div>
<p>and as soon as we have <span class="math notranslate nohighlight">\(\mathbf{y}\)</span> we can obtain <span class="math notranslate nohighlight">\(\mathbf{x}\)</span>
through <span class="math notranslate nohighlight">\(\mathbf{Ux}=\mathbf{y}\)</span>.</p>
<p>For our four-dimentional example this takes the form</p>
<div class="math notranslate nohighlight">
\[
y_1=w_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
l_{21}y_1 + y_2=w_2\nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
l_{31}y_1 + l_{32}y_2 + y_3 =w_3\nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
l_{41}y_1 + l_{42}y_2 + l_{43}y_3 + y_4=w_4. \nonumber
\]</div>
<p>and</p>
<div class="math notranslate nohighlight">
\[
u_{11}x_1 +u_{12}x_2 +u_{13}x_3 + u_{14}x_4=y_1 \nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
u_{22}x_2 + u_{23}x_3 + u_{24}x_4=y_2\nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
u_{33}x_3 + u_{34}x_4=y_3\nonumber
\]</div>
<div class="math notranslate nohighlight">
\[
u_{44}x_4=y_4 \nonumber
\]</div>
<p>This example shows the basis for the algorithm
needed to solve the set of <span class="math notranslate nohighlight">\(n\)</span> linear equations.</p>
<p>The algorithm goes as follows</p>
<ul class="simple">
<li><p>Set up the matrix <span class="math notranslate nohighlight">\(\bf A\)</span> and the vector <span class="math notranslate nohighlight">\(\bf w\)</span> with their correct dimensions. This determines the dimensionality of the unknown vector <span class="math notranslate nohighlight">\(\bf x\)</span>.</p></li>
<li><p>Then LU decompose the matrix <span class="math notranslate nohighlight">\(\bf A\)</span> through a call to the function <code class="docutils literal notranslate"><span class="pre">ludcmp(double</span> <span class="pre">a,</span> <span class="pre">int</span> <span class="pre">n,</span> <span class="pre">int</span> <span class="pre">indx,</span> <span class="pre">double</span> <span class="pre">&amp;d)</span></code>. This functions returns the LU decomposed matrix <span class="math notranslate nohighlight">\(\bf A\)</span>, 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.</p></li>
<li><p>Thereafter you call the function <code class="docutils literal notranslate"><span class="pre">lubksb(double</span> <span class="pre">a,</span> <span class="pre">int</span> <span class="pre">n,</span> <span class="pre">int</span> <span class="pre">indx,</span> <span class="pre">double</span> <span class="pre">w)</span></code> which uses the LU decomposed matrix <span class="math notranslate nohighlight">\(\bf A\)</span> and the vector <span class="math notranslate nohighlight">\(\bf w\)</span> and returns <span class="math notranslate nohighlight">\(\bf x\)</span> in the same place as <span class="math notranslate nohighlight">\(\bf w\)</span>. Upon exit the original content in <span class="math notranslate nohighlight">\(\bf w\)</span> is destroyed. If you wish to keep this information, you should make a backup of it in your calling function.</p></li>
</ul>
<div class="section" id="lu-decomposition-the-inverse-of-a-matrix">
<h3><span class="section-number">2.6.1. </span>LU Decomposition, the inverse of a matrix<a class="headerlink" href="#lu-decomposition-the-inverse-of-a-matrix" title="Permalink to this headline"></a></h3>
<p>If the inverse exists then</p>
<div class="math notranslate nohighlight">
\[
\mathbf{A}^{-1}\mathbf{A}=\mathbf{I},
\]</div>
<p>the identity matrix. With an LU decomposed matrix we can rewrite the last equation as</p>
<div class="math notranslate nohighlight">
\[
\mathbf{LU}\mathbf{A}^{-1}=\mathbf{I}.
\]</div>
<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</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\mathbf{A}_1^{-1}= \begin{bmatrix}
a_{11}^{-1} \\
a_{21}^{-1} \\
\dots \\
a_{n1}^{-1} \\
\end{bmatrix},
\end{split}\]</div>
<p>then we have a linear set of equations</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\mathbf{LU}\begin{bmatrix}
a_{11}^{-1} \\
a_{21}^{-1} \\
\dots \\
a_{n1}^{-1} \\
\end{bmatrix} =\begin{bmatrix}
1 \\
0 \\
\dots \\
0 \\
\end{bmatrix}.
\end{split}\]</div>
<p>In a similar way we can compute the unknow entries of the second column,</p>
<div class="math notranslate nohighlight">
\[\begin{split}
\mathbf{LU}\begin{bmatrix}
a_{12}^{-1} \\
a_{22}^{-1} \\
\dots \\
a_{n2}^{-1} \\
\end{bmatrix}=\begin{bmatrix}
0 \\
1 \\
\dots \\
0 \\
\end{bmatrix},
\end{split}\]</div>
<p>and continue till we have solved all <span class="math notranslate nohighlight">\(n\)</span> sets of linear equations.</p>
<p>The calculation of the inverse here assumes that it actually
exists. In many machine learning applications there may be strong
linear dependencies among the various columns and/or rows. In our
discussions of linear regression we will dive into the mathematics of
the singular value decomposition, an algorithm which will allow us to calculate the so-called pseudo-inverse.
These details will be presented in our linear regression chapter.</p>
</div>
</div>
</div>
<script type="text/x-thebe-config">
{
requestKernel: true,
binderOptions: {
repo: "binder-examples/jupyter-stacks-datascience",
ref: "master",
},
codeMirrorConfig: {
theme: "abcdef",
mode: "python"
},
kernelOptions: {
kernelName: "python3",
path: "./."
},
predefinedOutput: true
}
</script>
<script>kernelName = 'python3'</script>
</div>
<!-- Previous / next buttons -->
<div class='prev-next-area'>
<a class='left-prev' id="prev-link" href="statistics.html" title="previous page">
<i class="fas fa-angle-left"></i>
<div class="prev-next-info">
<p class="prev-next-subtitle">previous</p>
<p class="prev-next-title"><span class="section-number">1. </span>Elements of Probability Theory and Statistical Data Analysis</p>
</div>
</a>
<a class='right-next' id="next-link" href="chapter1.html" title="next page">
<div class="prev-next-info">
<p class="prev-next-subtitle">next</p>
<p class="prev-next-title"><span class="section-number">3. </span>Linear Regression</p>
</div>
<i class="fas fa-angle-right"></i>
</a>
</div>
</div>
</div>
<footer class="footer">
<p>
By Morten Hjorth-Jensen<br/>
&copy; Copyright 2021.<br/>
</p>
</footer>
</main>
</div>
</div>
<script src="_static/js/index.be7d3bbb2ef33a8344ce.js"></script>
</body>
</html>
Reference in New Issue View Git Blame Copy Permalink