Website build

Author: colinjcotter

Computationallinearalgebracourse.pdf

@@ -564,9 +564,7 @@ <h2><span class="section-number">1.5. </span>Inner products and orthogonality<a
 <p>Let <span class="math notranslate nohighlight">\(x,y\in \mathbb{C}^m\)</span>. Then the inner product of <span class="math notranslate nohighlight">\(x\)</span> and <span class="math notranslate nohighlight">\(y\)</span> is</p>
 <div class="math notranslate nohighlight">
 \[x^*y = \sum_{i=1}^m \bar{x}_iy_i.\]</div>
-</div></div><p>(Exercise: check that the inner product is bilinear, i.e. linear in
-both of the arguments.)</p>
-<p>We will frequently use the natural norm derived from the inner product
+</div></div><p>We will frequently use the natural norm derived from the inner product
 to define size of vectors.</p>
 <div class="proof proof-type-definition" id="id19">
 

L1_preliminaries.html

@@ -195,14 +195,14 @@ <h2><span class="section-number">2.2. </span>QR factorisation by classical Gram-
         <span class="proof-type">Exercise 2.4</span>
 
     </div><div class="proof-content">
-<p><span class="math notranslate nohighlight">\((\ddagger)\)</span> The <a class="reference internal" href="cla_utils.html#cla_utils.exercises2.GS_classical" title="cla_utils.exercises2.GS_classical"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises2.GS_classical()</span></code></a> function has been
-left unimplemented. It should implement the classical Gram-Schmidt
-algorithm above, using Numpy slice notation so that only one Python
-for loop is used. The function should work “in place” by making a
-copy of <span class="math notranslate nohighlight">\(A\)</span> and then changing the values in the copy, without
-introducing additional intermediate arrays (you will need to create
-a new array to store <span class="math notranslate nohighlight">\(R\)</span>). The test script <code class="docutils literal notranslate"><span class="pre">test_exercises2.py</span></code>
-in the <code class="docutils literal notranslate"><span class="pre">test</span></code> directory will test this function.</p>
+<p><span class="math notranslate nohighlight">\((\ddagger)\)</span> The <a class="reference internal" href="cla_utils.html#cla_utils.exercises2.GS_classical" title="cla_utils.exercises2.GS_classical"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises2.GS_classical()</span></code></a> function
+has been left unimplemented. It should implement the classical
+Gram-Schmidt algorithm above, using Numpy slice notation so that
+only one Python for loop is used. The function should work “in
+place” by changing the values in <span class="math notranslate nohighlight">\(A\)</span>, without introducing
+additional intermediate arrays (you will need to create a new array
+to store <span class="math notranslate nohighlight">\(R\)</span>). The test script <code class="docutils literal notranslate"><span class="pre">test_exercises2.py</span></code> in the
+<code class="docutils literal notranslate"><span class="pre">test</span></code> directory will test this function.</p>
 </div></div><div class="admonition hint">
 <p class="admonition-title">Hint</p>
 <p>The <span class="math notranslate nohighlight">\((\ddagger)\)</span> symbol in an exercise indicates that the code for
@@ -321,9 +321,9 @@ <h2><span class="section-number">2.4. </span>Modified Gram-Schmidt<a class="head
 algorithm above, using Numpy slice notation where possible.
 What is the minimal number of Python
 for loops possible?</p>
-<p>The function should work “in place” by making a
-copy of <span class="math notranslate nohighlight">\(A\)</span> and then changing those values, without introducing
-additional intermediate arrays. The test script
+<p>The function should work “in place” by changing the values in <span class="math notranslate nohighlight">\(A\)</span>,
+without introducing additional intermediate arrays (you will need
+to create a new array to store <span class="math notranslate nohighlight">\(R\)</span>). The test script
 <code class="docutils literal notranslate"><span class="pre">test_exercises2.py</span></code> in the <code class="docutils literal notranslate"><span class="pre">test</span></code> directory will test this
 function.</p>
 </div></div><div class="proof proof-type-exercise" id="id6">
