Opened 9 months ago
Closed 9 months ago
#27352 closed defect (fixed)
Add checks for matrix multiplication
Reported by:  slelievre  Owned by:  

Priority:  major  Milestone:  sage8.7 
Component:  linear algebra  Keywords:  
Cc:  slelievre  Merged in:  
Authors:  John Palmieri  Reviewers:  Vincent Delecroix 
Report Upstream:  N/A  Work issues:  
Branch:  a2bac95 (Commits)  Commit:  a2bac95e56e6855f3393038c04c2eb537929dc30 
Dependencies:  Stopgaps: 
Description
It seems the compatibility of dimensions is not always tested when multiplying matrices.
This ticket is opened following the report at
In particular, multiplying two 3 by 2 matrices should fail, but:
sage: A = matrix(QQ, [[1, 2], [1, 0], [1, 1]]) sage: B = matrix(QQ, [[0, 4], [1, 1], [1, 2]]) sage: A [ 1 2] [1 0] [ 1 1] sage: B [ 0 4] [ 1 1] [ 1 2] sage: A*B [ 1 0] [ 1 2] [ 1 3]
In this case A * B
amounts to A.__mul__(B)
which ends up calling A._multiply_flint(B)
.
For matrices over the integers, the multiplication above fails as it should:
sage: A = matrix(ZZ, [[1, 2], [1, 0], [1, 1]]) sage: B = matrix(ZZ, [[0, 4], [1, 1], [1, 2]]) sage: A * B Traceback (most recent call last) ... IndexError: Number of columns of self must equal number of rows of right.
but we get a similar surprise by calling
sage: A._multiply_linbox(B) [ 2 2] [ 0 4] [ 1 3]
In case _multiply_linbox
and _multiply_flint
skip dimension tests for speed,
tests should be performed before calling them, which seems to be the case for
_multiply_linbox
but not _multiply_flint
.
Change History (16)
comment:1 Changed 9 months ago by
comment:2 Changed 9 months ago by
Those random answers are likely because upstream uses whatever is in that next spot in memory. I think this also has a chance of causing a segfault or destroying some other computation for the same reason.
comment:3 Changed 9 months ago by
Regarding the ticket description, is _multiply_linbox
used anywhere? _multiply_flint
is funny because it has a doctest
sage: matrix(QQ, 2, 3) * matrix(QQ, 4, 5) Traceback (most recent call last): ... TypeError: unsupported operand parent(s) for *: 'Full MatrixSpace of 2 by 3 dense matrices over Rational Field' and 'Full MatrixSpace of 4 by 5 dense matrices over Rational Field'
but I think if the two matrices have the same parent, then the shapes aren't checked.
comment:4 followup: ↓ 5 Changed 9 months ago by
Can we just do something like this?

src/sage/structure/element.pyx
diff git a/src/sage/structure/element.pyx b/src/sage/structure/element.pyx index 1b167b7ab4..eef1798d54 100644
a b cdef class Matrix(ModuleElement): 3676 3676 """ 3677 3677 cdef int cl = classify_elements(left, right) 3678 3678 if HAVE_SAME_PARENT(cl): 3679 return (<Matrix>left)._matrix_times_matrix_(<Matrix>right) 3679 if left.ncols() == right.nrows(): 3680 return (<Matrix>left)._matrix_times_matrix_(<Matrix>right) 3681 else: 3682 raise ValueError 3680 3683 if BOTH_ARE_ELEMENT(cl): 3681 3684 return coercion_model.bin_op(left, right, mul) 3682 3685
(with an appropriate message for the ValueError
)? I don't know the element.pyx
code well. Or should the change be done in _matrix_times_matrix
in matrix_rational_dense.pyx
?
comment:5 in reply to: ↑ 4 ; followup: ↓ 6 Changed 9 months ago by
Replying to jhpalmieri:
Can we just do something like this?
src/sage/structure/element.pyx
diff git a/src/sage/structure/element.pyx b/src/sage/structure/element.pyx index 1b167b7ab4..eef1798d54 100644
a b cdef class Matrix(ModuleElement): 3676 3676 """ 3677 3677 cdef int cl = classify_elements(left, right) 3678 3678 if HAVE_SAME_PARENT(cl): 3679 return (<Matrix>left)._matrix_times_matrix_(<Matrix>right) 3679 if left.ncols() == right.nrows(): 3680 return (<Matrix>left)._matrix_times_matrix_(<Matrix>right) 3681 else: 3682 raise ValueError 3680 3683 if BOTH_ARE_ELEMENT(cl): 3681 3684 return coercion_model.bin_op(left, right, mul) 3682 3685 (with an appropriate message for the
ValueError
)? I don't know theelement.pyx
code well. Or should the change be done in_matrix_times_matrix
inmatrix_rational_dense.pyx
?
Yes but use self._ncols
(Py_ssize_t
attribute) instead of ncols()
(Python function).
Note that the coercion model raises TypeError
sage: Zmod(3).an_element() * Zmod(2).an_element()
and I think it is better to follow this convention.
comment:6 in reply to: ↑ 5 Changed 9 months ago by
Replying to vdelecroix:
Replying to jhpalmieri:
Yes but use
self._ncols
(Py_ssize_t
attribute) instead ofncols()
(Python function).
This won't work, use
(<Matrix> left)._nrows != (<Matrix> left)._ncols
(note that since left
and right
have the same parent we want them to be square matrices)
comment:7 Changed 9 months ago by
 Branch set to u/jhpalmieri/matrixmult
