pardiso is random and inaccurate

[Zohar]

First issue is that it produces a random solution each time. The difference is negligible, but it affects my program behavior, which makes it harder to debug.

I thought that it's due to parallelism, but limiting the process to a single core didn't help:

SetProcessAffinityMask( GetCurrentProcess(), 1 );

Second issue is that often the solution is inaccurate. Consider the attached example in matlab format. The error of the solution, |A*x-g0|, is large (~0.1).

In comparison, matlab's ldl gives the same solution each time with error 1e-14.

The matrix is badly conditioned (condest=1.4071e+07), but a solver should be robust enough to handle it.

[Zohar]

I'm using pardiso?

[Ruqiu_C_Intel]

Thanks for posting the issue. Please share out a simple reproducer and tell us your environment, for example, oneMKL version, OS, HW, reproduce steps, etc.

[Zohar]

In my previous example, I set iparm[12]=2, which is for the other pardiso lib. Yours accepts 0 or 1.

So, here is a new linear system and a python script to solve it.

It may help other users if you give them such a script to test a system--create a reproducer.

My OS is win11, and the latest pypardiso uses mkl-2025.0.1.

import scipy.io
import numpy as np
import pypardiso

#fl = 'pardiso_small.mat'
#fl = 'pardiso_small2.mat'
#fl = 'pardiso_inaccuracy.mat'
fl = 'pardiso_inaccuracy2.mat'

mat = scipy.io.loadmat( fl )

A = mat['A']
A = scipy.sparse.csr_matrix( A )
b = mat['g0']
x0 = mat['x']

print( 'Loaded solution: |Ax-b|=', scipy.linalg.norm( A*x0 - b ), '\n' )

if 0:
    #[x, info] = scipy.sparse.linalg.cg( A, b )
    #[x, info] = scipy.sparse.linalg.bicgstab( A, b )
    x = scipy.sparse.linalg.spsolve( A, b )
    #x = pypardiso.spsolve( A, b )

    x = np.atleast_2d( x ).T # convert to column; else no error, but wrong result
else:
    solver = pypardiso.PyPardisoSolver()
    x = solver.solve( A, b )
    iparm = solver.get_iparms()
    iparm = [ int(i) for i in iparm.values() ]
    print( 'iparm =', iparm, '\n' )

print( '|Ax-b|=', scipy.linalg.norm( A*x - b ) )
#print( '|x-x0|=', scipy.linalg.norm( x - x0 ) )

#print( f'A={A}\nb={b}\nx0={x0}\nx={x}' )

