The impedance matrix localization (IML) method replaces the usual full N¿N moment method matrix by a sparse N¿N matrix. Typically, this matrix has about 50N nonzero elements. Iterative methods require one (or a few) matrix vector product(s) at each iteration. Each product can be accomplished in 50N operations rather than N2 for the usual moment method. There is a correspondence between the clumps of large numbers in the IML matrix and terms with a ray description such as in geometrical optics and the geometrical theory of diffraction. It is known that using these methods approximate solutions are possible while taking into account only a limited number of multiple interactions. Using this analogy an approximate factored form of the IML matrix is generated by explicitly using only a limited number of direct interactions. However, all multiple interactions involving these direct interactions are (implicitly) taken into account. This effectively gives an approximate inverse which is used in an iterative approach. A reduction in the residual of orders of magnitude per iteration is achieved. The result is a more accurate solution than was possible before and which requires only a few iterations. This method is generally applicable since the factored form is based on the structure of the IML matrix and not on any specific geometry. ¿ American Geophysical Union 1995