Floating Point Numbers are Not Associative

On a recent exam, I asked students for suggestions on how to test for equality of two floating point numbers.

Loren Shure and Mike Croucher present simple examples which highlight the perils. For example, (0.1 + 0.2) + 0.3 is not equal to 0.1 + (0.2 + 0.3).

As pointed out in a comment, this also implies:

octave:2> x = 0.3 + 0.2 + 0.1;
octave:3> y = 0.1 + 0.2 + 0.3;
octave:4> x==y
ans = 0
octave:5> x-y
ans = -1.1102e-16

Tridiagonal Matrices

Tridiagonal matrices have nonzero element only in a narrow band around the diagonal \[A = \begin{bmatrix} b_1 & c_1 & 0 & ... & 0\\ a_2 & b_2 & c_2 & ... & 0 \\ 0 & & \ddots & & 0 \\ 0 & ... & a_{n-1} & b_{n-1} & c_{n-1} \\ 0 & & ... & a_n & b_n\\ \end{bmatrix}.\] Solving a system of linear equations \(Ax=b\) can be done very efficiently by the so called Thomas algorithm, named after the physicist and computational scientist Llewellyn Thomas.

Since the algorithm is very well documented on wikipedia, I won't repeat it here. It is essentially a specialization of the standard Gauss elimination method to the special structure of tridiagonal matrices.

It is an order \(\mathcal{O}(n)\) algorithm, which means it is very efficient compared to \(\mathcal{O}(n^3)\) for a general matrix.

The attribution to Llewellyn Thomas is somewhat interesting. Cleve Moler writes:

Llewellyn H. Thomas is a distinquished physicist who in the 50's held positions at Columbia University and at IBM's Watson Research Laboratory when it was located adjacent to the Columbia campus. He is probably best known in connection with the Thomas-Fermi electron gas model. 

The so-called Thomas Algorithm is indeed just a form of elimination for solving tridiagonal systems of linear equations. But it usually associated with the systems that arise from finite difference approximations to partial differential equations. The attribution to Thomas seems to be more common in some engineering disciplines than it is in numerical analysis. 

W.F. Ames writes in his book [1]: "The method we describe was discovered independently by many and has been called the Thomas algorithm (see [2]) by (David) Young. Its general description first appeared in widely distributed published form in an article by Bruce et al. [3]."

Nye-Ham Debate

So, I finally finished watching the Bill Nye-Ken Ham "evolution v/s creationism" debate (YouTube link), just to see what the fuss was all about. While the world wide internet is saturated with opinions on who "won", I found the debate itself to be quite disappointing - and would be surprised if it shifted any opinions.

