reason scipy factors only square sparse matrices Problem is basically in _superlu.gstrf <https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/sparse/linalg/_dsolve/linsolve.py#L437>it accepts only one dimension as input rather than both row and col. In c implementation <https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/sparse/linalg/_dsolve/_superlumodule.c#L202> we can see it uses N only <https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/sparse/linalg/_dsolve/_superlumodule.c#L246> and assumes A to be square matrix.
Currently I am going with the solution given by *Dima* which uses block matrices. On Wednesday, February 28, 2024 at 9:49:07 PM UTC+5:30 Animesh Shree wrote: > Yes > Because of same reason I tried to commented scipy code > <https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/sparse/linalg/_dsolve/linsolve.py#L423C1-L424C77> > > 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. > > > On Wednesday, February 28, 2024 at 8:59:36 PM UTC+5:30 Dima Pasechnik > wrote: > >> 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 >> >> >> >> >> >> >> >> >> On Wed, Feb 28, 2024 at 1:09 PM 'Animesh Shree' via sage-devel < >> sage-...@googlegroups.com> wrote: >> >>> 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 >>> <https://github.com/scikit-umfpack/scikit-umfpack/blob/ce77944bce003a29771ae07be182af348c3eadce/scikits/umfpack/interface.py#L199C1-L200C81> >>> On Wednesday, February 28, 2024 at 6:17:26 PM UTC+5:30 Max Alekseyev >>> wrote: >>> >>>> One more option would be umfack via scikits.umfpack: >>>> >>>> https://scikit-umfpack.github.io/scikit-umfpack/reference/scikits.umfpack.UmfpackLU.html >>>> >>>> Regards, >>>> Max >>>> On Wednesday, February 28, 2024 at 7:07:53 AM UTC-5 Animesh Shree wrote: >>>> >>>>> 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 >>>>> <https://github.com/sagemath/sage/blob/acbe15dcd87085d4fa177705ec01107b53a026ef/src/sage/matrix/matrix2.pyx#L13160> >>>>> (Note >>>>> - link >>>>> <https://github.com/sagemath/sage/blob/acbe15dcd87085d4fa177705ec01107b53a026ef/src/sage/matrix/matrix2.pyx#L13249C1-L13254C51> >>>>> ) >>>>> 3. convert to dense then factor >>>>> >>>>> It will be up to user to choose based on the complexity. >>>>> Is it fine? >>>>> >>>>> On Wednesday, February 28, 2024 at 4:30:51 PM UTC+5:30 Animesh Shree >>>>> wrote: >>>>> >>>>>> 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. >>>>>> >>>>>> On Tuesday, February 27, 2024 at 11:57:02 PM UTC+5:30 Dima Pasechnik >>>>>> wrote: >>>>>> >>>>>>> >>>>>>> >>>>>>> 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. >>>>>>> >>>>>>> >>>>>>> > >>>>>>> >On Tuesday, February 27, 2024 at 10:03:25 PM UTC+5:30 Dima >>>>>>> Pasechnik wrote: >>>>>>> > >>>>>>> >> >>>>>>> >> >>>>>>> >> 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. >>>>>>> >> >>>>>>> >> >>>>>>> >> >On Tuesday, February 6, 2024 at 11:29:07 PM UTC+5:30 Dima >>>>>>> Pasechnik wrote: >>>>>>> >> > >>>>>>> >> >> 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 >>>>>>> >> >> >>>>>>> >> >> >>>>>>> >> >> >>>>>>> >> >> >>>>>>> >> >> On Mon, Feb 5, 2024 at 6:00 PM Nils Bruin <nbr...@sfu.ca> >>>>>>> wrote: >>>>>>> >> >> > >>>>>>> >> >> > 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. >>>>>>> >> >> > >>>>>>> >> >> > -- >>>>>>> >> >> > 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 >>>>>>> >>>>>>> >> >> . >>>>>>> >> >> >>>>>>> >> > >>>>>>> >> >>>>>>> > >>>>>>> >>>>>> -- >>> 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/7a28f65b-a68d-4758-862e-b07d2e859d8bn%40googlegroups.com >>> >>> <https://groups.google.com/d/msgid/sage-devel/7a28f65b-a68d-4758-862e-b07d2e859d8bn%40googlegroups.com?utm_medium=email&utm_source=footer> >>> . >>> >> -- 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+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/sage-devel/d13d3977-7067-436c-a108-a3cfd13199ban%40googlegroups.com.