comment:8 Changed 9 months ago by
 Commit set to 4a59d8e5d817894193b8f9ebe1ac5168d5954536
 Status changed from new to needs_review
Here's a first draft. Feel free to modify it.
New commits:
4a59d8e  trac 27352: add one more check on matrix sizes for matrix multiplication

comment:9 Changed 9 months ago by
 Reviewers set to Vincent Delecroix
 Status changed from needs_review to needs_work
`trac`:27352:
should be :trac:`27352`
comment:10 Changed 9 months ago by
 Commit changed from 4a59d8e5d817894193b8f9ebe1ac5168d5954536 to c9c2109d2c2647facf1781d40e416e91cf54b9a7
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
c9c2109  trac 27352: add one more check on matrix sizes for matrix multiplication

comment:11 Changed 9 months ago by
Fixed. (I built the reference manual, but of course this doesn't appear in the reference manual, so I didn't catch it.)
comment:13 Changed 9 months ago by
Usually, no space between the type declaration <Matrix>
and the name left
or right
.
See other uses elsewhere in the source code as found by search_src("<matrix>")
:
matrix/action.pyx:248: cdef Matrix A = <Matrix>g matrix/action.pyx:249: cdef Matrix B = <Matrix>s matrix/action.pyx:305: cdef Matrix A = <Matrix>g matrix/action.pyx:356: cdef Matrix A = <Matrix>g matrix/action.pyx:367: return (<Matrix>A)._vector_times_matrix_(v) # v * A matrix/args.pyx:652: M = <Matrix>self.entries matrix/args.pyx:880: m = <Matrix>self.entries matrix/matrix1.pyx:1461: other = <Matrix>bottom matrix/matrix1.pyx:1503: cdef Matrix other = <Matrix>bottom matrix/matrix2.pyx:790: return (<Matrix>A)._elementwise_product(B) matrix/matrix2.pyx:2450: cdef Matrix M = <Matrix> self matrix/matrix2.pyx:2478: cdef Matrix a = <Matrix> matrix(R, n1, n) matroids/lean_matrix.pyx:1033: if int((<Matrix>M).get_unsafe(i, j)) & 1: matroids/lean_matrix.pyx:1655: s = int((<Matrix>M).get_unsafe(i, j)) % 3 matroids/lean_matrix.pyx:2208: self._gf4 = (<Matrix>M).base_ring() matroids/lean_matrix.pyx:2215: self.set(i, j, (<Matrix>M).get_unsafe(i, j)) structure/element.pyx:3679: return (<Matrix>left)._matrix_times_matrix_(<Matrix>right)
Not sure if it's still time to change this here or it should happen in a followup ticket or we don't care.
comment:14 Changed 9 months ago by
 Commit changed from c9c2109d2c2647facf1781d40e416e91cf54b9a7 to a2bac95e56e6855f3393038c04c2eb537929dc30
 Status changed from positive_review to needs_review
Branch pushed to git repo; I updated commit sha1 and set ticket back to needs_review. This was a forced push. New commits:
a2bac95  trac 27352: add one more check on matrix sizes for matrix multiplication

comment:15 Changed 9 months ago by
 Status changed from needs_review to positive_review
I don't care much, but it's easy to change. If Volker complains, we will know it was too late.
comment:16 Changed 9 months ago by
 Branch changed from u/jhpalmieri/matrixmult to a2bac95e56e6855f3393038c04c2eb537929dc30
 Resolution set to fixed
 Status changed from positive_review to closed
With 8.7.beta5, I actually get random (?) answers: