Efficient solution of huge sparse linear system

Question

I'm trying to improve efficiency of my code in which main task is a solution to huge (~$10^4\times 10^4$) but sparse linear system $$Ax=b$$ (In fact my aim is to solve nonlinear equations $F(x)=0$, but the problem can be reduced with Newton's method to iterative procedure, with $A=J=F'(x)$ as Jacobian matrix of a system, for which I have an explicit form. Since I have good guess for a starting point it takes only 1-2 steps for FindRoot to converge, but it consumes great ammount of RAM ~9GB, which increases with number of steps taken). The problem is both the time it takes for Mathematica to get the solution and the ammount of memory it uses e.g.

LinearSolve[A, b, Method-> "Multifrontal"]

takes around 140s (from the beginning it works on two cores) while the memory consumption is enormous, with a peak of 6GB (where it comes from?). Thanks to the question (82645) I've found the "Pardiso" method which is better in terms of memory (at the beginning it consumes only 1GB, at the end ~3.1GB) but is much slower 415s. Interesting thing about "Pardiso" is that for most of the time (6*60s) it uses only one core, then for the next 60s it uses two (see CPU load graph below; 2 physical cores with HT) and consumes additional 2.1GB of RAM. What it does at these phases (unpacking sparse array, data transfer to MKLs Pardiso solver)?

Is there any chance to improve efficiency for this problem? Both in terms of execution time and memory consumption? Matrix $A$ is quite sparse (density $\leq 1\%$), at the end of the day I would like to increase size of the the problem (dimension of $A$) wich with current timings will make it impossible.

I think that this is independent problem of any details of the system (its sparsity pattern). Below I've uploaded files with compressed (use Uncompress[]) $A$ matrix (large file 80MB) and $b$ vector (the $32769\times 32769$ test case).

I'm suspecting that FindRoot uses "Multifrontal" method internally in this case (consistent timings, but for 2 intermediate steps I see that it consumes ~1.5 times memory as LinearSolve). Iterative methods are not an option here (as far as I know) since these take forever (both for SparseArray`KrylovLinearSolve[] and LinearSolve[..., Method -> "Krylov"] as these seem to be equivalent).

Answer 1

It's been a while since this question was asked but perhaps a sketch of an answer is useful for future reference. In version 12.0 FindRoot got a new method option for an Affine Corvariant Newton solver that is fairly efficient for large scale sets of equations. It achieves this performance by making use of a few tricks.

1. It's optimized to use minimal Jacobian evaluations.
2. It makes use of Broyden updates.
3. It stores the Jacobian decomposition that is computed with LinearSolve with the method "Pardiso".

Since your post concentrated more on the LinearSolve aspect I can not test FindRoot with your application. For LinearSolve you should make use the Method option "Pardiso" as that is quite efficient.

You can also use Iterative solves - which may or may not be efficient in your case. Since you have a good starting value, this might pan out. If you feel like it you could replace the entire 'LinearSolve` command with something geared towards your specific problem.

For more information of the use of the Affine Covariant Newton solver I'd like to refer you to the documentation