KPCA Procedure

This example shows how to use the KPCA procedure 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 value 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 values 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;

Output 18.3.1: 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;

Output 18.3.2: Plot of Clean USPS Digit Data

In Output 18.3.1 and Output 18.3.2, noisy digits are blurry and clean digits are sharp and clear.

Next, run the KPCA procedure with the mapping pre-image method specified. In the following code, the RBF kernel is used in both KPCA (BW=8) and kernel ridge regression (BW=1), and 80 principal components are used for KPCA projection. When the mapping pre-image training is completed, the ASTORE procedure 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 mylib.digitTrain;
  set digitTrain;
run;

data mylib.digitNoiseTest;
  set digitNoiseTest;
run;

proc KPCA data=mylib.digitTrain method=exact;
  input COL1-COL256;
  kernel RBF/bw=8;
  preimage npc=80 method=map(L2=1 kernel=RBF(bw=1));
  savestate rstore=mylib.state;
  displayout modelinfo="TaskTiming";
run;


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

proc ASTORE;
    setoption scoring_mode 1;
    score data=mylib.digitNoiseTest rstore=mylib.state out=mylib.preimage_results
               copyvars=id;
run;
quit;

data preimage_results;
    set mylib.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;

/* 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;

Output 18.3.3: Plot of Denoised USPS Digits

In Output 18.3.3, you can see that the denoised digits are much clearer and easier to identify than the original noisy digits shown in Output 18.3.1.

The following code implements the approximate mapping pre-image method by specifying METHOD=APPROXIMATE. 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 KPCA data=mylib.digitTrain method=approximate;
  input COL1-COL256;
  kernel RBF/bw=8;
  lrapproximation clusmethod = KMPP;
  preimage npc=80 method=map(L2=1 kernel=RBF(bw=1));
  savestate rstore=mylib.state;
  displayout modelinfo="TaskTiming";
run;


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

proc ASTORE;
    setoption scoring_mode 1;
    score data=mylib.digitNoiseTest rstore=mylib.state
          out=mylib.preimage_results_approx copyvars=id;
run;
quit;


data preimage_results_approx;
    set mylib.preimage_results_approx;
run;

proc sort data=preimage_results_approx;
    by id;
run;

/* Read approximate pre-image results data into iml matrix */
proc iml;

use preimage_results_approx;
read all var _all_ into preimage_results_approx;

/* second appearance of the 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;

Output 18.3.4: Plot of Denoised USPS Digits with Approximate Mapping Method

Output 18.3.4 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 Output 18.3.3.

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 18.3.5 and Output 18.3.6, 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.