`vignettes/profiling.Rmd`

`profiling.Rmd`

This vignette demonstrates how to use the new profiling functionality introduced in CmdStan 2.26.0.

Profiling identifies which parts of a Stan program are taking the longest time to run and is therefore a useful guide when working on optimizing the performance of a model.

However, be aware that the statistical assumptions that go into a model are the most important factors in overall model performance. It is often not possible to make up for model problems with just brute force computation. For ideas on how to address performance of your model from a statistical perspective, see (Gelman et al. 2020).

Consider a simple logistic regression with parameters `alpha`

and `beta`

, covariates `X`

, and outcome `y`

.

```
data {
int<lower=1> k;
int<lower=0> n;
matrix[n, k] X;
int y[n];
}
parameters {
vector[k] beta;
real alpha;
}
model {
beta ~ std_normal();
alpha ~ std_normal();
y ~ bernoulli_logit(X * beta + alpha);
}
```

A simple question is how much time do the prior calculations take compared against the likelihood? To answer this we surround the prior and likelihood calculations with `profile`

statements.

```
profile("priors") {
target += std_normal_lpdf(beta);
target += std_normal_lpdf(alpha);
}
profile("likelihood") {
target += bernoulli_logit_lpmf(y | X * beta + alpha);
}
```

We can then run the model as usual and Stan will collect the profiling information for any sections with `profile`

statements.

library(cmdstanr) # Compile the model model <- cmdstan_model("profiling-files/profiling_bernoulli_logit.stan") # Generate some fake data n <- 1000 k <- 20 X <- matrix(rnorm(n * k), ncol = k) y <- 3 * X[,1] - 2 * X[,2] + 1 p <- runif(n) y <- ifelse(p < (1 / (1 + exp(-y))), 1, 0) mdata <- list(k = ncol(X), n = nrow(X), y = y, X = X) # Run one chain of the model fit <- model$sample(data = mdata, chains = 1)

The raw profiling information can then be accessed with the `$profiles()`

method, which returns one data frame per chain (profiles across multiple chains are not automatically aggregated). Details on the column names are available in the CmdStan documentation.

fit$profiles()

```
[[1]]
name thread_id total_time forward_time reverse_time chain_stack
1 likelihood 140075679106880 0.902559 0.76531100 0.137247000 51639
2 priors 140075679106880 0.004681 0.00381098 0.000870017 34426
no_chain_stack autodiff_calls no_autodiff_calls
1 34426000 17213 1
2 0 17213 1
```

The `total_time`

column is the total time spent inside a given profile statement. It is clear that the vast majority of time is spent in the likelihood function.

Stan’s specialized glm functions can be used to make models like this faster. In this case the likelihood can be replaced with

`target += bernoulli_logit_glm_lpmf(y | X, alpha, beta);`

The profiling information for the new model is collected automatically just like the last one.

model_glm <- cmdstan_model("profiling-files/profiling_bernoulli_logit_glm.stan") fit_glm <- model_glm$sample(data = mdata, chains = 1)

fit_glm$profiles()

```
[[1]]
name thread_id total_time forward_time reverse_time chain_stack
1 likelihood 140157106120512 0.48736500 0.48669500 0.000670490 17487
2 priors 140157106120512 0.00439922 0.00372599 0.000673234 34974
no_chain_stack autodiff_calls no_autodiff_calls
1 0 17487 1
2 0 17487 1
```

We can see from the `total_time`

column that this is much faster.

The other columns of the profiling output are documented in the CmdStan documentation.

The timing numbers are broken down by forward pass and reverse pass, and the `chain_stack`

and `no_chain_stack`

columns contain information about how many autodiff variables were saved in the process of performing a calculation.

These numbers are all totals – times are the total times over the whole calculation, and `chain_stack`

counts are similarly the total counts of autodiff variables used over the whole calculation. It is often convenient to have per-gradient calculations (which will be more stable across runs with different seeds). To compute these, use the `autodiff_calls`

column.

profile_chain_1 <- fit$profiles()[[1]] per_gradient_timing <- profile_chain_1["total_time"]/profile_chain_1["autodiff_calls"] print(per_gradient_timing)

```
total_time
1 5.243473e-05
2 2.719456e-07
```

After sampling (or optimization or variational inference) finishes, CmdStan stores the profiling data in CSV files in a temporary location. The paths of the profiling CSV files can be retrieved using `$profile_files()`

.

fit$profile_files()

`[1] "/tmp/RtmpzU9HVj/profiling_bernoulli_logit-profile-202101252058-1-4e28cc.csv"`

These can be saved to a more permanent location with the `$save_profile_files()`

method.

# see ?save_profile_files for info on optional arguments fit$save_profile_files(dir = "path/to/directory")

Gelman, Andrew, Aki Vehtari, Daniel Simpson, Charles C. Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, and Martin Modrák. 2020. “Bayesian Workflow.” http://arxiv.org/abs/2011.01808.