reduced_rank_LDA = function( XTrain, yTrain, XTest, yTest ){
  #
  # R code to implement classification using Regularized Discriminant Analysis
  #
  # See the section with the same name as this function in Chapter 4 from the book ESLII
  #
  # Inputs:
  #   XTrain = training data frame
  #   yTrain = training labels of true classifications with indices 1 - K (where K is the number of classes)
  #   xTest = testing data frame
  #   yTest = testing response
  #
  # Written by:
  # -- 
  # John L. Weatherwax                2009-04-21
  # 
  # email: wax@alum.mit.edu
  # 
  # Please send comments and especially bug reports to the
  # above email address.
  # 
  #-----

  K = length(unique(yTrain)) # the number of classes (expect the classes to be labeled 1, 2, 3, ..., K-1, K 
  N = dim( XTrain )[1] # the number of samples
  p = dim( XTrain )[2] # the number of features
  
  # Compute the class dependent probabilities and class dependent centroids: 
  #
  PiK = matrix( data=0, nrow=K, ncol=1 )
  M = matrix( data=0, nrow=K, ncol=p )
  ScatterMatrices = list()
  for( ci in 1:K ){
    inds = yTrain == ci
    Nci = sum(inds)
    PiK[ci] = Nci/N
    M[ci,] = t( as.matrix( mean( XTrain[ inds, ] ) ) )
  }

  # Compute W:
  #
  W = cov( XTrain ) 

  # Compute M^* = M W^{-1/2} using the eigen-decomposition of W :
  #
  e = eigen(W)
  V = e$vectors # W = V %*% diag(e$values) %*% t(V)
  W_Minus_One_Half = V %*% diag( 1./sqrt(e$values) ) %*% t(V) 
  MStar = M %*% W_Minus_One_Half 

  # Compute B^* the covariance matrix of M^* and its eigen-decomposition:
  #
  BStar = cov( MStar )
  e = eigen(BStar)
  VStar = - e$vectors # note taking the negative to match the results in the book (results are independent of this)

  V = W_Minus_One_Half %*% VStar # the full projection matrix
  
  # Project the data into the invariant subspaces:
  #
  XTrainProjected = t( t(V) %*% t(XTrain) )
  XTestProjected = t( t(V) %*% t(XTest) )
  MProjected = t( t(V) %*% t(M) ) # the centroids projected

  # Classify the training/testing data for each possible projection dimension:
  # 
  TrainClassification = matrix( data=0, nrow=N, ncol=p ) # number of samples x number of projection dimensions 

  discriminant = matrix( data=0, nrow=1, ncol=K )
  for( si in 1:N ){ # for each sample
      for( pi in 1:p ){ # for each projection dimension 
	  for( ci in 1:K ){ # for each class centroid 
	      discriminant[ci] = 0.5 * sum( ( XTrainProjected[si,1:pi] - MProjected[ci,1:pi] )^2 ) - log( PiK[ci] )
	  }
          TrainClassification[si,pi] = which.min( discriminant )
      }
  } 

  N = dim(XTest)[1]
  TestClassification = matrix( data=0, nrow=N, ncol=p ) # number of samples x number of projection dimensions 

  discriminant = matrix( data=0, nrow=1, ncol=K )
  for( si in 1:N ){ # for each sample 
      for( pi in 1:p ){ # for each projection dimension 
	  for( ci in 1:K ){ # for each class centroid
	      discriminant[ci] = 0.5 * sum( ( XTestProjected[si,1:pi] - MProjected[ci,1:pi] )^2 ) - log( PiK[ci] )
	  }
          TestClassification[si,pi] = which.min( discriminant )
      }
  } 
  
  return( list(XTrainProjected,XTestProjected,MProjected,TrainClassification,TestClassification) )

}