Opened 6 years ago

Closed 11 months ago

# Definition of LU decomposition of a matrix depends on the base ring

Reported by: Owned by: tmonteil major sage-9.1 linear algebra kcrisman, gh-ChamanAgrawal Chaman Agrawal Markus Wageringel N/A 96a9f92 (Commits) 96a9f92abda84d6dd863c8604f23081425043110

### 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.

### comment:1 Changed 6 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:3 follow-up: ↓ 6 Changed 2 years ago by kcrisman

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

### comment:4 Changed 23 months ago by gh-ChamanAgrawal

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

New commits:

 ​2512540 `Reverted transpose to maintain consistency with scipy, numpy`

### comment:5 Changed 23 months 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.

### comment:6 in reply to: ↑ 3 ; follow-up: ↓ 7 Changed 22 months ago by SimonKing

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 in reply to: ↑ 6 ; follow-up: ↓ 8 Changed 22 months ago by 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.

### comment:8 in reply to: ↑ 7 ; follow-up: ↓ 9 Changed 22 months ago by gh-ChamanAgrawal

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 in reply to: ↑ 8 Changed 22 months ago by SimonKing

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`.
```
Version 0, edited 22 months ago by SimonKing (next)

### comment:10 Changed 22 months ago by git

• Commit changed from 251254086d59c476e0b3db3b90c94c357063a357 to 6fdae81bf4fd40caff8e877a1b4caa8c58db4e66

Branch pushed to git repo; I updated commit sha1. New commits:

 ​6fdae81 `Added Note for behaviour change`

### comment:11 Changed 22 months ago by gh-ChamanAgrawal

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 2e-16
[ 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 13 months ago by chapoton

• Milestone changed from sage-6.7 to sage-9.1
• Summary changed from Definition of LU descomposition of a matrix depends on the base ring to Definition of LU decomposition of a matrix depends on the base ring

### comment:13 Changed 11 months ago by gh-mwageringel

I suggest that this gets merged now, despite the backwards incompatibility.

### comment:14 Changed 11 months ago by gh-mwageringel

• Branch changed from u/gh-ChamanAgrawal/18365_RDF_LU to u/gh-mwageringel/18365
• Commit changed from 6fdae81bf4fd40caff8e877a1b4caa8c58db4e66 to 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 11 months ago by gh-mwageringel

• Reviewers set to Markus Wageringel

### comment:16 Changed 11 months ago by gh-mwageringel

• Status changed from needs_review to positive_review

### comment:17 Changed 11 months ago by vbraun

• Branch changed from u/gh-mwageringel/18365 to 96a9f92abda84d6dd863c8604f23081425043110
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.