Improving the efficiency of the Gaussian and (Gaussian) copula methods

[LHBO]

In this pull request, we improve the efficiency of the Gaussian and (Gaussian) copula methods.

The efficiency has been improved by

1. Computing the conditional Gaussian distributions only once and then adding the explicand-specific parts to the precomputed quantities as outlined in Improve efficency for Gaussian method #231.
2. Sampling MC samples only once from a standard normal and then transforming them to the correct conditional normal distribution.
3. Rewriting the computationally heavy parts of the code into C++ code to speed up the computations.

In this pull request, the following C++ functions replace the equivalent R functions:

• The C++ function prepare_data_gaussian_cpp replaces the R function prepare_data_gaussian.
• The C++ function prepare_data_copula_cpp replaces the R function prepare_data_copula.
• The C++ function inv_gaussian_transform_cpp replaces the R function inv_gaussian_transform.
• The C++ function quantile_type7_cpp replaces the R function quantile(..., type = 7).

Closes #231

[LHBO]

Currently, I have just added an extra file where I have written six versions of the Gaussian method. We need to decide which one to use and replace the prepare_data.gaussian() function in approach_gaussian.R with the one we decide to use.

0. In the old version, we compute and test the conditional covariance matrix for each coalition and explicand. Generate new MC samples for each explicand and coalition.
1. Only compute and test the covariance matrix once for each coalition, but mvnfast computes the Cholensky decomposition for each explicand. So it is computed n_explain times. Generate new MC samples for each explicand and coalition.
2. Same as 1, but we use R's chol() function to compute the Cholensky decomposition once and provide that directly to rmvn().
3. We only sample once per coalition and only add the explicand dependent mean in a secondary call. I.e., we sample the n_samples MC samples from N(0, Sigma_{Sbar|S}) for each coalition, and then add the mean for the different explicands. Let rmvn() compute the Cholensky decomposition, but this is only done once now.
4. Same as above, but we use R's chol() function to compute the Cholensky decomposition once and provide that directly to rmvn().
5. Only generate the n_samples MC samples once. We generate n_samples from N(0, I) and then use Cholensky decomposition to transform to N(0, Sigma_{Sbar|S}), and then add the explicand specific means. That is, we use the same n_samples to create the MC samples for all n_explain explicands and n_combinations coalitions.
6. Same as 5, but generate n_samples * n_combinations MC samples from N(0, I), such that each coalition is based on different samples. Different explicands are still based on the same generated samples from N(0,I).

In general, we see in a small 8-dim example that the new versions are between two to four times faster than the old version. Version 3-6 are faster than 1-2. In general number 5 seems to be the fastest (in the setup I looked at).

One can probably write more efficient code, too.

[martinju]

This is great! I am most in favor v5, which I added a simplification to, using rnorm directly instead of rmvnorm. This is slightly faster. I didn't check, but I think one can then also go directly to the transpose of B so that may help even more.
Further, since we are improving code here, it would be nice to see whether the whole double-lapply think can be done with 2 basic for loops in C++. All of that code is essentially just matrix multiplication etc which is much fast in C++. It is quite much code to transform, but ChatGPT should be super-helpful here, and you can also see the other cpp-files in shapr for how to set it up. Note that it does not matter if the c++-code is not "efficient" using loops and so on, it will anyway be fast since it is being complied.

Feel free to try that out if you want. Otherwise, I can give it a go next week.

[LHBO]

This is great! I am most in favor v5, which I added a simplification to, using rnorm directly instead of rmvnorm. This is slightly faster. I didn't check, but I think one can then also go directly to the transpose of B so that may help even more. Further, since we are improving code here, it would be nice to see whether the whole double-lapply think can be done with 2 basic for loops in C++. All of that code is essentially just matrix multiplication etc which is much fast in C++. It is quite much code to transform, but ChatGPT should be super-helpful here, and you can also see the other cpp-files in shapr for how to set it up. Note that it does not matter if the c++-code is not "efficient" using loops and so on, it will anyway be fast since it is being complied.

Feel free to try that out if you want. Otherwise, I can give it a go next week.

I will look at it next week. I agree that using C++ will be beneficial

[martinju]

@Lars: All good except the manual update mentioned above, and the ties-issue mentioned on slack. I will approve the failing tests once that is fixed, and you have confirmed that the tests differences are due to sampling error only.

[LHBO]

@Lars: All good except the manual update mentioned above, and the ties-issue mentioned on slack. I will approve the failing tests once that is fixed, and you have confirmed that the tests differences are due to sampling error only.

• Fixed the manual update mentioned above.
• I fixed the ties issues by using the R rank function for that problem. We do not need C++ code there, as the code takes milliseconds to run.
• I added a code example in shapr/inst/scripts/compare_copula_in_R_and_C++.R where I compare the Shapley values produced by the R and C++ versions. The difference tends to zero when n_samples tends to infinity. For n_samples = 1000000, the max absolute difference is in the third decimal.