SuiteSparse and sage and sparse_matrix.LU()

[Animesh Shree]

I have been looking into sparse matrix LU decomposition issue. 

In that conversation SuiteSparse is proposed to handle sparse matrix decomposition.

What I could find that it is mostly C and Matlab. I don't know where it has been used before in sage. There is doc present on site.

[Dima Pasechnik]

I've just left a comment there - copying here just in case:

suitesparse is best reached via a Sage package cvxopt.
https://cvxopt.org/userguide/spsolvers.html#cvxopt.umfpack.numeric
says you can get LU factorisation of sparse matrices there.

Dima

[Animesh Shree]

Thank You I am looking into it

[Animesh Shree]

I tried using that but I am unable to access the factors. I asked them(CVXOPT) about this.

For now I am falling back to Scipy implementation. But it (scipy.sparse.linalg.splu) handles square sparse matrix only.

[Animesh Shree]

It returns PyCapsule object, I am looking if we can decapsulate it in python/cython. If it doesn't work then I will go for Scipy

On Friday, January 19, 2024 at 4:04:37 PM UTC+5:30 Animesh Shree wrote:

[Nils Bruin]

By the looks of it, the routines you'd be using would be coming from umfpack. cvxopt has chosen to package the details of LU factorization in opaque objects and instead offer routines to use these decompositions (via taking the opaque object as input). In scipy, umfpack is also used: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.spsolve.html#scipy.sparse.linalg.spsolve . Scipy also offers LU decomposition routines: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html#scipy.sparse.linalg.splu but that uses a different library. It looks suspicious that two packages offer umfpack for solving sparse problems, but don't give explicit access to LU factorizations produced in the process, when using UMFPACK. Perhaps UMFPACK isn't suited to provide the explicit factorization (but it may be very good at using its internally computed data).

[Dima Pasechnik]

IMHO opaque C objects may be accessed from Python, it just needs an implementation

[Animesh Shree]

I tried to use ctypes pythonapi to extract PyCapsule_New objects using PyCapsule_GetPointer. But got erorr :

              ArgumentError: argument 1: TypeError: Don't know how to convert parameter 1

Attaching file

---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Also, please tell me about the implementations for accessing opaque C objects I will go through it

[Animesh Shree]

The other library that scipy uses is SuperLU : https://portal.nersc.gov/project/sparse/superlu/

for the function scipy.sparse.linalg.splu : https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html#scipy-sparse-linalg-splu

Scipy supports only factorization for square matrices

the developer asks whether its a limitation or not : https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/sparse/linalg/_dsolve/linsolve.py#L424

In reference manual of SuperLU : https://portal.nersc.gov/project/sparse/superlu/ug.pdf

In Section 2.4, page 23 we can see code where "nrow = 5; ncol = 5; "  are defined separately. (  The example intended to solved a 5x5 system )

[Dima Pasechnik]

On 5 February 2024 05:36:45 GMT, 'Animesh Shree' via sage-devel <sage-...@googlegroups.com> wrote:
>I tried to use ctypes pythonapi to extract PyCapsule_New objects using
>PyCapsule_GetPointer. But got erorr :

> *ArgumentError: argument 1: TypeError: Don't know how to
>convert parameter 1*

>Attaching file
>
>---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>
>Also, please tell me about the implementations for accessing opaque C
>objects I will go through it

check out the source code of cvxopt.

[Dima Pasechnik]

On 5 February 2024 06:30:45 GMT, 'Animesh Shree' via sage-devel <sage-...@googlegroups.com> wrote:
>The other library that scipy uses is SuperLU :
>https://portal.nersc.gov/project/sparse/superlu/
>for the function scipy.sparse.linalg.splu :
>https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.splu.html#scipy-sparse-linalg-splu
>
>Scipy supports only factorization for square matrices
>the developer asks whether its a limitation or not :
>https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/sparse/linalg/_dsolve/linsolve.py#L424
>
>In reference manual of SuperLU :
>https://portal.nersc.gov/project/sparse/superlu/ug.pdf

>In Section 2.4, page 23 we can see code where *"nrow = 5; ncol = 5; *" are

>defined separately. ( The example intended to solved a 5x5 system )

it is the matter of adding extra zero rows or columns to the matrix you want to decompose. This could be a quick fix.

