Sequential monte carlo techniques for the solution of linear systems

Given alinear system Ax=b, wherex is anm-vector,direct numerical methods, such as Gaussian elimination, take timeO(m³) to findx. Iterative numerical methods, such as the Gauss-Seidel method or SOR, reduce the system to the formx=a+Hx, whence and then apply the iterationsx₀=a,x_s+1=a+Hx_s, until sufficient accuracy is achieved; this takes timeO(m²) per iteration. They generate the truncated sumsThe usualplain Monte Carlo approach uses independent random walks, to give an approximation to the truncated sumx_s, 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].