Repository: sympy/sympy
Pull request: Fixes diff and Derivative different results in matix expressions
<!-- 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://tinyurl.com/auto-closing for more information). Also, please
write a comment on that issue linking back to this pull request once it is
open. -->
Fixes #24859 
#### Brief description of what is fixed or changed


#### Other comments


#### Release Notes

<!-- Write the release notes for this release below between the BEGIN and END
statements. The basic format is a bulleted list with the name of the subpackage
and the release note for this PR. For example:

* solvers
  * Added a new solver for logarithmic equations.

* functions
  * Fixed a bug with log of integers. Formerly, `log(-x)` incorrectly gave `-log(x)`.

* physics.units
  * Corrected a semantical error in the conversion between volt and statvolt which
    reported the volt as being larger than the statvolt.

or if no release note(s) should be included use:

NO ENTRY

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 -->
NO ENTRY
<!-- END RELEASE NOTES -->


Problem statement:
diff and Derivative return different results
Here's a short code snippet to create this different results:

```
import sympy as sp
A=sp.MatrixSymbol('A', 3, 4)
B=sp.MatrixSymbol('B', 4, 3)
J=A*B
Jinv=sp.Matrix(J).adjugate()#/sp.Matrix(J).det()

alpha=sp.symbols('alpha')
u=sp.MatrixSymbol('u', 3, 4)
x=sp.symbols('x')
Jk=Jinv.subs(A,A+x*u)

print('Result from diff')
sp.pprint(Jk[0,0].diff(x))

print('Result from Derivative')
sp.Derivative(Jk[0,0], x).doit()
```

I'm anticipating both would return the same result. Is this intended?

