This works! However the main problem is a bit more tricky. In the feature_test function, I am using N in np.tile, so I would need N.eval(). I still thus get problems when performing output = output*pm.math.matmul(feature_test(X[:,d],self.N,self.M),factors[d]).
I could avoid this whole mess if I could write the whole feature_test in pytensor without ever calling .eval()?