Kernel Principal Component Analysis Action Set

This section contains PROC CAS code.

Note: Input data must be accessible in your CAS session, either as a CAS table or as a transient-scope table. A CAS table has a two-level name: the first level is your CAS engine libref, and the second level is the table name. You refer to this table in the CAS procedure by specifying only the second level. For more information about two-level names, see Chapter 2, Shared Concepts (SAS Visual Data Mining and Machine Learning: Procedures). A transient-scope table is called directly from the action and exists in memory for the duration of the action. For more information about accessing data, see SAS Viya: System Programming Guide. For more information about PROC CAS and programming in CASL, see SAS Cloud Analytic Services: CASL Programmer’s Guide and SAS Cloud Analytic Services: CASL Reference.

This example shows how to use the kPca action to get the pre-image of the noisy USPS digit data. It contains -pixel labeled images of handwritten digits scanned by the US Postal Service. There are a total of 9,298 digits, and each digit is represented by a 256-dimensional real value vector and the corresponding integer label. For each of the 10 digits, 300 examples were randomly selected for training and 50 examples for testing. The test data were rendered noisy by adding Gaussian noise, which is characterized by a zero mean and a standard deviation of 0.5.

The following code first uses the IML procedure to randomly select the training and test data sets and then plots the second appearance of each digit in the test data set, with and without Gaussian noise, respectively:

proc iml;
option source2;
ods graphics on/ reset=all;
title;
/* read variables from a clean USPS digit data set into a matrix */
use usps_clean;
read all var _all_ into digitClean;

/* flip pixel values so that the digit is black and background is white */

digitClean[,2:257]=-digitClean[,2:257];

/* read variables from a noisy USPS digit data set into a matrix */
use usps_noise;
read all var _all_ into digitNoise;

/* flip pixel value so that the digit is black and background is white */
digitNoise[,2:257]=-digitNoise[,2:257];