It was generally a diffuse debate (like US presidential debates); the only pointed moment, in my opinion, came when an audience member asked Ken Ham whether he could imagine the form of any new evidence that would compel him to abandon his theory. He said "no", which essentially means his theory is non-falsifiable, and hence mostly a fairy tale. (For a moment I wondered what if God himself came down and told him that this literal interpretation was stretching it! Wouldn't that be compelling, even within the restrictive set of assumptions he works with?)

Example: Sherman-Morrison Formula

Consider the following equation \(Ax=b\),\[\begin{bmatrix}
 1 &  2 &  3 \\
 1 &  4 &  9 \\
 4 &  5 &  7 \\
\end{bmatrix}x = \begin{bmatrix}
 2 \\
 4 \\
 6 \\
\end{bmatrix}.\] The solution of this equation is,\[x = \begin{bmatrix}
 0.75 \\
 0.25 \\
 0.25 \\
\end{bmatrix}.\] Now let us suppose we wanted to solve a similar problem ( \(A'x' = b\) )with a slightly perturbed matrix,\[A' =  \begin{bmatrix}
 1 &  2 &  3 \\
 1 &  6 &  9 \\
 4 &  5 &  7 \\
\end{bmatrix}.\] Notice that this matrix differs from \(A\) in the (2,2) element (6 instead of 4).

The Sherman-Morrison formula tells us to form \(u = [0~1~0]^T, v = [0~~2~~0]^T\), so that \[A' = A +  u v^T.\] Hence \[y = x = A^{-1}b = \begin{bmatrix}
 0.75 \\
 0.25 \\
 0.25 \\
\end{bmatrix},\text{ and } z = A^{-1} u = \begin{bmatrix}
 0.125 \\
 -0.625 \\
 0.375 \\
\end{bmatrix}.\] Thus,\[ x' = y - \left({v^T y \over 1 + v^T z}\right) z\]  This implies \[x' = \begin{bmatrix}
 0.75 \\
 0.25 \\
 0.25 \\
\end{bmatrix} \frac{0.25}{1-0.625} \begin{bmatrix}
 0.125 \\
 -0.625 \\
 0.375 \\
\end{bmatrix} =\begin{bmatrix}
   1 \\
  -1 \\
   1 \\
\end{bmatrix},\] which is indeed the correct solution.

Retirement of Scientific Ideas

The Edge has an interesting collection of scientific ideas that are ready for retirement (H/T Skeptic's Guide to the Universe).

As expected, it includes its fair share of straw men (ideas which are already dead or were never seriously considered in the scientific community), opinions with a bias induced by the author's competing theory, or some potential conflict of interest, and some truly thought provoking ones (especially in areas far removed from one's own scientific work/reading).

Anyway, the list has been around for a while (I bumped into it only recently). It is worth a read.

The Sherman-Morrison Formula

In a recent post, I mentioned how the upfront cost of LU decomposition pays for itself when solving a system of linear equations \[Ax=b_i, \quad \quad i = 1, 2, ... N.\] Now suppose the (presumably large) matrix \(A\) was modified slightly to form a new matrix \(A'\). Let us assume that we had a LU factorization of the original unmodified matrix, which would have allowed us to solve \(Ax=b\) very efficiently.

If we now wanted to solve \(A'x=b\), can we somehow avoid having to perform a brand new LU factorization on the modified matrix \(A'\)?

Recall that the factorization itself doesn't come very cheap.

The Sherman-Morrison formula (SMF) says "sure can!" as long as the modification can be expressed as an outer-product of two vectors. That is, \[A' = A + uv^T.\] As an example, suppose that the element in the \((j,k)\) position \(a_{jk}\) was changed to \(a'_{kj}\). This modification can be expressed as,\[A' = A + (a'_{jk} - a_{jk}) \mathbf{e}_{j} \mathbf{e}_{k}^T,\] where \(\mathbf{e}_{j}\) is a vector with "1" in the row \(j\) and zeroes everywhere else.

The Sherman-Morrison formula says \[(A')^{-1} = (A+uv^T)^{-1} = A^{-1} - {A^{-1}uv^T A^{-1} \over 1 + v^T A^{-1}u}.\] The formula can be readily verified, as shown on wikipedia.

Note: The scalar analog to question the SMF seeks to address is: given \(1/a\), how can you find \(1/(a+b)\)?

The SMF answer would be: \[1/a \times \left(1- {b/a \over 1 + b/a} \right).\] Returning to the original problem \(A'x = b \implies x = (A')^{-1} b\). Thus \[x = A^{-1}b - {A^{-1}uv^T A^{-1}b \over 1 + v^T A^{-1}u}.\] We can set \(y = A^{-1}b\), and \(z=A^{-1}u\). Since a LU factorization is available for \(A\), we can solve for \(y\) and \(z\) in \(\mathcal{O}(m^2)\) operations, where \(m\) is the size of the square matrix \(A\).

This gives us \[x = y - \left({v^T y \over 1 + v^T z}\right) z.\] In general, this can be performed in \(\mathcal{O}(3 m^2)\) operations, as opposed to to \(\mathcal{O}(m^3)\) for a fresh factorization.

Issues with MathJax and Blogger Preview

In the past few weeks, I noticed that I was having trouble with MathJax (which I have come to love) and Blogger. While there is a plethora of problems being reported on the internet, my symptoms were

1. MathJax and Blogger seemed to work fine for the most part
2. However, while composing a new post (with LaTeX commands in it), and previewing it using Blogger's Preview feature, the LaTeX would not be rendered.
3. Once I published the blog, the rendering would work fine.

After some googling, I found that problem has a somewhat easy - albeit manual fix. The "preview" page normally has an "https://" prefix. Changing that to "http://" in the browser seems to solve the problem.

For example, upon normal preview:

After the change, the preview changed to:

I hope this is a temporary band-aid, because it is still a little annoying!

EDIT: Changing the source of the MathJax engine to https as pointed out here does fix it for good. Just swap the old "src=http" part with the new "src=https" part, as described in the link.

LU Decomposition: Analogy

Systems of linear equations \(Ax = b\) in which only the "right hand side" \(b\) changes between successive problems occur very frequently in science and engineering. In such cases, it is prudent to factorize or decompose the matrix \(A = LU\).

This factorization is often viewed as the matrix encoding of Gauss elimination (GE).

It front-loads the cost of having to solve the series of N matrix equations, \(Ax = b_1, Ax = b_2, ..., Ax = b_N\).

An imperfect analogy might be a subscription to Amazon Prime, which costs (let's say) $100/year upfront (the factorization), and gives subscribers free shipping on all products. If you buy a lot of stuff from Amazon (large number of matrix equations N), the up front cost pays for itself very quickly.

The only difference is that the deal offered by LU decomposition beats anything on Amazon.

For a square \(m \times m\) matrix \(A\), the computational cost of LU decomposition is \(\mathcal{O}(m^3)\). Once the factorization is available, the matrix equation can be solved for a particular right hand side using forward and back substitution steps in \(\mathcal{O}(m^2)\) operations.

Thus, the total cost of solving the series of N equations is \(\mathcal{O}(m^3 + N m^2)\). If GE is used to solve all the equations, the corresponding cost is \(\mathcal{O}(N m^3)\), which is much worse.

Links: Things I did not know edition

1. All Froot Loops taste the same! The Reality Check had a segment on it in a recent episode.

2. The earth year was not always 365 days long. No, we are not quibbling about the 1/4, plus or minus. Long time ago it was more than 400 days long, and the moon appeared 10 times larger in our skies. For some reason, I felt I should have known this earlier. Thanks to Radiolab.

Moler: Complex Step Differentiation

This is an old post by Cleve Moler, but given the amount of attention I've devoted to numerical differentiation, I ought to include a reference to this technique.

It works for analytic functions (infinitely differentiable), and consists of taking a small step in the imaginary axis. Thus an \(\mathcal{O}(h^2)\) method is given by,
\[f(x_0) = Im(f(x_0 + i h))/h)\]
A particularly useful feature of this algorithm is that as the step-size "h" decreases, it is not susceptible to round-off error like most other finite difference based schemes.