2
$\begingroup$

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
$\endgroup$

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – amoeba
    Commented Mar 1, 2019 at 13:59
  • $\begingroup$ @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. $\endgroup$
    – shani
    Commented Mar 1, 2019 at 21:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.