Softmax with cross-entropy

Posted on June 25, 2017

backpropogation, matrix calculus, softmax, cross-entropy, neural networks, deep learning

A matrix-calculus approach to deriving the sensitivity of cross-entropy cost to the weighted input to a softmax output layer. We use row vectors and row gradients, since typical neural network formulations let columns correspond to features, and rows correspond to examples. This means that the input to our softmax layer is a row vector with a column for each class.


Softmax

Softmax is a vector-to-vector transformation that turns a row vector

x=[x1 x2 ... xn]\displaystyle \mathbf x = \begin{bmatrix} x_1 \ \, x_2 \ \, ... \ \, x_n \end{bmatrix}

into a normalized row vector

s(x)=[s(x)1 s(x)2 ... s(x)n]\displaystyle \mathbf s(\mathbf x) = \begin{bmatrix} s(\mathbf x)_1 \ \, s(\mathbf x)_2 \ \, ... \ \, s(\mathbf x)_n \end{bmatrix}

The transformation is easiest to describe element-wise. The iith output s(x)is(\mathbf x)_i is a function of the entire input x\mathbf x, and is given by

s(x)i=exi∑exj\displaystyle s(\mathbf x)_i = \frac{e^{x_i}}{\sum e^{x_j}}

Softmax is nice because it turns x\mathbf x into a probability distribution.


Invariance to scaling

Softmax is invariant to additively scaling x\mathbf x by a constant cc.

s(x+c)i=exi+c∑exj+c=exiec∑exjec=exi∑exj\displaystyle \begin{aligned} s(\mathbf x + c)_i &= \frac{e^{x_i + c}}{\sum e^{x_j + c}} \\[1.6em] &= \frac{e^{x_i}e^c}{\sum e^{x_j}e^c} \\[1.6em] &= \frac{e^{x_i}}{\sum e^{x_j}} \end{aligned}

In other words, softmax only cares about the relative differences in the elements of x\mathbf x. That means we can protect softmax from overflow by subtracting the maximum element of x\mathbf x from every element of x\mathbf x. This will also protect against underflow because the denominator will contain a sum of non-negative terms, one of which is exmax−xmax=1e^{x_\text{max} - x_\text{max}} = 1.


Softmax in Python

def softmax(x):
    """Return the softmax of a vector x.
    
    :type x: ndarray
    :param x: vector input
    
    :returns: ndarray of same length as x
    """
    x = x - np.max(x)
    row_sum = np.sum(np.exp(x))
    return np.array([np.exp(x_i) / row_sum for x_i in x])


Jacobian of softmax

Since softmax is a vector-to-vector transformation, its derivative is a Jacobian matrix. The Jacobian has a row for each output element sis_i, and a column for each input element xjx_j.

The entries of the Jacobian take two forms, one for the main diagonal entry, and one for every off-diagonal entry. We’ll show how to compute these entries for an arbitrary row ii of the Jacobian.


Diagonal row entry

First compute the diagonal entry of row ii. That is, compute the derivative of the ii‘th output, sis_i, with respect to its ii‘th input, xix_i. All we need is the division rule from calculus.

∂s(x)i∂xi=∑jexjexi−exiexi(∑jexj)2=s(x)i−s(x)i2\displaystyle \begin{aligned} \frac{\partial s(\mathbf x)_i}{\partial x_i} &= \frac{\sum_j e^{x_j} e^{x_i} - e^{x_i} e^{x_i}}{(\sum_j e^{x_j})^2} \\[1.6em] &= s(\mathbf x)_i - s(\mathbf x)_i^2 \end{aligned}


Off-diagonal row entries

Now compute every off-diagonal entry of row ii. That is, compute the derivative of the ii‘th output, sis_i, with respect to its jj‘th input, xjx_j, where j≠ij \neq i. Again we use the division rule, but in this case the derivative of the numerator, exie^{x_i} with respect to xjx_j is zero, because j≠ij \neq i means the numerator is constant with respect to xjx_j.

