Back Propagation

In a recent post on neural networks, using R I described neural networks and presented the following visualization from R:

 I have also described a multilayer perceptron as a weighted average or ensemble of logits. But how are the weights in each hidden layer logistic activation function (or any activation function for other network architectures) estimated? How are the weights in the combination functions estimated? Neural networks can be estimated using back propagation, described in Hastie as 'a generic approach to minimizing R(θ) (the cost function) by gradient descent.'

Given a neural network with inputs X with hidden layers comprised of hidden units Z used to predict some target T,  we can represent a neural network schematically (simplifying the notation in Hastie by omitting key subscripts and summations)

X -> Z -> T

Z = σ( α₀ +  α^Tx)
T =  β₀ + βZ
f(X) = g(T)   [1]

where σ = the activation function

Given weights {α₀,α₀ ,  β₀ , β} find the values that minimize the specified error function:

R(θ) =∑∑ ( y-f(x)² )      [2] (note a number of possible error functions may be used)

Backpropogation equations:

s = σ'( α^Tx )βδ   [3]

Errors can be re-specified as:

∂R/ ∂β = δZ   [4]

∂R/ ∂α = sx    [5]

Gradient Descent Update:

β^r+1 =  β^r - γ ∂R/ ∂β    [6]

α^r+1 =  α^r - γ ∂R/ ∂α   [7]

Algorithm:

Forward Pass: use initial or current weights (guesses) and calculate f(X), and errors δ from the output layer [2]

Backward Pass:  'back propagate' via back propagation equation [3] to obtain s. Both sets of errors (δ) and (s) are used to derive the derivative terms in [4] and [5]  which are then used in the gradient descent update weight estimates via equations [6]& [7].

In Predictive modeling with SAS Enterprise Miner by Sarma,  the following basic description of back propagation is given:

Specify an error function E.

1) 1st iteration- set initial weights, use to evaluate E

2) 2nd iteration- weights are changed by a small amount such that the error is redced

-repeat until convergence

As Sarma explains, with each iteration a number of weights are produced, so if it takes 100 iterations to converge, 100 possible models are specified, giving 100 sets of weights. Using validation data, the best iteration can be chosen calculating E via the validation data.

R Code Example for Neural Networks

See also NEURAL NETWORKS.

In this past June's issue of R journal, the 'neuralnet' package was introduced. I had recently been familiar with utilizing neural networks via the 'nnet' package (see my post on Data Mining in A Nutshell) but I find the neuralnet package more useful because it will allow you to actually plot the network nodes and connections. (it may be possible to do this with nnet, but I'm not aware of how).The neuralnet package was written primarily for multilayer perceptron architectures, which may be a limitation if you are interested in other architectures.

The data set used was a default data set found in the package 'datasets' and consisted of 248 observations and 8 variables:

"education"      "age"            "parity"         "induced"        "case"           "spontaneous"    "stratum"        "pooled.stratum"

The following code runs the network (with 2 hidden layers) classifying 'case'  (a binary variable) as a function of several independent varaibles. The neural network is estimated, and the results are stored in the data frame 'nn.'

nn <- neuralnet(

 case~age+parity+induced+spontaneous,

 data=infert, hidden=2, err.fct="ce",

 linear.output=FALSE)

The weight estimates can be obtained with the following command:

nn$result.matrix

And, the network can be plotted or visualized with the simple command:

plot(nn)

As can be seen below, the output weights correspond directly with the visualized network. While this tool may be more difficult to interpret with more complex network architectures, I find this simplified version very useful for demonstration/instructional purposes.

Example: The weight for the path from input 'age' to the first hidden layer is  -3.0689 (age.to.1layhid1) which can easily be found in the network diagram. After all inputs feed into hidden layer 1, the weight associated with the path from hidden layer 1(1layhid.1.to.case) to the output layer (which along with information from the other layers of the network will give us the classification of 'case') is -1001.15.


Estimated Network Weights

