miércoles, 31 de julio de 2013

Montecarlo method for matrix inversion

One of the countless applications of stochastic concepts to real life problems is the approximate solution of large systems of algebraic equations.
After some research, I have accumulated some references on the topic and tried to implement their described algorithms in order to use them in my thesis. All the works are based on a 40's paper by Ulam and Neumann, later re-visited by Leibler and Forsythe in 1950 (Matrix inversion by a Montecarlo method). The most prominent contributions afterwards are (in my opinion) those by Halton, Dimov, Alexandrov and Mascagni, all of them available in the web.
Unfortunately, these descriptions are not really suitable for the layman and in all cases come wrapped in the (allegedly) more precise "mathematical language". This has caused me an enormous amount of trouble just to get to the conclusion that Montecarlo matrix inversion methods are not suitable for real time simulation or engineering purposes.
This assert is based on the fact that those methods are only directly applicable to diagonally dominant matrices. This means matrices whose diagonal values are bigger than the sum of all the other values in their respective rows. I made the mistake of assuming that, in FEM, as the diagonal always has the highest values, it should be equivalent, but no...

Of course, different authors have developed and adapted the method to more general cases, but I believe that then the gain in computational velocity one achieves by simply approximating is lost in all the preparatives.

As a conclusion, a much more simple and better documented algorithm to be used instead is Conjugate Gradient.