∂s(x)i∂xj=∑jexj⋅0−exiexj(∑jexj)2=−s(x)is(x)j\displaystyle \begin{aligned} \frac{\partial s(\mathbf x)_i}{\partial x_j} &= \frac{\sum_j e^{x_j} \cdot 0 - e^{x_i} e^{x_j}}{(\sum_j e^{x_j})^2} \\[1.6em] &= - s(\mathbf x)_i s(\mathbf x)_j \end{aligned}

This is nice! The derivative of softmax is always phrased in terms of softmax.

From now on, to keep things clear, we won’t write dependence on x\mathbf x. Instead we’ll write s(x)\mathbf s(\mathbf x) as s\mathbf s and s(x)is(\mathbf x)_i as sis_i, understanding that s\mathbf s and sis_i are each a function of the entire vector x\mathbf x.


Full Jacobian

The form of the off-diagonals tells us that the Jacobian of softmax is a symmetric matrix. This is nice because symmetric matrices have great numeric and analytic properties. We expand it below. Each row is a gradient of one output element sis_i with respect to each of its input elements xjx_j.

Jx(s)=[s0−s02−s0s1...−s0sn−s1s0s1−s12...−s1sn............−sns0−sns1...sn−sn2]\displaystyle \mathbf J_{\mathbf x}(\mathbf s) = \begin{bmatrix} s_0 - s_0^2 & -s_0 s_1 & ... & -s_0 s_n \\[0.6em] -s_1 s_0 & s_1 - s_1^2 & ... & -s_1 s_n \\[0.6em] ... & ... & ... & ... \\[0.6em] -s_n s_0 & -s_n s_1 & ... & s_n - s_n^2 \end{bmatrix}

Notice that we can express this matrix as

Jx(s)=diag(s)−s⊤s\displaystyle \mathbf J_{\mathbf x}(\mathbf s) = \text{diag}(\mathbf s) - \mathbf s^\top \mathbf s

where the second term is the n×nn \times n outer product, because we defined s\mathbf s as a row vector.


Jacobian of softmax in Python

def jacobian_softmax(s):
    """Return the Jacobian matrix of softmax vector s.

    :type s: ndarray
    :param s: vector input

    :returns: ndarray of shape (len(s), len(s))
    """
    return np.diag(s) - np.outer(s, s)


Cross-entropy

Cross-entropy measures the difference between two probability distributions. We saw that s\mathbf s is a distribution. The correct class is also a distribution if we encode it as a one-hot vector:

y=[y1 y2 ... yn]=[0  0  ...  1  ...  0]\displaystyle \begin{aligned} \mathbf y &= \begin{bmatrix} y_1 \ \, y_2 \ \, ... \ \, y_n \end{bmatrix} \\[1.1em] &= \begin{bmatrix} 0 \ \ \, 0 \ \ \, ... \ \ \, 1 \ \ \, ... \ \ \, 0 \end{bmatrix} \end{aligned}

where the 11 appears at the index of the correct class of this single example.

The cross-entropy between our predicted distribution over classes, s\mathbf s, and the true distribution over classes, y\mathbf y, is a scalar measure of their difference, which is perfect for a cost function. It’ll drive our softmax distribution toward the one-hot distribution. We can write this cost function as

H(y,s)=−∑i=1nyilogsi=−ylogs⊤\displaystyle \begin{aligned} H(\mathbf y, \mathbf s) &= -\sum_{i=1}^n y_i \log s_i \\[1.6em] &= -\mathbf y \log \mathbf s^\top \end{aligned}

which is the dot product since we’re using row vectors. This formula comes from information theory. It measures the information gained about our softmax distribution when we sample from our one-hot distribution.


Cross-entropy in Python

def cross_entropy(y, s):
    """Return the cross-entropy of vectors y and s.

    :type y: ndarray
    :param y: one-hot vector encoding correct class

    :type s: ndarray
    :param s: softmax vector

    :returns: scalar cost
    """
    # Naively computes log(s_i) even when y_i = 0
    # return -y.dot(np.log(s))
    
    # Efficient, but assumes y is one-hot
    return -np.log(s[np.where(y)])


Gradient of cross-entropy

