Ticket #6301: trac_6301_elementwise_matrix_product.patch

File trac_6301_elementwise_matrix_product.patch, 11.5 KB (added by Rob Beezer, 14 years ago)
  • sage/matrix/matrix2.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1247609671 25200
    # Node ID 910a3bd6262f5ded569c185bf313ab0a067ecfe4
    # Parent  48b53a103bc444ea3ce12b822dcdd7e04969b8f9
    Implement elementwise (Hadamard) product of matrices
    
    diff -r 48b53a103bc4 -r 910a3bd6262f sage/matrix/matrix2.pyx
    a b  
    436436            pr = pr * sum(tmp)
    437437        return pr
    438438
     439    def elementwise_product(self, right):
     440        r"""
     441        Returns the elementwise product of two matrices
     442        of the same size (also known as the Hadamard product).
     443
     444        INPUT:
     445
     446        - ``right`` - the right operand of the product.  A matrix
     447          of the same size as ``self`` such that multiplication
     448          of elements of the base rings of ``self`` and ``right``
     449          is defined, once Sage's coercion model is applied.  If
     450          the matrices have different sizes, or if multiplication
     451          of individual entries cannot be achieved, a ``TypeError``
     452          will result.
     453
     454        OUTPUT:
     455
     456        A matrix of the same size as ``self`` and ``right``.  The
     457        entry in location `(i,j)` of the output is the product of
     458        the two entries in location `(i,j)` of ``self`` and
     459        ``right`` (in that order).
     460
     461        The parent of the result is determined by Sage's coercion
     462        model.  If the base rings are identical, then the result
     463        is dense or sparse according to this property for
     464        the left operand.  If the base rings must be adjusted
     465        for one, or both, matrices then the result will be sparse
     466        only if both operands are sparse.  No subdivisions are
     467        present in the result.
     468
     469        If the type of the result is not to your liking, or
     470        the ring could be "tighter," adjust the operands with
     471        :meth:`~sage.matrix.matrix0.Matrix.change_ring`.
     472        Adjust sparse versus dense inputs with the methods
     473        :meth:`~sage.matrix.matrix1.Matrix.sparse_matrix` and
     474        :meth:`~sage.matrix.matrix1.Matrix.dense_matrix`.
     475
     476        EXAMPLES::
     477
     478            sage: A = matrix(ZZ, 2, range(6))
     479            sage: B = matrix(QQ, 2, [5, 1/3, 2/7, 11/2, -3/2, 8])
     480            sage: C = A.elementwise_product(B)
     481            sage: C
     482            [   0  1/3  4/7]
     483            [33/2   -6   40]
     484            sage: C.parent()
     485            Full MatrixSpace of 2 by 3 dense matrices over Rational Field
     486
     487
     488        Notice the base ring of the results in the next two examples.  ::
     489
     490            sage: D = matrix(ZZ[x],2,[1+x^2,2,3,4-x])
     491            sage: E = matrix(QQ,2,[1,2,3,4])
     492            sage: F = D.elementwise_product(E)
     493            sage: F
     494            [  x^2 + 1         4]
     495            [        9 -4*x + 16]
     496            sage: F.parent()
     497            Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in x over Rational Field
     498
     499        ::
     500
     501            sage: G = matrix(GF(3),2,[0,1,2,2])
     502            sage: H = matrix(ZZ,2,[1,2,3,4])
     503            sage: J = G.elementwise_product(H)
     504            sage: J
     505            [0 2]
     506            [0 2]
     507            sage: J.parent()
     508            Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 3
     509
     510        Non-commutative rings behave as expected.  These are the usual quaternions. ::
     511
     512            sage: R.<i,j,k> = QuaternionAlgebra(-1, -1)
     513            sage: A = matrix(R, 2, [1,i,j,k])
     514            sage: B = matrix(R, 2, [i,i,i,i])
     515            sage: A.elementwise_product(B)
     516            [ i -1]
     517            [-k  j]
     518            sage: B.elementwise_product(A)
     519            [ i -1]
     520            [ k -j]
     521
     522        Input that is not a matrix will raise an error.  ::
     523
     524            sage: A = random_matrix(ZZ,5,10,x=20)
     525            sage: A.elementwise_product(vector(ZZ, [1,2,3,4]))
     526            Traceback (most recent call last):
     527            ...
     528            TypeError: operand must be a matrix, not an element of Ambient free module of rank 4 over the principal ideal domain Integer Ring
     529
     530        Matrices of different sizes for operands will raise an error. ::
     531
     532            sage: A = random_matrix(ZZ,5,10,x=20)
     533            sage: B = random_matrix(ZZ,10,5,x=40)
     534            sage: A.elementwise_product(B)
     535            Traceback (most recent call last):
     536            ...
     537            TypeError: incompatible sizes for matrices from: Full MatrixSpace of 5 by 10 dense matrices over Integer Ring and Full MatrixSpace of 10 by 5 dense matrices over Integer Ring
     538
     539        Some pairs of rings do not have a common parent where
     540        multiplication makes sense.  This will raise an error. ::
     541
     542            sage: A = matrix(QQ, 3, range(6))
     543            sage: B = matrix(GF(3), 3, [2]*6)
     544            sage: A.elementwise_product(B)
     545            Traceback (most recent call last):
     546            ...
     547            TypeError: no common canonical parent for objects with parents: 'Full MatrixSpace of 3 by 2 dense matrices over Rational Field' and 'Full MatrixSpace of 3 by 2 dense matrices over Finite Field of size 3'
     548
     549        We illustrate various combinations of sparse and dense matrices.
     550        Notice how if base rings are unequal, both operands must be sparse
     551        to get a sparse result.  When the base rings are equal, the left
     552        operand dictates the sparse/dense property of the result. This
     553        behavior is entirely a consequence of the coercion model.  ::
     554
     555            sage: A = matrix(ZZ, 5, range(30), sparse=False)
     556            sage: B = matrix(ZZ, 5, range(30), sparse=True)
     557            sage: C = matrix(QQ, 5, range(30), sparse=True)
     558            sage: A.elementwise_product(C).is_dense()
     559            True
     560            sage: B.elementwise_product(C).is_sparse()
     561            True
     562            sage: A.elementwise_product(B).is_dense()
     563            True
     564            sage: B.elementwise_product(A).is_sparse()
     565            True
     566
     567        TESTS:
     568
     569        Implementation for dense and sparse matrices are
     570        different, this will provide a trivial test that
     571        they are working identically. ::
     572
     573            sage: A = random_matrix(ZZ, 10, x=1000, sparse=False)
     574            sage: B = random_matrix(ZZ, 10, x=1000, sparse=False)
     575            sage: C = A.sparse_matrix()
     576            sage: D = B.sparse_matrix()
     577            sage: E = A.elementwise_product(B)
     578            sage: F = C.elementwise_product(D)
     579            sage: E.is_dense() and F.is_sparse() and (E == F)
     580            True
     581
     582        If the ring has zero divisors, the routines for setting
     583        entries of a sparse matrix should intercept zero results
     584        and not create an entry. ::
     585
     586            sage: R = Integers(6)
     587            sage: A = matrix(R, 2, [3, 2, 0, 0], sparse=True)
     588            sage: B = matrix(R, 2, [2, 3, 1, 0], sparse=True)
     589            sage: C = A.elementwise_product(B)
     590            sage: len(C.nonzero_positions()) == 0
     591            True
     592
     593        AUTHOR:
     594
     595        - Rob Beezer (2009-07-13)
     596        """
     597        # Optimized routines for specialized classes would likely be faster
     598        # See the "pairwise_product" of vectors for some guidance on doing this
     599        from sage.structure.element import canonical_coercion
     600        if not PY_TYPE_CHECK(right, Matrix):
     601            raise TypeError('operand must be a matrix, not an element of %s' % right.parent())
     602        if (self.nrows() != right.nrows()) or (self.ncols() != right.ncols()):
     603            raise TypeError('incompatible sizes for matrices from: %s and %s'%(self.parent(), right.parent()))
     604        if self._parent is not (<Matrix>right)._parent:
     605            self, right = canonical_coercion(self, right)
     606        return self._elementwise_product(right)
     607
    439608    def permanent(self):
    440609        r"""
    441610        Calculate and return the permanent of this `m \times n`
  • sage/matrix/matrix_dense.pyx

    diff -r 48b53a103bc4 -r 910a3bd6262f sage/matrix/matrix_dense.pyx
    a b  
    248248                             [nr - t for t in reversed(row_divs)])
    249249        return atrans
    250250
     251    def _elementwise_product(self, right):
     252        r"""
     253        Returns the elementwise product of two dense
     254        matrices with identical base rings.
     255
     256        This routine assumes that ``self`` and ``right``
     257        are both matrices, both dense, with identical
     258        sizes and with identical base rings.  It is
     259        "unsafe" in the sense that these conditions
     260        are not checked and no sensible errors are
     261        raised.
     262
     263        This routine is meant to be called from the
     264        meth:`~sage.matrix.matrix2.Matrix.elementwise_product`
     265        method, which will ensure that this routine receives
     266        proper input.  More thorough documentation is provided
     267        there.
     268
     269        EXAMPLE::
     270
     271            sage: A = matrix(ZZ, 2, range(6), sparse=False)
     272            sage: B = matrix(ZZ, 2, [1,0,2,0,3,0], sparse=False)
     273            sage: A._elementwise_product(B)
     274            [ 0  0  4]
     275            [ 0 12  0]
     276
     277        AUTHOR:
     278
     279        - Rob Beezer (2009-07-14)
     280        """
     281        cdef Py_ssize_t r, c
     282        cdef Matrix_dense other, prod
     283
     284        nc, nr = self.ncols(), self.nrows()
     285        other = right
     286        prod = self.new_matrix(nr, nc, copy=False, coerce=False)
     287        for r in range(nr):
     288            for c in range(nc):
     289                entry = self.get_unsafe(r,c)*other.get_unsafe(r,c)
     290                prod.set_unsafe(r,c,entry)
     291        return prod
     292
    251293    def apply_morphism(self, phi):
    252294        """
    253295        Apply the morphism phi to the coefficients of this dense matrix.
  • sage/matrix/matrix_sparse.pyx

    diff -r 48b53a103bc4 -r 910a3bd6262f sage/matrix/matrix_sparse.pyx
    a b  
    384384        self.cache('charpoly', f)
    385385        return f
    386386
     387    def _elementwise_product(self, right):
     388        r"""
     389        Returns the elementwise product of two sparse
     390        matrices with identical base rings.
     391
     392        This routine assumes that ``self`` and ``right``
     393        are both matrices, both sparse, with identical
     394        sizes and with identical base rings.  It is
     395        "unsafe" in the sense that these conditions
     396        are not checked and no sensible errors are
     397        raised.
     398
     399        This routine is meant to be called from the
     400        meth:`~sage.matrix.matrix2.Matrix.elementwise_product`
     401        method, which will ensure that this routine receives
     402        proper input.  More thorough documentation is provided
     403        there.
     404
     405        EXAMPLE::
     406
     407            sage: A = matrix(ZZ, 2, range(6), sparse=True)
     408            sage: B = matrix(ZZ, 2, [1,0,2,0,3,0], sparse=True)
     409            sage: A._elementwise_product(B)
     410            [ 0  0  4]
     411            [ 0 12  0]
     412
     413        AUTHOR:
     414
     415        - Rob Beezer (2009-07-14)
     416        """
     417        cdef Py_ssize_t k, r, c
     418        cdef Matrix_sparse other, prod
     419
     420        nc, nr = self.ncols(), self.nrows()
     421        other = right
     422        prod = self.new_matrix(nr, nc, copy=False, coerce=False)
     423        nzleft = self.nonzero_positions(copy=False)
     424        nzright = other.nonzero_positions(copy=False)
     425        for k from 0 <= k < len(nzleft):
     426            r = get_ij(nzleft, k, 0)
     427            c = get_ij(nzleft, k, 1)
     428            if (r,c) in nzright:
     429                entry = self.get_unsafe(r,c)*other.get_unsafe(r,c)
     430                prod.set_unsafe(r,c,entry)
     431        return prod
     432
    387433    def apply_morphism(self, phi):
    388434        """
    389435        Apply the morphism phi to the coefficients of this sparse matrix.