calling all Python gurus...

[Ether]

Originally Posted by Jared Russell

I used the cvxopt example implementation of L1 norm minimization on this page.

Hi Jared,

As I mentioned to you privately, I was able to duplicate your Python/cvxopt results exactly on a Debian Linux installation of Python/cvxopt.

Since then I've had some free time and came back to the challenge of getting it to work in Octave (MatLab).

The solution turned out to be to convert the original nonlinear/nondifferentiable problem -- min(sum(abs(Ax-b))) -- into a linear program and solve it with Octave's glpk (or MatLab's linprog) LP solver.

According to tic/toc, Octave/glpk took 0.125 seconds to solve on my 8-year-old desktop running XP.

The resulting solution, however, is interestingly different from yours. Although the L1 norm of the residuals of both solutions are virtually identical (to the 4th decimal place), the L2 norm of the glpk LAD residuals is 0.2 smaller than cvxopt's.

Although min L2 norm of residuals is mathematically unique, L1 is not. I assume that accounts for the difference.

From all the data that I've looked at so far, min L1 norm of residuals (Least Absolute Deviation or LAD) appears to be a slightly better metric than L2 (Least Squares). Now that efficient methods of computing LAD have been identified perhaps that opens a whole new world of fun number crunching for the stats mavens here on CD.

Are there any R, Mathematica, SciLab, etc gurus out there who would be interested in doing the LAD computation and posting results?

I also got an implementation using fmin_cobyla working, but it is much slower (and/or needs some parameter tuning). Minutes instead of sub-second.

That doesn't surprise me. I think any approach that doesn't first convert the LAD problem to an LP problem is going to be much more computation-intensive.