Since our y\mathbf y is given and fixed, cross-entropy is a vector-to-scalar function of only our softmax distribution. That means it will have a gradient with respect to our softmax distribution. This vector-to-scalar cost function is actually made up of two steps: (1) a vector-to-vector element-wise log\log and (2) a vector-to-scalar dot product. The vector-to-vector logarithm will have a Jacobian, but since it’s applied element-wise, the Jacobian will be diagonal, holding each elementwise derivative. The gradient of a dot product, being a linear operation, is just the vector y\mathbf y.

∇sH=−∇sylogs⊤=−y∇slogs=−y diag(1s)=−ys\displaystyle \begin{aligned} \nabla_{\mathbf s} H &= -\nabla_{\mathbf s} \mathbf y \log \mathbf s^\top \\[1.6em] &= -\mathbf y \nabla_{\mathbf s} \log \mathbf s \\[1.1em] &= -\mathbf y \ \text{diag}\left(\frac{\mathbf 1}{\mathbf s}\right) \\[1.6em] &= -\frac{\mathbf y}{\mathbf s} \end{aligned}

where we used equation (69) of the matrix cookbook for the derivative of the dot product.


Gradient of cross-entropy in Python

def gradient_cross_entropy(y, s):
    """Return the gradient of cross-entropy of vectors y and s.

    :type y: ndarray
    :param y: one-hot vector encoding correct class

    :type s: ndarray
    :param s: softmax vector

    :returns: ndarray of size len(s)
    """
    return -y / s


Error at input to softmax layer

By the chain rule, the sensitivity of cost HH to the input to the softmax layer, x\mathbf x, is given by a gradient-Jacobian product, each of which we’ve already computed:

∇xH=∇sH Jx(s)\displaystyle \nabla_{\mathbf x} H = \nabla_{\mathbf s} H \ \mathbf J_{\mathbf x}(\mathbf s)

The first term is the gradient of cross-entropy to softmax activation. The second term is the Jacobian of softmax activation to softmax input. Remember that we’re using row gradients - so this is a row vector times a matrix, resulting in a row vector. Expanding and simplifying, we get

∇xH=−ys[diag(s)−s⊤s]=yss⊤s−ys diag(s)=y Srepeated row−y=s−y\displaystyle \begin{aligned} \nabla_{\mathbf x} H &= -\frac{\mathbf y}{\mathbf s} \bigg[\text{diag}(\mathbf s) -\mathbf s^\top \mathbf s \bigg] \\[1.6em] &= \frac{\mathbf y}{\mathbf s} \mathbf s^\top \mathbf s - \frac{\mathbf y}{\mathbf s} \ \text{diag}(\mathbf s) \\[1.6em] &= \mathbf y \ \mathbf S^{\text{repeated row}} - \mathbf y \\[1.6em] &= \mathbf s - \mathbf y \end{aligned}

The last line follows from the fact that y\mathbf y was one-hot and applied to a matrix whose rows are identically our softmax distribution. But actually, any y\mathbf y whose elements sum to 11 would satisfy the same property. To be more specific, the equation above would hold not just for one-hot y\mathbf y, but for any y\mathbf y specifying a distribution over classes.


Error at input to softmax layer in Python

def error_softmax_input(y, s):
    """Return the sensitivity of cross-entropy cost to input of softmax.

    :type y: ndarray
    :param y: one-hot vector encoding correct class

    :type s: ndarray
    :param s: softmax vector

    :returns: ndarray of size len(s)
    """
    return s - y


Now with a batch of examples

Our work thus far considered a single example. Hence x\mathbf x, our input to the softmax layer, was a row vector. Alternatively, if we feed forward a batch of mm examples, then X\mathbf X contains a row vector for each example in the minibatch.

X=[x1x2...xm]∼m×n\displaystyle \mathbf X = \begin{bmatrix} \mathbf x_1 \\[0.6em] \mathbf x_2 \\[0.6em] ... \\[0.6em] \mathbf x_m \end{bmatrix} \sim m \times n


Batch softmax

Softmax is still a vector-to-vector transformation, but it’s applied independently to each row of X\mathbf X.

S=[s(x1)s(x2)...s(xm)]∼m×n\displaystyle \mathbf S = \begin{bmatrix} \mathbf s(\mathbf x_1) \\[0.6em] \mathbf s(\mathbf x_2) \\[0.6em] ... \\[0.6em] \mathbf s(\mathbf x_m) \end{bmatrix} \sim m \times n


