Fix bug in gradient of Elemwise containing multi-output scalars

Author: ricardoV94

pytensor/tensor/elemwise.py

@@ -638,6 +638,9 @@ def transform(r):
                 return DimShuffle((), ["x"] * nd)(res)
 
             new_r = Elemwise(node.op, {})(*[transform(ipt) for ipt in node.inputs])
+            if isinstance(new_r, (list, tuple)):
+                # Scalar Op with multiple outputs
+                new_r = new_r[r.owner.outputs.index(r)]
             return new_r
 
         ret = []