Incorrect results for some problems when using cupy in CholeskySolverF

[din14970]

Thanks for publishing this project! I want to see if I can use it for solving least squares problems directly with A.T.dot(A.dot(x)) = A.T.dot(b). To see if this works I played around with some simple test with some made up values

A = array(
    [[1., 2., 3.],
     [4., 5., 6.],
     [7., 8., 9.],
     [7., 8., 9.]]
)
b = array([0.1, 0.5, -0.2, -0.1])

I then pass the details from a sparse COO representation of A.T.dot(A) to instantiate the CholeskySolverF:

solver = CholeskySolverF(
   3,
   array([0, 0, 0, 1, 1, 1, 2, 2, 2], dtype=int32),
   array([0, 1, 2, 0, 1, 2, 0, 1, 2], dtype=int32),
   array([115., 134., 153., 134., 157., 180., 153., 180., 207.]),
   MatrixType.COO,
)
x = zeros_like(b)
# b = A.T.dot(b)
solver.solve(array([-8.32667268e-17,  3.00000000e-01,  6.00000000e-01]), x)

When I use numpy for this entire procedure it seems to work fine and the test passes. The answer I get is [-0.29087737, 0.11811837, 0.11518325]. For cupy I get a clearly incorrect answer which looks a lot like some overflow issue is going on: [ 65558.29 , -131117.05 , 65558.695]. By using CholeskySolverD both numpy and cupy give me a reasonable answer.

I thought I would report this here. Perhaps if there runaway numerical errors occurring the code could give a warning/error?

[wjakob]

Your linear system is extremely ill-conditioned:

>>> np.linalg.cond(A.T @ A)
4.786320491726312e+17

In other words, this is not a reasonable input for Cholespy (in particular, the single precision variant). You can read more on conditioning of linear systems here: https://en.wikipedia.org/wiki/Condition_number#Matrices