@@ -362,7 +362,7 @@ <h2><span class="section-number">2.5. </span>Modified Gram-Schmidt as triangular
 \underbrace{
 \begin{pmatrix}
 \frac{1}{r_{11}} &amp; -\frac{r_{12}}{r_{11}} &amp; \ldots &amp;
-\ldots &amp; -\frac{r_{11}}{r_{11}} \\
+\ldots &amp; -\frac{r_{1n}}{r_{11}} \\
 0 &amp; 1 &amp; 0 &amp; \ldots &amp; 0 \\
 0 &amp; 0 &amp; 1 &amp; \ldots &amp; 0 \\
 \vdots &amp; \ddots &amp; \ddots &amp; \ldots &amp; \vdots \\
@@ -392,7 +392,7 @@ <h2><span class="section-number">2.5. </span>Modified Gram-Schmidt as triangular
 \begin{pmatrix}
 1 &amp; 0 &amp; 0 &amp;
 \ldots &amp; 0 \\
-0 &amp; \frac{1}{r_{22}} &amp; -\frac{r_{23}}{r_{22}} &amp; \ldots &amp; -\frac{r_{2n}}{r_{nn}} \\      0 &amp; 0 &amp; 1 &amp; \ldots &amp; 0 \\
+0 &amp; \frac{1}{r_{22}} &amp; -\frac{r_{23}}{r_{22}} &amp; \ldots &amp; -\frac{r_{2n}}{r_{22}} \\      0 &amp; 0 &amp; 1 &amp; \ldots &amp; 0 \\
 \vdots &amp; \ddots &amp; \ddots &amp; \ldots &amp; \vdots \\
 0 &amp; 0 &amp; 0 &amp; \ldots &amp; 1 \\
 \end{pmatrix}}_{R_2}
@@ -660,7 +660,30 @@ <h2><span class="section-number">2.6. </span>Householder triangulation<a class="
 <summary>
 Supplementary video</summary><div class="video_wrapper" style="">
 <iframe allowfullscreen="true" src="https://player.vimeo.com/video/450202242" style="border: 0; height: 345px; width: 560px">
-</iframe></div></details><div class="proof proof-type-exercise" id="id9">
+</iframe></div></details><p>In this section we will frequently encounter systems of the form</p>
+<div class="math notranslate nohighlight">
+\[\hat{R}x = y.\]</div>
+<p>This is an upper triangular system that can be solved efficiently
+using back-substitution.</p>
+<p>Written in components, this equation is</p>
+<blockquote>
+<div><div class="math notranslate nohighlight">
+\[ \begin{align}\begin{aligned}R_{11}x_1 + R_{12}x_2 + \ldots + R_{1(m-1)}x_{m-1} + R_{1m}x_m = y_1,\\0x_1 + R_{22}x_2 + \ldots + R_{2(m-1)}x_{m-1} + R_{2m}x_m = y_2,\\\vdots\\0x_1 + 0x_2 + \ldots + R_{(m-1)(m-1)}x_{m-1} + R_{(m-1)m}x_m = y_{m-1},\\ 0x_1 + 0x_2 + \ldots + 0x_{m-1} + R_{mm}x_m = y_{m}.\end{aligned}\end{align} \]</div>
+</div></blockquote>
+<p>The last equation yields <span class="math notranslate nohighlight">\(x_m\)</span> directly by dividing by <span class="math notranslate nohighlight">\(R_{mm}\)</span>, then
+we can use this value to directly compute <span class="math notranslate nohighlight">\(x_{m-1}\)</span>. This is repeated
+for all of the entries of <span class="math notranslate nohighlight">\(x\)</span> from <span class="math notranslate nohighlight">\(m\)</span> down to 1. This procedure is
+called back substitution, which we summarise in the following
+pseudo-code.</p>
+<ul class="simple">
+<li><p><span class="math notranslate nohighlight">\(x_m  \gets y_m/R_{mm}\)</span></p></li>
+<li><p>FOR <span class="math notranslate nohighlight">\(i= m-1\)</span> TO 1 (BACKWARDS)</p>
+<ul>
+<li><p><span class="math notranslate nohighlight">\(x_i \gets (y_i - \sum_{k=i+1}^mR_{ik}x_k)/R_{ii}\)</span></p></li>
+</ul>
+</li>
+</ul>
+<div class="proof proof-type-exercise" id="id9">
 
     <div class="proof-title">
         <span class="proof-type">Exercise 2.9</span>