error	123.811
reached.threshold	0.009853
steps	16822
Intercept.to.1layhid1	0.535237
age.to.1layhid1	-3.0689
parity.to.1layhid1	2.262792
induced.to.1layhid1	30.50824
spontaneous.to.1layhid1	0.000961
Intercept.to.1layhid2	-5.52491
age.to.1layhid2	0.11538
parity.to.1layhid2	-1.98372
induced.to.1layhid2	2.550932
spontaneous.to.1layhid2	3.829613
Intercept.to.case	-1.7714
1layhid.1.to.case	-1001.15
1layhid.2.to.case	4.570324

   
Neural Network Visualization (click to enlarge)

There are many other options avaiable with the neural net package, I encourage you to read the article referenced below for more details.

Reference: neuralnet: Training of Neural Networks  by Frauke G¸nther and Stefan Fritsch The R Journal Vol. 2/1, June 2010

R-Code:

#  ------------------------------------------------------------------
#  |PROGRAM NAME: NEURALNET_PKG_R
#  |DATE: 12/3/10
#  |CREATED BY: MATT BOGARD 
#  |PROJECT FILE:   P:\R  Code References\Data Mining_R           
#  |----------------------------------------------------------------
#  | PURPOSE: DEMO OF THE 'neuralnet' PACKAGE AND OUTPUT INTERPRETATION              
#  | 
#  |  ADAPTED FROM:  neuralnet: Training of Neural Networks
#  |    by Frauke Günther and Stefan Fritsch The R Journal Vol. 2/1, June 2010 
#  |    ISSN 2073-4859  (LOCATED: P:\TOOLS AND REFERENCES (Copy)\R References\Neural Networks
#  | 
#  | 
#  |------------------------------------------------------------------
#  |DATA USED: 'infert' FROM THE 'datasets' LIBRARY              
#  |------------------------------------------------------------------
#  |CONTENTS:               
#  |
#  |  PART 1: get the data 
#  |  PART 2: train the network
#  |  PART 3: 
#  |  PART 4: 
#  |  PART 5: 
#  |------------------------------------------------------------------
#  |COMMENTS:               
#  |
#  |-----------------------------------------------------------------
#  |UPDATES:               
#  |
#  |
#  |------------------------------------------------------------------
 
#  *------------------------------------------------------------------*
#  | get the data
#  *------------------------------------------------------------------* 
 
library(datasets)
 
names(infert)
 
#  *------------------------------------------------------------------*
#  | train the network
#  *------------------------------------------------------------------* 
 
library(neuralnet)
 
nn <- neuralnet(
 case~age+parity+induced+spontaneous,
 data=infert, hidden=2, err.fct="ce",
 linear.output=FALSE)
 
#  *------------------------------------------------------------------*
#  | output training results
#  *------------------------------------------------------------------*  
 
# basic
nn
 
# reults options
names(nn)
 
 
# result matrix
 
nn$result.matrix
 
# The given data is saved in nn$covariate and
# nn$response as well as in nn$data for the whole data
# set inclusive non-used variables. The output of the
# neural network, i.e. the fitted values o(x), is provided
# by nn$net.result:
 
out <- cbind(nn$covariate,nn$net.result[[1]])
 
dimnames(out) <- list(NULL, c("age", "parity","induced","spontaneous","nn-output"))
 
head(out)
 
# generalized weights
 
# The generalized weight expresses the effect of each
# ovariate xi and thus has an analogous interpretation
# as the ith regression parameter in regression models.
# However, the generalized weight depends on all
# other covariates. Its distribution indicates whether
# the effect of the covariate is linear since a small variance
# suggests a linear effect
 
# The columns refer to the four covariates age (j =
# 1), parity (j = 2), induced (j = 3), and spontaneous (j=4)
 
head(nn$generalized.weights[[1]])
 
# visualization
 
plot(nn)

Created by Pretty R at inside-R.org

Neural Networks

NEURAL NETWORK: A nonlinear model of complex relationships composed of multiple 'hidden' layers (similar to composite functions)