Hints from issue discussion:
This is the output:
```python
In [3]: Jk[0, 0]
Out[3]: 
(x⋅u + A)[1, 0]⋅(x⋅u + A)[2, 1]⋅B₀₁⋅B₁₂ - (x⋅u + A)[1, 0]⋅(x⋅u + A)[2, 1]⋅B₀₂⋅B₁₁ + (x⋅u + A)[1, 0]
⋅(x⋅u + A)[2, 2]⋅B₀₁⋅B₂₂ - (x⋅u + A)[1, 0]⋅(x⋅u + A)[2, 2]⋅B₀₂⋅B₂₁ + (x⋅u + A)[1, 0]⋅(x⋅u + A)[2, 3
]⋅B₀₁⋅B₃₂ - (x⋅u + A)[1, 0]⋅(x⋅u + A)[2, 3]⋅B₀₂⋅B₃₁ - (x⋅u + A)[1, 1]⋅(x⋅u + A)[2, 0]⋅B₀₁⋅B₁₂ + (x⋅
u + A)[1, 1]⋅(x⋅u + A)[2, 0]⋅B₀₂⋅B₁₁ + (x⋅u + A)[1, 1]⋅(x⋅u + A)[2, 2]⋅B₁₁⋅B₂₂ - (x⋅u + A)[1, 1]⋅(x
⋅u + A)[2, 2]⋅B₁₂⋅B₂₁ + (x⋅u + A)[1, 1]⋅(x⋅u + A)[2, 3]⋅B₁₁⋅B₃₂ - (x⋅u + A)[1, 1]⋅(x⋅u + A)[2, 3]⋅B
₁₂⋅B₃₁ - (x⋅u + A)[1, 2]⋅(x⋅u + A)[2, 0]⋅B₀₁⋅B₂₂ + (x⋅u + A)[1, 2]⋅(x⋅u + A)[2, 0]⋅B₀₂⋅B₂₁ - (x⋅u +
 A)[1, 2]⋅(x⋅u + A)[2, 1]⋅B₁₁⋅B₂₂ + (x⋅u + A)[1, 2]⋅(x⋅u + A)[2, 1]⋅B₁₂⋅B₂₁ + (x⋅u + A)[1, 2]⋅(x⋅u 
+ A)[2, 3]⋅B₂₁⋅B₃₂ - (x⋅u + A)[1, 2]⋅(x⋅u + A)[2, 3]⋅B₂₂⋅B₃₁ - (x⋅u + A)[1, 3]⋅(x⋅u + A)[2, 0]⋅B₀₁⋅
B₃₂ + (x⋅u + A)[1, 3]⋅(x⋅u + A)[2, 0]⋅B₀₂⋅B₃₁ - (x⋅u + A)[1, 3]⋅(x⋅u + A)[2, 1]⋅B₁₁⋅B₃₂ + (x⋅u + A)
[1, 3]⋅(x⋅u + A)[2, 1]⋅B₁₂⋅B₃₁ - (x⋅u + A)[1, 3]⋅(x⋅u + A)[2, 2]⋅B₂₁⋅B₃₂ + (x⋅u + A)[1, 3]⋅(x⋅u + A
)[2, 2]⋅B₂₂⋅B₃₁

In [4]: Jk[0, 0].diff(x)
Out[4]: 0

In [5]: diff(Jk[0, 0], x)
Out[5]: 0

In [6]: Derivative(Jk[0, 0], x).doit()
Out[6]: 
(x⋅u₁₀ + A₁₀)⋅B₀₁⋅B₁₂⋅u₂₁ + (x⋅u₁₀ + A₁₀)⋅B₀₁⋅B₂₂⋅u₂₂ + (x⋅u₁₀ + A₁₀)⋅B₀₁⋅B₃₂⋅u₂₃ - (x⋅u₁₀ + A₁₀)⋅B
₀₂⋅B₁₁⋅u₂₁ - (x⋅u₁₀ + A₁₀)⋅B₀₂⋅B₂₁⋅u₂₂ - (x⋅u₁₀ + A₁₀)⋅B₀₂⋅B₃₁⋅u₂₃ - (x⋅u₁₁ + A₁₁)⋅B₀₁⋅B₁₂⋅u₂₀ + (x
⋅u₁₁ + A₁₁)⋅B₀₂⋅B₁₁⋅u₂₀ + (x⋅u₁₁ + A₁₁)⋅B₁₁⋅B₂₂⋅u₂₂ + (x⋅u₁₁ + A₁₁)⋅B₁₁⋅B₃₂⋅u₂₃ - (x⋅u₁₁ + A₁₁)⋅B₁₂
⋅B₂₁⋅u₂₂ - (x⋅u₁₁ + A₁₁)⋅B₁₂⋅B₃₁⋅u₂₃ - (x⋅u₁₂ + A₁₂)⋅B₀₁⋅B₂₂⋅u₂₀ + (x⋅u₁₂ + A₁₂)⋅B₀₂⋅B₂₁⋅u₂₀ - (x⋅u
₁₂ + A₁₂)⋅B₁₁⋅B₂₂⋅u₂₁ + (x⋅u₁₂ + A₁₂)⋅B₁₂⋅B₂₁⋅u₂₁ + (x⋅u₁₂ + A₁₂)⋅B₂₁⋅B₃₂⋅u₂₃ - (x⋅u₁₂ + A₁₂)⋅B₂₂⋅B
₃₁⋅u₂₃ - (x⋅u₁₃ + A₁₃)⋅B₀₁⋅B₃₂⋅u₂₀ + (x⋅u₁₃ + A₁₃)⋅B₀₂⋅B₃₁⋅u₂₀ - (x⋅u₁₃ + A₁₃)⋅B₁₁⋅B₃₂⋅u₂₁ + (x⋅u₁₃
 + A₁₃)⋅B₁₂⋅B₃₁⋅u₂₁ - (x⋅u₁₃ + A₁₃)⋅B₂₁⋅B₃₂⋅u₂₂ + (x⋅u₁₃ + A₁₃)⋅B₂₂⋅B₃₁⋅u₂₂ - (x⋅u₂₀ + A₂₀)⋅B₀₁⋅B₁₂
⋅u₁₁ - (x⋅u₂₀ + A₂₀)⋅B₀₁⋅B₂₂⋅u₁₂ - (x⋅u₂₀ + A₂₀)⋅B₀₁⋅B₃₂⋅u₁₃ + (x⋅u₂₀ + A₂₀)⋅B₀₂⋅B₁₁⋅u₁₁ + (x⋅u₂₀ +
 A₂₀)⋅B₀₂⋅B₂₁⋅u₁₂ + (x⋅u₂₀ + A₂₀)⋅B₀₂⋅B₃₁⋅u₁₃ + (x⋅u₂₁ + A₂₁)⋅B₀₁⋅B₁₂⋅u₁₀ - (x⋅u₂₁ + A₂₁)⋅B₀₂⋅B₁₁⋅u
₁₀ - (x⋅u₂₁ + A₂₁)⋅B₁₁⋅B₂₂⋅u₁₂ - (x⋅u₂₁ + A₂₁)⋅B₁₁⋅B₃₂⋅u₁₃ + (x⋅u₂₁ + A₂₁)⋅B₁₂⋅B₂₁⋅u₁₂ + (x⋅u₂₁ + A
₂₁)⋅B₁₂⋅B₃₁⋅u₁₃ + (x⋅u₂₂ + A₂₂)⋅B₀₁⋅B₂₂⋅u₁₀ - (x⋅u₂₂ + A₂₂)⋅B₀₂⋅B₂₁⋅u₁₀ + (x⋅u₂₂ + A₂₂)⋅B₁₁⋅B₂₂⋅u₁₁
 - (x⋅u₂₂ + A₂₂)⋅B₁₂⋅B₂₁⋅u₁₁ - (x⋅u₂₂ + A₂₂)⋅B₂₁⋅B₃₂⋅u₁₃ + (x⋅u₂₂ + A₂₂)⋅B₂₂⋅B₃₁⋅u₁₃ + (x⋅u₂₃ + A₂₃
)⋅B₀₁⋅B₃₂⋅u₁₀ - (x⋅u₂₃ + A₂₃)⋅B₀₂⋅B₃₁⋅u₁₀ + (x⋅u₂₃ + A₂₃)⋅B₁₁⋅B₃₂⋅u₁₁ - (x⋅u₂₃ + A₂₃)⋅B₁₂⋅B₃₁⋅u₁₁ +
 (x⋅u₂₃ + A₂₃)⋅B₂₁⋅B₃₂⋅u₁₂ - (x⋅u₂₃ + A₂₃)⋅B₂₂⋅B₃₁⋅u₁₂
```

