On Wed, Feb 28, 2024 at 5:42 PM 'Animesh Shree' via sage-devel < sage-devel@googlegroups.com> wrote:
> 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. > 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... > > 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 > <https://groups.google.com/d/msgid/sage-devel/d13d3977-7067-436c-a108-a3cfd13199ban%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/CAAWYfq27onkM5qo9zuS_3eLNZzGgmjGhBqX850s%2B5q1uF2raDQ%40mail.gmail.com.
