Links:

1. The Psychology of Small Packages (H/T Farnam Street)

2.  Wilkinson's Matrices

3. The Grammar of Anarchy

Linear Least Squares: What technique should I use? (part 1)

The quintessential least-squares problem looks like \[A x = b,\] where \(A_{m \times n}\) is a "tall" matrix with \(m > n\), \(x_{n \times 1}\) is the vector we are trying to "regress", and \(b_{m \times 1}\) is a vector as tall as \(A\).

Since we have too many equations or constraints (m) compared with unknowns (n), we can't solve the problem in the usual sense. Hence we seek a least-squares solution, which can be obtained by hitting both sides of the original equation with the transpose \(A'\), to obtain the so called normal equation, \[A' A x = A' b.\]Now we have a linear system \(n\) equations and \(n\) unknowns, which can be solved for \(x\).

There are three methods one can use:

1. Direct solution of Normal Equations: The idea is to factorize \(A' A\) using Cholesky decomposition, and solve the resulting system using forward and back substitution. This is straightforward, fast (asymptotic cost \(\sim m n^2 + n^3/3\)), but susceptible to round-off. We will explore this last part in a little bit.
2. QR Decomposition: The idea is to first factorize \(A = QR\). Then, we can write the normal equations as \[A' A x = b \implies R' Q'  Q R x = R' Q' b.\] Since \(Q\) is orthogonal, we can simplify the equations to the triangular system \[R x = Q' b.\] The asymptotic cost is \(2mn^2 - 2n^3/3\). This method is less susceptible to round-off.
3. SVD: We can always bring out the ultimate thermonuclear weapon. We first factorize \(A = U S V'\). The normal equations can then be simplified as above into \[S V' x = U' b.\]The asymptotic cost is \(2mn^2 + 11 n^3\). This method is also less susceptible to round-off. An advantage is that, if required, we can throw in some regularization at a problem to make it less nasty.

In the usual case when \(m \gg n\), the cost of SVD and QR is comparable. However as \(m \sim n\), SVD becomes more than an order of magnitude slower.

Let us explore the relative merits of the three methods by considering a Matlab example.

Chris Jordan

So I was at a dinner this evening, where this guy Chris Jordan gave an amazing talk.

You really should check out his artwork. He specializes in "garbage" - using photographs to give us a sense of the maddening material consumption in our society. It is funny how pictures can make big mind-numbing numbers (of discarded plastic bottles, cellphones, cars etc.) spring to life; how they replace apathy with grief.

His remarkable ability to make us "feel" these big numbers is captured in this TED talk.

A particularly poignant story that he told revolves around the albatrosses on Midway, and the amount of plastic found in the guts of dead baby albatrosses. As a friend apparently said to him, "if I start crying, I am worried I will never stop."

Conditional Column Switching in Octave or Matlab

Consider a simple Matlab/Octave exercise.

I have an N by 2 vector x. Thus, x(i,1) and x(i,2) are the entries of the i-th row. Say I want to rearrange x so that the first entry in any row is larger than the second. That is, I want to switch the entries of a row, only if a certain condition is met.

It is easy to write this in a element-by-element fashion

for i = 1:N
  if x(i,2) > x(i,1)
    c = x(i,1);
    x(i,1) = x(i,2);
    x(i,2) = c;
  endif
endfor

But we all know that element-by-element operations in Matlab or Octave suck! Here is another vectorized method:


SwitchCond = x(:,2) < x(:,1);
x(SwitchCond,[1,2]) = x(SwitchCond,[2,1]);

Besides being faster, it is also more succinct.

Links

1. The Art of Splash: Beautiful pictures of fluids in motion!

2. Nice back and forth between Deborah Meier and Eric Hanushek at Bridging Differences.

3. Yet another writing guide, with some nice references.