Ok, I will do the same. On Thursday, February 29, 2024 at 12:11:56 AM UTC+5:30 Dima Pasechnik wrote:
> 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 >> <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+...@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/5ba39e26-c75c-44f4-803a-9be41b2cb9e2n%40googlegroups.com.
