DirectSolver yields incorrect result for random and close to random matrices

[TimoEckstein]

I was considering swapping cupy sparse solvers for nvmath Python/cuDSS DirectSolver, but I am unable to obtain correct results even for small problem instances.

Minimal example

In example01_cupy.py replace

n = 8
...
a += sp.diags([2.0] * n, format="csr", dtype="float64")

with

n = 100
...
# a += sp.diags([2.0] * n, format="csr", dtype="float64")

The same remains true, if I used small, but not tiny values for the diagonal like 0.01.

Example observation 1

||A||: 	40.65133430135082
||b||: 	14.142135623730951

det(A): 7.262125991606542e+28

02-10 17:35:07 userlogger   INFO     = SPECIFICATION PHASE =
02-10 17:35:07 userlogger   INFO     The LHS package is cupyx.
02-10 17:35:07 userlogger   INFO     The RHS package is cupy.
02-10 17:35:07 userlogger   INFO     The device_id=0, dtype = float64, index type = int32.
02-10 17:35:07 userlogger   INFO     The number of equations = 100.
02-10 17:35:07 userlogger   INFO     The operands' memory space is cuda, and the execution space is on device 0.
02-10 17:35:07 userlogger   INFO     The specified stream for the DirectSolver ctor is <cuda.core.experimental._stream.Stream object at 0x14b6709b87c0>.
02-10 17:35:07 userlogger   INFO     The library handle has been created: 94793135484384.
02-10 17:35:07 userlogger   INFO     The sparse direct solver operation has been created.
02-10 17:35:07 userlogger   INFO     Starting solver phase ANALYSIS...
02-10 17:35:07 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:35:07 userlogger   INFO     Starting solver phase FACTORIZATION...
02-10 17:35:07 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:35:07 userlogger   INFO     Starting solver phase SOLVE...
02-10 17:35:07 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:35:07 userlogger   INFO     The DirectSolver object's resources have been released.

 CuDSS:
||x||: 		29391823863.993862
||Ax - b||: 	92575181911.9934

CuPy:
||x||: 		10.647411227999223
||Ax - b||: 	4.289313285599112e-14

NumPy:
||x||: 		10.647411227999223
||Ax - b||: 	4.29904389535207e-14

Note that I computed the determinant to check whether A is invertible.

Example observation 2

Vice versa I made the observation that DirectSolver sometimes returns not nan even though the determinant is 0. Note that here b = cp.ones((n, 1), order="F") was only a vector

||A||: 	5.6922584694712395
||b||: 	10.0

det(A): 0.0

02-10 17:30:48 userlogger   INFO     = SPECIFICATION PHASE =
02-10 17:30:48 userlogger   INFO     The LHS package is cupyx.
02-10 17:30:48 userlogger   INFO     The RHS package is cupy.
02-10 17:30:48 userlogger   INFO     The device_id=0, dtype = float64, index type = int32.
02-10 17:30:48 userlogger   INFO     The number of equations = 100.
02-10 17:30:48 userlogger   INFO     The operands' memory space is cuda, and the execution space is on device 0.
02-10 17:30:48 userlogger   INFO     The specified stream for the DirectSolver ctor is <cuda.core.experimental._stream.Stream object at 0x14c6b6504580>.
02-10 17:30:48 userlogger   INFO     The library handle has been created: 94064601263200.
02-10 17:30:48 userlogger   INFO     The sparse direct solver operation has been created.
02-10 17:30:48 userlogger   INFO     Starting solver phase ANALYSIS...
02-10 17:30:48 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:30:48 userlogger   INFO     Starting solver phase FACTORIZATION...
02-10 17:30:48 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:30:48 userlogger   INFO     Starting solver phase SOLVE...
02-10 17:30:48 userlogger   INFO     This call is non-blocking and will return immediately after the operation is launched on the device.
02-10 17:30:48 userlogger   INFO     The DirectSolver object's resources have been released.

 CuDSS:
||x||: 		1.5379474527090538e+102
||Ax - b||: 	1.5379474527090537e+89

CuPy:
||x||: 		nan
||Ax - b||: 	nan

NumPy:
||x||: 		nan
||Ax - b||: 	nan

[yangcal]

Hello @TimoEckstein ,

Thanks for sharing your observation.

RE: incorrect results for invertible A

Can you try either setting reordering_algorithm to nvmath.sparse.advanced.DirectSolverAlgType.ALG_1, or use_matching to True (not both), e,g

algo = nvmath.sparse.advanced.DirectSolverAlgType.ALG_1
with Solver(a, b) as solver:
    solver.plan_config.reordering_algorithm = algo
    # alternatively, solver.plan_config.use_matching = True
    plan_info = solver.plan()
    fac_info = solver.factorize()
    x = solver.solve()
    print(f"{fac_info.diag=}")
    print(f"{abs(a @ x - b).max()=}")

On my side, when the cudss ALG_DEFAULT was on, factorization seems to encounter numerical instability with diagonals ranging over orders of magnitude. When ALG_1 was used, I was able to get converged results with reasonable diagonals for my calculation. Meanwhile we will check with internal dev team to get advice on algorithm selection for general cases.

RE: Not returning nan values

We will report to the cuDSS dev team.

[yangcal]

To follow up on my previous comment, here are some words of wisdom shared from cuDSS team:

use_matching

1. Matching often reduces the number of perturbed pivots and improves accuracy of the solution (especially for ill-conditioned and badly scaled matrices). Nevertheless, this is not guaranteed.
2. For symmetric positive definite (SPD) matrices, matching is unnecessary. For symmetric indefinite matrices, the benefits of matching is unknown.

algorithm

1. PlanConfig allow setting reordering/factorization/matching algorithms. These algorithms are all independent - and often have incompatibilities, e.g, matching is not supported for CUDSS_ALG_1 and CUDSS_ALG_2 reordering algorithms.
2. ALG_1  and ALG_2 forces global pivoting factorization, so it is the most accurate factorization, while performance is inferior (that is why it is not default)

[yangcal]

RE: nan value returns, sharing insights from internal dev team:

LU factorization (either dense or sparse) generally is not a rank-revealing decomposition.  This means that it might not necessarily be able to detect the case where a matrix is singular (and thus has zero determinant).  What it does in that case is dependent on the algorithm and the problem.  Furthermore, a matrix's determinant might be zero in finite-precision arithmetic, yet not be numerically singular.  An example is a 400 x 400 diagonal matrix whose diagonal entries are all 0.1.