glmnet package: "mgaussian" vs "gaussian" for $\alpha = 0$

In multiresponse Gaussian family the objective function when $$\alpha = 0$$:

\begin{align} \frac{1}{2n}||Y-XB||_F^2 + \frac{\lambda}{2}||B||_F^2. \end{align}

This can also mathematically solved as \begin{align*} \frac{1}{2nk}||vec(Y)- (I \otimes X)vec(B)||_2^2 + \frac{\lambda}{2}||vec(B)||_2^2. \end{align*} yielding the same solution. However, when I tried to do it using glmnet package the answers I got are different. Can anyone explain why it is?

require(glmnet)
set.seed(7)
n <- 100
p <- 5
k <- 2

X <- matrix(rnorm(n * p), ncol=p)
beta <- matrix(rnorm(p * k , 0, 1), ncol = k)
Y <- X %*% beta + matrix(rnorm(n * k), ncol=k)

Xnw <- diag(k) %x% X
Ynw <- c(Y)

sd_ynw <- sqrt(var(Ynw)*(200-1)/200)
fit_glmnet1 <- glmnet(X, Y/sd_ynw, alpha = 0, standardize = F,
standardize.response = FALSE,
intercept = FALSE,
thresh = 1e-20, lambda = c(1, 2, 3), family = "mgaussian")

fit_glmnet2 <- glmnet(Xnw, Ynw/sd_ynw, alpha=0, standardize = F, intercept = FALSE,
thresh = 1e-20, lambda = c(1, 2, 3))
do.call('cbind', coef(fit_glmnet1, s = 1))[-1, ]
matrix(coef(fit_glmnet2, s = 1)[-1], 5, 2)

> do.call('cbind', coef(fit_glmnet1, s = 1))[-1, ]
5 x 2 sparse Matrix of class "dgCMatrix"
1          1
V1 -0.01826165 -0.3290432
V2 -0.09559630 -0.2715274
V3 -0.13507974 -0.1161845
V4  0.18951137 -0.3014415
V5 -0.09837942 -0.2702261
> matrix(coef(fit_glmnet2, s = 1)[-1], 5, 2)
[,1]        [,2]
[1,] -0.005284682 -0.21435828
[2,] -0.063831599 -0.19521015
[3,] -0.091713598 -0.08387311
[4,]  0.124559675 -0.21271566
[5,] -0.068698728 -0.17432439


I have made a mistake in the comparison. We can't directly compare the same $$\lambda$$.

set.seed(7)
n <- 100
p <- 5
k <- 2

X <- matrix(rnorm(n * p), ncol=p)
beta <- matrix(rnorm(p * k , 0, 1), ncol = k)
Y <- X %*% beta + matrix(rnorm(n * k), ncol=k)

Xnw <- diag(k) %x% X
Ynw <- c(Y)

sd_ynw <- sqrt(var(Ynw)*(200-1)/200)
fit_glmnet1 <- glmnet(X, Y/sd_ynw, alpha = 0, standardize = F, standardize.response = FALSE,
intercept = FALSE,
thresh = 1e-20, lambda = c(1, 2, 3, 4), family = "mgaussian")
solve(t(X) %*% X + 4*n*diag(p), t(X)%*%(Y/sd_ynw))
do.call('cbind', coef(fit_glmnet1, s = 4))[-1, ]

fit_glmnet2 <- glmnet(Xnw, Ynw/sd_ynw, alpha=0, standardize = F, intercept = FALSE,
thresh = 1e-20, lambda = c(1, 2, 3), family = "gaussian")
matrix(solve(t(Xnw) %*%Xnw + 2*200 * diag(p*k), t(Xnw)%*%(Ynw/sd_ynw)), 5, 2)
matrix(coef(fit_glmnet2, s = 2)[-1], 5, 2)


Then the solutions are identical.

• I don't understand this answer. Why did you insert the solve lines? They don't seem to influence anything. If you can't compare the same lambdas, then how to you change the lambdas? Both glmnet commands still use 1,2,3. Commented Mar 1, 2019 at 13:59
• @amoeba I just included manual calculations of $B$. Both glmnet commands don't use the same set of lambda. There is an additional 4 in fit_glmnet1. The sample sizes are different in the two objective functions: one has $n$ and other has $n*k$. This needs to be taken into account to compare lambda. Commented Mar 1, 2019 at 21:38