Why do Deep Learning libraries force the cost function to output a scalar? Let's say we have a neural net with:


*

*5 input neurons

*some arbitrary amount of hidden layers

*3 output neurons


Let's say we're using minibatches of size 32. So,


*

*if we input a 5x32 matrix into the neural net,

*we will then get out a 3x32 matrix of output activations. 


Assume we are using a simple MSE loss function. We take the difference of our 3x32 target matrix and 3x32 output matrix, and do elementwise squaring of each entry in the resulting new 3x32 matrix, which we'll call M.
What I'm confused about, is that I've seen lots of code in Theano and other deep learning libraries that takes the mean of this matrix for its cost function (i.e. outputs a scalar--not a vector). 
For example,     T.mean(T.pow(T-Y, 2))    . Why is this? In my example, shouldn't Theano backprop M, not a scalar value? Even when I'm using minibatches, are these libraries just backpropping scalars?
 A: All machine learning is about minimizing cost of some model. The most elementary thing when you try to find minimum value is ability to compare two values. You can do it only with scalar values. For example, given two vectors [0,2], [2,2] how would you compare those tuples? You have to define some norm function. Euclidean, max, Manhattan, or your own fancy one. Whatever you use it must output scalar values.
A: 
  
*
  
*3 output neurons
  

The loss function in most applications is chosen such that it calculates a combined loss for these three neurons (e.g. cross entropy loss). This defines the tradeoff between better matching of the target value of neuron 1 at the expense of worse matching of their respective target values of the other two neurons. You can of course define a different loss function, e.g. one that reflects that matching the target value is twice as important for neuron 1 than for the other two output neurons etc.

For example, T.mean(T.pow(T-Y, 2)) . Why is this? 

This is the average error over your (training) sample. Usually all samples of the training sample are treated equal so matching the network outputs to the target values for sample 1 is as important as matching the target values for sample 2 etc. There are situations where some elements of the training sample are treated as more important than others, this can be modeled by adding weight terms in the loss function.

In my example, shouldn't Theano backprop M, not a scalar value?

This would correspond to doing backpropagation M times for a single sample where the network weights are adjusted M times, potentially in opposite directions at each of the M backpropagation operations. Plain stochastic gradient descent works like this for example.
What you usually want to minimize (i.e. the original loss function) is the average over the entire training sample. You would get one weight update after calculating the loss over the entire training sample. 
It turns out that in practice you can get faster close to the optimum by using stochastic gradient descent. On the other hand, modern computing hardware (CPUs or GPUs) is vectorized (i.e. support running the same code on multiple values) so minibatching allows to take into account more than one sample per upgrade at similar execution speed. 

when I'm using minibatches, are these libraries just backpropping
  scalars?

They backpropagate the derivatives with respect to the average loss over the entire minibatch (which is a scalar value). 
If you have a network with three output neurons contributing to a single loss function, each neuron gets its share of the correction according how badly it contributed to the loss, averaged over the M rows in the minibatch.
A: You ask in a comment:

Let's say we backprop the matrix. If we do that we will get 32 gradient updates for each parameter in the neural net. For each parameter, can't we just take the average of the 32 gradients and use that in gradient descent?

Well, suppose you have 1,000,000 parameters. You're suggesting that we calculate 32,000,000 partial derivatives and then average them, 32 at a time, in order to get 1,000,000 partial derivatives that we can then apply to the parameters.
If the cost is a scalar, on the other hand, then we only need to calculate 1,000,000 partial derivatives in the first place, and then we can apply them to the parameters immediately, without needing to do any averaging first.
So, you're asking if we can't "just" do 32 times as much work. And the answer is yes, we can... but it's 32 times as much work.

Isn't doing that different than just taking a mean of the output activations and backpropping the scalar?

No, it's the same. Suppose you have several loss functions $L_1$, $L_2$, $L_3$, $L_4$ that you want to optimize by changing some parameter $p$. Then the average of the derivatives,
$$\frac14 \left (\frac{\partial L_1}{\partial p} + \frac{\partial L_2}{\partial p} + \frac{\partial L_3}{\partial p} + \frac{\partial L_4}{\partial p} \right),$$
is simply the derivative of the average,
$$\frac14 \cdot \frac{\partial}{\partial p} (L_1 + L_2 + L_3 + L_4).$$
So, taking the average first and then doing backpropagation is equivalent to first doing backpropagation and then taking the average. It's also much faster.