Batch softmax in Python

def batch_softmax(x):
    """Return matrix of row-wise softmax of x.

    :type x: ndarray
    :param x: row per example and column per feature

    :returns: ndarray of x.shape after row-wise softmax
    """
    # Stabilize by subtracting row max from each row
    row_maxes = np.max(x, axis=1)
    row_maxes = row_maxes[:, np.newaxis]  # for broadcasting
    x = x - row_maxes

    return np.array([softmax(row) for row in x])


Jacobian of batch softmax

Because rows are independently mapped, the Jacobian of row ii of S\mathbf S with respect to row j≠ij \neq i of X\mathbf X is a zero matrix.

Jxj≠i(si)=0\displaystyle \mathbf J_{\mathbf x_{j \neq i}}(\mathbf s_i) = \mathbf 0

and the Jacobian of row ii of S\mathbf S with respect to row ii of X\mathbf X is our familiar matrix from before

Jxi(si)=diag(si)−si⊤si\displaystyle \mathbf J_{\mathbf x_{i}}(\mathbf s_i) = \text{diag}(\mathbf s_i) - \mathbf s_i^\top \mathbf s_i

That means our grand Jacobian of S\mathbf S with respect to X\mathbf X is a diagonal m×mm \times m matrix of n×nn \times n matrices, most of which are zero matrices:

JX(S)=[Jx1(s1)0...00Jx2(s2)...0............00...Jxm(sm)]\displaystyle \mathbf J_{\mathbf X}(\mathbf S) = \begin{bmatrix} \mathbf J_{\mathbf x_1}(\mathbf s_1) & \mathbf 0 & ... & \mathbf 0 \\[1.1em] \mathbf 0 & \mathbf J_{\mathbf x_2}(\mathbf s_2) & ... & \mathbf 0 \\[1.1em] ... & ... & ... & ... \\[1.1em] \mathbf 0 & \mathbf 0 & ... & \mathbf J_{\mathbf x_m}(\mathbf s_m) \end{bmatrix}


Jacobian of batch softmax in Python

def jacobian_batch_softmax(s):
    """Return array of row-wise Jacobians of s.

    :type s: ndarray
    :param s: matrix whose rows are softmaxed

    :returns: ndarray of shape
              (s.shape[0], s.shape[0], s.shape[1], s.shape[1])
    """
    # Array of nonzero Jacobians lying along tensor diagonal
    return np.array([jacobian_softmax(row) for row in s])


Mean cross-entropy of batch

Let each row of Y\mathbf Y be a one-hot label for an example:

Y=[y1y2...ym]∼m×n\displaystyle \mathbf Y = \begin{bmatrix} \mathbf y_1 \\[0.6em] \mathbf y_2 \\[0.6em] ... \\[0.6em] \mathbf y_m \end{bmatrix} \sim m \times n

Then we compute the mean cross-entropy by averaging the cross-entropy of every matching pair of rows of Y\mathbf Y and S\mathbf S. That is, we average over examples, the cross-entropy of each example:

H(Y,S)=−1m∑i=1myilogsi⊤=−1mTr(YlogS⊤)\displaystyle \begin{aligned} H(\mathbf Y, \mathbf S) &= -\frac{1}{m} \sum_{i=1}^m \mathbf y_i \log \mathbf s_i^\top \\[1.6em] &= -\frac{1}{m} \text{Tr}(\mathbf Y \log \mathbf S^\top) \end{aligned}

The above simplification works because each row of S\mathbf S is si\mathbf s_i. So each column of S⊤\mathbf S^\top is si\mathbf s_i. So the matrix product YlogS⊤\mathbf Y \log \mathbf S^\top dots rows of Y\mathbf Y with columns of (logS⊤)(\log \mathbf S^\top), which is exactly what we want for cross-entropy. Now, we only care about entries where the row index equals the column index. That’s because cross-entropy sums the dot products of matching rows of Y\mathbf Y and S\mathbf S. We can sum over matching dot products by using a trace.

