src/test/scripts/applications/validation/CV_MultiClassSVM.dml

#-------------------------------------------------------------
#
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.
#
#-------------------------------------------------------------

# 50k Dataset w/ 4-fold CV w/ 2 classes
#    Total time:		3.7473180000000004 sec.
#    Number of executed MR Jobs:	0.
# hadoop jar SystemDS.jar -f CV_MultiClassSVM.dml -args itau/svm/X_50k_10 itau/svm/y_50k 4 0 2 0.001 1.0 100

# 100k Dataset w/ 4-fold CV w/ 2 classes
#    Total time:		206.549618 sec.
#    Number of executed MR Jobs:	16.
# hadoop jar SystemDS.jar -f CV_MultiClassSVM.dml -args itau/svm/X_100k_500 itau/svm/y_100k 4 0 2 0.001 1.0 100


X = read( $1 );
y = read( $2 );
m = nrow( X );
n = ncol( X );

k = $3;

#parameters for model training
intercept = $4
num_classes = $5
epsilon = $6
lambda = $7 
maxiter = $8

#CV
P = Rand(rows=m, cols=1, min=0.0, max=1.0, pdf = "uniform");
P = round(0.5+P*k);

ones = matrix(1, rows=1, cols=n);
onem = matrix(1, rows=1, cols=1);
stats = matrix(0, rows=k, cols=1); #k-folds x 1-stats
   
parfor( i in 1:k )
{
   #prepare train/test fold projections
   vPxi = (P == i);   #  Select 1/k fraction of the rows
   mPxi = (vPxi %*% ones);       #  for the i-th fold TEST set
   #nvPxi = (P != i);
   #nmPxi = (nvPxi %*% ones);  #note: inefficient for sparse data  

   #create train/test folds
   Xi = X * mPxi;  #  Create the TEST set with 1/k of all the rows
   yi = y * vPxi;  #  Create the labels for the TEST set
   nXi = X - Xi;   #  Create the TRAINING set with (k-1)/k of the rows
   nyi = y - yi;   #  Create the labels for the TRAINING set
   Xyi = cbind(Xi,yi); #keep alignment on removeEmpty
   Xyi = removeEmpty( target=Xyi, margin="rows" );
   Xi = Xyi[ , 1:n];
   yi = Xyi[ , n+1];   
   nXyi = cbind(nXi,nyi); #keep alignment on removeEmpty
   nXyi = removeEmpty( target=nXyi, margin="rows" );
   nXi = nXyi[ , 1:n];
   nyi = nXyi[ , n+1];

   #train multiclass SVM model per fold, use the TRAINING set
   wi = multiClassSVM( nXi, nyi, intercept, num_classes, epsilon, lambda, maxiter)

   #score multiclass SVM model per fold, use the TEST set
   out_correct_pct = scoreMultiClassSVM( Xi, yi, wi, intercept);
   
   stats[i,1] = out_correct_pct * onem;
}

# print output of stats
z = printFoldStatistics( stats );


################################################################################

printFoldStatistics = function( Matrix[double] stats)
   return( Integer err)
{
   mean_correct_pct = mean( stats[,1])

   print (" Mean Correct Percentage of the " + nrow( stats) + " Folds: " + mean_correct_pct);

   err = 0
}

################################################################################

scoreMultiClassSVM = function( Matrix[double] X, Matrix[double] y, Matrix[double] W, Integer intercept) 
   return (Double out_correct_pct)
{
   Nt = nrow(X);
   num_classes = ncol(W)
   b = Rand( rows=1, cols=num_classes, min=0, max=0, pdf="uniform")
   n = ncol(X);

   if (intercept == 1) 
   {
      b = W[n+1,]
   }
   
   ones = Rand( rows=Nt, cols=1, min=1, max=1, pdf="uniform")

   scores = X %*% W[1:n,] + ones %*% b;
   
   predicted_y = rowIndexMax( scores);
   
   correct_percentage = sum((predicted_y - y) == 0) / Nt * 100;

   out_correct_pct = correct_percentage;

}


################################################################################

