How to Perform Cross-Validation for LASSO in GAMLSS to Find Optimal $\lambda$?

I am working with a Generalized Additive Model for Location, Scale, and Shape (GAMLSS) and trying to determine the optimal $\lambda$ values for LASSO-penalized regression using cross-validation. However, I am struggling to understand how to properly set up the cross-validation procedure in this context.

For standard Generalized Linear Models (GLMs), this can be done using cv.glmnet(), which follows these steps:

1. The dataset is split into $k$ folds.
2. For each fold, a model is trained on the remaining data using a grid of 100 $\lambda$ values.
3. These $k$ models per $\lambda$ are then evaluated on the corresponding $k$ test sets.
4. The average performance across the $k$ folds is computed for each $\lambda$.
5. The $\lambda$ that minimizes the average test error is selected.

I could not figure out how to implement LASSO regression using the gamlss2 package, but I found gamlss.lasso, which seems to serve the same purpose. However, my main issue is not the LASSO implementation itself but rather how to correctly set up the cross-validation procedure to select the best $\lambda$.

For a GLM, we only tune one $\lambda$, but for GAMLSS, we need to tune two: one for the location parameter $\mu$ and one for the scale/dispersion parameter $\theta$. This suggests that we should use a two-dimensional grid of $\lambda$ values, meaning that if we use 100 values for $\lambda_\mu$ and 100 for $\lambda_\theta$, we would have $100 \times 100 = 10,000$ different combinations. Given $k = 10$ folds, this would result in training 100,000 models in total, which seems computationally expensive.

I attempted to implement LASSO for GAMLSS without cross-validation as follows:

set.seed(123) n <- 500 d <- 50 X <- matrix(rnorm(n*d), n, d) BETA <- cbind("mu" = rbinom(d, 1, .1), "sigma" = rbinom(d, 1, .1) * .3) ysd <- exp(1 + tcrossprod(BETA[,2], X)) data <- cbind(y = as.numeric(rnorm(n, sd = ysd)) + t(tcrossprod(BETA[,1], X)), as.data.frame(X)) # Estimating the model with gnet default setting mod <- gamlss(y ~ gnet(x.vars = names(data)[-1]), sigma.fo = ~ gnet(x.vars = names(data)[-1]), data = data, family = NO, i.control = glim.control(cyc = 1, bf.cyc = 1)) mu_betas <- as.matrix(mu[[1]]$beta) mu_lambdas <- mu[[1]]$lambda sigma_betas <- as.matrix(sigma[[1]]$beta) sigma_lambdas <- sigma[[1]]$lambda 

From the estimated coefficients for each $\lambda$ value, one could then evaluate performance on test data and construct a cross-validation procedure. However, I am unclear about the following points:

1. Since we model two parameters in GAMLSS, do we need a two-dimensional grid of $\lambda$ values? That is, should we test all $100 \times 100$ combinations $(\lambda_\mu, \lambda_\theta)$ independently? This would allow for cases where, for example, $\lambda_\mu$ is large while $\lambda_\theta$ is small if $\theta$ depends on many variables while $\mu$ depends on few.
2. In the above implementation, it appears that only 100 combinations are tested, following the sequence $(\lambda_{\mu}[1], \lambda_{\theta}[1]), ..., (\lambda_{\mu}[100], \lambda_{\theta}[100])$ rather than evaluating all possible pairwise combinations from two independent grids. Is this the correct approach, or should all $10,000$ combinations be tested explicitly?
3. I read that it is possible to combine LASSO and cross-validation using gamlss2, but I couldn’t find any clear documentation on how to do this. If anyone knows how to use gamlss2 for cross-validation with LASSO, I would appreciate any guidance!