(Note to self: previously, I didn't see the expansion "..." in the toolbar since the page was scaled.)

[Zohar]

I'd like to emphasize that in this example, it's not just accuracy; Pardiso fails, giving a difference of |Ax-b|=~1e+7. Standard solvers had no trouble. Each one uses default param.

I provided a convenient python script, but I get the same result using c++ Tuxfamily Eigen wrapper.

Right now, it appears that Pardiso is completely unreliable. I was hoping that someone would look into this.

[morskaya_svinka_1 · Accepted solution]

Hi. When you send a python script, you kinda force people reading this topic to search for subscript files and figure out what is written there. I think not many people have experience with these python libraries. If you send a c/cpp reader to read these csr files, I'll try to solve the systems and find out whether PARDISO is so bad with them. What mtype are these matrices, by the way?
Among other shared memory sparse solvers (I can't tell much on cluster solvers), MKL PARDISO is the best solver on memory and time so far for symmetric positive definite general sparse systems and well-conditioned symmetric indefinite ones. As for ill-conditioned symmetric matrices and unsymmetric matrices, MKL PARDISO can have troubles with accuracy, as it does not do dynamic pivoting (row interchanges in case of small pivots). Often such problems still can be solved, but you have to tune parameters (most helpful are iparm[7] and iparm[9] in my experience) for each matrix or class of matrices.

[Zohar]

The matrix file (.mat) is in matlab format.

If I implement this in c++, then I'll use Eigen, which again may be something that you don't like. Can you perhaps provide your favorite code that solves a system along with a sample file of a matrix? I'll convert my matrix file to that format and make sure it runs on your example. I can accommodate any language. Again, I'm not doing anything special; just solving a (sparse, symmetric, indefinite) system with Pardiso using default parameters.

[Zohar]

I updated iparm 7 & 9 as you suggested, and it solved my problem above. Now, I'll need to run tests again, and find another problem.

About the code. It would still be nice if you could supply something, which may help future users to create minimal examples. You could simply tell me to use pardiso_sym.c from MKL examples, but the matrix there is hardcoded.

Here is my updated script:

import scipy.io
import numpy as np
import pypardiso

#fl = 'pardiso_small.mat'
#fl = 'pardiso_small2.mat'
#fl = 'pardiso_inaccuracy.mat'
fl = 'pardiso_inaccuracy2.mat'

mat = scipy.io.loadmat( fl )

A = mat['A']
A = scipy.sparse.csr_matrix( A )
b = mat['g0']
x0 = mat['x']

print( 'Loaded solution: |Ax-b| =', scipy.linalg.norm( A*x0 - b ), '\n' )

if 0: # other solvers
    #[x, info] = scipy.sparse.linalg.cg( A, b )
    #[x, info] = scipy.sparse.linalg.bicgstab( A, b )
    x = scipy.sparse.linalg.spsolve( A, b )
    #x = pypardiso.spsolve( A, b )

    # convert
    x = np.atleast_2d( x ).T # convert to column; else no error, but wrong result
    
else: # pardiso
    solver = pypardiso.PyPardisoSolver()
    
    # mtype
    solver.set_matrix_type = -2 # symmetric indefinite
    
    # iparm
    if 0:
        # from mkl example pardiso_sym.c with my tweaks for 7 and 9
        ip = [ 1, 2, 0, 0, 0, 0, 0, 5, 0, 5, 1, 0, 0, 0, 0, 0, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ]
        
        for i, v in enumerate( ip ):
            solver.set_iparm( i +1, v ) # 1-based indexing
            
    if 1:
        solver.set_iparm( 0 +1, 1 )
        solver.set_iparm( 1 +1, 2 )
        
        # https://community.intel.com/t5/Intel-oneAPI-Math-Kernel-Library/pardiso-is-random-and-inaccurate/m-p/1658985
        solver.set_iparm( 7 +1, 3 )
        solver.set_iparm( 9 +1, 7 )
        
        solver.set_iparm( 10 +1, 1 )
        solver.set_iparm( 17 +1, -1 )
        solver.set_iparm( 18 +1, -1 )
        
    iparm = solver.get_iparms()
    iparm = [ int( i ) for i in iparm.values() ]
    print( 'Set iparm =', iparm, '\n' )

    # solve
    x = solver.solve( A, b )
    
    # iparm after
    iparm = solver.get_iparms()
    iparm = [ int( i ) for i in iparm.values() ]
    print( '(after) iparm =', iparm, '\n' )

print( '|Ax-b| =', scipy.linalg.norm( A*x - b ) )
#print( '|x-x0| =', scipy.linalg.norm( x - x0 ) )

#print( f'A={A}\nb={b}\nx0={x0}\nx={x}' )

[Zohar]

It seems fine so far. I'll keep testing it.

[morskaya_svinka_1]

Here is working code for int32 spd systems. I caught myself I don't actually have a good int32 IC single-file launcher, so it can contain some deprecated comments and unnecessary actions. Still, I think it explains file format that I like (4 separate binary files and 1 text one). I can't really tell what iparms are good for indefinite systems as for different systems I found different sets of iparms. I used such strategy: set some fixed steps of refinement, say iparm[7]=-5, and then run through iparm[9], like 5, 7, 9, ..., 13, 15. After choosing most reliable iparm[9], I proceeded to set iparm[7] and then other iparms. But oher ones like iparm[10] and iparm[12] did not really affect the accuracy, as I guess I had zero-free diagonals mostly. I did not meet a matrix from practical application that I would not be able to solve at least with 1e-3 relative residual using MKL PARDISO, but I think many of such examples exist and I had to spend plenty of time adjusting iparms for indefinite and unsymmetric matrices. This is the reason I prefer other solvers for solving indefinite or unsymmetric problems (MKL PARDISO cannot be automatized to provide reliable results in general case on the contrary with solvers that can do dynamic pivoting). Though it is fine when you need to factor only a few classes of similar matrices, so you have to adjust iparms just once.
If anyone knows another strategies and another experience I am glad to hear about it.

[Zohar]

1. I slightly changed your code and attached an example (if someone is interested).

It appears that the wrappers aren't flexible enough to get peak performance/accuracy. For the attached example, pypardiso and Eigen/PardisoSupport produced at best a residual of 0.6. Going low level with your code, I got 1e-6. The right mtype and phases are important...

2. "I prefer other solvers for solving indefinite or unsymmetric problems"
Which ones?
Normally, I just use matlab's mldivide (and suffer the overhead of c++ to matlab and back), which uses ldl (MA57).
I tried suitesparse, but there's nothing specialized for this type of problems.

By the way, have you tried the solvers on this list:

https://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html#title3

If they made it to eigen, they might be good.

3. "Pardiso can't be automatized as opposed to solvers that can do dynamic pivoting."
Pardiso has dynamic pivoting:

http://www.icm.tu-bs.de/~bolle/Publicat/siamreview.pdf

[morskaya_svinka_1]

1. I was able to solve the system with MKL PARDISO (iparm [1]=3, [7]=-2, other defaults) relative residual 6.38e-13, which is a decent result against other solvers umfpack 8.06e-13, SuperLU 5.74e-11, mumps 9.385e-15, pastix 3.105e-10. So this matrix is not the case to complain about MKL PARDISO incaccuracy.
2. I don't know if this is appropriate place to discuss other solvers, you may drop me your e-mail.
3. Nope, what made you think that? I remember there was some article that was somewhat related to some version of PARDISO solver, and different pivoting techniques were discussed including dynamic ones, but it probably remained on paper. As far as I understood from PARDISO solver manual (https://panua.ch/pardiso/), it can do more advanced pivoting techniques compared to MKL PARDISO, but still does not support dynamic pivoting.

[Zohar]

1. I stopped complaining about accuracy once I realized the wrappers aren't good enough. My example isn't a problem, I was fine with 1e-6, and your improvement is even better.

2. I sent you a private message.

3. Right, I misread the paragraph on my link:

"With the introduction of symmetric maximum weighted matchings [24] as an alternative to complete pivoting... This approach has recently been successfully implemented in the sparse direct solver Pardiso [60, 59]."

They don't use pivoting but an alternative.

[Zohar]

1. Just noting (for future reference) that you also used iparm[12]=1.

nnz of iparm: 0:1, 1:3, 7:-2, 12:1, 17:-1, 34:1,

[Zohar]

Another tip. Prune the matrix, removing all exact zeros.

I had a case where the residual was 1e9 because of that. It took me 3 hours to get to the bottom of it: I had also a pruned (apparently) matrix from matlab, which worked fine, and when I compared both using norm, they were the same... Saving and loading triplets in Eigen didn't make a difference since it doesn't prune.

Pruning compromises reusing an analyzed pattern, where an indication is a large residual.

[Zohar]

I'd like to share my current conclusion.

I'm solving an optimization problem using Newton method. It's an iterative process, where in each iteration I'm solve a linear system for the search direction (involves factoring a constrained Hessian:  https://community.intel.com/t5/Intel-oneAPI-Math-Kernel-Library/Pardiso-Schur-complement-for-an-SPD-matrix/m-p/1661844 ).

Pardiso is fast but not fully reliable. I couldn't find a golden set of parameters that will work for all iterations. Trying to automate fidgeting them while running isn't bulletproof.

What I simply do is use matlab (ldl, MA57) for a second solve in case of failure. This method is slower since pattern analysis isn't reused, but it's reliable. This is important since it's not practical to stop in the middle of the iterations and mess with a specific system.

About the earlier pruning suggestion. It requires re-analyzing the pattern, which is costly and disables pardiso's main advantage.

(My chosen iparm: 0:1, 1:3, 7:-5, 9:5, 12:1, 17:-1, 18:-1, 34:1.)