title: Theano gradient with function on tensors
body: <p>I have a function that calculates a value of a scalar field on a 3D space, so I feed it 3D tensors for x, y and z coordinates (obtained by numpy.meshgrid) and use elementwise operations everywhere. This works as expected.</p>
  
  <p>Now I need to calculate a gradient of the scalar field. I've been playing around with <code>theano.tensor.grad</code> and <code>theano.tensor.jacobian</code> and I don't understand how a derivative of elementwise operation is supposed to work.</p>
  
  <p>This is a MWE that I don't understand:</p>
  
  <pre><code>import theano.tensor as T 
  
  x, y = T.matrices("xy")
  
  expr = x**2 + y
  grad = T.grad(expr[0, 0], x)
  print(grad.eval({x: [[1, 2], [1, 2]], y: [[1, 1], [2, 2]]}))
  </code></pre>
  
  <p>It prints</p>
  
  <pre><code>[[ 2.  0.]
   [ 0.  0.]]
  </code></pre>
  
  <p>while I would expect </p>
  
  <pre><code>[[ 2.  4.]
   [ 2.  4.]]
  </code></pre>
  
  <p>I also tried with jacobian:</p>
  
  <pre><code>import theano.tensor as T
  
  x, y = T.matrices("xy")
  
  expr = x**2 + y
  grad = T.jacobian(expr.flatten(), x)
  print(grad.eval({x: [[1, 2], [1, 2]], y: [[1, 1], [2, 2]]}))
  </code></pre>
  
  <p>which returns</p>
  
  <pre><code>[[[ 2.  0.]
    [ 0.  0.]]
  
   [[ 0.  4.]
    [ 0.  0.]]
  
   [[ 0.  0.]
    [ 2.  0.]]
  
   [[ 0.  0.]
    [ 0.  4.]]]
  </code></pre>
  
  <p>(the nonzero elements together would give me my expected matrix from the previous example)</p>
  
  <p>Is there some way to get the elmentwise gradients I need?</p>
  
  <p>Can I for example somehow define the function as scalar (three scalars into a scalar) apply it elementwise over the coordinate tensors? This way the derivative would also be just a simple scalar and everything would work smoothly.</p>