I haven't seen an answer from a trusted source, but I'll try to answer this myself, with a simple example (with my current knowledge).
In general, note that training an MLP using back-propagation is usually implemented with matrices.
Time complexity of matrix multiplication
The time complexity of matrix multiplication for $M_{ij} * M_{jk}$ is simply $\mathcal{O}(i*j*k)$.
Notice that we are assuming the simplest multiplication algorithm here: there exist some other algorithms with somewhat better time complexity.
Feedforward pass algorithm
The feedforward propagation algorithm is as follows.
First, to go from layer $i$ to $j$, you do
$$S_j = W_{ji}*Z_i$$
Then you apply the activation function
$$Z_j = f(S_j)$$
If we have $N$ layers (including input and output layer), this will run $N-1$ times.
Example
As an example, let's compute the time complexity for the forward pass algorithm for an MLP with $4$ layers, where $i$ denotes the number of nodes of the input layer, $j$ the number of nodes in the second layer, $k$ the number of nodes in the third layer and $l$ the number of nodes in the output layer.
Since there are $4$ layers, you need $3$ matrices to represent weights between these layers. Let's denote them by $W_{ji}$, $W_{kj}$ and $W_{lk}$, where $W_{ji}$ is a matrix with $j$ rows and $i$ columns ($W_{ji}$ thus contains the weights going from layer $i$ to layer $j$).
Assume you have $t$ training examples. For propagating from layer $i$ to $j$, we have first
$$S_{jt} = W_{ji} * Z_{it}$$
and this operation (i.e. matrix multiplication) has $\mathcal{O}(j*i*t)$ time complexity. Then we apply the activation function
$$
Z_{jt} = f(S_{jt})
$$
and this has $\mathcal{O}(j*t)$ time complexity, because it is an element-wise operation.
So, in total, we have
$$\mathcal{O}(j*i*t + j*t) = \mathcal{O}(j*t*(i + 1)) = \mathcal{O}(j*i*t)$$
Using same logic, for going $j \to k$, we have $\mathcal{O}(k*j*t)$, and, for $k \to l$, we have $\mathcal{O}(l*k*t)$.
In total, the time complexity for feedforward propagation will be
$$\mathcal{O}(j*i*t + k*j*t + l*k*t) = \mathcal{O}(t*(ij + jk + kl))$$
I'm not sure if this can be simplified further or not. Maybe it's just $\mathcal{O}(t*i*j*k*l)$, but I'm not sure.
Back-propagation algorithm
The back-propagation algorithm proceeds as follows. Starting from the output layer $l \to k$, we compute the error signal, $E_{lt}$, a matrix containing the error signals for nodes at layer $l$
$$
E_{lt} = f'(S_{lt}) \odot {(Z_{lt} - O_{lt})}
$$
where $\odot$ means element-wise multiplication. Note that $E_{lt}$ has $l$ rows and $t$ columns: it simply means each column is the error signal for training example $t$.
We then compute the "delta weights", $D_{lk} \in \mathbb{R}^{l \times k}$ (between layer $l$ and layer $k$)
$$
D_{lk} = E_{lt} * Z_{tk}
$$
where $Z_{tk}$ is the transpose of $Z_{kt}$.
We then adjust the weights
$$
W_{lk} = W_{lk} - D_{lk}
$$
For $l \to k$, we thus have the time complexity $\mathcal{O}(lt + lt + ltk + lk) = \mathcal{O}(l*t*k)$.
Now, going back from $k \to j$. We first have
$$
E_{kt} = f'(S_{kt}) \odot (W_{kl} * E_{lt})
$$
Then
$$
D_{kj} = E_{kt} * Z_{tj}
$$
And then
$$W_{kj} = W_{kj} - D_{kj}$$
where $W_{kl}$ is the transpose of $W_{lk}$. For $k \to j$, we have the time complexity $\mathcal{O}(kt + klt + ktj + kj) = \mathcal{O}(k*t(l+j))$.
And finally, for $j \to i$, we have $\mathcal{O}(j*t(k+i))$. In total, we have
$$\mathcal{O}(ltk + tk(l + j) + tj (k + i)) = \mathcal{O}(t*(lk + kj + ji))$$
which is the same as the feedforward pass algorithm. Since they are the same, the total time complexity for one epoch will be $$O(t*(ij + jk + kl)).$$
This time complexity is then multiplied by the number of iterations (epochs). So, we have $$O(n*t*(ij + jk + kl)),$$ where $n$ is number of iterations.
Notes
Note that these matrix operations can greatly be parallelized by GPUs.
Conclusion
We tried to find the time complexity for training a neural network that has 4 layers with respectively $i$, $j$, $k$ and $l$ nodes, with $t$ training examples and $n$ epochs. The result was $\mathcal{O}(nt*(ij + jk + kl))$.
We assumed the simplest form of matrix multiplication that has cubic time complexity. We used the batch gradient descent algorithm. The results for stochastic and mini-batch gradient descent should be the same. (Let me know if you think the otherwise: note that batch gradient descent is the general form, with little modification, it becomes stochastic or mini-batch)
Also, if you use momentum optimization, you will have the same time complexity, because the extra matrix operations required are all element-wise operations, hence they will not affect the time complexity of the algorithm.
I'm not sure what the results would be using other optimizers such as RMSprop.
Sources
The following article http://briandolhansky.com/blog/2014/10/30/artificial-neural-networks-matrix-form-part-5 describes an implementation using matrices. Although this implementation is using "row major", the time complexity is not affected by this.
If you're not familiar with back-propagation, check this article:
http://briandolhansky.com/blog/2013/9/27/artificial-neural-networks-backpropagation-part-4