Note: this formulation is computationally wasteful. We shouldn’t implement batch cross-entropy this way in a computer. We’re only using it for its analytic simplicity to work out the backpropogating error. However, the end analytic result is actually computationally efficient.


Mean cross-entropy of batch in Python

def mean_cross_entropy(y, s):
    """Return the mean row-wise cross-entropy of y and s.

    :type y: ndarray
    :param y: matrix whose rows are one-hot vectors encoding
              the correct class of each example.

    :type s: ndarray
    :param s: matrix whose every row is a softmax distribution over
              class predictions for a given example.

    :returns: scalar, mean row-wise cross-entropy cost
    """
    return np.mean([cross_entropy(y_row, s_row)
                    for y_row, s_row in zip(y, s)])


Jacobian of mean cross-entropy of batch

Since mean cross-entropy maps a matrix to a scalar, its Jacobian with respect to S\mathbf S will be a matrix.

JS(H)=−JS(1mTr(YlogS⊤))=−1mJSTr(YlogS⊤)=−1mYJSlogS=−1mY⊙1S=−1mYS\displaystyle \begin{aligned} \mathbf J_{\mathbf S}(H) &= -\mathbf J_{\mathbf S} \bigg(\frac{1}{m} \text{Tr}(\mathbf Y \log \mathbf S^\top) \bigg) \\[1.6em] &= -\frac{1}{m} \mathbf J_{\mathbf S} \text{Tr}(\mathbf Y \log \mathbf S^\top) \\[1.6em] &= -\frac{1}{m} \mathbf Y \mathbf J_{\mathbf S} \log \mathbf S \\[1.6em] &= -\frac{1}{m} \mathbf Y \odot \frac{1}{\mathbf S} \\[1.6em] &= -\frac{1}{m} \frac{\mathbf Y}{\mathbf S} \end{aligned}

Since logS\log \mathbf S is an element-wise operation mapping a matrix to a matrix, its Jacobian is a matrix of element-wise derivatives which we chain rule by a Hadamard product, rather than by a dot product.


Why this works

This procedure is always true for any element-wise operations. We can see this by concatenating the rows of S\mathbf S.

s=[s1 s2 ... sm]\displaystyle \mathbf s = \begin{bmatrix} \mathbf s_1 \ \, \mathbf s_2 \ \, ... \ \, \mathbf s_m \end{bmatrix}

such that s\mathbf s is a row vector of length m⋅nm \cdot n. Then log\log is an element-wise vector-to-vector transformation again. So it has an m⋅n×m⋅nm \cdot n \times m \cdot n diagonal Jacobian matrix.

Js(logs)=diag(1s)\displaystyle \mathbf J_{\mathbf s}(\log \mathbf s) = \text{diag}\left(\frac{\mathbf 1}{\mathbf s}\right)

If we flatten Y\mathbf Y in the same way

y=[y1 y2 ... ym]\displaystyle \mathbf y = \begin{bmatrix} \mathbf y_1 \ \, \mathbf y_2 \ \, ... \ \, \mathbf y_m \end{bmatrix}

then we get

H(y,s)=−1mylogs⊤\displaystyle H(\mathbf y, \mathbf s) = - \frac{1}{m} \mathbf y \log \mathbf s^\top

and so

Js(H)=−1my∂∂s(logs)=−1mydiag(1s)=−1mys\displaystyle \begin{aligned} \mathbf J_{\mathbf s}(H) &= -\frac{1}{m} \mathbf y \frac{\partial}{\partial \mathbf s} \Big(\log \mathbf s \Big) \\[1.6em] &= -\frac{1}{m} \mathbf y \text{diag}\left(\frac{\mathbf 1}{\mathbf s}\right) \\[1.6em] &= -\frac{1}{m} \frac{\mathbf y}{\mathbf s} \end{aligned}

Now since y\mathbf y and s\mathbf s are each of length m⋅nm \cdot n, we can reshape this formulation back into matrices, understanding that in both cases the division is element-wise:

Js(H)=−1mYS\displaystyle \mathbf J_{\mathbf s}(H) = -\frac{1}{m} \frac{\mathbf Y}{\mathbf S}

and we have our result. □\square


Jacobian of mean cross-entropy of batch in Python