@@ -709,11 +732,7 @@ <h2><span class="section-number">2.6. </span>Householder triangulation<a class="
 <a class="reference internal" href="cla_utils.html#cla_utils.exercises3.householder" title="cla_utils.exercises3.householder"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises3.householder()</span></code></a>.  If you have not already
 done so, you will need to modified
 <a class="reference internal" href="cla_utils.html#cla_utils.exercises3.householder" title="cla_utils.exercises3.householder"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises3.householder()</span></code></a> to use the <code class="docutils literal notranslate"><span class="pre">kmax</span></code>
-argument. You may make use of the built-in triangular solve
-algorithm <a class="reference external" href="http://scipy.github.io/devdocs/reference/generated/scipy.linalg.solve_triangular.html#scipy.linalg.solve_triangular" title="(in SciPy v1.10.0.dev0+1717.ebbd597)"><code class="xref py py-func docutils literal notranslate"><span class="pre">scipy.linalg.solve_triangular()</span></code></a> (we shall consider
-triangular matrix algorithms briefly later). The test script
-<code class="docutils literal notranslate"><span class="pre">test_exercises3.py</span></code> in the <code class="docutils literal notranslate"><span class="pre">test</span></code> directory will also test this
-function.</p>
+argument. You will also need <a class="reference internal" href="cla_utils.html#cla_utils.exercises3.solve_U" title="cla_utils.exercises3.solve_U"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises3.solve_U()</span></code></a>.</p>
 </div></div><p>If we really need <span class="math notranslate nohighlight">\(Q\)</span>, we can get it by matrix-vector products with
 each element of the canonical basis <span class="math notranslate nohighlight">\((e_1,e_2,\ldots,e_n)\)</span>.  This
 means that first we need to compute a matrix-vector product <span class="math notranslate nohighlight">\(Qx\)</span> with
@@ -809,24 +828,7 @@ <h2><span class="section-number">2.7. </span>Application: Least squares problems
 \[\hat{R}x = \hat{Q}^*b.\]</div>
 </div></blockquote>
 <p>This is an upper triangular system that can be solved efficiently
-using back-substitution. Written in components, this equation is</p>
-<blockquote>
-<div><div class="math notranslate nohighlight">
-\[ \begin{align}\begin{aligned}R_{11}x_1 + R_{12}x_2 + \ldots + R_{1(m-1)}x_{m-1} + R_{1m}x_m = y_1,\\0x_1 + R_{22}x_2 + \ldots + R_{2(m-1)}x_{m-1} + R_{2m}x_m = y_2,\\\vdots\\0x_1 + 0x_2 + \ldots + R_{(m-1)(m-1)}x_{m-1} + R_{(m-1)m}x_m = y_{m-1},\\ 0x_1 + 0x_2 + \ldots + 0x_{m-1} + R_{mm}x_m = y_{m}.\end{aligned}\end{align} \]</div>
-</div></blockquote>
-<p>The last equation yields <span class="math notranslate nohighlight">\(x_m\)</span> directly by dividing by <span class="math notranslate nohighlight">\(R_{mm}\)</span>, then
-we can use this value to directly compute <span class="math notranslate nohighlight">\(x_{m-1}\)</span>. This is repeated
-for all of the entries of <span class="math notranslate nohighlight">\(x\)</span> from <span class="math notranslate nohighlight">\(m\)</span> down to 1. This procedure is
-called back substitution, which we summarise in the following
-pseudo-code.</p>
-<ul class="simple">
-<li><p><span class="math notranslate nohighlight">\(x_m  \gets y_m/R_{mm}\)</span></p></li>
-<li><p>FOR <span class="math notranslate nohighlight">\(i= m-1\)</span> TO 1 (BACKWARDS)</p>
-<ul>
-<li><p><span class="math notranslate nohighlight">\(x_i \gets (y_i - \sum_{k=i+1}^mR_{ik}x_k)/R_{ii}\)</span></p></li>
-</ul>
-</li>
-</ul>
+using back-substitution.</p>
 <div class="proof proof-type-exercise" id="id12">
 
     <div class="proof-title">
