Example of backpropagation for neural network with softmax and sigmoid activation I am trying to produce a NN algorithm to classify the species of Iris into three species (versicolor, virginica, setosa) - preferably in R. 
The scaffolding / source is this code in R with ReLU activation of the hidden unit and softmax. The question is code-neutral, and an alternative source is this post in Python, probably by the same authors.
I want to solve the backpropagation algorithm with sigmoid activation (as opposed to ReLU) of a 6-neuron single hidden layer without using packaged functions (just to gain insight into backpropagation). As in the linked posts the architecture is as follows:

The compiled code for the links is here for R and here for Python. 

Here is the ready to paste (now working) code reproduced in the answer.
 A: Here is the code with some substantial modifications, and adapted to the idea of both sigmoid activation of the intermediate (6-node hidden) layer, and softmax output:
First some housekeeping mainly to get the dataset split into training and testing:
set.seed(9)                                                      # reproducible example
data = iris[order(runif(nrow(iris))), ] ; rownames(data) <- NULL # shuffling the dataset
training_rows = sample(1:nrow(data), round(0.60 * nrow(data)))   # 60% training
data_train = data[training_rows, ]                               # final training set
data_test  = data[-training_rows, ]                              # 40% testing set

The main function is train, which takes as input the number of columns in the dataset that will act as explanatory variables (in this iris example, x = 1:4); the column containing the multi-class labels (in this case y = 5); the dataset to be trained; the number of nodes in the hidden layer; the number of iterations; the lr (learning rate); and the regularization rate:
train =  function( x, y, data = data, hidden = 6, iter = 2000, lr = .1, reg = .001 )  
{
  N = nrow(data)                          # total number of examples in the data set 
  X = unname(data.matrix(data[ ,x]))      # model matrix without headers
  X = cbind(X, rep(1, nrow(X)))           # add a bias column to model matrix
  Y = data[ ,y]                           # extract dependent variable
  if(is.factor(Y)) {Y <- as.integer(Y)}   # numerical coding of species
  Y.len = length(unique(Y))               # number of unique species = 3
  Y.set = sort(unique(Y))                 # in order - [1] 1 2 3
  Y.index = cbind(1:N, match(Y, Y.set))   # matching example and "answer" (1, 2, 3)
  D = ncol(X)                             # number columns in model matrix

  # Initiating weights:  
  W1 = 0.01 * matrix(rnorm(ncol(X) * hidden), nrow = ncol(X)) 
  W2 = 0.01 * matrix(rnorm((hidden + 1) * Y.len ), nrow = hidden + 1)

  # Training:
  i = 0
  while (i < iter){
    i = i + 1

  # Forward pass:
    hidden.in       = X %*% W1                  # input of the hidden layer
    hidden.out      = 1 / (1 + exp(-hidden.in)) # sigmoid activation 
    # Introducing bias as a column of 1's:
    hidden.out.bias = cbind(hidden.out,rep(1, nrow(hidden.out))) 
    outer.in        = hidden.out.bias %*% W2    # input of outer layer
    # Softmax:
    outer.net.exp   = exp(outer.in)             # exp inputs of outer layer
    # Normalizing into probabilities:
    outer.out       = sweep(outer.net.exp, 1, rowSums(outer.net.exp), '/') 

  # Backpropagation:
    dscores = outer.out                         # probability scores
    # Softmax derivative wrt the input of outer layer (applied to "true" column only):
    dscores[Y.index] = dscores[Y.index] - 1     
    loss.wrt.outer.in = dscores / N             # dividing no. of examples.
    # Partial derivative of error with respect to W2:
    dW2.temp = t(hidden.out.bias) %*% loss.wrt.outer.in
    dW2 = dW2.temp  + reg * W2                  # regularization
    W2 = W2 - lr * dW2                          # updating W2 weight matrix
    # Updating first matrix of weights (W1):
    loss.wrt.hidden.out = loss.wrt.outer.in %*% t(W2[-nrow(W2),]) 
    # excluding last row of W2 because we can't attribute error 
    # to the hidden layer when the bias was added later.
    # Derivative of the sigmoid function in the hidden layer:
    sigmoid = hidden.out * (1 - hidden.out) 
    # Partial derivative of error wrt to W1:
    dW1.temp = t(X) %*% (loss.wrt.hidden.out * sigmoid)
    dW1 = dW1.temp + reg * W1                   # regularization
    W1 = W1 - lr * dW1                          # updating W1 weights
  }
  # Results:
  model <- list( W1= W1, W2= W2)                # final W1 and W2 matrices
  return(model)
}        

Continuing with the "manual" approach to this toy example, we can get the weight matrices in model$W1 and model$W2, and then apply them to the test data:
model = train(x = 1:4, y = 5, data = data_train, hidden = 6)

# inserting a column of 1's to the test data as the bias:
new.data = as.matrix(cbind(data_test[-5], rep(1, nrow(data_test))))
# applying W1 matrix of weights to produce the input of the hidden layer:
hidden.layer = new.data %*% model$W1
# activating the hidden layer:
hidden.layer = 1 / (1 + exp(- hidden.layer))
# introducing the bias term:
hidden.layer.bias = cbind(hidden.layer, rep(1, nrow(hidden.layer)))
# producing the input for the output layer:
score = hidden.layer.bias %*% model$W2
# exponentiating and normalizing:
score.exp = exp(score)
probs <-sweep(score.exp, 1, rowSums(score.exp), '/') 
# selecting "winners" (highest probability values):
labels.predicted <- max.col(probs)

