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/0ef9f470-b42a-42f8-8e2d-18338770e8e0n%40googlegroups.com.
