Gradient of a Theano custom operator returning a matrix

In your code above, you dont need to wrap it in a custom operator, as your operation is reshaping the 3 scaler into a matrix, the easiest way to do is something like:

with pm.Model() as m:
    c11 = ...
    c12 = ...
    c22 = ...
    x = tt.zeros((2, 2))
    x = tt.inc_subtensor(x[0, 0], c11)
    x = tt.inc_subtensor(x[0, 1], c12)
    x = tt.inc_subtensor(x[1, 0], c12)
    x = tt.inc_subtensor(x[1, 1], c22)

And theano will take care of the gradient for each scaler.

I will reply to you in Multivariate Gaussian with Custom Operator Mean shortly.