@@ -840,7 +842,7 @@ <h2><span class="section-number">2.7. </span>Application: Least squares problems
 by forming an appropriate augmented matrix <span class="math notranslate nohighlight">\(\hat{A}\)</span>, calling
 <a class="reference internal" href="cla_utils.html#cla_utils.exercises3.householder" title="cla_utils.exercises3.householder"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises3.householder()</span></code></a> and extracting appropriate
 subarrays using slice notation, before using
-<a class="reference external" href="http://scipy.github.io/devdocs/reference/generated/scipy.linalg.solve_triangular.html#scipy.linalg.solve_triangular" title="(in SciPy v1.10.0.dev0+1717.ebbd597)"><code class="xref py py-func docutils literal notranslate"><span class="pre">scipy.linalg.solve_triangular()</span></code></a> to solve the resulting upper
+<a class="reference internal" href="cla_utils.html#cla_utils.exercises3.solve_U" title="cla_utils.exercises3.solve_U"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises3.solve_U()</span></code></a> to solve the resulting upper
 triangular system, before returning the solution <span class="math notranslate nohighlight">\(x\)</span>. The test
 script <code class="docutils literal notranslate"><span class="pre">test_exercises3.py</span></code> in the <code class="docutils literal notranslate"><span class="pre">test</span></code> directory will also
 test this function.</p>

L2_QR_factorisation.html

@@ -731,8 +731,10 @@ <h2><span class="section-number">3.9. </span>Stability<a class="headerlink" href
 <div><div class="math notranslate nohighlight">
 \[\frac{\|\tilde{f}(x)-f(x)\|}{\|f(x)\|} \leq C\varepsilon,\]</div>
 </div></blockquote>
-<p>for sufficiently small <span class="math notranslate nohighlight">\(\varepsilon\)</span>. We shall see below that
-we have to lower our aspirations depending on the condition number of <span class="math notranslate nohighlight">\(A\)</span>.</p>
+<p>for sufficiently small <span class="math notranslate nohighlight">\(\varepsilon\)</span> (assuming, albeit
+unrealistically, that we have a sequence of computers with smaller and
+smaller $varepsilon$). We shall see below that we have to lower our
+aspirations depending on the condition number of <span class="math notranslate nohighlight">\(A\)</span>.</p>
 <div class="proof proof-type-definition" id="id21">
 
     <div class="proof-title">
@@ -741,14 +743,12 @@ <h2><span class="section-number">3.9. </span>Stability<a class="headerlink" href
             <span class="proof-title-name">(Stability)</span>
 
     </div><div class="proof-content">
-<p>An algorithm <span class="math notranslate nohighlight">\(\tilde{f}\)</span> for <span class="math notranslate nohighlight">\(f\)</span> is stable if for each <span class="math notranslate nohighlight">\(x\in X\)</span>,</p>
-<blockquote>
-<div><div class="math notranslate nohighlight">
-\[\frac{\|\tilde{f}(x)-f(\tilde{x})\|}{\|f(\tilde{x})\|} = \mathcal{O}(\varepsilon),\]</div>
-</div></blockquote>
-<p>there exists <span class="math notranslate nohighlight">\(\tilde{x}\)</span> with</p>
+<p>An algorithm <span class="math notranslate nohighlight">\(\tilde{f}\)</span> for <span class="math notranslate nohighlight">\(f\)</span> is stable if for each <span class="math notranslate nohighlight">\(x\in X\)</span>,
+there exists <span class="math notranslate nohighlight">\(\tilde{x}\)</span> with</p>
 <blockquote>
 <div><div class="math notranslate nohighlight">
