Opened 7 years ago
Closed 3 years ago
#18365 closed defect (fixed)
Definition of LU decomposition of a matrix depends on the base ring
Reported by:  Thierry Monteil  Owned by:  

Priority:  major  Milestone:  sage9.1 
Component:  linear algebra  Keywords:  
Cc:  KarlDieter Crisman, Chaman Agrawal  Merged in:  
Authors:  Chaman Agrawal  Reviewers:  Markus Wageringel 
Report Upstream:  N/A  Work issues:  
Branch:  96a9f92 (Commits, GitHub, GitLab)  Commit:  96a9f92abda84d6dd863c8604f23081425043110 
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 (17)
comment:1 Changed 7 years ago by
comment:2 Changed 4 years ago by
Cc:  KarlDieter Crisman added 

comment:3 followup: 6 Changed 4 years ago by
Probably in principle the old behavior should be deprecated ... but the inconsistency is annoying, agreed.
comment:4 Changed 4 years ago by
Authors:  → Chaman Agrawal 

Branch:  → u/ghChamanAgrawal/18365_RDF_LU 
Cc:  Chaman Agrawal added 
Commit:  → 251254086d59c476e0b3db3b90c94c357063a357 
Status:  new → needs_review 
New commits:
2512540  Reverted transpose to maintain consistency with scipy, numpy

comment:5 Changed 4 years ago by
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.
comment:6 followup: 7 Changed 3 years ago by
Replying to kcrisman:
Probably in principle the old behavior should be deprecated ... but the inconsistency is annoying, agreed.
The proposed change does not change the syntax of using the .LU()
method, but changes the semantic of the output. How would it be possible to deprecate it?
Do you say that if the .LU()
method is used for a matrix over RDF then a deprecation warning should be used? That warning would appear also in the case that the user expects the new (not the old) behaviour, i.e., when a warning is inappropriate.
comment:7 followup: 8 Changed 3 years ago by
The proposed change does not change the syntax of using the
.LU()
method, but changes the semantic of the output. How would it be possible to deprecate it?
Good point. Perhaps we need to just put a visible note in the documentation that might help clarify if someone is perplexed by the output changing  I guess that was my point, that it is a "silent" change. But they should agree over all rings, yes.
comment:8 followup: 9 Changed 3 years ago by
Replying to kcrisman:
The proposed change does not change the syntax of using the
.LU()
method, but changes the semantic of the output. How would it be possible to deprecate it?Good point. Perhaps we need to just put a visible note in the documentation that might help clarify if someone is perplexed by the output changing  I guess that was my point, that it is a "silent" change. But they should agree over all rings, yes.
I think along with a note in the documentation LU()
should also give a warning on using for this behavior change, what do you think?
comment:9 Changed 3 years ago by
Replying to ghChamanAgrawal:
I think along with a note in the documentation
LU()
should also give a warning on using for this behavior change, what do you think?
I don't think so. A warning (in the sense of: "If M.LU() is called then a warning is printed") is disturbing, because the mere fact of calling .LU()
is not a problem. What is a potential problem is whether or not the user expects the output in a particular form. And we cannot guess on the expectation of the user. Hence, certainly it should be written in the documentation that there has been a change, give an example for which that change happens, and insert a link to this ticket. The syntax for such link in the docstrings is
See :trac:`18365`.
But there should be no warning inserted into the code.
comment:10 Changed 3 years ago by
Commit:  251254086d59c476e0b3db3b90c94c357063a357 → 6fdae81bf4fd40caff8e877a1b4caa8c58db4e66 

Branch pushed to git repo; I updated commit sha1. New commits:
6fdae81  Added Note for behaviour change

comment:11 Changed 3 years ago by
I have added a note for the change as below, any suggestion for change?
.. NOTE:: The behaviour of ``LU()`` has been changed. Earlier ``LU()`` returned ``P,L,U`` such that ``P*A=L*U``, where ``P`` represents the permutation and is the matrix inverse of the ``P`` returned by this method. The computation of this matrix inverse can be accomplished quickly with just a transpose as the matrix is orthogonal/unitary. For Details See :trac:`18365`. EXAMPLES:: sage: m = matrix(RDF,4,range(16)) sage: P,L,U = m.LU() sage: P*L*U # rel tol 2e16 [ 0.0 1.0 2.0 3.0] [ 4.0 5.0 6.0 7.0] [ 8.0 9.0 10.0 11.0] [12.0 13.0 14.0 15.0] Below example illustrate the change in behaviour of ``LU()``. :: sage: m == P*L*U True sage: P*m == L*U False
comment:12 Changed 3 years ago by
Milestone:  sage6.7 → sage9.1 

Summary:  Definition of LU descomposition of a matrix depends on the base ring → Definition of LU decomposition of a matrix depends on the base ring 
comment:13 Changed 3 years ago by
I suggest that this gets merged now, despite the backwards incompatibility.
comment:14 Changed 3 years ago by
Branch:  u/ghChamanAgrawal/18365_RDF_LU → u/ghmwageringel/18365 

Commit:  6fdae81bf4fd40caff8e877a1b4caa8c58db4e66 → 96a9f92abda84d6dd863c8604f23081425043110 
I have made small formatting changes to the documentation and removed trailing whitespace. I will let the patchbot run on this ticket once more and then set it to positive.
New commits:
96a9f92  18365: reviewer changes

comment:15 Changed 3 years ago by
Reviewers:  → Markus Wageringel 

comment:16 Changed 3 years ago by
Status:  needs_review → positive_review 

comment:17 Changed 3 years ago by
Branch:  u/ghmwageringel/18365 → 96a9f92abda84d6dd863c8604f23081425043110 

Resolution:  → fixed 
Status:  positive_review → closed 
It is worth noticing that the source code for matrices over
RDF
callsscipy.linalg.lu
which returns a triple(P,L,U)
such thatA=PLU
and then transposeP
with the documentation:So, i guess this extra transposition should be reverted, so that it becomes consistent with both
scipy
and the implementation for exact rings.