Opened 4 years ago

Last modified 2 weeks ago

#18365 needs_review defect

Definition of LU descomposition of a matrix depends on the base ring

Reported by: tmonteil Owned by:
Priority: major Milestone: sage-6.7
Component: linear algebra Keywords:
Cc: kcrisman, gh-ChamanAgrawal Merged in:
Authors: Chaman Agrawal Reviewers:
Report Upstream: N/A Work issues:
Branch: u/gh-ChamanAgrawal/18365_RDF_LU (Commits) Commit: 251254086d59c476e0b3db3b90c94c357063a357
Dependencies: Stopgaps:

Description

As reported on this ask question, there is an inconsistency in the definition of the LU decomposition of a matrix depending on its base ring. If a matrix A belongs to ZZ, QQ, AA, QQbar, A.LU() returns a triple (P,L,U) such that A=PLU. If a matrix A belongs to RDF, A.LU() returns a triple (P,L,U) such that PA=LU. For example:

sage: A = random_matrix(ZZ,4)
sage: A.LU()
(
[0 1 0 0]  [    1     0     0     0]  [   72    -4   -38     0]
[0 0 1 0]  [    0     1     0     0]  [    0    -1    -2    10]
[1 0 0 0]  [ 1/36   8/9     1     0]  [    0     0  17/6 -80/9]
[0 0 0 1], [-1/36   1/9    -1     1], [    0     0     0   -17]
)
sage: B = A.change_ring(RDF)
sage: B.LU()
(
[0.0 0.0 1.0 0.0]
[1.0 0.0 0.0 0.0]
[0.0 1.0 0.0 0.0]
[0.0 0.0 0.0 1.0],

[             1.0              0.0              0.0              0.0]
[             0.0              1.0              0.0              0.0]
[ 0.0277777777778   0.888888888889              1.0              0.0]
[-0.0277777777778   0.111111111111             -1.0              1.0],

[          72.0           -4.0          -38.0            0.0]
[           0.0           -1.0           -2.0           10.0]
[           0.0            0.0  2.83333333333 -8.88888888889]
[           0.0            0.0            0.0          -17.0]
)


sage: B.LU()[0] == A.LU()[0]
False
sage: B.LU()[0] == A.LU()[0].transpose()
True
sage: 

The aim of this ticket is to fix this inconsistency and choose a common definition for all rings.

Change History (5)

comment:1 Changed 4 years ago by tmonteil

It is worth noticing that the source code for matrices over RDF calls scipy.linalg.lu which returns a triple (P,L,U) such that A=PLU and then transpose P with the documentation:

# Numpy has a different convention than we had with GSL
# So we invert (transpose) the P to match our prior behavior
# TODO: It's an awful waste to store a huge matrix for P, which
# is just a simple permutation, really.

So, i guess this extra transposition should be reverted, so that it becomes consistent with both scipy and the implementation for exact rings.

comment:2 Changed 3 weeks ago by kcrisman

  • Cc kcrisman added

comment:3 Changed 3 weeks ago by kcrisman

Probably in principle the old behavior should be deprecated ... but the inconsistency is annoying, agreed.

comment:4 Changed 2 weeks ago by gh-ChamanAgrawal

  • Authors set to Chaman Agrawal
  • Branch set to u/gh-ChamanAgrawal/18365_RDF_LU
  • Cc gh-ChamanAgrawal added
  • Commit set to 251254086d59c476e0b3db3b90c94c357063a357
  • Status changed from new to needs_review

New commits:

2512540Reverted transpose to maintain consistency with scipy, numpy

comment:5 Changed 2 weeks ago by gh-ChamanAgrawal

I have removed to the extra copy of the matrix and also reverted the transposing step to maintain consistent behavior with scipy and also with LU() in matrices of other class in sagemath.

Note: See TracTickets for help on using tickets.