The results:
> table(data_test[,5], labels.predicted)
            labels.predicted
              1  2  3
  setosa     16  0  0
  versicolor  0 21  1
  virginica   0  1 21
> mean(as.integer(data_test[, 5]) == labels.predicted)
[1] 0.9666667

Here is a graphic illustration of multiple simulations including different random selection of training and testing sets, showing how among the $60$ examples in the testing set, there typically are from $0 - 2$ misclassified examples, which consistently lie in areas of data cloud overlap between species:


Points to emphasize:


*

*We want to update $W_2$, backpropagating the error $(L)$ through the NN tree. So we need to see how much the loss $(L)$ changes with changes in $W_2$. This is the $\Delta_o$ (output delta):


$$\frac{\partial L}{\partial W_2}=\Delta_0=\color{blue}{\large\frac{\partial \text{outer}_{input}}{\partial W_2}}\,\color{red}{\large \frac{\partial L}{\partial \text{outer}_{input}}}$$
$$\color{red}{\delta_0}=\color{red}{\large \frac{\partial L}{\partial \text{outer}_{input}}} = \color{orange}{E_0} \circ \color{brown}{D_0}= \color{orange}{\frac{\partial L}{\partial \text{outer}_{output}}}\,\circ\color{brown}{\frac{\partial \text{outer}_{output}}{\partial \text{outer}_{input}}}$$
where $E_0$ is the error calculated of the output layer and $D_0$ is the derivative of the activation function in the outer layer. $\circ$ is the Hadamard product.
However, in the softmax case there is no real activation function of the output layer, and $\color{red}{\delta_0} = p_k-\mathbf{1} (y_i=k)$, where $\mathbf{1} (y_i=k)$ is the indicator variable that denotes that the calculated probability matches the correct class. From this post:

Lets introduce the intermediate variable $p$, which is a vector of the
  (normalized) probabilities. The loss for one example is:
$$p_k = \frac{e^{f_k}}{ \sum_j e^{f_j} } \hspace{1in} L_i
 =-\log\left(p_{y_i}\right)$$
We now wish to understand how the computed scores inside $f$ should
  change to decrease the loss $L_i$ that this example contributes to the
  full objective. In other words, we want to derive the gradient
  $\partial L_i/\partial f_k$. The loss $L_i$ is computed from $p$,
  which in turn depends on $f$. It’s a fun exercise to the reader to use
  the chain rule to derive the gradient, but it turns out to be
  extremely simple and interpretible [sic] in the end, after a lot of
  things cancel out:
$$\frac{\partial L_i }{ \partial f_k } = p_k - \mathbb{1}(y_i = k)$$
Notice how elegant and simple this expression is. Suppose the
  probabilities we computed were $p = [0.2, 0.3, 0.5]$, and that the
  correct class was the middle one (with probability $0.3$). According
  to this derivation the gradient on the scores would be $\text{df} =
 [0.2, -0.7, 0.5]$.



*The error in the output layer can be propagated back to the output of the hidden layer:


$$\begin{align}E_H=E_h=E_{\text{hidden}}&=\frac{\partial L}{\partial \text{hidden}_{output}}\\[2ex]
&=\color{orange}{\frac{\partial L}{\partial \text{outer}_{output}}}\,\circ\color{brown}{\frac{\partial \text{outer}_{output}}{\partial \text{outer}_{input}}}\,\frac{\partial\text{outer}_{input}}{\partial \text{hidden}_{output}}\\[2ex]
&=\color{red}{\delta_0}\cdot W_2^\top= \color{orange}{E_0}\circ \color{brown}{D_0} \cdot W_2^T
\end{align}$$
However, it is important to note that it is necessary to exclude the last row in $W_2$ before performing the matrix multiplication in $\color{orange}{E_0}\circ \color{brown}{D_0} \cdot W_2^T$: the last row of $W_2$ is the product of having added a bias $1$ column to the output of the hidden layer; and since this was done after activating the hidden layer, no backpropagation of the eventual error to the hidden layer can be channeled through the bias term (or its corresponding row in $W_2$.


*Calculating delta for the hidden layer (change in error wrt input in the hidden layer) requires calculating the derivative of the activation function (sigmoid), and then the change in the error wrt the input of the hidden layer (small $\delta$):


$$\begin{align}\color{red}{\delta_H}&=\color{orange}{\frac{\partial L}{\partial \text{hidden}_{output}}}\,\circ\color{brown}{\frac{\partial \text{hidden}_{output}}{\partial\text{hidden}_{input}}}\\[2ex]
&=\color{orange}{E_H}\circ \color{brown}{D_H}\\[2ex]
&=E_H \circ \left(\text{hidden}_{out} \circ (1 - \text{hidden}_{out})\right)
\end{align}$$
and finally computing the the $\Delta$ of the weight matrix (change in error wrt $W_1$):
$$\begin{align}\Delta_H&=\color{orange}{\frac{\partial L}{\partial \text{hidden}_{output}}}\,\circ\,\color{brown}{\frac{\partial \text{hidden}_{output}}{\partial\text{hidden}_{input}}}\frac{\partial \text{hidden}_{input}}{W_1}\\[2ex]
&=X^\top \cdot \delta_H
\end{align}$$
