Wednesday, January 30, 2013

A Randomized Parallel Algorithm with Run Time $O(n^2)$ for Solving an $n \times n$ System of Linear Equations - implementation -



Around the time Curiosity landed on Mars, an Around the blogs in 80 summer hours pointed to Dick Lipton's blog's entry entitled A New Way To Solve Linear Equations where he mentioned a stunning preprint by Prasad Raghavendra. Yesterday, Riccardo Murri commented the following: 

Benjamin Jonen and I have implemented some example code in Python for this solver; you can read the source files at:
https://code.google.com/p/gc3pie/source/browse/#svn%2Ftrunk%2Fgc3pie%2Fexamples%2Foptimizer%2Flinear-systems-solver
There’s also an implementation of the algorithm as modified by Fliege (arXiv 1209.3995v1), which works over the reals, and some variants with which we are experimenting.

In this note, following suggestions by Tao, we extend the randomized algorithm for linear equations over prime fields by Raghavendra to a randomized algorithm for linear equations over the reals. We also show that the algorithm can be parallelized to solve a system of linear equations $A x = b$ with a regular $n \times n$ matrix $A$ in time $O(n^2)$, with probability one. Note that we do not assume that $A$ is symmetric.
Several comments here. This is the second time I see a preprint being borne out of a comment on a blog. I wonder how many there are. Second, while the randomized algorithm works for square matrices, I wonder if it still works out for short and fat matrices (underdetermined linear systems of equations) and how easy it is to have a solution from a certain family of vectors. In other words, does it still work if one chooses only vectors that are 5-sparse, that have a bernouilli distribution, I am sure y'all get my drift. 



Image Credit: NASA/JPL-Caltech
This image was taken by Navcam: Right A (NAV_RIGHT_A) onboard NASA's Mars rover Curiosity on Sol 172 (2013-01-29 21:21:57 UTC) .
Full Resolution

3 comments:

Anonymous said...

For undetermined system it should find *some* nontrivial solution.

The question about sparse solution and non-linearly constrains is quite interesting.
Algorithm based on the producing new partial solution form existing partial solutions, so for exactly solvable constraint which can be mapped to convex set it can work. But I don't see how it can be modified for sparsity constraint - combination of n-sparse solution to make new n-sparse....

Igor said...

I am sure that a graph based approach could handle this constraint issue.

Antti Lipponen said...

Hi all. I implemented the Fliege's algorithm for real valued linear systems in C++ for Matlab (based on Fliege's paper and also took some ideas from Riccardo's and Benjamin's code). However, I cannot get the performance even close to Matlab's backslash operator (probably I am using too small matrices or my code is not optimal or...). If someone is interested in my code, it can be found at: http://pastebin.com/XpSaTqPQ
All suggestions and comments are welcome!

Printfriendly