The first time you see matrices, if someone asked you how you multiply two matrices together, your first idea might be to multiply every element of the first matrix by the element in the same position of the corresponding matrix, analogous to the way you add matrices.

But that’s not usually how we multiply matrices. That notion of multiplication hardly involves the matrix structure; it treats the matrix as an ordered container of numbers, but not as a way of representing a linear transformation. Once you have a little experience with linear algebra, the customary way of multiplying matrices seems natural, and the way that may have seemed natural at first glance seems kinda strange.

The componentwise product of matrices is called the **Hadamard product** or sometimes the **Schur product**. Given two *m* by *n* matrices *A* and *B*, the Hadamard product of *A* and *B*, written *A* ∘ *B*, is the *m* by *n* matrix *C* with elements given by

*c*_{ij} = *a*_{ij} *b*_{ij}.

Because the Hadamard product hardly uses the linear structure of a matrix, you wouldn’t expect it to interact nicely with operations that depend critically on the linear structure. And yet we can give a couple theorems that do show a nice interaction, at least when *A* and *B* are positive semi-definite matrices.

The first is the **Schur product theorem**. It says that if *A* and *B* are positive semi-definite *n* by *n* matrices, then

det(A ∘ *B*) ≥ det(*A*) det(*B*)

where det stands for determinant.

Also, there is the following theorem of Pólya and Szegö. Assume *A* and *B* are symmetric positive semi-definite *n* by *n* matrices. If the eigenvalues of *A* and *B*, listed in increasing order, are α_{i} and β_{i} respectively, then for every eigenvalue λ of A ∘ *B*, we have

α_{1} β_{1} ≤ λ ≤ α_{n} β_{n}.

## Python implementation

If you multiply two (multidimensional) arrays in NumPy, you’ll get the componentwise product. So if you multiply two matrices *as arrays* you’ll get the Hadamard product, but if you multiply them *as matrices* you’ll get the usual matrix product. We’ll illustrate that below. Note that the function `eigvalsh`

returns the eigenvalues of a matrix. The name may look a little strange, but the “h” on the end stands for “Hermitian.” We’re telling NumPy that the matrix is Hermitian so it can run software specialized for that case [1].

from numpy import array, matrix, array_equal, all from numpy.linalg import det, eigvalsh A = array([ [3, 1], [1, 3] ]) B = array([ [5, -1], [-1, 5] ]) H = array([ [15, -1], [-1, 15] ]) AB = array([ [14, 2], [ 2, 14] ]) # Hadamard product assert(array_equal(A*B, H)) # Ordinary matrix product assert(array_equal(A@B, AB)) # Schur product theorem assert(det(H) >= det(A)*det(B)) # Eigenvalues eigA = eigvalsh(A) eigB = eigvalsh(B) eigH = eigvalsh(A*B) lower = eigA[0]*eigB[0] upper = eigA[1]*eigB[1] assert(all(eigH >= lower)) assert(all(eigH <= upper))

The code above shows that the eigenvalues of *A* are [2, 4], the eigenvalues of *B* are [4, 6], and the eigenvalues of A ∘ *B* are [14, 16].

## Related posts

- How fast can you multiply matrices?
- Computing SVD and pseudoinverse
- Circular law for random matrices

[1] For complex matrices, Hermitian means conjugate symmetric, which in the real case reduces to simply symmetric. The theorem of Pólya and Szegö is actually valid for Hermitian matrices, but I simplified the statement for the case of real-valued matrices.