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: sage-8.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 jhpalmieri

With 8.7.beta5, I actually get random (?) answers:

sage: A = matrix(QQ, [[1, 2], [-1, 0], [1, 1]])
sage: B = matrix(QQ, [[0, 4], [1, -1], [1, 2]])
sage: A * B
[          1           0]
[          1          -2]
[13204699393 26409398787]
sage: A * B
[          1           0]
[          1          -2]
[13226416241 26452832483]
sage: A * B
[          1           0]
[          1          -2]
[13226415713 26452831427]
sage: A * B
[          1           0]
[          1          -2]
[13226429585 26452859171]
sage: A * B
[          1           0]
[          1          -2]
[13226415905 26452831811]
sage: A * B
[                    1                     0]
[                    1                    -2]
[ -6917529027641081855 -13835058055282163709]
sage: A * B
[                  1                   0]
[                  1                  -2]
[-411335435850571967 -822670871701143933]
sage: A * B
[              1               0]
[              1              -2]
[140331841062305 280663682124611]
sage: A * B
[ 1  0]
[ 1 -2]
[ 1  3]

comment:2 Changed 9 months ago by tscrim

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 jhpalmieri

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 follow-up: Changed 9 months ago by 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): 
    36763676        """
    36773677        cdef int cl = classify_elements(left, right)
    36783678        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
    36803683        if BOTH_ARE_ELEMENT(cl):
    36813684            return coercion_model.bin_op(left, right, mul)
    36823685

(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 ; follow-up: Changed 9 months ago by vdelecroix

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): 
    36763676        """
    36773677        cdef int cl = classify_elements(left, right)
    36783678        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
    36803683        if BOTH_ARE_ELEMENT(cl):
    36813684            return coercion_model.bin_op(left, right, mul)
    36823685

(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?

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 vdelecroix

Replying to vdelecroix:

Replying to jhpalmieri:

Yes but use self._ncols (Py_ssize_t attribute) instead of ncols() (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 jhpalmieri

  • Branch set to u/jhpalmieri/matrix-mult

comment:8 Changed 9 months ago by jhpalmieri

  • Authors set to John Palmieri
  • Commit set to 4a59d8e5d817894193b8f9ebe1ac5168d5954536
  • Status changed from new to needs_review

Here's a first draft. Feel free to modify it.


New commits:

4a59d8etrac 27352: add one more check on matrix sizes for matrix multiplication

comment:9 Changed 9 months ago by vdelecroix

  • 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 git

  • Commit changed from 4a59d8e5d817894193b8f9ebe1ac5168d5954536 to c9c2109d2c2647facf1781d40e416e91cf54b9a7

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

c9c2109trac 27352: add one more check on matrix sizes for matrix multiplication

comment:11 Changed 9 months ago by jhpalmieri

Fixed. (I built the reference manual, but of course this doesn't appear in the reference manual, so I didn't catch it.)

comment:12 Changed 9 months ago by vdelecroix

  • Status changed from needs_work to positive_review

Indeed...

comment:13 Changed 9 months ago by slelievre

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, n-1, 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 follow-up ticket or we don't care.

comment:14 Changed 9 months ago by git

  • 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:

a2bac95trac 27352: add one more check on matrix sizes for matrix multiplication

comment:15 Changed 9 months ago by jhpalmieri

  • 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 vbraun

  • Branch changed from u/jhpalmieri/matrix-mult to a2bac95e56e6855f3393038c04c2eb537929dc30
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.