> I'm anticipating both would return the same result. Is this intended?

No it is not intended. They should give the same result.
The problem is here:
https://github.com/sympy/sympy/blob/940f1af711d0148d44f6418c3d92ec781f701341/sympy/matrices/expressions/matexpr.py#L635-L639
In the debugger:
```python
635  	        if not isinstance(v, MatrixElement):
636  	            from sympy.matrices.matrices import MatrixBase
637  	            if isinstance(self.parent, MatrixBase):
638  	                return self.parent.diff(v)[self.i, self.j]
639  ->	            return S.Zero
640  	
641  	        M = self.args[0]
642  	
643  	        m, n = self.parent.shape
644  	
(Pdb) p self
(x*u + A)[1, 0]
(Pdb) p self.parent
x*u + A
(Pdb) p type(self.parent)
<class 'sympy.matrices.expressions.matadd.MatAdd'>
(Pdb) p isinstance(self.parent, MatrixBase)
False
(Pdb) p self.parent.diff(v)
u + 0
(Pdb) p self.parent.diff(v)[self.i, self.j]
u[1, 0]
```
I'm not sure why the check is for `MatrixBase` which is the superclass for explicit matrix types (as opposed to matrix expressions). It should be guaranteed that the parent is some kind of matrix.

The difference with doit is that it is called recursively to change other parts of the expression before computing the derivative:
```python
In [3]: from sympy.matrices.expressions.matexpr import MatrixElement

In [4]: MatrixElement((x*u + A), 1, 0)
Out[4]: (x⋅u + A)[1, 0]

In [5]: MatrixElement((x*u + A), 1, 0).doit()
Out[5]: x⋅u₁₀ + A₁₀

In [6]: (x*u + A)[1, 0]
Out[6]: x⋅u₁₀ + A₁₀

In [7]: MatrixElement((x*u + A), 1, 0).diff(x)
Out[7]: 0

In [8]: MatrixElement((x*u + A), 1, 0).doit().diff(x)
Out[8]: u₁₀

In [9]: (x*u + A)[1, 0].diff(x)
Out[9]: u₁₀
```

As far as I can tell this fixes it:
```diff
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index a433c8b..75ee617 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -633,10 +633,7 @@ def indices(self):
     def _eval_derivative(self, v):
 
         if not isinstance(v, MatrixElement):
-            from sympy.matrices.matrices import MatrixBase
-            if isinstance(self.parent, MatrixBase):
-                return self.parent.diff(v)[self.i, self.j]
-            return S.Zero
+            return self.parent.diff(v)[self.i, self.j]
 
         M = self.args[0]
```

