Repository: sympy/sympy

Issue:
HadamardPower(A, B).as_explicit() fails
Matrix to matrix power fails:

```python
>>> A = MatrixSymbol('A', 2, 3)
>>> HadamardPower(A, A).as_explicit()
---------------------------------------------------------------------------

NotImplementedError: Matrix Power not defined
```

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v145). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* matrices
  - `HadamardProduct.as_explicit` now supports matrix to matrix exponential. ([#16508](https://github.com/sympy/sympy/pull/16508) by [@sylee957](https://github.com/sylee957))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.5.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234". See
    https://github.com/blog/1506-closing-issues-via-pull-requests . Please also
    write a comment on that issue linking back to this pull request once it is
    open. -->

    Fixes #16505

    #### Brief description of what is fixed or changed

    If both base and exponent are matrix, hadamard power can be defined as 
    ![image](https://user-images.githubusercontent.com/34944973/55284042-e3ec3f80-53a9-11e9-82a5-f16b29465304.png)

    And since hadamard exponential and hadamard log are elementwise exponential and logarithm, it can be reduced to
    ![image](https://user-images.githubusercontent.com/34944973/55284069-53fac580-53aa-11e9-9560-cdc6ebdd1579.png)
    for each matrix elements.


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. See
    https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
    on how to write release notes. The bot will check your release notes
    automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    - matrices
      - `HadamardProduct.as_explicit` now supports matrix to matrix exponential.
    <!-- END RELEASE NOTES -->


</details><p>

codecov[bot]: # [Codecov](https://codecov.io/gh/sympy/sympy/pull/16508?src=pr&el=h1) Report
> Merging [#16508](https://codecov.io/gh/sympy/sympy/pull/16508?src=pr&el=desc) into [master](https://codecov.io/gh/sympy/sympy/commit/a7e6f093c98a3c4783848a19fce646e32b6e0161?src=pr&el=desc) will **decrease** coverage by `0.022%`.
> The diff coverage is `66.666%`.

```diff
@@             Coverage Diff              @@
##           master    #16508       +/-   ##
============================================
- Coverage   73.69%   73.667%   -0.023%     
============================================
  Files         618       618               
  Lines      158567    158575        +8     
  Branches    37224     37227        +3     
============================================
- Hits       116849    116819       -30     
- Misses      36302     36334       +32     
- Partials     5416      5422        +6
```

oscarbenjamin: This has merge conflicts now...
sylee957: Whatever, I couldn't find any good article for the hadamard product notation to be defined between a matrix and scalar, or the hadamard power notation to be defined between a matrix and a matrix, but at least having an option to raise a matrix base to a matrix exponent can be convenient if defined.

If other properties related to the hadamard product or hadamard power do not break with such new definition, 
I would define a hadamard product between a matrix and a scalar as a usual matrix-scalar multiplication, and
hadamard power between matrices as entrywise power between each matrix elements, as used in most of the numeric libraries, unless some other people can tackle in.

Patch:
diff --git a/doc/src/modules/matrices/expressions.rst b/doc/src/modules/matrices/expressions.rst
index efe6837eb6..64897b8865 100644
--- a/doc/src/modules/matrices/expressions.rst
+++ b/doc/src/modules/matrices/expressions.rst
@@ -34,6 +34,10 @@ Matrix Expressions Core Reference
    :members:
 .. autoclass:: MatPow
    :members:
+.. autoclass:: HadamardProduct
+   :members:
+.. autoclass:: HadamardPower
+   :members:
 .. autoclass:: Inverse
    :members:
 .. autoclass:: Transpose
diff --git a/sympy/matrices/expressions/hadamard.py b/sympy/matrices/expressions/hadamard.py
index 36e1ebe5b9..8657f3c5e8 100644
--- a/sympy/matrices/expressions/hadamard.py
+++ b/sympy/matrices/expressions/hadamard.py
@@ -294,13 +294,75 @@ def hadamard_power(base, exp):
 
 
 class HadamardPower(MatrixExpr):
-    """
+    r"""
     Elementwise power of matrix expressions
+
+    Parameters
+    ==========
+
+    base : scalar or matrix
+
+    exp : scalar or matrix
+
+    Notes
+    =====
+
+    There are four definitions for the hadamard power which can be used.
+    Let's consider `A, B` as `(m, n)` matrices, and `a, b` as scalars.
+
+    Matrix raised to a scalar exponent:
+
+    .. math::
+        A^{\circ b} = \begin{bmatrix}
+        A_{0, 0}^b   & A_{0, 1}^b   & \cdots & A_{0, n-1}^b   \\
+        A_{1, 0}^b   & A_{1, 1}^b   & \cdots & A_{1, n-1}^b   \\
+        \vdots       & \vdots       & \ddots & \vdots         \\
+        A_{m-1, 0}^b & A_{m-1, 1}^b & \cdots & A_{m-1, n-1}^b
+        \end{bmatrix}
+
+    Scalar raised to a matrix exponent:
+
+    .. math::
+        a^{\circ B} = \begin{bmatrix}
+        a^{B_{0, 0}}   & a^{B_{0, 1}}   & \cdots & a^{B_{0, n-1}}   \\
+        a^{B_{1, 0}}   & a^{B_{1, 1}}   & \cdots & a^{B_{1, n-1}}   \\
+        \vdots         & \vdots         & \ddots & \vdots           \\
+        a^{B_{m-1, 0}} & a^{B_{m-1, 1}} & \cdots & a^{B_{m-1, n-1}}
+        \end{bmatrix}
+
+    Matrix raised to a matrix exponent:
+
+    .. math::
+        A^{\circ B} = \begin{bmatrix}
+        A_{0, 0}^{B_{0, 0}}     & A_{0, 1}^{B_{0, 1}}     &
+        \cdots & A_{0, n-1}^{B_{0, n-1}}     \\
+        A_{1, 0}^{B_{1, 0}}     & A_{1, 1}^{B_{1, 1}}     &
+        \cdots & A_{1, n-1}^{B_{1, n-1}}     \\
+        \vdots                  & \vdots                  &
+        \ddots & \vdots                      \\
+        A_{m-1, 0}^{B_{m-1, 0}} & A_{m-1, 1}^{B_{m-1, 1}} &
+        \cdots & A_{m-1, n-1}^{B_{m-1, n-1}}
+        \end{bmatrix}
+
+    Scalar raised to a scalar exponent:
+
+    .. math::
+        a^{\circ b} = a^b
     """
 
     def __new__(cls, base, exp):
         base = sympify(base)
         exp = sympify(exp)
+        if base.is_scalar and exp.is_scalar:
+            return base ** exp
+
+        if base.is_Matrix and exp.is_Matrix and base.shape != exp.shape:
+            raise ValueError(
+                'The shape of the base {} and '
+                'the shape of the exponent {} do not match.'
+                .format(base.shape, exp.shape)
+                )
+
         obj = super(HadamardPower, cls).__new__(cls, base, exp)
         return obj
 
@@ -314,10 +376,31 @@ def exp(self):
 
     @property
     def shape(self):
-        return self.base.shape
+        if self.base.is_Matrix:
+            return self.base.shape
+        return self.exp.shape
 
     def _entry(self, i, j, **kwargs):
-        return self.base._entry(i, j, **kwargs)**self.exp
+        base = self.base
+        exp = self.exp
+
+        if base.is_Matrix:
+            a = base._entry(i, j, **kwargs)
+        elif base.is_scalar:
+            a = base
+        else:
+            raise ValueError(
+                'The base {} must be a scalar or a matrix.'.format(base))
+
+        if exp.is_Matrix:
+            b = exp._entry(i, j, **kwargs)
+        elif exp.is_scalar:
+            b = exp
+        else:
+            raise ValueError(
+                'The exponent {} must be a scalar or a matrix.'.format(exp))
+
+        return a ** b
 
     def _eval_transpose(self):
         from sympy.matrices.expressions.transpose import transpose


Test patch:
diff --git a/sympy/matrices/expressions/tests/test_hadamard.py b/sympy/matrices/expressions/tests/test_hadamard.py
index 2f2b368882..08ca8ce46b 100644
--- a/sympy/matrices/expressions/tests/test_hadamard.py
+++ b/sympy/matrices/expressions/tests/test_hadamard.py
@@ -11,6 +11,7 @@
 B = MatrixSymbol('B', n, m)
 C = MatrixSymbol('C', m, k)
 
+
 def test_HadamardProduct():
     assert HadamardProduct(A, B, A).shape == A.shape
 
@@ -26,9 +27,11 @@ def test_HadamardProduct():
 
     assert set(HadamardProduct(A, B, A).T.args) == set((A.T, A.T, B.T))
 
+
 def test_HadamardProduct_isnt_commutative():
     assert HadamardProduct(A, B) != HadamardProduct(B, A)
 
+
 def test_mixed_indexing():
     X = MatrixSymbol('X', 2, 2)
     Y = MatrixSymbol('Y', 2, 2)
@@ -53,6 +56,7 @@ def test_canonicalize():
     assert HadamardProduct(X, Z, U, Y).doit() == Z
 
 
+
 def test_hadamard():
     m, n, p = symbols('m, n, p', integer=True)
     A = MatrixSymbol('A', m, n)
@@ -82,3 +86,27 @@ def test_hadamard_power():
     assert hadamard_power(A, n)[0, 0] == A[0, 0]**n
     assert hadamard_power(m, n) == m**n
     raises(ValueError, lambda: hadamard_power(A, A))
+
+
+def test_hadamard_power_explicit():
+    from sympy.matrices import Matrix
+    A = MatrixSymbol('A', 2, 2)
+    B = MatrixSymbol('B', 2, 2)
+    a, b = symbols('a b')
+
+    assert HadamardPower(a, b) == a**b
+
+    assert HadamardPower(a, B).as_explicit() == \
+        Matrix([
+            [a**B[0, 0], a**B[0, 1]],
+            [a**B[1, 0], a**B[1, 1]]])
+
+    assert HadamardPower(A, b).as_explicit() == \
+        Matrix([
+            [A[0, 0]**b, A[0, 1]**b],
+            [A[1, 0]**b, A[1, 1]**b]])
+
+    assert HadamardPower(A, B).as_explicit() == \
+        Matrix([
+            [A[0, 0]**B[0, 0], A[0, 1]**B[0, 1]],
+            [A[1, 0]**B[1, 0], A[1, 1]**B[1, 1]]])