Morten Hjorth-Jensen
1337 lines
67 KiB
HTML
1337 lines
67 KiB
HTML
|
||
<!DOCTYPE html>
|
||
|
||
|
||
<html lang="en" data-content_root="./" >
|
||
|
||
<head>
|
||
<meta charset="utf-8" />
|
||
<meta name="viewport" content="width=device-width, initial-scale=1.0" /><meta name="viewport" content="width=device-width, initial-scale=1" />
|
||
|
||
<title>2. Linear Algebra, Handling of Arrays and more Python Features — Applied Data Analysis and Machine Learning</title>
|
||
|
||
|
||
|
||
<script data-cfasync="false">
|
||
document.documentElement.dataset.mode = localStorage.getItem("mode") || "";
|
||
document.documentElement.dataset.theme = localStorage.getItem("theme") || "";
|
||
</script>
|
||
|
||
<!-- Loaded before other Sphinx assets -->
|
||
<link href="_static/styles/theme.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
|
||
<link href="_static/styles/bootstrap.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
|
||
<link href="_static/styles/pydata-sphinx-theme.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
|
||
|
||
|
||
<link href="_static/vendor/fontawesome/6.5.2/css/all.min.css?digest=dfe6caa3a7d634c4db9b" rel="stylesheet" />
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-solid-900.woff2" />
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-brands-400.woff2" />
|
||
<link rel="preload" as="font" type="font/woff2" crossorigin href="_static/vendor/fontawesome/6.5.2/webfonts/fa-regular-400.woff2" />
|
||
|
||
<link rel="stylesheet" type="text/css" href="_static/pygments.css?v=fa44fd50" />
|
||
<link rel="stylesheet" type="text/css" href="_static/styles/sphinx-book-theme.css?v=eba8b062" />
|
||
<link rel="stylesheet" type="text/css" href="_static/togglebutton.css?v=13237357" />
|
||
<link rel="stylesheet" type="text/css" href="_static/copybutton.css?v=76b2166b" />
|
||
<link rel="stylesheet" type="text/css" href="_static/mystnb.8ecb98da25f57f5357bf6f572d296f466b2cfe2517ffebfabe82451661e28f02.css?v=6644e6bb" />
|
||
<link rel="stylesheet" type="text/css" href="_static/sphinx-thebe.css?v=4fa983c6" />
|
||
<link rel="stylesheet" type="text/css" href="_static/sphinx-design.min.css?v=95c83b7e" />
|
||
|
||
<!-- Pre-loaded scripts that we'll load fully later -->
|
||
<link rel="preload" as="script" href="_static/scripts/bootstrap.js?digest=dfe6caa3a7d634c4db9b" />
|
||
<link rel="preload" as="script" href="_static/scripts/pydata-sphinx-theme.js?digest=dfe6caa3a7d634c4db9b" />
|
||
<script src="_static/vendor/fontawesome/6.5.2/js/all.min.js?digest=dfe6caa3a7d634c4db9b"></script>
|
||
|
||
<script src="_static/documentation_options.js?v=9eb32ce0"></script>
|
||
<script src="_static/doctools.js?v=9a2dae69"></script>
|
||
<script src="_static/sphinx_highlight.js?v=dc90522c"></script>
|
||
<script src="_static/clipboard.min.js?v=a7894cd8"></script>
|
||
<script src="_static/copybutton.js?v=f281be69"></script>
|
||
<script src="_static/scripts/sphinx-book-theme.js?v=887ef09a"></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?v=4a39c7ea"></script>
|
||
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown';</script>
|
||
<script src="_static/design-tabs.js?v=f930bc37"></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?v=c100c467"></script>
|
||
<script>var togglebuttonSelector = '.toggle, .admonition.dropdown';</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>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>
|
||
<script>DOCUMENTATION_OPTIONS.pagename = 'linalg';</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="en"/>
|
||
</head>
|
||
|
||
|
||
<body data-bs-spy="scroll" data-bs-target=".bd-toc-nav" data-offset="180" data-bs-root-margin="0px 0px -60%" data-default-mode="">
|
||
|
||
|
||
|
||
<div id="pst-skip-link" class="skip-link d-print-none"><a href="#main-content">Skip to main content</a></div>
|
||
|
||
<div id="pst-scroll-pixel-helper"></div>
|
||
|
||
<button type="button" class="btn rounded-pill" id="pst-back-to-top">
|
||
<i class="fa-solid fa-arrow-up"></i>Back to top</button>
|
||
|
||
|
||
<input type="checkbox"
|
||
class="sidebar-toggle"
|
||
id="pst-primary-sidebar-checkbox"/>
|
||
<label class="overlay overlay-primary" for="pst-primary-sidebar-checkbox"></label>
|
||
|
||
<input type="checkbox"
|
||
class="sidebar-toggle"
|
||
id="pst-secondary-sidebar-checkbox"/>
|
||
<label class="overlay overlay-secondary" for="pst-secondary-sidebar-checkbox"></label>
|
||
|
||
<div class="search-button__wrapper">
|
||
<div class="search-button__overlay"></div>
|
||
<div class="search-button__search-container">
|
||
<form class="bd-search d-flex align-items-center"
|
||
action="search.html"
|
||
method="get">
|
||
<i class="fa-solid fa-magnifying-glass"></i>
|
||
<input type="search"
|
||
class="form-control"
|
||
name="q"
|
||
id="search-input"
|
||
placeholder="Search this book..."
|
||
aria-label="Search this book..."
|
||
autocomplete="off"
|
||
autocorrect="off"
|
||
autocapitalize="off"
|
||
spellcheck="false"/>
|
||
<span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd>K</kbd></span>
|
||
</form></div>
|
||
</div>
|
||
|
||
<div class="pst-async-banner-revealer d-none">
|
||
<aside id="bd-header-version-warning" class="d-none d-print-none" aria-label="Version warning"></aside>
|
||
</div>
|
||
|
||
|
||
<header class="bd-header navbar navbar-expand-lg bd-navbar d-print-none">
|
||
</header>
|
||
|
||
|
||
<div class="bd-container">
|
||
<div class="bd-container__inner bd-page-width">
|
||
|
||
|
||
|
||
<div class="bd-sidebar-primary bd-sidebar">
|
||
|
||
|
||
|
||
<div class="sidebar-header-items sidebar-primary__section">
|
||
|
||
|
||
|
||
|
||
</div>
|
||
|
||
<div class="sidebar-primary-items__start sidebar-primary__section">
|
||
<div class="sidebar-primary-item">
|
||
|
||
|
||
|
||
|
||
|
||
<a class="navbar-brand logo" href="intro.html">
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<img src="_static/logo.png" class="logo__image only-light" alt="Applied Data Analysis and Machine Learning - Home"/>
|
||
<script>document.write(`<img src="_static/logo.png" class="logo__image only-dark" alt="Applied Data Analysis and Machine Learning - Home"/>`);</script>
|
||
|
||
|
||
</a></div>
|
||
<div class="sidebar-primary-item">
|
||
|
||
<script>
|
||
document.write(`
|
||
<button class="btn search-button-field search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
|
||
<i class="fa-solid fa-magnifying-glass"></i>
|
||
<span class="search-button__default-text">Search</span>
|
||
<span class="search-button__kbd-shortcut"><kbd class="kbd-shortcut__modifier">Ctrl</kbd>+<kbd class="kbd-shortcut__modifier">K</kbd></span>
|
||
</button>
|
||
`);
|
||
</script></div>
|
||
<div class="sidebar-primary-item"><nav class="bd-links bd-docs-nav" aria-label="Main">
|
||
<div class="bd-toc-item navbar-nav active">
|
||
|
||
<ul class="nav bd-sidenav bd-sidenav__home-link">
|
||
<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">Course setting</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: Linear Regression and Gradient descent</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: Gradient descent methods</a></li>
|
||
<li class="toctree-l1"><a class="reference internal" href="exercisesweek38.html">Exercises week 38</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 6 (midnight), 2025</a></li>
|
||
</ul>
|
||
|
||
</div>
|
||
</nav></div>
|
||
</div>
|
||
|
||
|
||
<div class="sidebar-primary-items__end sidebar-primary__section">
|
||
</div>
|
||
|
||
<div id="rtd-footer-container"></div>
|
||
|
||
|
||
</div>
|
||
|
||
<main id="main-content" class="bd-main" role="main">
|
||
|
||
|
||
|
||
<div class="sbt-scroll-pixel-helper"></div>
|
||
|
||
<div class="bd-content">
|
||
<div class="bd-article-container">
|
||
|
||
<div class="bd-header-article d-print-none">
|
||
<div class="header-article-items header-article__inner">
|
||
|
||
<div class="header-article-items__start">
|
||
|
||
<div class="header-article-item"><button class="sidebar-toggle primary-toggle btn btn-sm" title="Toggle primary sidebar" data-bs-placement="bottom" data-bs-toggle="tooltip">
|
||
<span class="fa-solid fa-bars"></span>
|
||
</button></div>
|
||
|
||
</div>
|
||
|
||
|
||
<div class="header-article-items__end">
|
||
|
||
<div class="header-article-item">
|
||
|
||
<div class="article-header-buttons">
|
||
|
||
|
||
|
||
|
||
|
||
<div class="dropdown dropdown-download-buttons">
|
||
<button class="btn dropdown-toggle" type="button" data-bs-toggle="dropdown" aria-expanded="false" aria-label="Download this page">
|
||
<i class="fas fa-download"></i>
|
||
</button>
|
||
<ul class="dropdown-menu">
|
||
|
||
|
||
|
||
<li><a href="_sources/linalg.ipynb" target="_blank"
|
||
class="btn btn-sm btn-download-source-button dropdown-item"
|
||
title="Download source file"
|
||
data-bs-placement="left" data-bs-toggle="tooltip"
|
||
>
|
||
|
||
|
||
<span class="btn__icon-container">
|
||
<i class="fas fa-file"></i>
|
||
</span>
|
||
<span class="btn__text-container">.ipynb</span>
|
||
</a>
|
||
</li>
|
||
|
||
|
||
|
||
|
||
<li>
|
||
<button onclick="window.print()"
|
||
class="btn btn-sm btn-download-pdf-button dropdown-item"
|
||
title="Print to PDF"
|
||
data-bs-placement="left" data-bs-toggle="tooltip"
|
||
>
|
||
|
||
|
||
<span class="btn__icon-container">
|
||
<i class="fas fa-file-pdf"></i>
|
||
</span>
|
||
<span class="btn__text-container">.pdf</span>
|
||
</button>
|
||
</li>
|
||
|
||
</ul>
|
||
</div>
|
||
|
||
|
||
|
||
|
||
<button onclick="toggleFullScreen()"
|
||
class="btn btn-sm btn-fullscreen-button"
|
||
title="Fullscreen mode"
|
||
data-bs-placement="bottom" data-bs-toggle="tooltip"
|
||
>
|
||
|
||
|
||
<span class="btn__icon-container">
|
||
<i class="fas fa-expand"></i>
|
||
</span>
|
||
|
||
</button>
|
||
|
||
|
||
|
||
<script>
|
||
document.write(`
|
||
<button class="btn btn-sm nav-link pst-navbar-icon theme-switch-button" title="light/dark" aria-label="light/dark" data-bs-placement="bottom" data-bs-toggle="tooltip">
|
||
<i class="theme-switch fa-solid fa-sun fa-lg" data-mode="light"></i>
|
||
<i class="theme-switch fa-solid fa-moon fa-lg" data-mode="dark"></i>
|
||
<i class="theme-switch fa-solid fa-circle-half-stroke fa-lg" data-mode="auto"></i>
|
||
</button>
|
||
`);
|
||
</script>
|
||
|
||
|
||
<script>
|
||
document.write(`
|
||
<button class="btn btn-sm pst-navbar-icon search-button search-button__button" title="Search" aria-label="Search" data-bs-placement="bottom" data-bs-toggle="tooltip">
|
||
<i class="fa-solid fa-magnifying-glass fa-lg"></i>
|
||
</button>
|
||
`);
|
||
</script>
|
||
<button class="sidebar-toggle secondary-toggle btn btn-sm" title="Toggle secondary sidebar" data-bs-placement="bottom" data-bs-toggle="tooltip">
|
||
<span class="fa-solid fa-list"></span>
|
||
</button>
|
||
</div></div>
|
||
|
||
</div>
|
||
|
||
</div>
|
||
</div>
|
||
|
||
|
||
|
||
<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 id="searchbox"></div>
|
||
<article class="bd-article">
|
||
|
||
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
|
||
doconce format html linalg.do.txt --><section class="tex2jax_ignore mathjax_ignore" 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="Link to this heading">#</a></h1>
|
||
<section id="introduction">
|
||
<h2><span class="section-number">2.1. </span>Introduction<a class="headerlink" href="#introduction" title="Link to this heading">#</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>
|
||
</section>
|
||
<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="Link to this heading">#</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>
|
||
</section>
|
||
<section id="basic-matrix-features">
|
||
<h2><span class="section-number">2.3. </span>Basic Matrix Features<a class="headerlink" href="#basic-matrix-features" title="Link to this heading">#</a></h2>
|
||
<p>Matrix properties reminder</p>
|
||
<div class="math notranslate nohighlight">
|
||
\[\begin{split}
|
||
\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}
|
||
\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><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="Link to this heading">#</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 > 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 < 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 > 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 < 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| > 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 > 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 < 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>
|
||
</section>
|
||
</section>
|
||
<section id="numpy-and-arrays">
|
||
<h2><span class="section-number">2.4. </span>Numpy and arrays<a class="headerlink" href="#numpy-and-arrays" title="Link to this heading">#</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
n = 10
|
||
x = np.random.normal(size=n)
|
||
print(x)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
x = np.array([1, 2, 3])
|
||
print(x)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
x = np.log(np.array([4, 7, 8]))
|
||
print(x)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
from math import log
|
||
x = np.array([4, 7, 8])
|
||
for i in range(0, len(x)):
|
||
x[i] = log(x[i])
|
||
print(x)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
x = np.log(np.array([4, 7, 8], dtype = np.float64))
|
||
print(x)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
x = np.log(np.array([4.0, 7.0, 8.0]))
|
||
print(x)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
x = np.log(np.array([4.0, 7.0, 8.0]))
|
||
print(x.itemsize)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
|
||
print(A)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
|
||
# print the first column, row-major order and elements start with 0
|
||
print(A[:,0])
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
A = np.log(np.array([ [4.0, 7.0, 8.0], [3.0, 10.0, 11.0], [4.0, 5.0, 7.0] ]))
|
||
# print the first column, row-major order and elements start with 0
|
||
print(A[1,:])
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
n = 10
|
||
# define a matrix of dimension 10 x 10 and set all elements to zero
|
||
A = np.zeros( (n, n) )
|
||
print(A)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
n = 10
|
||
# define a matrix of dimension 10 x 10 and set all elements to one
|
||
A = np.ones( (n, n) )
|
||
print(A)
|
||
</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-none notranslate"><div class="highlight"><pre><span></span>import numpy as np
|
||
n = 10
|
||
# define a matrix of dimension 10 x 10 and set all elements to random numbers with x \in [0, 1]
|
||
A = np.random.rand(n, n)
|
||
print(A)
|
||
</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} & \sigma_{xy} & \sigma_{xz} \\
|
||
\sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\
|
||
\sigma_{zx} & \sigma_{zy} & \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 & y_0 & z_0 \\
|
||
x_1 & y_1 & z_1 \\
|
||
x_2 & y_2 & z_2 \\
|
||
\dots & \dots & \dots \\
|
||
x_{n-2} & y_{n-2} & z_{n-2} \\
|
||
x_{n-1} & y_{n-1} & z_{n-1}
|
||
\end{bmatrix},
|
||
\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-none notranslate"><div class="highlight"><pre><span></span># Importing various packages
|
||
import numpy as np
|
||
|
||
n = 100
|
||
x = np.random.normal(size=n)
|
||
print(np.mean(x))
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
print(np.mean(y))
|
||
z = x**3+np.random.normal(size=n)
|
||
print(np.mean(z))
|
||
W = np.vstack((x, y, z))
|
||
Sigma = np.cov(W)
|
||
print(Sigma)
|
||
Eigvals, Eigvecs = np.linalg.eig(Sigma)
|
||
print(Eigvals)
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
<div class="cell docutils container">
|
||
<div class="cell_input docutils container">
|
||
<div class="highlight-none notranslate"><div class="highlight"><pre><span></span>%matplotlib inline
|
||
|
||
import numpy as np
|
||
import matplotlib.pyplot as plt
|
||
from scipy import sparse
|
||
eye = np.eye(4)
|
||
print(eye)
|
||
sparse_mtx = sparse.csr_matrix(eye)
|
||
print(sparse_mtx)
|
||
x = np.linspace(-10,10,100)
|
||
y = np.sin(x)
|
||
plt.plot(x,y,marker='x')
|
||
plt.show()
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</section>
|
||
<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="Link to this heading">#</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-none notranslate"><div class="highlight"><pre><span></span>"""
|
||
Simple code that tests various numpy functions
|
||
"""
|
||
|
||
import numpy as np
|
||
# Simple test-matrix of dim 3 x 4
|
||
a = np.array([ [1, 2, 3], [4, 5, 6], [7, 8, 9],[10, 11, 12]],dtype=np.float64)
|
||
print(f"The test matrix:{a}")
|
||
# This is the total mean summed over all elements, which here has to be 6.5
|
||
print(f"This is the total mean summed over all elements:{np.mean(a,dtype=np.float64)}")
|
||
# This is the mean for each column, it returns an array with the mean values for each column. It returns a row-like vector
|
||
print(f"This is the mean for each column:{np.mean(a, axis=0, keepdims=True,dtype=np.float64)}")
|
||
# This is the mean value for each row, it returns an array via the keepdims option which is a column-like vector if
|
||
# keepdims=True. Else it return a row-like vector
|
||
# Try setting keepdims=False
|
||
print(f"This is the mean value for each row:{np.mean(a, axis=1, keepdims=True,dtype=np.float64)}")
|
||
# We print then the mean value for each row by setting keepdims=False
|
||
print(f"This is the mean value for each row with keepdims false:{np.mean(a, axis=1, keepdims=False,dtype=np.float64)}")
|
||
</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-none notranslate"><div class="highlight"><pre><span></span># Ravel return a contiguous flattened array.
|
||
print(f"Flatten the matrix:{np.ravel(a)}")
|
||
# It is the same as reshaping the matrix into a one-dimensional array
|
||
print(f"Reshape the matrix to a one-dim array:{a.reshape(-1)}")
|
||
# ‘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.
|
||
# ‘F’ means to index the elements in column-major, Fortran-style order, with the first index changing fastest, and the last index changing slowest
|
||
print(np.ravel(a, order='F'))
|
||
# When order is ‘A’, it will preserve the array’s ‘C’ or ‘F’ ordering
|
||
# ‘A’ means to read the elements in Fortran-like index order if a is Fortran contiguous in memory, C-like order otherwise.
|
||
# ‘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.
|
||
# Transposing it
|
||
print(np.ravel(a.T))
|
||
print(np.ravel(a.T, order='A'))
|
||
</pre></div>
|
||
</div>
|
||
</div>
|
||
</div>
|
||
</section>
|
||
<section id="gaussian-elimination">
|
||
<h2><span class="section-number">2.6. </span>Gaussian Elimination<a class="headerlink" href="#gaussian-elimination" title="Link to this heading">#</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}& 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}.
|
||
\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}& 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},
|
||
\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} & 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}.
|
||
\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">&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>
|
||
<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="Link to this heading">#</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>
|
||
</section>
|
||
</section>
|
||
</section>
|
||
|
||
<script type="text/x-thebe-config">
|
||
{
|
||
requestKernel: true,
|
||
binderOptions: {
|
||
repo: "binder-examples/jupyter-stacks-datascience",
|
||
ref: "master",
|
||
},
|
||
codeMirrorConfig: {
|
||
theme: "abcdef",
|
||
mode: "python"
|
||
},
|
||
kernelOptions: {
|
||
name: "python3",
|
||
path: "./."
|
||
},
|
||
predefinedOutput: true
|
||
}
|
||
</script>
|
||
<script>kernelName = 'python3'</script>
|
||
|
||
</article>
|
||
|
||
|
||
|
||
|
||
|
||
|
||
<footer class="prev-next-footer d-print-none">
|
||
|
||
<div class="prev-next-area">
|
||
<a class="left-prev"
|
||
href="statistics.html"
|
||
title="previous page">
|
||
<i class="fa-solid 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"
|
||
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="fa-solid fa-angle-right"></i>
|
||
</a>
|
||
</div>
|
||
</footer>
|
||
|
||
</div>
|
||
|
||
|
||
|
||
<div class="bd-sidebar-secondary bd-toc"><div class="sidebar-secondary-items sidebar-secondary__inner">
|
||
|
||
|
||
<div class="sidebar-secondary-item">
|
||
<div class="page-toc tocsection onthispage">
|
||
<i class="fa-solid fa-list"></i> Contents
|
||
</div>
|
||
<nav class="bd-toc-nav page-toc">
|
||
<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>
|
||
<footer class="bd-footer-content">
|
||
|
||
<div class="bd-footer-content__inner container">
|
||
|
||
<div class="footer-item">
|
||
|
||
<p class="component-author">
|
||
By Morten Hjorth-Jensen
|
||
</p>
|
||
|
||
</div>
|
||
|
||
<div class="footer-item">
|
||
|
||
|
||
<p class="copyright">
|
||
|
||
© Copyright 2023.
|
||
<br/>
|
||
|
||
</p>
|
||
|
||
</div>
|
||
|
||
<div class="footer-item">
|
||
|
||
</div>
|
||
|
||
<div class="footer-item">
|
||
|
||
</div>
|
||
|
||
</div>
|
||
</footer>
|
||
|
||
|
||
</main>
|
||
</div>
|
||
</div>
|
||
|
||
<!-- Scripts loaded after <body> so the DOM is not blocked -->
|
||
<script src="_static/scripts/bootstrap.js?digest=dfe6caa3a7d634c4db9b"></script>
|
||
<script src="_static/scripts/pydata-sphinx-theme.js?digest=dfe6caa3a7d634c4db9b"></script>
|
||
|
||
<footer class="bd-footer">
|
||
</footer>
|
||
</body>
|
||
</html> |