SuperLU: advanced usage schemes

[sschnug]

Interest

I'm interested in SuperLUs support for:

• A: Fact='SamePattern'
• B: Fact='SamePattern SameRowPerm'
• C: ColPerm='MY_PERMC'

Additionally i'm interested in:

• D: Error Bounds, like BERR

Previous scipy-related discussions

Now i don't think any of those are supported or successfully tested by anyone although SuperLU supports of all this (at least with the expert driver). The only relevant scipy-user discussion i found is this:

I'm sure you are the first to try to use the Python wrapper for this
purpose, so the answers to your two questions here are not known.

If the matter in using SuperLU is simply to pass in this flag, then it
will probably work, but if you also need to pass in additional data
telling SuperLU about the previous permutation, then it probably won't
work out of the box.

Quick glance at the code

• _superluobject.c seems not support getters for any error-information
• gstrf has no info-return, opposed to gssv
• I don't think there is a way to pass perm_c (as indicated in above citation)

Failed test

As there (imho) is no way of passing perm_c or any other internal states, i dropped approaches based on splu and looked for the next inner-level. Futhermore i don't see any kind of usage where i can hold my factorization-object to do non-complete (from scratch) work.

I tried to address A+B with using the low-level access to _superlu.gstrf with the assumption, internally perm c will be reused (indicated by code).

This sometimes gives a segfault (not debugged more detailed) or runs for a quite some time (which is not expected!).

Edit: on a second look, my approach does not make any sense as i do not reuse any superLU-object (internal ones get probably destroyed). In this case it might be possible that this is just calling superLU in any mode expecting some given perm_c which is not available.

Code:

import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as splin
np.random.seed(0)

A = sp.random(6, 6, density=0.5)

print('Classic API')
lhs_fact = splin.splu(A)
print(lhs_fact.L[:10])
print('---')

from scipy.sparse.linalg.dsolve import _superlu
A_ = sp.csc_matrix(A)
A_.sort_indices()
A_ = A_.asfptype()  # upcast to a floating point format
N = A.shape[0]
res = _superlu.gstrf(N, A_.nnz, A_.data, A_.indices, A_.indptr,
                          ilu=False)
print('Low-level usage')
print(res.L[:10])
print('---')

print('Low-level usage with alternative ColPrerm option')
_options = dict(ColPerm='NATURAL')
res = _superlu.gstrf(N, A_.nnz, A_.data, A_.indices, A_.indptr,
                          ilu=False, options=_options)
print(res.L[:10])
print('---')

print('Low-level usage with Fact=SamePattern')
_options = dict(Fact='SamePattern')
res = _superlu.gstrf(N, A_.nnz, A_.data, A_.indices, A_.indptr,
                          ilu=False, options=_options)
print(res.L[:10])
print('---')

Output:

Classic API
C:\Miniconda3\lib\site-packages\scipy\sparse\linalg\dsolve\linsolve.py:295: SparseEfficiencyWarning: splu requires CSC matrix format
  warn('splu requires CSC matrix format', SparseEfficiencyWarning)
  (0, 0)        1.0
  (2, 0)        0.773083658044
  (4, 0)        0.18196564555
  (1, 1)        1.0
  (2, 1)        -0.793686219299
  (4, 1)        0.437604079313
  (2, 2)        1.0
  (3, 2)        0.382942859147
  (4, 2)        -0.551356529409
  (3, 3)        1.0
  (4, 3)        -0.351780778845
  (4, 4)        1.0
  (5, 5)        1.0
---
Low-level usage
  (0, 0)        1.0
  (2, 0)        0.773083658044
  (4, 0)        0.18196564555
  (1, 1)        1.0
  (2, 1)        -0.793686219299
  (4, 1)        0.437604079313
  (2, 2)        1.0
  (3, 2)        0.382942859147
  (4, 2)        -0.551356529409
  (3, 3)        1.0
  (4, 3)        -0.351780778845
  (4, 4)        1.0
  (5, 5)        1.0
---
Low-level usage with alternative ColPrerm option
  (0, 0)        1.0
  (1, 0)        0.303936470084
  (1, 1)        1.0
  (2, 1)        0.697885278513
  (2, 2)        1.0
  (3, 2)        0.773083658044
  (4, 2)        0.18196564555
  (3, 3)        1.0
  (4, 3)        0.235376396405
  (4, 4)        1.0
  (5, 5)        1.0
---
Low-level usage with Fact=SamePattern
# ...BREAKS...
# Here: Anaconda 64-bit Scipy 1.0.0 (but i don't think it matters much)

General use-case

The general idea is to save computation in some algorithmic environments by reusing previous information, which A, B, C might achieve. This is for the moment only a theoretical-idea, but i think it's worthwile to check out!

For some more aggressive reuse, one should read out D: BERR (errors) and act according to these (SuperLU manual 1.3.8).

In my case, although other things are of higher priority right now, i'm implementing an Interior-point method for Convex Quadratic Programming (with the goal of adding it to scipy).

Without detailed checking, i expect supporting the above to have consequences at least for Matt's IPM for LPs. He also asked for C in the past (issue #7700).

Support for perm_c could also be used to incorporate METIS (without some complex intertwined core) when license-usage is not a problem.