Multinomial Logistic Loss vs (Cross Entropy vs Square Error) I observed that Caffe (a deep learning framework) used the Softmax Loss Layer SoftmaxWithLoss as output layer for most of the model samples.
As far as I know, Softmax Loss layer is the combination of Multinomial Logistic Loss Layer and Softmax Layer. 
From Caffe, they said that

Softmax Loss Layer gradient computation is more numerically stable

However, this explanation is not the answer that I want, the explanation is just compare the combination of Multinomial Logistic Loss Layer and Softmax Loss layer instead of layer by layer. But not compare with other type of loss function.
However, I would like to know more what is the differences/advantages/disadvantages of these 3 error function which is Multinomial Logistic Loss, Cross Entropy (CE) and Square Error (SE) in supervised learning perspective? Any supportive articles?
 A: SHORT ANSWER
According to other answers Multinomial Logistic Loss and Cross Entropy Loss are the same.
Cross Entropy Loss is an alternative cost function for NN with sigmoids activation function introduced artificially to eliminate the dependency on $\sigma'$ on the update equations. Some times this term slows down the learning process. Alternative methods are regularised cost function.
In these type of networks one might want to have probabilities as output but this does not happen with the sigmoids in a multinomial network. The softmax function normalizes the outputs and force them in the range $[0,1]$. This can be useful for example in MNIST classification.
LONG ANSWER WITH SOME INSIGHTS
The answer is quite long but I'll try to summarise.
The first modern artificial neurons that have been used are the sigmoids whose function is:
$$\sigma(x) = \frac{1}{1+e^{-x}}$$
which has the following shape:

