Skip to content

Commit

Permalink
Trac #28402: incorrect inverse of sparse matrix over inexact rings
Browse files Browse the repository at this point in the history
As reported on [https://ask.sagemath.org/question/47587/inverse-of-real-
sparse-matrix/ this ask question], we have:

{{{
sage: B=matrix(RR, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0],
[-1/24,1/60,1/60, 1/420, -1/24, 0, 0, 0, 0], [-1/30,1/60, 2/105, 1/140,
-1/20, 0, 0, 0, 0], [1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0
....: ], [1/12, -1/24, -1/20, -1/40, 1/3, -1/24, -1/30, 1/120,1/12], [0,
0, 0, 0, -1/24,1/60,1/60, 1/420, -1/24], [0, 0, 0, 0, -1/30,1/60, 2/105,
1/140, -1/20], [0, 0, 0, 0, 1/120, 1/420, 1/140, 13/1260, -
....: 1/40], [0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]],sparse=true)
sage: (B.inverse()*B).norm(1)
138.4999999999923
}}}

Note that reverting #24122, by removing the call to
`build_inverse_from_augmented_sparse` in the `__invert__` method of
`sage/matrix/matrix0.pyx` (line 5410), leads to a correct answer.

URL: https://trac.sagemath.org/28402
Reported by: tmonteil
Ticket author(s): Travis Scrimshaw
Reviewer(s): Thierry Monteil
  • Loading branch information
Release Manager committed Sep 4, 2019
2 parents e049a00 + df1cbce commit d2eae07
Showing 1 changed file with 29 additions and 3 deletions.
32 changes: 29 additions & 3 deletions src/sage/matrix/matrix0.pyx
Expand Up @@ -5356,6 +5356,33 @@ cdef class Matrix(sage.structure.element.Matrix):
sage: N = ~M
sage: (N*M).norm()
0.9999999999999999
Check that :trac:`28402` is fixed::
sage: B = matrix(RR, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0],
....: [-1/24,1/60,1/60, 1/420, -1/24, 0, 0, 0, 0],
....: [-1/30,1/60, 2/105, 1/140, -1/20, 0, 0, 0, 0],
....: [1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0],
....: [1/12, -1/24, -1/20, -1/40, 1/3, -1/24, -1/30, 1/120,1/12],
....: [0, 0, 0, 0, -1/24,1/60,1/60, 1/420, -1/24],
....: [0, 0, 0, 0, -1/30,1/60, 2/105, 1/140, -1/20],
....: [0, 0, 0, 0, 1/120, 1/420, 1/140, 13/1260, -1/40],
....: [0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]],
....: sparse=True)
sage: (B.inverse()*B).norm(1) # rel tol 2e-12
1.0
sage: B = matrix(QQ, [[1/6, -1/24, -1/30, 1/120,1/12, 0, 0, 0, 0],
....: [-1/24,1/60,1/60, 1/420, -1/24, 0, 0, 0, 0],
....: [-1/30,1/60, 2/105, 1/140, -1/20, 0, 0, 0, 0],
....: [1/120, 1/420, 1/140, 13/1260, -1/40, 0, 0, 0, 0],
....: [1/12, -1/24, -1/20, -1/40, 1/3, -1/24, -1/30, 1/120,1/12],
....: [0, 0, 0, 0, -1/24,1/60,1/60, 1/420, -1/24],
....: [0, 0, 0, 0, -1/30,1/60, 2/105, 1/140, -1/20],
....: [0, 0, 0, 0, 1/120, 1/420, 1/140, 13/1260, -1/40],
....: [0, 0, 0, 0,1/12, -1/24, -1/20, -1/40, 1/6]],
....: sparse=True)
sage: (B.inverse()*B).norm(1)
1.0
"""
if not self.base_ring().is_field():
try:
Expand Down Expand Up @@ -5403,13 +5430,12 @@ cdef class Matrix(sage.structure.element.Matrix):
if self.base_ring().is_exact():
if A[self._nrows-1, self._ncols-1] != 1:
raise ZeroDivisionError("input matrix must be nonsingular")
if self.is_sparse():
return self.build_inverse_from_augmented_sparse(A)
else:
if not A[self._nrows-1, self._ncols-1]:
raise ZeroDivisionError("input matrix must be nonsingular")

if self.is_sparse():
return self.build_inverse_from_augmented_sparse(A)

return A.matrix_from_columns(list(xrange(self._ncols, 2 * self._ncols)))

cdef build_inverse_from_augmented_sparse(self, A):
Expand Down

0 comments on commit d2eae07

Please sign in to comment.