Is there a package providing Neumann series expansion for solving linear problem?
I usually solve Ax = b by calculating x = b + Sb + SSb + SSSb ... where S = 1- A, and I’d like to convert my matlab code into Julia.
So far I’m considering implementing the part by myself using IterativeSolve.jl as a last option.
If you have better idea, please give your advice.
Thank you for your suggestion in advance.
Well, that method is not always applicable, and even when applicable, could also be slowly converging. I guess that’s why it didin’t make it into a package. There are general-purpose linear solvers like GMRES you can use, as well, that are very fast. For example, GMRES uses the same information you are computing (powers of Ab), and then computes the least squares solution using that data.
That’s essentially a variant of a Jacobi iteration for D=I. IterativeSolvers.jl has a Jacobi method, but nowadays Jacobi methods and other “stationary” iterative solvers are generally considered obsolete (except in very specialized cases, e.g. as smoothers for multigrid).
You should probably use a Krylov method (e.g. CG, GMRES, BiCGSTAB, …) instead, as @martin.d.maas suggested above. (If you have a matrix A that is close enough to I for Jacobi to work well with I-A, then Krylov solvers should converge extremely quickly.)
(Of course if you really want to use this Jacobi-style iteration, e.g. for some reason it works amazingly well for your specific problem, then it is trivial to implement yourself. I suppose IterativeSolvers.jl could add an option to its jacobi function to accept an arbitrary diagonal matrix D instead of using Diagonal(A), but it’s the sort of thing that may be of limited general interest these days.)
Thank you for your suggestions! Yes, the matrix I use has a very specific form, so I think I’ve been using it without any troubles so far.
It would be interesting to compare the performance with general-purpose solvers. Thank you for the kind answers!