Sequential monte carlo techniques for the solution of linear systems

AbstractGiven alinear system Ax=b, wherex is anm-vector,direct numerical methods, such as Gaussian elimination, take timeO(m3) to findx. Iterative numerical methods, such as the Gauss-Seidel method or SOR, reduce the system to the formx=a+Hx, whence $$x = \sum\nolimits_{r = 0}^\infty {H^r a} ;$$ and then apply the iterationsx0=a,xs+1=a+Hxs, until sufficient accuracy is achieved; this takes timeO(m2) per iteration. They generate the truncated sums $$x_s = \sum\nolimits_{r = 0}^\infty {H^r a} .$$ The usualplain Monte Carlo approach uses independent “random walks”, to give an approximation to the truncated sumxs, taking timeO(m) per random step. Unfortunately, millions of random steps are typically needed to achieve reasonable accuracy (say, 1% r.m.s. error). Nevertheless, this is what has had to be done, ifm is itself of the order of a million or more.The alternative presented here, is to apply a sequential Monte Carlo method, in which the sampling scheme is iteratively improved. Simply put, ifx=y+z, wherey is a current estimate ofx, then its correction,z, satisfiesz=d+Hz, whered=a+Hy−y. At each stage, one uses plain Monte Carlo to estimatez, and so, the new estimatey. If the sequential computation ofd is itself approximated, numerically or stochastically, then the expected time for this process to reach a given accuracy is againO(m) per random step; but the number of steps is dramatically reduced [improvement factors of about 5,000, 26,000, 550, and 1,500 have been obtained in preliminary tests].