Y = f(g(h(x)) or
x -> hidden layers ->Y


Example 1 With a logistic/sigmoidal activation function, a neural network can be visualized as a sum of weighted logits:



Y = α Σ w_i e ^θ_i/1 + e ^θ_i + ε
w_i = weights θ = linear function Xβ
Y= 2 + w₁ Logit A + w₂ Logit B + w₃ Logit C + w₄ Logit D
( Adapted from ‘A Guide to Econometrics, Kennedy, 2003)

Example 2

 Where:  Y= W₀ + W₁ Logit H₁ + W₂ Logit H₂ + W₃ Logit H₃ + W₄ Logit H₄   and

                        H₁= logit(w₁₀ +w₁₁ x₁ + w₁₂ x₂ )

H₂ = logit(w₂₀ +w₂₁ x₁ + w₂₂ x₂ )

H₃ = logit(w₃₀ +w₃₁ x₁ + w₃₂ x₂ )

H₄ = logit(w₄₀ +w₄₁ x₁ + w₄₂ x₂ )



 

The links between each layer in the diagram correspond to the weights (w’s) in each equation. The weights can be estimated via back propagation.

( Adapted from ‘A Guide to Econometrics, Kennedy, 2003 and Applied Analytics Using SAS Enterprise Miner 6.1)

 MULTILAYER PERCEPTRON: a neural network architecture that has one or more hidden layers, specifically having linear combination functions in the hidden and output layers, and sigmoidal activation functions in the hidden layers.  (note: a basic logistic regression function can be visualized as a single layer perceptron)

RADIAL BASIS FUNCTION (architecture):  a neural network architecture with exponential or softmax (generalized multinomial logistic) activation functions and radial basis combination functions in the hidden layers and linear combination functions in the output layers.

 RADIAL BASIS FUNCTION: A combination function that is based on the Euclidean distance between inputs and weights

ACTIVATION FUNCTION: formula used for transforming values from inputs and the outputs in a neural network.

COMBINATION FUNCTION: formula used for combining transformed values from activation functions in neural networks.
HIDDEN LAYER: The layer between input and output layers in a neural network.

R CODE EXAMPLE HERE

Data Mining in A Nutshell

# The following code may look rough, but simply paste into R or

# a text editor (especially Notepad++) and it will look

# much better.

# PROGRAM NAME: MACHINE_LEARNING_R

# DATE: 4/19/2010

# AUTHOR : MATT BOGARD

# PURPOSE: BASIC EXAMPLES OF MACHINE LEARNING IMPLEMENTATIONS IN R

# DATA USED: GENERATED VIA SIMULATION IN PROGRAM

# COMMENTS: CODE ADAPTED FROM : Joshua Reich (josh@i2pi.com) April 2, 2009

# SUPPORT VECTOR MACHINE CODE ALSO BASED ON :

# ëSupport Vector Machines in Rî Journal of Statistical Software April 2006 Vol 15 Issue 9

#

# CONTENTS: SUPPORT VECTOR MACHINES

# DECISION TREES

# NEURAL NETWORK

# load packages

library(rpart)
library(MASS)
library(class)
library(e1071)

# get data

# A simple function for producing n random samples

# from a multivariate normal distribution with mean mu

# and covariance matrix sigma

rmulnorm <- function (n, mu, sigma) {

    M <- t(chol(sigma))
    d <- nrow(sigma)
    Z <- matrix(rnorm(d*n),d,n)
    t(M %*% Z + mu)
}

# Produce a confusion matrix

# there is a potential bug in here if columns are tied for ordering

cm <- function (actual, predicted)

{
   
    t<-table(predicted,actual) 
    t[apply(t,2,function(c) order(-c)[1]),]
}

# Total number of observations

N <- 1000 * 3

# Number of training observations

Ntrain <- N * 0.7

# The data that we will be using for the demonstration consists

# of a mixture of 3 multivariate normal distributions. The goal

# is to come up with a classification system that can tell us,

# given a pair of coordinates, from which distribution the data

# arises.

A <- rmulnorm (N/3, c(1,1), matrix(c(4,-6,-6,18), 2,2))
B <- rmulnorm (N/3, c(8,1), matrix(c(1,0,0,1), 2,2))
C <- rmulnorm (N/3, c(3,8), matrix(c(4,0.5,0.5,2), 2,2))

data <- data.frame(rbind (A,B,C))
colnames(data) <- c('x', 'y')
data$class <- c(rep('A', N/3), rep('B', N/3), rep('C', N/3))

# Lets have a look

plot_it <- function () {
    plot (data[,1:2], type='n')
    points(A, pch='A', col='red')
    points(B, pch='B', col='blue')
    points(C, pch='C', col='orange')
}
plot_it()

# Randomly arrange the data and divide it into a training

# and test set.

data <- data[sample(1:N),]
train <- data[1:Ntrain,]
test <- data[(Ntrain+1):N,]

# SVM

# Support vector machines take the next step from LDA/QDA. However

# instead of making linear voronoi boundaries between the cluster

# means, we concern ourselves primarily with the points on the

# boundaries between the clusters. These boundary points define

# the 'support vector'. Between two completely separable clusters

# there are two support vectors and a margin of empty space

# between them. The SVM optimization technique seeks to maximize

# the margin by choosing a hyperplane between the support vectors

# of the opposing clusters. For non-separable clusters, a slack

# constraint is added to allow for a small number of points to

# lie inside the margin space. The Cost parameter defines how

# to choose the optimal classifier given the presence of points

# inside the margin. Using the kernel trick (see Mercer's theorem)

# we can get around the requirement for linear separation

# by representing the mapping from the linear feature space to

# some other non-linear space that maximizes separation. Normally

# a kernel would be used to define this mapping, but with the

# kernel trick, we can represent this kernel as a dot product.

# In the end, we don't even have to define the transform between

# spaces, only the dot product distance metric. This leaves

# this algorithm much the same, but with the addition of

# parameters that define this metric. The default kernel used

# is a radial kernel, similar to the kernel defined in my

# kernel method example. The addition is a term, gamma, to

# add a regularization term to weight the importance of distance.

s <- svm( I(factor(class)) ~ x + y, data = train, cost = 100, gama = 1)

s # print model results

#plot model and classification-my code not originally part of this

plot(s,test)

(m <- cm(train$class, predict(s)))
1 - sum(diag(m)) / sum(m)

(m <- cm(test$class, predict(s, test[,1:2])))
1 - sum(diag(m)) / sum(m)

# Recursive Partitioning / Regression Trees

# rpart() implements an algorithm that attempts to recursively split

# the data such that each split best partitions the space according

# to the classification. In a simple one-dimensional case with binary

# classification, the first split will occur at the point on the line

# where there is the biggest difference between the proportion of

# cases on either side of that point. The algorithm continues to

# split the space until a stopping condition is reached. Once the

# tree of splits is produced it can be pruned using regularization

# parameters that seek to ameliorate overfitting.

names(train)
names(test)

(r <- rpart(class ~ x + y, data = train)) plot(r) text(r) summary(r) plotcp(r) printcp(r) rsq.rpart(r)

cat("\nTEST DATA Error Matrix - Counts\n\n")
print(table(predict(r, test, type="class"),test$class, dnn=c("Predicted", "Actual")))

# Here we look at the confusion matrix and overall error rate from applying

# the tree rules to the training data.

predicted <- as.numeric(apply(predict(r), 1, function(r) order(-r)[1]))
(m <- cm (train$class, predicted))
1 - sum(diag(m)) / sum(m)


# Neural Network

# Recall from the heuristic on data mining and machine learning , a neural network is # a nonlinear model of complex relationships composed of multiple 'hidden' layers

# (similar to composite functions)

# Build the NNet model.

require(nnet, quietly=TRUE)

crs_nnet <- nnet(as.factor(class) ~ ., data=train, size=10, skip=TRUE, trace=FALSE, maxit=1000)

# Print the results of the modelling.

print(crs_nnet)
print("

Network Weights:

")
print(summary(crs_nnet))

# Evaluate Model Performance

# Generate an Error Matrix for the Neural Net model.

# Obtain the response from the Neural Net model.

crs_pr <- predict(crs_nnet, train, type="class")

# Now generate the error matrix.

table(crs_pr, train$class, dnn=c("Predicted", "Actual"))

# Generate error matrix showing percentages.

round(100*table(crs_pr, train$class, dnn=c("Predicted", "Actual"))/length(crs_pr))

Calucate overall error percentage.

print( "Overall Error Rate")

(function(x){ if (nrow(x) == 2) cat((x[1,2]+x[2,1])/sum(x)) else cat(1-(x[1,rownames(x)])/sum(x))}) (table(crs_pr, train$class, dnn=c("Predicted", "Actual")))