+\[ \begin{align}\begin{aligned}\frac{\|\tilde{f}(x)-f(\tilde{x})\|}{\|f(\tilde{x})\|} = \mathcal{O}(\varepsilon),\\and\end{aligned}\end{align} \]</div>
+<div class="math notranslate nohighlight">
 \[\frac{\|\tilde{x}-x\|}{\|x\|} = \mathcal{O}(\varepsilon).\]</div>
 </div></blockquote>
 </div></div><p>We say that a stable algorithm gives nearly the right answer to nearly the

L3_analysing_algorithms.html

@@ -499,7 +499,7 @@ <h2><span class="section-number">5.4. </span>Rayleigh quotient<a class="headerli
 <li><p>Finding an eigenvector <span class="math notranslate nohighlight">\(v\)</span> of <span class="math notranslate nohighlight">\(A\)</span> with eigenvalue <span class="math notranslate nohighlight">\(\lambda\)</span> (you can use <a class="reference external" href="https://numpy.org/doc/stable/reference/generated/numpy.linalg.eig.html#numpy.linalg.eig" title="(in NumPy v1.23)"><code class="xref py py-func docutils literal notranslate"><span class="pre">numpy.linalg.eig()</span></code></a> for this),</p></li>
 <li><p>Choosing a perturbation vector <span class="math notranslate nohighlight">\(r\)</span>, and perturbation parameter <span class="math notranslate nohighlight">\(\epsilon&gt;0\)</span>,</p></li>
 <li><p>Comparing the Rayleigh quotient of <span class="math notranslate nohighlight">\(v + \epsilon r\)</span> with <span class="math notranslate nohighlight">\(\lambda\)</span>,</p></li>
-<li><p>Plotting (on a log-log graph, use <a class="reference external" href="https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.loglog.html#matplotlib.pyplot.loglog" title="(in Matplotlib v3.5.3)"><code class="xref py py-func docutils literal notranslate"><span class="pre">matplotlib.pyplot.loglog()</span></code></a>) the error in estimating the eigenvalue as a function of <span class="math notranslate nohighlight">\(\epsilon\)</span>.</p></li>
+<li><p>Plotting (on a log-log graph, use <a class="reference external" href="https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.loglog.html#matplotlib.pyplot.loglog" title="(in Matplotlib v3.6.0)"><code class="xref py py-func docutils literal notranslate"><span class="pre">matplotlib.pyplot.loglog()</span></code></a>) the error in estimating the eigenvalue as a function of <span class="math notranslate nohighlight">\(\epsilon\)</span>.</p></li>
 </ol>
 <p>The best way to do this is to plot the computed data values as points,
 and then superpose a line plot of <span class="math notranslate nohighlight">\(a\epsilon^k\)</span> for appropriate
@@ -773,7 +773,7 @@ <h2><span class="section-number">5.8. </span>The pure QR algorithm<a class="head
 <a class="reference internal" href="cla_utils.html#cla_utils.exercises9.get_B100" title="cla_utils.exercises9.get_B100"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises9.get_B100()</span></code></a>,
 <a class="reference internal" href="cla_utils.html#cla_utils.exercises9.get_C100" title="cla_utils.exercises9.get_C100"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises9.get_C100()</span></code></a>, and
 <a class="reference internal" href="cla_utils.html#cla_utils.exercises9.get_D100" title="cla_utils.exercises9.get_D100"><code class="xref py py-func docutils literal notranslate"><span class="pre">cla_utils.exercises9.get_D100()</span></code></a>. You can use
-<a class="reference external" href="https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.pcolor.html#matplotlib.pyplot.pcolor" title="(in Matplotlib v3.5.3)"><code class="xref py py-func docutils literal notranslate"><span class="pre">matplotlib.pyplot.pcolor()</span></code></a> to visualise the entries,
+<a class="reference external" href="https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.pcolor.html#matplotlib.pyplot.pcolor" title="(in Matplotlib v3.6.0)"><code class="xref py py-func docutils literal notranslate"><span class="pre">matplotlib.pyplot.pcolor()</span></code></a> to visualise the entries,
 or compute norms of the components of the matrices below the diagonal,
 for example. What do you observe? How does this relate to the structure
 of the four matrices?</p>