def jacobian_mean_cross_entropy(y, s):
    """Return the Jacobian matrix for mean cross-entropy.

    :type y: ndarray
    :param y: matrix whose rows are one-hot vectors encoding
              the correct class of each example.

    :type s: ndarray
    :param s: matrix whose every row is a softmax distribution over
              class predictions for a given example.

    :returns: ndarray of shape y.shape holding gradients as rows
    """
    return -(1 / y.shape[0]) * (y / s)


Error at input to softmax layer for batch

We apply the chain rule just as before. The only difference is that our gradient-Jacobian product is now a matrix-tensor product. Multiplying a matrix against a tensor is difficult. One approach is to flatten everything, do a vector-matrix product as before, and then reshape everything, but this is not elegant or intuitive. Instead, we dot rows of JS(H)\mathbf J_{\mathbf S}(H), each a gradient of a row-wise cross-entropy, against diagonal elements of JX(S)\mathbf J_{\mathbf X}(\mathbf S), each a Jacobian matrix of a row-wise softmax.

We are able to do this because of the fact that JX(S)\mathbf J_{\mathbf X}(\mathbf S) is diagonal, which breaks the matrix-tensor product into an element-wise dot product of gradients and Jacobians. We owe this entirely to the fact that softmax is a row-to-row transformation, such that its Jacobian tensor is diagonal.

JX(H)=JS(H) JX(S)=(−1mYS)JX(S)=1m[−y1/s1−y2/s2...−ym/sm]JX(S)=1m[−y1s1Jx1(s1)−y2s2Jx2(s2)...−ymsmJxm(sm)]=1m[s1−y1s2−y2...sm−ym]=1m(S−Y)\displaystyle \begin{aligned} \mathbf J_{\mathbf X}(H) &= \mathbf J_{\mathbf S}(H) \ \mathbf J_{\mathbf X}(\mathbf S) \\[1.6em] &= \bigg(-\frac{1}{m} \frac{\mathbf Y}{\mathbf S}\bigg) \mathbf J_{\mathbf X}(\mathbf S) \\[2.4em] &= \frac{1}{m} \begin{bmatrix} -\mathbf y_1 / \mathbf s_1 \\[0.4em] -\mathbf y_2 / \mathbf s_2 \\[0.4em] ... \\[0.4em] -\mathbf y_m / \mathbf s_m \end{bmatrix} \mathbf J_{\mathbf X}(\mathbf S) \\[2.4em] \\[0.6em] &= \frac{1}{m} \begin{bmatrix} -\frac{\mathbf y_1}{\mathbf s_1} \mathbf J_{\mathbf x_1}(\mathbf s_1) \\[1.4em] -\frac{\mathbf y_2}{\mathbf s_2} \mathbf J_{\mathbf x_2}(\mathbf s_2) \\[1.4em] ... \\[1.4em] -\frac{\mathbf y_m}{\mathbf s_m} \mathbf J_{\mathbf x_m}(\mathbf s_m) \end{bmatrix} \\[1.4em] \\[0.6em] &= \frac{1}{m} \begin{bmatrix} \mathbf s_1 - \mathbf y_1 \\[0.4em] \mathbf s_2 - \mathbf y_2 \\[0.4em] ... \\[0.4em] \mathbf s_m - \mathbf y_m \end{bmatrix} \\[1.4em] \\[0.6em] &= \frac{1}{m} \Big(\mathbf S - \mathbf Y\Big) \end{aligned}

Where the third step followed by the fact that JX(S)J_{\mathbf X}(\mathbf S) is diagonal. So the sensitivity of cost to the weighted input to our softmax layer is just the difference of our softmax matrix and our matrix of one-hot labels, where every element is divided by the number of examples in the batch.


Error at input to softmax layer for batch in Python

def batch_error_softmax_input(y, s):
    """Return the sensitivity of cross-entropy cost to input of softmax.

    :type y: ndarray
    :param y: matrix whose rows are one-hot vectors encoding
              the correct class of each example.

    :type s: ndarray
    :param s: matrix whose every row is a softmax distribution over
              class predictions for a given example.

    :returns: ndarray of shape y.shape
    """
    return (1 / y.shape[0]) * (S - Y)