In R2016 (not sure if it was a or b. I think b, but not really relevant, so will not check.) they introduced the capability to expand singleton dimensions for operations like this.
So now things like this work:
M = magic(3)
M =
8 1 6
3 5 7
4 9 2
M - mean(M)
ans =
3 -4 1
-2 0 2
-1 4 -3
Which is as you describe. Prior to the change, you needed to use bsxfun to do that, and before bsxfun existed, you had to expand the vector yourself, using a tool like repmat.
