DYNAMICLINEAR Procedure

Example 17.2 Comparing CPU and GPU Workflows in Simultaneous Graphical Dynamic Linear Models

This example demonstrates how to use the GPU statement in PROC DYNAMICLINEAR to accelerate the filtering and forecasting of a large-scale simultaneous graphical dynamic linear model (SGDLM). It compares the Kullback-Leibler (KL) divergence and effective sample size of the CPU and graphics processing unit (GPU) workflows, as well as the mean, median, and variance predictions for one of the variables of an SGDLM that has 300 time series () and 300 observations (), where each variable has 10 elements in its parental set. The example is similar to Tracking Volatility with Simultaneous Graphical Dynamic Linear Models. Note: Because of the difference in implementation, the GPU workflow can generate results that differ from the CPU workflow results, even when using the same random seed.

Using a data generation process similar to the one in Example 17.1, the following statements simulate 300 observations for 300 time series. As in Example 17.1, the covariance matrix is selected as the initialization covariance matrix for all variables in the SGDLM:

*macro variables;
%let seed = 12345;
%let m = 300;
%let T = 300;
%let nparent = 10;
%let sigma_phi = 0.0001;
%let sigma_gamma = 0.01;
%let nSim = 10000;
%let lead = 1;

proc iml;
   call randseed(&seed.);
   *define the parental sets;
   nparent_j = j(&m, 1, 0);
   do h = 1 to &m.;
      nparent_j[h] = &nparent.;
   end;

   parents = j(&m., &m., 0);
   all_ind = 1:&m.;
   do h = 1 to &m.;
      p_gam = nparent_j[h];
      if p_gam > 0 then do;
         all_ind_j = all_ind[loc(all_ind ^= h)];
         parent_ju = sample(all_ind_j, p_gam, "NoReplace");
         parent_ju_c = parent_ju`;
         call sort(parent_ju_c);
         parent_j = parent_ju_c`;
         parents[h, parent_j] = 1;
      end;
   end;

   parentsc = {};
   all_ind = 1:&m.;
   do h = 1 to &m.;
      parents_h = loc(parents[h,] = 1);
      npar_h = ncol(parents_h);
      if npar_h > 0 then do;
         do k = 1 to npar_h;
            pjk = h || parents_h[k];
            parentsc = parentsc // pjk;
         end;
      end;
   end;

   *generate data;
   y = j(&T, &m. + 1);
   do t = 1 to &T;
      gamma_t = j(&m, &m, 0);
      all_ind = 1:&m.;
      do h = 1 to &m.;
         p_gam = nparent_j[h];
         if p_gam > 0 then do;
            parent_j = all_ind[loc(parents[h,] = 1)];
            gamj = j(p_gam, 1);
            call randgen(gamj, "normal");
            gamma_t[h, parent_j`] = &sigma_gamma. * gamj`;
         end;
      end;

      *generate mut;
      mu_ts = j(&m, 1);
      call randgen(mu_ts, "Normal");
      mu_t = &sigma_phi. * mu_ts;

      *generate Lambdat;
      lambda_t = j(&m, 1);
      call randgen(lambda_t, "Gamma", 1);
      lambda_t = inv(diag(lambda_t));

      *compute A_t, Sigma_t;
      A_t = inv(I(&m) - gamma_t);
      mean_t = A_t * mu_t;
      Sigma_t = A_t * inv(lambda_t) * A_t`;
      /* the Cholesky root of the full covariance matrix */
      L_t = root(Sigma_t);
      z_t = j(&m, 1);
      call randgen(z_t, "Normal"); /* z is MVN(0,I(&m)) */
      y_t = mean_t + L_t` * z_t;
      y[t, 1:&m.] = y_t`;
      y[t, &m. + 1] = t;
   end;

   yr = y[loc(y[, &m. + 1] <= &T - 1), ];

   yvarNames = "y1":"y&m.";
   varNames = yvarNames || "t";
   create data_correlated from y [colname=varNames];
   append from y;
   close data_correlated;

   create data_correlatedr from yr [colname=varNames];
   append from yr;
   close data_correlatedr;

   create parents from parents [colname=yvarNames];
   append from parents;
   close parents;

   varNamesc = {'ParentID', 'ChildID'};
   create adj_table from parentsc [colname=varNamesc];
   append from parentsc;
   close adj_table;

   *generate initial covariance data programmatically;
   nParmPerVar = 1 + &nparent.;
   nRows = &m. * nParmPerVar;
   init_cov = j(nRows, 3, 0);
   row = 0;
   do v = 1 to &m.;
      do p = 1 to nParmPerVar;
         row = row + 1;
         init_cov[row, 1] = v;
         init_cov[row, 2] = p;
      if p = 1 then
         init_cov[row, 3] = 0.0001;
      else
         init_cov[row, 3] = 0.01;
      end;
   end;

   covNames = {'VarID', 'ParmID', 'ParmValue'};
   create initial_cov_with_adj from init_cov [colname=covNames];
   append from init_cov;
   close initial_cov_with_adj;
run;

proc sort data=data_correlatedr;
   by t;
run;

The following statements load the generated data set data_correlatedr, the defined adjacency matrix data set adj_table, and the initial covariance data set initial_cov_with_adj into data tables in your session that is associated with the mylib libref. The DATA step assumes that your libref is named mylib, but you can substitute any appropriately defined libref.

data mylib.test_data;
   set data_correlatedr;
run;

data mylib.adj_table;
   set adj_table;
run;

data mylib.initial_cov_with_adj;
   set initial_cov_with_adj;
run;

The following statements run the CPU workflow of PROC DYNAMICLINEAR on the simulated data set:

proc dynamicLinear data=mylib.test_data nSimulations=10000 seed=12345;
   id t;
   model y1 - y300;
   parentalSet inAdjacency=mylib.adj_table;
   filter out(percentiles=(2.5, 25, 50, 75, 97.5))=mylib.out_filter_cpu
      outFilterInfo=(klDivergence=mylib.kldiv_cpu)
      back=0;
   forecast out(percentiles=(2.5, 25, 50, 75, 97.5))=mylib.out_forecast_cpu;
   initialDistribution inCovariance(diagonal)=mylib.initial_cov_with_adj
      shape=(default=5) rate=(default=1.5);
   discountFactor beta=(default=0.98) deltaPhi=(default=0.98)
      deltaGamma=(default=0.99);
run;

The following statements run the GPU workflow. The only difference from the CPU workflow is the addition of the GPU statement, which directs PROC DYNAMICLINEAR to perform the computations on GPU hardware. By default, the GPU calculations use single-precision floating-point arithmetic, which typically increases the speed significantly while maintaining sufficient accuracy for the SGDLM:

proc dynamicLinear data=mylib.test_data nSimulations=10000 seed=12345;
   id t;
   model y1 - y300;
   parentalSet inAdjacency=mylib.adj_table;
   filter out(percentiles=(2.5, 25, 50, 75, 97.5))=mylib.out_filter_gpu
      outFilterInfo=(klDivergence=mylib.kldiv_gpu)
      lead=1 back=0;
   forecast out(percentiles=(2.5, 25, 50, 75, 97.5))=mylib.out_forecast_gpu
      lead=1;
   initialDistribution inCovariance(diagonal)=mylib.initial_cov_with_adj
      shape=(default=5) rate=(default=1.5);
   discountFactor beta=(default=0.98) deltaPhi=(default=0.98)
      deltaGamma=(default=0.99);
   gpu;
run;

Comparing KL Divergence

The following statements compare the KL divergence and effective sample size of the CPU and GPU workflows. The plot is shown in Output 17.2.1.

data kldiv_cpu;
   set mylib.kldiv_cpu;
run;
data kldiv_gpu;
   set mylib.kldiv_gpu;
run;

proc sort data=kldiv_cpu; by t; run;
proc sort data=kldiv_gpu; by t; run;

data kldiv_comb;
   merge kldiv_cpu(rename=(KLDivergence=KLDiv_CPU
                           EffectiveSampleSize=ESS_CPU))
         kldiv_gpu(rename=(KLDivergence=KLDiv_GPU
                           EffectiveSampleSize=ESS_GPU));
   by t;
run;

proc sgplot data=kldiv_comb;
   series x=t y=KLDiv_CPU / legendlabel='KL Divergence (CPU)'
      lineattrs=(color=blue pattern=solid thickness=2);
   series x=t y=KLDiv_GPU / legendlabel='KL Divergence (GPU)'
      lineattrs=(color=red pattern=shortdash thickness=2);
   series x=t y=ESS_CPU / legendlabel='Effective Sample Size (CPU)'
      lineattrs=(color=cyan pattern=solid) Y2Axis;
   series x=t y=ESS_GPU / legendlabel='Effective Sample Size (GPU)'
      lineattrs=(color=orange pattern=shortdash) Y2Axis;
   yaxis label="KL Divergence";
   xaxis label="Step";
run;

Output 17.2.1: KL Divergence and Effective Sample Size: CPU vs. GPU

Output 17.2.1 shows that the GPU workflow produces better results than the CPU workflow in terms of KL divergence and effective sample size for the problem that has the specified parameters. The KL divergence for the GPU workflow is consistently lower than that of the CPU workflow, indicating that the GPU-accelerated computations are providing a closer approximation to the true posterior distribution. Additionally, the effective sample size of the GPU workflow is higher than that of the CPU workflow, suggesting that the GPU computations yield more efficient sampling and better convergence properties than the CPU computations.

Comparing Mean and Median Predictions

The following statements compare the 60-step moving average of the actual values of y1 to the mean and median predictions from the CPU and GPU workflows. The plots are shown in Output 17.2.2 and Output 17.2.3.

%let maWindow = 60;
%let varName = y1;

data out_filter_cpu;
   set mylib.out_filter_cpu;
run;

data out_filter_gpu;
   set mylib.out_filter_gpu;
run;

proc expand data=data_correlatedr out=test_data_ma method=none;
   id t;
   convert &varName. = var_ma / transout=(movave &maWindow.);
run;

data mean_comb;
   format t var_ma var_cpu var_gpu Best12.;
   set test_data_ma(keep=t var_ma where=(t >= &maWindow.));
   set out_filter_cpu(keep=t Mean_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Mean_&varName.=var_cpu));
   set out_filter_gpu(keep=t Mean_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Mean_&varName.=var_gpu));
   t = t + 1;
run;

proc sort data=mean_comb;
   by t;
run;

proc sgplot data=mean_comb;
   series x=t y=var_ma / legendlabel="Moving Average Values of &varName."
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=var_cpu / legendlabel="Mean_&varName. (CPU)"
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=var_gpu / legendlabel="Mean_&varName. (GPU)"
      markers markerattrs=(color=green size=5 symbol=diamond)
      lineattrs=(color=green pattern=longdash);
   label t="Step" var_ma="Values of the Means";
run;

Output 17.2.2: Comparison of the Mean Predictions: CPU vs. GPU

proc expand data=data_correlatedr out=test_data_med method=none;
   id t;
   convert &varName. = var_med / transout=(MOVMED &maWindow.);
run;

data median_comb;
   format t var_med var_cpu var_gpu Best12.;
   set test_data_med(keep=t var_med where=(t >= &maWindow.));
   set out_filter_cpu(keep=t P50_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(P50_&varName.=var_cpu));
   set out_filter_gpu(keep=t P50_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(P50_&varName.=var_gpu));
   t = t + 1;
run;

proc sort data=median_comb;
   by t;
run;

proc sgplot data=median_comb;
   series x=t y=var_med / legendlabel="Moving Median of &varName."
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=var_cpu / legendlabel="Median_&varName. (CPU)"
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=var_gpu / legendlabel="Median_&varName. (GPU)"
      markers markerattrs=(color=green size=5 symbol=diamond)
      lineattrs=(color=green pattern=longdash);
   label t="Step" var_med="Values of the Medians";
run;

Output 17.2.3: Comparison of the Median Predictions (y1): CPU vs. GPU

Output 17.2.2 and Output 17.2.3 show that the CPU and GPU workflows produce mean and median predictions that share similar patterns for the variable y1. Both workflows track the moving average and moving median trends of the actual data. However, the CPU workflow has a longer burn-in period, which is reflected by the larger deviations from the moving average and moving median values in the early steps than those of the GPU workflow. The GPU workflow converges more quickly to the underlying trends, demonstrating that the GPU-accelerated computations provide more efficient sampling and better convergence properties than the CPU computations.

Comparing Variance Predictions

The following statements compare the 60-step moving average variance of the actual values of y1 to the variance predictions from the CPU and GPU workflows. The plot is shown in Output 17.2.4.

proc expand data=data_correlatedr out=test_data_ma_var method=none;
   id t;
   convert &varName. = var_ma_var / transout=(MOVVAR &maWindow.);
run;

data var_comb;
   format t var_ma_var var_cpu var_gpu Best12.;
   set test_data_ma_var(keep=t var_ma_var where=(t >= &maWindow.));
   set out_filter_cpu(keep=t Variance_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Variance_&varName.=var_cpu));
   set out_filter_gpu(keep=t Variance_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Variance_&varName.=var_gpu));
   t = t + 1;
run;

proc sort data=var_comb;
   by t;
run;

proc sgplot data=var_comb;
   series x=t y=var_ma_var / legendlabel="Moving Average Variances of &varName."
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=var_cpu / legendlabel="Variance_&varName. (CPU)"
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=var_gpu / legendlabel="Variance_&varName. (GPU)"
      markers markerattrs=(color=green size=5 symbol=diamond)
      lineattrs=(color=green pattern=longdash);
   label t="Step" var_ma_var="Values of the Variances";
run;

Output 17.2.4: Comparison of the Variance Predictions (y1): CPU vs. GPU

Output 17.2.4 shows that the GPU workflow provides variance predictions that closely agree with those of the CPU workflow. Both workflows track the underlying moving variance of y1. However, the results of the CPU workflow are slightly more volatile than those of the GPU workflow, especially in the early steps. This is consistent with the observations in Output 17.2.2 and Output 17.2.3 that the CPU workflow has a longer burn-in period than the GPU workflow.

Comparing Performance

In addition to providing greater accuracy, the GPU workflow also provides significantly better performance than the CPU workflow. In this example, the GPU workflow reduces the computation time from about four hours to about five minutes, demonstrating the substantial efficiency gains (more than 40 times faster) that can be achieved by GPU acceleration for large-scale SGDLMs. The example ran on a machine with an NVIDIA L40 48 GB VRAM GPU and a machine with an Intel Xeon Platinum 8358 CPU at 2.60 GHz.

Output 17.2.5: Total Run Time: CPU vs. GPU (Speedup: 42.1x)

If you take a closer look at the step-by-step run time comparison of the CPU and GPU workflows, you can see that the GPU workflow provides a significant increase in speed for each step of the computations, especially the prediction generation step, which is typically the most computationally intensive step of the SGDLM algorithm. The GPU workflow runs more than 60 times faster for the posterior updating step than the CPU workflow, substantially contributing to the overall performance improvement.

Output 17.2.6: Step-by-Step Run Time Comparison: CPU vs. GPU

Step	CPU Avg (s)	GPU Avg (s)	Speedup
computePrior	0.0044	0.00012	38.1
computeForecast	41.2756	0.64029	64.5
computePosterior	0.0017	0.00052	3.2
computeVBPosterior	4.2127	0.17667	23.8
Total	46.1055	1.09045	42.3

Comparing Precisions: FP64 vs. FP32

You can run the following statements to execute the GPU workflow with FP64 precision. The only difference from the previous GPU workflow is the addition of the PRECISION=FP64 option in the GPU statement, which directs PROC DYNAMICLINEAR to use double-precision floating-point arithmetic for the computations on GPUs.

proc dynamicLinear data=mylib.test_data nSimulations=10000 seed=12345;
   id t;
   model y1 - y300;
   parentalSet inAdjacency=mylib.adj_table;
   filter out(percentiles=(2.5, 25, 50, 75, 97.5))=mylib.out_filter_gpu2
      outFilterInfo=(klDivergence=mylib.kldiv_gpu2)
      lead=1 back=0;
   forecast out(percentiles=(2.5, 25, 50, 75, 97.5))=mylib.out_forecast_gpu2
      lead=1;
   initialDistribution inCovariance(diagonal)=mylib.initial_cov_with_adj
      shape=(default=5) rate=(default=1.5);
   discountFactor beta=(default=0.98) deltaPhi=(default=0.98)
      deltaGamma=(default=0.99);
   gpu precision=FP64;
run;

Along with the memory savings, FP32 precision can deliver significantly greater speed than FP64 precision on the NVIDIA L40S with 48 GB VRAM.

Output 17.2.7: Total Run Time: GPU FP32 vs. GPU FP64 (Speedup: 1.6x)

Meanwhile, FP32 precision provides accuracy comparable to FP64 precision for the SGDLM in this example. The following statements produce plots that compare the KL divergence and effective sample size (Output 17.2.8), means (Output 17.2.9), and variances (Output 17.2.10) of the GPU workflow predictions for each precision setting:

data kldiv_gpu1;
   set mylib.kldiv_gpu;
run;

proc sort data=kldiv_gpu1; by t; run;
proc sort data=kldiv_gpu2; by t; run;

data kldiv_comb;
   merge kldiv_gpu1(rename=(KLDivergence=KLDiv_GPU1
                           EffectiveSampleSize=ESS_GPU1))
         kldiv_gpu2(rename=(KLDivergence=KLDiv_GPU2
                           EffectiveSampleSize=ESS_GPU2));
   by t;
run;

proc sgplot data=kldiv_comb;
   series x=t y=KLDiv_GPU1 / legendlabel="KL Divergence (GPU FP32)"
      lineattrs=(color=blue pattern=solid thickness=2);
   series x=t y=KLDiv_GPU2 / legendlabel="KL Divergence (GPU FP64)"
      lineattrs=(color=red pattern=shortdash thickness=2);
   series x=t y=ESS_GPU1 / legendlabel="Effective Sample Size (GPU FP32)"
      lineattrs=(color=cyan pattern=solid) Y2Axis;
   series x=t y=ESS_GPU2 / legendlabel="Effective Sample Size (GPU FP64)"
      lineattrs=(color=orange pattern=shortdash) Y2Axis;
   yaxis label="KL Divergence";
   y2axis label="Effective Sample Size";
   xaxis label="Step";
run;

Output 17.2.8: KL Divergence and Effective Sample Size: GPU FP32 vs. GPU FP64

data out_filter_gpu1;
   set mylib.out_filter_gpu;
run;

data out_filter_gpu2;
   set mylib.out_filter_gpu_fp64;
run;

data mean_comb;
   format t var_gpu1 var_gpu2 Best12.;
   set test_data_ma(keep=t var_ma where=(t >= &maWindow.));
   set out_filter_gpu1(keep=t Mean_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Mean_&varName.=var_gpu1));
   set out_filter_gpu2(keep=t Mean_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Mean_&varName.=var_gpu2));
   t = t + 1;
run;

proc sort data=mean_comb;
   by t;
run;

proc sgplot data=mean_comb;
   series x=t y=var_ma / legendlabel="Moving Average Values of &varName."
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=var_gpu1 / legendlabel="Mean_&varName. (GPU FP32)"
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=var_gpu2 / legendlabel="Mean_&varName. (GPU FP64)"
      markers markerattrs=(color=green size=5 symbol=diamond)
      lineattrs=(color=green pattern=longdash);
   label t="Step" var_ma="Values of the Means";
run;

Output 17.2.9: Comparison of Mean Predictions: GPU FP32 vs. GPU FP64

data var_comb;
   format t var_ma_var var_gpu1 var_gpu2 Best12.;
   set test_data_ma_var(keep=t var_ma_var where=(t >= &maWindow.));
   set out_filter_gpu1(keep=t Variance_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Variance_&varName.=var_gpu1));
   set out_filter_gpu2(keep=t Variance_&varName. where=(t >= %eval(&maWindow. - 1))
      rename=(Variance_&varName.=var_gpu2));
   t = t + 1;
run;

proc sort data=var_comb;
   by t;
run;

proc sgplot data=var_comb;
   series x=t y=var_ma_var / legendlabel="Moving Average Variances of &varName."
      markers markerattrs=(color=red size=5 symbol=circlefilled)
      lineattrs=(color=red pattern=shortdash);
   series x=t y=var_gpu1 / legendlabel="Variance_&varName. (GPU FP32)"
      markers markerattrs=(color=blue size=5 symbol=triangle)
      lineattrs=(color=blue pattern=solid);
   series x=t y=var_gpu2 / legendlabel="Variance_&varName. (GPU FP64)"
      markers markerattrs=(color=green size=5 symbol=diamond)
      lineattrs=(color=green pattern=longdash);
   label t="Step" var_ma_var="Values of the Variances";
run;

Output 17.2.10: Comparison of Variance Predictions: GPU FP32 vs. GPU FP64

Summary

Overall, this example demonstrates that the GPU workflow in PROC DYNAMICLINEAR can provide accuracy comparable to that of the CPU workflow while significantly reducing computation time for large SGDLMs. The GPU-accelerated computations yield better KL divergence and effective sample size than the CPU computations, as well as more efficient sampling and better convergence properties. The performance improvements from GPU acceleration make it a powerful option for analyzing large-scale problems with SGDLMs.

Last updated: July 09, 2026