A good implementation of LU decomposition ought actually to take non-square matrix as input, and have the indices adjusted appropriately in the algorithm,
so it's indeed a bit strange that superLU only takes square matrices (?). Perhaps it's a good idea to look at its docs and the source code.

[Thierry Dumont]

Le 05/02/2024 à 11:30, Dima Pasechnik a écrit :
>
>
>
>
> A good implementation of LU decomposition ought actually to take non-square matrix as input, and have the indices adjusted appropriately in the algorithm,
> so it's indeed a bit strange that superLU only takes square matrices (?). Perhaps it's a good idea to look at its docs and the source code.
>
>
>

Actually, SuperLU whas developed for people (like me) who solve
discretized Partial Differential Equations. In that case, matrices are
quite always square matrices. It is an old program. The idea whas that
given the fact that new non zero terms appear during the
LU-factorization, but that we can predict where they appear and
minimize the amount of these news terms, there exists a class of rather
large problems (but not too large! -say less than 10^6 unknowns-) for
which a direct method, efficiently coded, could be faster than an
iterative one (and this is true).

t.d.

[Nils Bruin]

On Monday 5 February 2024 at 02:31:04 UTC-8 Dima Pasechnik wrote:

it is the matter of adding extra zero rows or columns to the matrix you want to decompose. This could be a quick fix.

(in reference to computing LU decompositions of non-square matrices) -- in a numerical setting, adding extra zero rows/columns may not be such an attractive option: if previously you know you had a maximal rank matrix, you have now ruined it by the padding. It's worth checking the documentation and literature if padding is appropriate/desirable for the target algorithm/implementation.

[Dima Pasechnik]

Non-square case for LU is in fact easy. Note that if you have A=LU as
a block matrix
A11 A12
A21 A22

then its LU-factors L and U are
L11 0 and U11 U12
L21 L22 0 U22

and A11=L11 U11, A12=L11 U12, A21=L21 U11, A22=L21 U12+L22 U22

Assume that A11 is square and full rank (else one may apply
permutations of rows and columns in the usual way). while A21=0 and
A22=0. Then one can take L21=0, L22=U22=0, while A12=L11 U12
implies U12=L11^-1 A12.
That is, we can first compute LU-decomposition of a square matrix A11,
and then compute U12 from it and A.

Similarly, if instead A12=0 and A22=0, then we can take U12=0,
L22=U22=0, and A21=L21 U11,
i.e. L21=A21 U11^-1, and again we compute LU-decomposition of A11, and
then L21=A21 U11^-1.

----------------

Note that in some cases one cannot get LU, but instead must go for an
PLU,with P a permutation matrix.
For non-square matrices this seems a bit more complicated, but, well,
still doable.

HTH
Dima

> --
> You received this message because you are subscribed to the Google Groups "sage-devel" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sage-devel+...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/622a01e0-9197-40c5-beda-92729c4e4a32n%40googlegroups.com.

[Animesh Shree]

I tried scipy which uses superLU. We get the result but there is little bit of issue.

--For Dense--

The dense matrix factorization gives this output using permutation matrix

sage: a = Matrix(RDF, [[1, 0],[2, 1]], sparse=True)
sage: a
[1.0 0.0]
[2.0 1.0]
sage: p,l,u = a.dense_matrix().LU()
sage: p
[0.0 1.0]
[1.0 0.0]
sage: l
[1.0 0.0]
[0.5 1.0]
sage: u
[ 2.0  1.0]
[ 0.0 -0.5]

--For Sparse--

But the scipy LU decomposition uses permutations which involves taking transpose, also the output permutations are represented as array.

sage: p,l,u = a.LU(force=True)
sage: p
{'perm_r': [1, 0], 'perm_c': [1, 0]}
sage: l
[1.0 0.0]
[0.0 1.0]
sage: u
[1.0 2.0]
[0.0 1.0]

Shall I continue with this?

[Dima Pasechnik]

On 27 February 2024 15:34:20 GMT, 'Animesh Shree' via sage-devel <sage-...@googlegroups.com> wrote:
>I tried scipy which uses superLU. We get the result but there is little bit
>of issue.
>
>
>--For Dense--
>The dense matrix factorization gives this output using permutation matrix
>sage: a = Matrix(RDF, [[1, 0],[2, 1]], sparse=True)
>sage: a
>[1.0 0.0]
>[2.0 1.0]
>sage: p,l,u = a.dense_matrix().LU()
>sage: p
>[0.0 1.0]
>[1.0 0.0]
>sage: l
>[1.0 0.0]
>[0.5 1.0]
>sage: u
>[ 2.0 1.0]
>[ 0.0 -0.5]
>

