# Efficient solution of huge sparse linear system

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 SparseArrayKrylovLinearSolve[] and LinearSolve[..., Method -> "Krylov"] as these seem to be equivalent).

• Is A constructed as an explicit SparseArray? If not, try making it one to see if that improves on memory use and/or LinearSolve speed. Nov 15, 2015 at 22:36
• Have you seen this? Nov 16, 2015 at 3:11
• @DanielLichtblau Yes, $A$ is constructed as SparseArray object. @J.M. Yes I've seen that, I'm calling FindRoot with specified Jacobian function FindRoot[F[X], {X, X0}, Method -> {"Newton", "StepControl" -> "LineSearch", "UpdateJacobian" -> 1}, Jacobian :> J[X], AccuracyGoal -> 10, PrecisionGoal -> 10] with F and J a 'black box' functions.
– mmal
Nov 16, 2015 at 9:02
• I suppose that one solution would be to explicitly call Pardiso solver from Fortran and use LibraryLink. I'll try to experiment with that later today.
– mmal
Nov 16, 2015 at 10:46

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.
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.