The curve is nice because it guarantees the output is in the range $[0,1]$.
Regarding the choice of a cost function, a natural choice is the quadratic cost function, whose derivative is guaranteed to exist and we know it has a minimum.
Now consider a NN with sigmoids trained with the quadratic cost function, with $L$ layers.
We define the cost function as the sum of the squared errors in the output layer for a set of inputs $X$:
$$C = \frac{1}{2N}\sum_x^N\sum_{j=1}^K (y_j(x) - a_j^L(x))^2$$
where $a_j^L$ is the j-th neuron in the output layer $L$, $y_j$ the desired output and $N$ is the number of training examples.
For simplicity let's consider the error for a single input:
$$C = \sum_{j=1}^K (y_j(x) - a_j^L(x))^2$$
Now an activation output of for the $j$ neuron in the $\ell$ layer, $a_j^\ell$ is:
$$a_j^\ell = \sum_k w_{jk}^\ell \cdot a_j^{\ell-1}+b_j^\ell =  \mathbf{w}_{j}^\ell \cdot \mathbf{a}_j^{\ell-1}+b_j^\ell$$
Most of the times (if not always) NN are trained with one of the gradient descent techniques, which basically consists updating the weights $w$ and biases $b$ by small steps towards the direction of minimization. The goal is to a apply a small change in the weights and biases towards the direction that minimizes the cost function.
For small steps the following holds:
$$\Delta C \approx \frac{\partial C}{\partial v_i}\Delta v_i$$
Our $v_i$ are the weights and biases.
Being it a cost function we want to minimise, i.e., find the proper value $\Delta v_i$. Suppose we choose $$\Delta v_i = -\eta \frac{\partial C}{\partial v_i}$$, then:
$$\Delta C \approx -\eta \left(\frac{\partial C}{\partial v_i}\right)$$
which means the change $\Delta v_i$ in the parameter decreased the cost function by $\Delta C$.
Consider the $j$-th output neuron:
$$C = \frac{1}{2}(y(x)-a_j^L(x)^2$$
$$a_j^L =\sigma = \frac{1}{ 1+e^{ -(\mathbf{w}_j^\ell \cdot \mathbf{a}_j^{\ell-1}+b_j^\ell)}}$$
Suppose we want to update the weight $w_{jk}^\ell$ which is the weight from the neuron $k$ in the $\ell-1$ layer to the $j$-th neuron in the \ell layer. Then we have:
$$w_{jk}^\ell \Rightarrow w_{jk}^\ell -\eta \frac{\partial C}{\partial w_{jk}^\ell}$$
$$b_{j}^\ell \Rightarrow b_{j}^\ell -\eta \frac{\partial C}{\partial b_{j}^\ell}$$
Taking the derivatives using the chain rule:
$$ \frac{\partial C}{\partial w_{jk}^\ell} = \left(a_j^L(x)-y(x)\right) \sigma' a_k^{\ell-1}$$
$$ \frac{\partial C}{\partial b_{j}^\ell} = \left(a_j^L(x)-y(x)\right) \sigma'$$
You see the dependency on the derivative of the sigmoid (in the first is w.r.t. $w$ in the second w.r.t $b$ actually, but it does not change a lot since both are exponents).
Now the derivative for a generic single variable sigmoid $z$ is:
$$\frac{d \sigma(z)}{d z} = \sigma(z)(1-\sigma(z))$$
Now consider a single output neuron and suppose you that neuron should output $0$ instead it is outputting a value close to $1$: you'll see both from the graph that the sigmoid for values close to $1$ is flat, i.e. its derivative is close to $0$, i.e. updates of the parameter are very slow (since the update equations depend on $\sigma'$.
Motivation of the cross-entropy function
To see how cross-entropy has been originally derived, suppose one has just found out that the term $\sigma'$ is slowing down the learning process. We might wonder if it is possible to choose a cost function to make the term  $\sigma'$ disappear. Basically one might want:
\begin{equation}
\begin{aligned}
\frac{\partial C}{\partial w} & =x \left( a - y\right)\\
\frac{\partial C}{\partial b} =\left( a - y\right)
\end{aligned}
\end{equation}
From the chain-rule we have:
\begin{equation}
\frac{\partial C}{\partial b} =\frac{\partial C}{\partial a} \frac{\partial a}{\partial b }  =\frac{\partial C}{\partial a}\sigma'(z) = \frac{\partial C}{\partial a} \sigma(1-\sigma)
\end{equation}
Comparing the desired equation with the one of the chain rule, one gets
\begin{equation}
\frac{\partial C}{\partial a} = \frac{a-y}{a(1-a)} 
\end{equation}
Using the cover-up method :
\begin{equation}
\frac{\partial C}{\partial a} = -\left[ y\ln a + (1-y)\ln(1-a)\right]+const
\end{equation}
To get the full cost function, we must average over all the training samples
\begin{equation}
\frac{\partial C}{\partial a} = -\frac{1}{n}\sum_x\left[y\ln a + (1-y)\ln(1-a)\right]+const
\end{equation}
where the constant here is the average of the individual constants for each training example.
There is a standard way of interpreting the cross-entropy that comes from the field of information theory. Roughly speaking, the idea is that the cross-entropy is a measure of surprise. We get low surprise if the output $a$ is what we expect ($y$), and high surprise if the output is unexpected.
Softmax
For a binary classification cross-entropy resembles the definition in information theory and the values can still be interpreted as probabilities.
With multinomial classification this does not hold true anymore: the outputs do note sum up to $1$.
If you want them to sum up to $1$ you use softmax function, which normalize the outputs so that the sum is $1$.
Also if the output layer is made up of softmax functions, the slowing down term is not present. If you use log-likelihood cost function with a softmax output layer, the result you will obtain a form of the partial derivatives, and in turn of the update equations, similar to the one found for a cross-entropy function with sigmoid neurons
However
A: 
I would like to know more what is the
  differences/advantages/disadvantages of these 3 error function which
  is Multinomial Logistic Loss, Cross Entropy (CE) and Square Error (SE)
  in supervised learning perspective?

The multinomial logistic loss is acturally the same as cross entropy. Look at this function(the cost function in softmax): 
$$ J( \theta ) = - \sum^m_{i=1} \sum^K_{k=1} 1 \{ y^{(i)} = k \} \log p(y^{(i)} =  k \mid x^{(i)} ; \theta) $$
where m is the sample number, K is the class number. 
The indicator function($1 \{ y^{(i)} = k \}$) determines whether the $p(x)$ bellow is 0 or 1 in the cross entropy definition, which is labelled as one hot in the training data, and $ p(y^{(i)} =  k \mid x^{(i)} ; \theta) $ is the conditional likelihood of the softmax(q(x) as shown bellow). 
$$-\sum_x p(x) \log q(x)$$
And MSE is mostly for the situation where the link function is the unity function(the response distribution follows a normal distribution), the standard linear regression, while the cross entropy is normally for where the link function is the logit function. Here is an awesome comparison you can refer to. 

Any supportive articles?

Except those in the links, recommend you this illustrating one: https://github.com/rasbt/python-machine-learning-book/blob/master/faq/softmax_regression.md
A: In my opinion, loss function is the objective function that we want our neural networks to optimize its weights according to it. Therefore, it is task-specific and also somehow empirical. Just to be clear, Multinomial Logistic Loss and Cross Entropy Loss are the same (please look at http://ufldl.stanford.edu/tutorial/supervised/SoftmaxRegression/). The cost function of Multinomial Logistic Loss is like this
$J(\theta) = -\frac{1}{m} \left[ \sum_{i=1}^m y^{(i)} \log h_\theta(x^{(i)}) + (1-y^{(i)}) \log (1-h_\theta(x^{(i)})) \right].$
It is usually used for classification problem. The Square Error has equation like
$\frac 1 {2N} \sum_{i=1}^N \| x^1_i - x^2_i \|_2^2.$
Therefore, it is usually used for minimize using some construction errors.
EDIT: @MartinThoma The above formula of multinomial logistics loss is just for binary case, for general case, it should be $J(\theta) = -\left[ \sum_{i=1}^{m} \sum_{k=1}^{K} 1\left\{y^{(i)} = k\right\} \log P(y^{(i)} = k | x^{(i)} ; \theta) \right]$, where K is number of categories.