you'd probably want to convert the permutation matrix into a permutation.

>--For Sparse--
>But the scipy LU decomposition uses permutations which involves taking
>transpose, also the output permutations are represented as array.

It is very normal to represent permutations as arrays.
One can reconstruct the permutation matrix from such an array trivially (IIRC, Sage even has a function for it)

I am not sure what you mean by "taking transpose".

>sage: p,l,u = a.LU(force=True)
>sage: p
>{'perm_r': [1, 0], 'perm_c': [1, 0]}
>sage: l
>[1.0 0.0]
>[0.0 1.0]
>sage: u
>[1.0 2.0]
>[0.0 1.0]
>
>
>Shall I continue with this?

sure, you are quite close to getting it all done it seems.

[Animesh Shree]

This works good if input is square and I also checked on your idea of padding zeros for non square matrices.

I am currently concerned about the permutation matrix and L, U in case of padded 0s. Because if we pad then how will they affect the outputs, so that we can extract p,l,u for unpadded matrix.

[Animesh Shree]

For transpose :

In  the example we can see permutations are provided as arrays for rows and cols.

The permutation is equivalent of taking transpose of matrix.

But we cant represent transpose as a permutation matrix.

>>> a = np.matrix([[1,2],[3,5]])
>>> # a * perm = a.T            
>>> # perm = a.I * a.T
>>> a.I*a.T
matrix([[-1., -5.],
        [ 1.,  4.]])
>>>

the output is not permutation matrix. 

[Animesh Shree]

Sorry for multiple messages

I just want to say

>sage: p,l,u = a.LU(force=True)
>sage: p
>{'perm_r': [1, 0], 'perm_c': [1, 0]}

It  (  {'perm_r': [1, 0], 'perm_c': [1, 0]}  )  represents transpose and it cannot be represented as permutation matrix.

Similar cases may arise for other matrices.

[Nils Bruin]

See the documentation:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.SuperLU.html

As you can see there it computes permutation matrices Pr and Pc and lower and upper triangular matrices L,U such that

Pr*A*Pc = L*U

so it's not computing a normal LU decomposition of a PLU decomposition: it's doing both row- and column permutations on A; presumably to condition it to get  a more numerically stable L,U.

It's a different decomposition. You may be able to solve similar problems with it, but you'll have to expose the computed data to the user. I don't think it's suitable for wrapping as a normal LU-decomposition, as it won't be a drop-in substitute for it.

You do get

A = Pr^(-1) *L*U * Pc^(-1)

[Dima Pasechnik]

On 27 February 2024 17:25:51 GMT, 'Animesh Shree' via sage-devel <sage-...@googlegroups.com> wrote:
>This works good if input is square and I also checked on your idea of
>padding zeros for non square matrices.
>I am currently concerned about the permutation matrix and L, U in case of
>padded 0s. Because if we pad then how will they affect the outputs, so that
>we can extract p,l,u for unpadded matrix.

please read details I wrote on how to deal with the non-square case. There is no padding needed.

[Animesh Shree]

Thank you for reminding

I went through. 

We need to Decompose  A11 only and rest can be calculated via taking inverse of L11 or U11.

Here A11 is square matrix and we can use scipy to decompose square matrices.

Am I correct?

New and only problem that I see is the returned LU decomposition of scipy's splu is calculated by full permutation of row and column as pointed out by Nils Bruin. We will be returning row and col permutation array/matrix separately instead of single row permutation which sage usage generally for plu decomposition.

User will have to manage row and col permutations. 

or else

We can return handler function for reconstruction of matrix from  L, U & p={perm_r, perm_c}

or

We can leave that to user

User will have to permute its input data according to perm_c (like : perm_c * input) before using the perm_r^(-1) * L * U

as perm_r^(-1) * L * U is PLU decomposition of Original_matrix*perm_c

https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.SuperLU.html

>>> A = Pr^(-1) *L*U * Pc^(-1) # as told by Nils Bruin

or

scipy's splu will not do.

[Animesh Shree]

One thing I would like to suggest.