CC @Upabjojr the relevant code here is all from your commits. I don't understand why this is checking if the parent is MatrixBase. It is clearly some kind of matrix.

Commit message: added simple tests

Patch to review:
diff --git a/sympy/matrices/expressions/matexpr.py b/sympy/matrices/expressions/matexpr.py
index a433c8bbb49d..75ee6177c0f6 100644
--- a/sympy/matrices/expressions/matexpr.py
+++ b/sympy/matrices/expressions/matexpr.py
@@ -633,10 +633,7 @@ def indices(self):
     def _eval_derivative(self, v):
 
         if not isinstance(v, MatrixElement):
-            from sympy.matrices.matrices import MatrixBase
-            if isinstance(self.parent, MatrixBase):
-                return self.parent.diff(v)[self.i, self.j]
-            return S.Zero
+            return self.parent.diff(v)[self.i, self.j]
 
         M = self.args[0]
 
diff --git a/sympy/matrices/expressions/tests/test_matexpr.py b/sympy/matrices/expressions/tests/test_matexpr.py
index 36b3846c2ec5..53604c6c74ea 100644
--- a/sympy/matrices/expressions/tests/test_matexpr.py
+++ b/sympy/matrices/expressions/tests/test_matexpr.py
@@ -404,6 +404,18 @@ def test_issue_21195():
     assert A.diff(x) == B
 
 
+def test_issue_24859():
+    A = MatrixSymbol('A', 2, 3)
+    B = MatrixSymbol('B', 3, 2)
+    J = A*B
+    Jinv = Matrix(J).adjugate()
+    u = MatrixSymbol('u', 2, 3)
+    Jk = Jinv.subs(A, A + x*u)
+
+    assert Jk[0,0].diff(x).simplify() == B[0, 1]*u[1, 0] + B[1, 1]*u[1, 1] + B[2, 1]*u[1, 2]
+    assert diff(Jk[0,0], x).doit().simplify() == B[0, 1]*u[1, 0] + B[1, 1]*u[1, 1] + B[2, 1]*u[1, 2]
+
+
 def test_MatMul_postprocessor():
     z = zeros(2)
     z1 = ZeroMatrix(2, 2)


Code review comments:
diff_hunk: @@ -404,6 +404,18 @@ def test_issue_21195():
     assert A.diff(x) == B
 
 
+def test_issue_24859():
+    A = MatrixSymbol('A', 2, 3)
+    B = MatrixSymbol('B', 3, 2)
+    J = A*B
+    Jinv = Matrix(J).adjugate()
+    u = MatrixSymbol('u', 2, 3)
+    Jk = Jinv.subs(A, A + x*u)
+
+    assert Jk[0,0].diff(x).simplify() == B[0, 1]*u[1, 0] + B[1, 1]*u[1, 1] + B[2, 1]*u[1, 2]
+    assert diff(Jk[0,0], x).doit().simplify() == B[0, 1]*u[1, 0] + B[1, 1]*u[1, 1] + B[2, 1]*u[1, 2]
path: sympy/matrices/expressions/tests/test_matexpr.py
text: @user1:
Three things:

1. These lines are quite long. Can you refactor the RHSs of the equality checks to an assigned variable that can be reused, something like `expected = B[0, 1]*u[1, 0] + B[1, 1]*u[1, 1] + B[2, 1]*u[1, 2]`.
2. Can we please format the matrix indexing consistently, i.e. change the `Jk[0,0]` to `Jk[0, 0]`.
3. Calls to `simplify` should be avoided if possible because it uses a heuristic approach. If you know the type of simplification needed to achieve the desired result then use one of the [targeted simplification functions](https://docs.sympy.org/latest/tutorials/intro-tutorial/simplification.html). Running locally it looks like the `simplify` isn't needed.

```suggestion
    expected = B[0, 1]*u[1, 0] + B[1, 1]*u[1, 1] + B[2, 1]*u[1, 2]
    assert Jk[0, 0].diff(x) == expected
    assert diff(Jk[0, 0], x).doit() == expected
```

@user1:
Addressed by https://github.com/sympy/sympy/pull/24867/commits/0fd5f03d09d8be1728da7c3c815cdd39350a07e1.