search

haasad / PyPardiso

Python interface to the Intel MKL Pardiso library to solve large sparse linear systems of equations
BSD 3-Clause "New" or "Revised" License
133 stars 20 forks source link

Wrong results for ill-conditioned CSC-matrices #7

Closed haasad closed 7 years ago

haasad commented 7 years ago

The Pardiso library uses the CSR-matrix format internally. It is possible to pass a CSC-matrix to the solver and solve the transposed system to get the same result. This is the way it is handled in the Julia-Pardiso-Wrapper, because Julia only knows the CSC-format. The Umfack library uses the same relation the other way around, because it uses the CSC format internally.

This approach works fine, as long as matrix A is well-conditioned. If matrix A is singular, the solver will find different solutions, depending on the input format of A (csr or csc). It is also possible that no solution is found. Details in this notebook.

PyPardiso is intended to be used with brightway2, where the technosphere matrix is generally rather ill-conditioned. PyPardiso performs fine if the technosphere matrix is passed in CSR-format, but fails to find a valid solution for CSC-matrices and produces wrong results (see notebook).

Brightway uses the CSR-format by default, where this issue doesn't appear. There is a case however (ie. lca.lci(factorize=True)), where the technosphere matrix is passed in CSC-format to the solver. LCIA results obtained in this way are likely to be wrong. This is a very serious bug.

The same issue (but less extreme) can be observed with scipy's SuperLU library (uses CSC internally), where the CSR-format produces wrong results (even though computation time is 40x longer for CSC). Comparison of SuperLU, Umfpack and PyPardiso can be found in this notebook. I suspect that some pattern (eg. almost identical pairs of columns because of RoW datasets), makes the transposed system more error-prone.

The obvious immediate fix for PyPardiso is to always transform the matrix to CSR instead of solving the transposed system in case of CSC.

haasad commented 7 years ago

This issue has been fixed back in November with 7a60c843a1caa2f67a61629ffeeb9b6de47dedb3. But it's still unclear why this happens for CSC-matrices.