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-devel@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+unsubscr...@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/CAAWYfq2JFwCDZYmFXb1V7pwAX%2B79FcBkxCeYC7HcG-X7UQOx7A%40mail.gmail.com.
