# 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 pr = pr * sum(tmp) return pr def elementwise_product(self, right): r""" Returns the elementwise product of two matrices of the same size (also known as the Hadamard product). INPUT: - right - the right operand of the product.  A matrix of the same size as self such that multiplication of elements of the base rings of self and right is defined, once Sage's coercion model is applied.  If the matrices have different sizes, or if multiplication of individual entries cannot be achieved, a TypeError will result. OUTPUT: A matrix of the same size as self and right.  The entry in location (i,j) of the output is the product of the two entries in location (i,j) of self and right (in that order). The parent of the result is determined by Sage's coercion model.  If the base rings are identical, then the result is dense or sparse according to this property for the left operand.  If the base rings must be adjusted for one, or both, matrices then the result will be sparse only if both operands are sparse.  No subdivisions are present in the result. If the type of the result is not to your liking, or the ring could be "tighter," adjust the operands with :meth:~sage.matrix.matrix0.Matrix.change_ring. Adjust sparse versus dense inputs with the methods :meth:~sage.matrix.matrix1.Matrix.sparse_matrix and :meth:~sage.matrix.matrix1.Matrix.dense_matrix. EXAMPLES:: sage: A = matrix(ZZ, 2, range(6)) sage: B = matrix(QQ, 2, [5, 1/3, 2/7, 11/2, -3/2, 8]) sage: C = A.elementwise_product(B) sage: C [   0  1/3  4/7] [33/2   -6   40] sage: C.parent() Full MatrixSpace of 2 by 3 dense matrices over Rational Field Notice the base ring of the results in the next two examples.  :: sage: D = matrix(ZZ[x],2,[1+x^2,2,3,4-x]) sage: E = matrix(QQ,2,[1,2,3,4]) sage: F = D.elementwise_product(E) sage: F [  x^2 + 1         4] [        9 -4*x + 16] sage: F.parent() Full MatrixSpace of 2 by 2 dense matrices over Univariate Polynomial Ring in x over Rational Field :: sage: G = matrix(GF(3),2,[0,1,2,2]) sage: H = matrix(ZZ,2,[1,2,3,4]) sage: J = G.elementwise_product(H) sage: J [0 2] [0 2] sage: J.parent() Full MatrixSpace of 2 by 2 dense matrices over Finite Field of size 3 Non-commutative rings behave as expected.  These are the usual quaternions. :: sage: R. = QuaternionAlgebra(-1, -1) sage: A = matrix(R, 2, [1,i,j,k]) sage: B = matrix(R, 2, [i,i,i,i]) sage: A.elementwise_product(B) [ i -1] [-k  j] sage: B.elementwise_product(A) [ i -1] [ k -j] Input that is not a matrix will raise an error.  :: sage: A = random_matrix(ZZ,5,10,x=20) sage: A.elementwise_product(vector(ZZ, [1,2,3,4])) Traceback (most recent call last): ... TypeError: operand must be a matrix, not an element of Ambient free module of rank 4 over the principal ideal domain Integer Ring Matrices of different sizes for operands will raise an error. :: sage: A = random_matrix(ZZ,5,10,x=20) sage: B = random_matrix(ZZ,10,5,x=40) sage: A.elementwise_product(B) Traceback (most recent call last): ... 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 Some pairs of rings do not have a common parent where multiplication makes sense.  This will raise an error. :: sage: A = matrix(QQ, 3, range(6)) sage: B = matrix(GF(3), 3, [2]*6) sage: A.elementwise_product(B) Traceback (most recent call last): ... 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' We illustrate various combinations of sparse and dense matrices. Notice how if base rings are unequal, both operands must be sparse to get a sparse result.  When the base rings are equal, the left operand dictates the sparse/dense property of the result. This behavior is entirely a consequence of the coercion model.  :: sage: A = matrix(ZZ, 5, range(30), sparse=False) sage: B = matrix(ZZ, 5, range(30), sparse=True) sage: C = matrix(QQ, 5, range(30), sparse=True) sage: A.elementwise_product(C).is_dense() True sage: B.elementwise_product(C).is_sparse() True sage: A.elementwise_product(B).is_dense() True sage: B.elementwise_product(A).is_sparse() True TESTS: Implementation for dense and sparse matrices are different, this will provide a trivial test that they are working identically. :: sage: A = random_matrix(ZZ, 10, x=1000, sparse=False) sage: B = random_matrix(ZZ, 10, x=1000, sparse=False) sage: C = A.sparse_matrix() sage: D = B.sparse_matrix() sage: E = A.elementwise_product(B) sage: F = C.elementwise_product(D) sage: E.is_dense() and F.is_sparse() and (E == F) True If the ring has zero divisors, the routines for setting entries of a sparse matrix should intercept zero results and not create an entry. :: sage: R = Integers(6) sage: A = matrix(R, 2, [3, 2, 0, 0], sparse=True) sage: B = matrix(R, 2, [2, 3, 1, 0], sparse=True) sage: C = A.elementwise_product(B) sage: len(C.nonzero_positions()) == 0 True AUTHOR: - Rob Beezer (2009-07-13) """ # Optimized routines for specialized classes would likely be faster # See the "pairwise_product" of vectors for some guidance on doing this from sage.structure.element import canonical_coercion if not PY_TYPE_CHECK(right, Matrix): raise TypeError('operand must be a matrix, not an element of %s' % right.parent()) if (self.nrows() != right.nrows()) or (self.ncols() != right.ncols()): raise TypeError('incompatible sizes for matrices from: %s and %s'%(self.parent(), right.parent())) if self._parent is not (right)._parent: self, right = canonical_coercion(self, right) return self._elementwise_product(right) def permanent(self): r""" 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 [nr - t for t in reversed(row_divs)]) return atrans def _elementwise_product(self, right): r""" Returns the elementwise product of two dense matrices with identical base rings. This routine assumes that self and right are both matrices, both dense, with identical sizes and with identical base rings.  It is "unsafe" in the sense that these conditions are not checked and no sensible errors are raised. This routine is meant to be called from the meth:~sage.matrix.matrix2.Matrix.elementwise_product method, which will ensure that this routine receives proper input.  More thorough documentation is provided there. EXAMPLE:: sage: A = matrix(ZZ, 2, range(6), sparse=False) sage: B = matrix(ZZ, 2, [1,0,2,0,3,0], sparse=False) sage: A._elementwise_product(B) [ 0  0  4] [ 0 12  0] AUTHOR: - Rob Beezer (2009-07-14) """ cdef Py_ssize_t r, c cdef Matrix_dense other, prod nc, nr = self.ncols(), self.nrows() other = right prod = self.new_matrix(nr, nc, copy=False, coerce=False) for r in range(nr): for c in range(nc): entry = self.get_unsafe(r,c)*other.get_unsafe(r,c) prod.set_unsafe(r,c,entry) return prod def apply_morphism(self, phi): """ 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 self.cache('charpoly', f) return f def _elementwise_product(self, right): r""" Returns the elementwise product of two sparse matrices with identical base rings. This routine assumes that self and right are both matrices, both sparse, with identical sizes and with identical base rings.  It is "unsafe" in the sense that these conditions are not checked and no sensible errors are raised. This routine is meant to be called from the meth:~sage.matrix.matrix2.Matrix.elementwise_product method, which will ensure that this routine receives proper input.  More thorough documentation is provided there. EXAMPLE:: sage: A = matrix(ZZ, 2, range(6), sparse=True) sage: B = matrix(ZZ, 2, [1,0,2,0,3,0], sparse=True) sage: A._elementwise_product(B) [ 0  0  4] [ 0 12  0] AUTHOR: - Rob Beezer (2009-07-14) """ cdef Py_ssize_t k, r, c cdef Matrix_sparse other, prod nc, nr = self.ncols(), self.nrows() other = right prod = self.new_matrix(nr, nc, copy=False, coerce=False) nzleft = self.nonzero_positions(copy=False) nzright = other.nonzero_positions(copy=False) for k from 0 <= k < len(nzleft): r = get_ij(nzleft, k, 0) c = get_ij(nzleft, k, 1) if (r,c) in nzright: entry = self.get_unsafe(r,c)*other.get_unsafe(r,c) prod.set_unsafe(r,c,entry) return prod def apply_morphism(self, phi): """ Apply the morphism phi to the coefficients of this sparse matrix.