The PCA Procedure

This example illustrates that the iterative method based on random projection (RANDOM) of Halko, Martinsson, and Tropp (2011) provides a good approximation to the true principal components by using simulated data. The data for this example are the result of multiplying the low-rank matrices and and adding random noise. Matrix is 1,000,000 rows by 50 columns, and matrix is 50 rows by 5,000 columns, so the final data table consists of 1,000,000 observations and 5,000 variables.

The following DATA step generates the data:

data mycas.testdata / sessref=mysess single=no;
   title;
   keep x:;
   drop rank number_of_obs number_of_var sigma ii idum;
   drop rv1 rv2 rsq fac row col;
   drop nobs_per_thread nextras start_obs obs j k;
   array B[50,5000];          /* dimensions: rank, number_of_var */
   array A[50];               /* dimension: rank */
   array x[5000] x1-x5000;    /* dimension: number_of_var */

   target_nthreads = min(_nthreads_,96);
   if (_threadid_ = 1) then
      put "Number of threads = " target_nthreads;
   rank=50;
   number_of_obs=1000000;
   number_of_var=5000;
   sigma=0.1;

   call streaminit(1);

   if (_threadid_ <= target_nthreads) then do;
      ii = 0;
      idum = 0;
      do while (ii < rank * number_of_var);
         idum = mod(mod(1664525*idum,4294967296)+1013904223,4294967296);
         rv1 = 2.0*(idum/4294967296)-1.0;
         idum = mod(mod(1664525*idum,4294967296)+1013904223,4294967296);
         rv2 = 2.0*(idum/4294967296)-1.0;
         rsq = rv1*rv1+rv2*rv2;
         if ((rsq < 1.0) and (rsq ^= 0.0)) then do;
            fac = sqrt(-2.0*log(rsq)/rsq);
            row = int(ii/number_of_var)+1;
            col = mod(ii,number_of_var)+1;
            B[row,col] = rv1*fac;
            ii = ii + 1;
            if (ii < rank * number_of_var) then do;
               row = int(ii/number_of_var)+1;
               col = mod(ii,number_of_var)+1;
               B[row,col] = rv2*fac;
               ii = ii + 1;
            end;
         end;
      end;

      nobs_per_thread = int(number_of_obs /target_nthreads);
      nextras = number_of_obs - nobs_per_thread*target_nthreads;
      if (_threadid_ <= nextras) then do;
         nobs_per_thread = nobs_per_thread+1;
         start_obs = (_threadid_-1)*nobs_per_thread+1;
      end;
      else
         start_obs = nextras+(_threadid_-1)*nobs_per_thread+1;

      do obs = start_obs to (start_obs+nobs_per_thread-1);
         do j = 1 to rank;
            A[j] = rand('Normal');
         end;
         do k = 1 to number_of_var;
            x[k] = sigma*rand('Normal');
            do j = 1 to rank;
               x[k] = x[k] + A[j]*B[j,k];
            end;
         end;
         output;
      end;
   end;
run;

The following statements use PROC PCA to extract principal components by using the RANDOM and EIG methods. The suboptions for the RANDOM method specify 1, 5, or 10 iterations. Both the EIG and RANDOM methods extract 25 principal components. The DISPLAY statement suppresses the default output except for eigenvalues, and the DISPLAYOUT statement saves the computed eigenvalues to specific CAS output tables.

proc pca data=mycas.testdata n=25 method=random(niter=1);
   var x:;
   display Eigenvalues;
   displayout Eigenvalues=oneiter;
run;

proc pca data=mycas.testdata n=25 method=random(niter=5);
   var x:;
   display Eigenvalues;
   displayout Eigenvalues=fiveiter;
run;

proc pca data=mycas.testdata n=25 method=random(niter=10);
   var x:;
   display Eigenvalues;
   displayout Eigenvalues=teniter;
run;

proc pca data=mycas.testdata n=25 method=eig;
   var x:;
   display Eigenvalues;
   displayout Eigenvalues=trueeig;
run;

To assess the accuracy of the RANDOM method, compare the estimated eigenvalues to those computed by the EIG method. Using the DATA step and PROC SGPLOT (code not shown), you can combine the estimated eigenvalue data sets with the true eigenvalue data set and plot the results. Output 15.3.1 shows this composite plot and demonstrates the following heuristic rules for the RANDOM method:

For additional examples that show the impact of the number of iterations for both simulated and real data, see Halko, Martinsson, and Tropp (2011).

Output 15.3.1: RANDOM Method Performance with Varying Number of Iterations

Note: The heuristic rules provide guidance for selecting the number of iterations for data tables that are (approximately) low-rank. You might need to increase the number of iterations beyond what is recommended for data tables that are not low-rank.