/* define a utility function that plots the digit data*/
start implot_digit(image_m, data_m, coln);
   n=nrow(data_m);
   do i=1 to n;
     row=shape(data_m[i,1:coln-1],16,16);
      m_sub=j(16##2,4);
       do j=1 to 16;
         do k=1 to 16;
           m_sub[k+16*(j-1),1]=(17-j);
           m_sub[k+16*(j-1),2]=k;
           m_sub[k+16*(j-1),3]=row[j,k];
           m_sub[k+16*(j-1),4]=i;
         end;
       end;
     image_m=image_m//m_sub;
  end;
finish;



/* get the features for clean digits (first column is observation ID) */
digitClean=digitClean[,2:ncol(digitClean)];

/* get the features for noisy digits (first column is observation ID) */
digitNoise=digitNoise[,2:ncol(digitNoise)];

/* row and column number of the digit data
   (noise and clean digit data sets have same dimension) */
coln=ncol(digitClean);
rown=nrow(digitClean);

/* allocate space for noisy test data */
digitNoiseTest={};

/* allocate space for clean test data*/
digitCleanTest={};

/* allocate space for clean training data*/
digitTrain={};


/* random seed for generating the training and the test data set*/
call randseed(54323,1);

/* for each digit, choose 300 samples for the training data
   and 50 samples for the test data */
do i=0 to 9;
   s = sample(loc(digitClean[,coln]=i), 350);
   /* last column is digit label */
   digitCleanTest=digitCleanTest//digitClean[s[301:350],1:coln-1];
   digitNoiseTest=digitNoiseTest//digitNoise[s[301:350],1:coln-1];
   digitTrain=digitTrain//digitClean[s[1:300],];
end;


/* digit labels */
label={0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

/* second appearance of noisy test digits */
digitNoiseTest_one=digitNoiseTest[{2,52,102,152,202,252,302,352,402,452},];
digitNoiseTest_one=digitNoiseTest_one||label;

/* second appearance of clean test digits */
digitCleanTest_one=digitCleanTest[{2,52,102,152,202,252,302,352,402,452},];
digitCleanTest_one=digitCleanTest_one||label;


/* Convert the digit pixel matrix to SAS data table of the form that
   heatmap statement can accept */
run implot_digit(digitNoiseTest_trans, digitNoiseTest_one, coln);
create digitNoiseTest_trans from digitNoiseTest_trans;
append from digitNoiseTest_trans;
close digitNoiseTest_trans;


run implot_digit(digitCleanTest_trans, digitCleanTest_one, coln);
create digitCleanTest_trans from digitCleanTest_trans;
append from digitCleanTest_trans;
close digitCleanTest_trans;

/* save digitTrain to SAS data table */
create digitTrain from digitTrain;
append from digitTrain;
close digitTrain;

/* save digitNoiseTest to SAS data table */
create digitNoiseTest from digitNoiseTest;
append from digitNoiseTest;
close digitNoiseTest;

/* save digitCleanTest to SAS data table */
create digitCleanTest from digitCleanTest;
append from digitCleanTest;
close digitCleanTest;

/* plot test digits with Gaussian noise*/
ods graphics /width=8in height=1.3in;
proc sgpanel data=digitNoiseTest_trans;
   panelby COL4/noborder columns=10 proportional ;
   heatmap x=COL2 y=COL1 /colorresponse=COL3 colormodel=(black white)
           showybins showxbins xbinstart=0 ybinstart=0;
run;

Figure 5: Plot of Noisy USPS Digit Data

/* plot clean test digits without Gaussian noise */
ods graphics /width=8in height=1.3in;
proc sgpanel data=digitCleanTest_trans;
   panelby COL4/noborder columns=10 proportional ;
   heatmap x=COL2 y=COL1 /colorresponse=COL3 colormodel=(black white)
           showybins showxbins xbinstart=0 ybinstart=0;
run;

Figure 6: Plot of Clean USPS Digit Data

In Figure 5 and Figure 6, noisy digits are blurry and clean digits are sharp and clear. Next, run the kPca action with the mapping pre-image method specified. In the following code, the RBF kernel is used in both KPCA (kerParam parameter value is 8) and kernel ridge regression (kerParam parameter value is 1), and 80 principal components are used for KPCA projection. When the mapping pre-image training is completed, the aStore action is called to compute the pre-image score of the noisy test digits. To get a sense of how well the mapping pre-image method performs, you can plot the pre-image of the second appearance of all 10 digits.

data mycas.digitTrain;
  set digitTrain;
run;

data mycas.digitNoiseTest;
  set digitNoiseTest;
run;

proc cas;
    loadactionset "KernelPCA";
    run;
    action kPca result = r/
                 table  = {name='digitTrain'},
                 saveState={name="STATE", replace=TRUE},
                 input=${COL1-COL256},
                 method = "EXACT",
                 kerType="RBF",
                 kerParam=8,
                 preimageNPC=80,
                 mapCoeffs={name="COEFFS", replace=TRUE},
                 outputTables={names="TaskTiming", repeated=TRUE},
                 preimage=True,
                 preimageMethod="MAP",
                 mapParam={kerType='RBF',kerParam=1,lambda=1};
   run;
   print r;
quit;

data mycas.digitNoiseTest;
   set digitNoiseTest;
   id =_N_;
run;

proc cas  ;
    action aStore.score / table={name='digitNoiseTest'},
    options={{name="scoring_mode",value=1}},
    out={name='preimage_results', replace=true},
    rstore={name='state'},
    copyVars={'id'};
run;
quit;

data preimage_results;
    set mycas.preimage_results;
run;

proc sort data=preimage_results;
    by id;
run;

proc iml;

/* read pre-image results into an iml matrix */
use preimage_results;
read all var _all_ into preimage_results;

/* The second appearance of the pre-image of noisy test digits */
preimage_results=preimage_results[{2,52,102,152,202,252,302,352,402,452},];
preimage_results=preimage_results||label;

coln=ncol(preimage_results);
run implot_digit(preimage_results_trans, preimage_results, coln);
create preimage_results_trans from preimage_results_trans;
append from preimage_results_trans;
close preimage_results_trans;


/* plot pre-image of noisy digits (denoising) */

ods graphics /width=8in height=1.3in;
proc sgpanel data=preimage_results_trans;
   panelby COL4/noborder columns=10 proportional ;
   heatmap x=COL2 y=COL1 /colorresponse=COL3 colormodel=(black white)
           showybins showxbins xbinstart=0 ybinstart=0;
run;

Figure 7: Plot of Denoised USPS Digits

In Figure 7, you can see that the denoised digits are much clearer and easier to identify than the original noisy digits shown in Figure 5.

The following code implements the approximate mapping pre-image method by specifying the APPROXIMATE value in the method parameter. Here the k-means clustering uses KMPP (k-means plus plus) as the clustering centroids initialization method and uses the default number of centroids (100).

/* Run mapping pre-image method using Nystrom approximation with 100 centroids */
proc cas;
     loadactionset "KernelPCA";
    run;
    action kPca result = r/
                table  = {name='digitTrain'},
                saveState={name="STATE", replace=TRUE},
                input=${COL1-COL256},
                method = "APPROXIMATE",
                clusMethod="KMPP",
                kerType="RBF",
                kerParam=8,
                preimageNPC=80,
                mapCoeffs={name="COEFFS", replace=TRUE},
                outputTables={names="TaskTiming", repeated=TRUE},
                preimage=True,
                preimageMethod="MAP",
                mapParam={kerType='RBF',kerParam=1,lambda=1};
   run;
   print r;
quit;


data mycas.digitNoiseTest;
   set digitNoiseTest;
   id =_N_;
run;


proc cas;
    action aStore.score / table={name='digitNoiseTest'},
    options={{name="scoring_mode",value=1}},
    out={name='preimage_results_approx', replace=true},
    rstore={name='state'},
    copyVars={'id'};
run;
quit;

data preimage_results_approx;
    set mycas.preimage_results_approx;
run;

proc sort data=preimage_results_approx;
    by id;
run;

/* read the approximated pre-image results data into an iml matrix */
proc iml;

use preimage_results_approx;
read all var _all_ into preimage_results_approx;

/* second appearance of pre-image of noisy test digits */
preimage_results_approx=preimage_results_approx[{2,52,102,152,202,252,302,352,402,452},];
preimage_results_approx=preimage_results_approx||label;


coln=ncol(preimage_results_approx);
run implot_digit(preimage_results_approx_trans, preimage_results_approx, coln);
create preimage_results_approx_trans from preimage_results_approx_trans;
append from preimage_results_approx_trans;
close preimage_results_approx_trans;


/* plot the approximated pre-image of noisy digits (denoising) */


ods graphics /width=8in height=1.3in;
proc sgpanel data=preimage_results_approx_trans;
  panelby COL4/noborder columns=10 proportional ;
  heatmap x=COL2 y=COL1 /colorresponse=COL3 colormodel=(black white)
          showybins showxbins xbinstart=0 ybinstart=0;
run;


quit;

Figure 8: Plot of Denoised USPS Digits with Approximate Mapping Method

Figure 8 shows that the denoised digits that are obtained by the approximate mapping pre-image method are similar to those obtained by the exact method and shown in Figure 7.

In terms of the amount of time the two methods take, the training time for the approximate mapping method is much less than the training time for the exact mapping method, as shown in the "Task Timing" tables in Output 21.4.1 and Output 21.4.2, respectively. This demonstrates the efficiency of the approximate mapping method in reducing the training time while not compromising much in denoising performance. In addition to that improvement, in the exact mapping pre-image method all the training data must be stored in the state file to get the pre-image scores of the test data, whereas in the approximate mapping pre-image method only the k-means centroids of the KPCA projection of the training data must be stored in the state file. This significantly reduces the time required for creating the scoring state from the state file.