multiClassSVM = function (Matrix[double] X, Matrix[double] Y, Integer intercept, Integer num_classes, Double epsilon, Double lambda, Integer max_iterations) 
   return (Matrix[double] ret_W) 
{
   check_X = sum(X)
   if(check_X == 0){

     print("X has no non-zeros")

   } else {

      num_samples = nrow(X)
      num_features = ncol(X)
      
      if (intercept == 1) {
        ones = Rand( rows=num_samples, cols=1, min=1, max=1, pdf="uniform");
        X = cbind( X, ones);
      }
      
      iter_class = 1
      
      Y_local = 2 * (Y == iter_class) - 1
      w_class = Rand( rows=num_features, cols=1, min=0, max=0, pdf="uniform")
   
      if (intercept == 1) {
         zero_matrix = Rand( rows=1, cols=1, min=0.0, max=0.0);
         w_class = t( cbind( t( w_class), zero_matrix));
      }
      
      g_old = t(X) %*% Y_local
      s = g_old
      iter = 0
      continue = 1
   
      while(continue == 1) {
        # minimizing primal obj along direction s
        step_sz = 0
        Xd = X %*% s
        wd = lambda * sum(w_class * s)
        dd = lambda * sum(s * s)
        continue1 = 1
        while(continue1 == 1){
         tmp_w = w_class + step_sz*s
         out = 1 - Y_local * (X %*% tmp_w)
         sv = (out > 0)
         out = out * sv
         g = wd + step_sz*dd - sum(out * Y_local * Xd)
         h = dd + sum(Xd * sv * Xd)
         step_sz = step_sz - g/h
         if (g*g/h < 0.0000000001){
          continue1 = 0
         }
        }
       
        #update weights
        w_class = w_class + step_sz*s
       
        out = 1 - Y_local * (X %*% w_class)
        sv = (out > 0)
        out = sv * out
        obj = 0.5 * sum(out * out) + lambda/2 * sum(w_class * w_class)
        g_new = t(X) %*% (out * Y_local) - lambda * w_class
      
        tmp = sum(s * g_old)
        
        train_acc = sum((Y_local*(X%*%w_class)) >= 0)/num_samples*100
        print("For class " + iter_class + " iteration " + iter + " training accuracy: " + train_acc)
         
        if((step_sz*tmp < epsilon*obj) | (iter >= max_iterations-1)){
         continue = 0
        }
       
        #non-linear CG step
        be = sum(g_new * g_new)/sum(g_old * g_old)
        s = be * s + g_new
        g_old = g_new
      
        iter = iter + 1
       }
      
      
      w = w_class
      iter_class = iter_class + 1
      
      while(iter_class <= num_classes){
       Y_local = 2 * (Y == iter_class) - 1
      # w_class = Rand(rows=num_features, cols=1, min=0, max=0, pdf="uniform")
       w_class = Rand(rows=ncol(X), cols=1, min=0, max=0, pdf="uniform")
       if (intercept == 1) {
       	zero_matrix = Rand(rows=1, cols=1, min=0.0, max=0.0);
       	w_class = t(cbind(t(w_class), zero_matrix));
       }
       
       g_old = t(X) %*% Y_local
       s = g_old
      
       iter = 0
       continue = 1
       while(continue == 1)  {
        # minimizing primal obj along direction s
        step_sz = 0
        Xd = X %*% s
        wd = lambda * sum(w_class * s)
        dd = lambda * sum(s * s)
        continue1 = 1
        while(continue1 == 1){
         tmp_w = w_class + step_sz*s
         out = 1 - Y_local * (X %*% tmp_w)
         sv = (out > 0)
         out = out * sv
         g = wd + step_sz*dd - sum(out * Y_local * Xd)
         h = dd + sum(Xd * sv * Xd)
         step_sz = step_sz - g/h
         if (g*g/h < 0.0000000001){
          continue1 = 0
         }
        }
       
        #update weights
        w_class = w_class + step_sz*s
       
        out = 1 - Y_local * (X %*% w_class)
        sv = (out > 0)
        out = sv * out
        obj = 0.5 * sum(out * out) + lambda/2 * sum(w_class * w_class)
        g_new = t(X) %*% (out * Y_local) - lambda * w_class
      
        tmp = sum(s * g_old)
        
        train_acc = sum((Y_local*(X%*%w_class)) >= 0)/num_samples*100
        print("For class " + iter_class + " iteration " + iter + " training accuracy: " + train_acc)
         
        if((step_sz*tmp < epsilon*obj) | (iter >= max_iterations-1)){
         continue = 0
        }
       
        #non-linear CG step
        be = sum(g_new * g_new)/sum(g_old * g_old)
        s = be * s + g_new
        g_old = g_new
      
        iter = iter + 1
       }
      
       w = cbind(w, w_class) 
       iter_class = iter_class + 1
      }
      ret_W = w
   }
}