We can provide multiple ways to compute the sparse LU

1. scipy 

2. sage original implementation in src.sage.matrix.matrix2.LU (Note - link)

3. convert to dense then factor

It will be up to user to choose based on the complexity.

Is it fine?

[Max Alekseyev]

One more option would be umfack via scikits.umfpack:

https://scikit-umfpack.github.io/scikit-umfpack/reference/scikits.umfpack.UmfpackLU.html

Regards,

Max

[Animesh Shree]

I went through the link.

It also returns perm_c and perm_r and the solution is represented as

Pr * (R^-1) * A * Pc = L * U

It is similar to one returned by scipy

>>> lu.perm_r

       array([0, 2, 1, 3], dtype=int32)

>>> lu.perm_c

       array([2, 0, 1, 3], dtype=int32)

I think it doesn't support square matrix too. Link

[Dima Pasechnik]

There is a good reason for numerics people to adopt "SuperLU" factorisations over

the classical PLU sparse decomposition - namely, it's more stable.

Perhaps it should be made the Sage's default for sparse RDF matrices, too.

By the way, https://portal.nersc.gov/project/sparse/superlu/superlu_ug.pdf

(the manual for the upstream superlu) says it can handle non-square matrices.

Dima

To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/7a28f65b-a68d-4758-862e-b07d2e859d8bn%40googlegroups.com.

[Animesh Shree]

Yes

Because of same reason I tried to commented scipy code to test this.

I got some error saying RuntimeError: Factor is exactly singular

But same worked for sage LU factorization in dense matrix for same matrix.

-----Scipy (Modified)------

>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import splu
>>> import numpy as np
>>> A = csc_matrix([[1,0,0],[5,0,2]], dtype=np.float64)
>>> print(A)
  (0, 0) 1.0
  (1, 0) 5.0
  (1, 2) 2.0
>>> splu(A)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/shay/miniconda3/envs/py39/lib/python3.9/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py", line 440, in splu
    return _superlu.gstrf(N, A.nnz, A.data, indices, indptr,
RuntimeError: Factor is exactly singular
>>>

-----Sage-------

sage: A = Matrix(RDF, 2,3, [[1,0,0],[5,0,2]])
sage: A
[1.0 0.0 0.0]
[5.0 0.0 2.0]
sage: A.LU()
(
[0.0 1.0]  [1.0 0.0]  [ 5.0  0.0  2.0]
[1.0 0.0], [0.2 1.0], [ 0.0  0.0 -0.4]
)

I am looking into it too.

[Nils Bruin]

As Dima points out, SPLU decomposition is probably more useful in practice, so it's great to expose it. It should probably be exposed through a different method than PLU decomposition, though.

For the return format: permutations are much more compact and efficient to describe permutations than permutation matrices. Also because of the parent: a permutation matrix is a {0,1}-valued matrix so it can be well represented over ZZ, rather than floats. So primarily the routine should expose that scipy exposes. Whether other bells and whistles are required is another matter.

Is this only for system floats/complex numbers? or will you also be implementing multiprecision versions?

[Animesh Shree]

I am currently working for RDF and CDF mainly. I was planning to overload LU function of sage/matrix/matrix2.pyx  by defining new LU in sage/matrix/matrix_double_sparce.pyx .

I couldn't understand what is  multiprecision version but I guess you are asking for RR or other inexact field.

Please correct me if I am wrong.

[Animesh Shree]

reason scipy factors only square sparse matrices

Problem is basically in _superlu.gstrf it accepts only one dimension as input rather than both row and col.

In c implementation we can see it uses N only and assumes A to be square matrix.

Currently I am going with the solution given by Dima which uses block matrices.

[Dima Pasechnik]

On Wed, Feb 28, 2024 at 5:42 PM 'Animesh Shree' via sage-devel <sage-...@googlegroups.com> wrote:

reason scipy factors only square sparse matrices

Problem is basically in _superlu.gstrf it accepts only one dimension as input rather than both row and col.

In c implementation we can see it uses N only and assumes A to be square matrix.

We can probably change scipy's code (and then do a PR) to accept a non-square matrix.

The C code they have basically sets up a call to superlu, and gets the results...

 

To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/d13d3977-7067-436c-a108-a3cfd13199ban%40googlegroups.com.

[Animesh Shree]

Ok, I will do the same.