See how correlation influence.

#library(MASS)
library(glmnet)
## Loading required package: Matrix
## Loading required package: foreach
## Warning: package 'foreach' was built under R version 3.4.3
## Loaded glmnet 2.0-13

n: the number of samples (“rows” of data) required. mu: a vector giving the means of the variables. Sigma: positive-definite symmetric matrix specifying the covariance matrix of the variables.

The problem: n x p -> 3 x 4. N(0, 1)

Define the global variable X, y1, y2, y3.

set.seed(65587)
# n = 2500, p = 4 
X <- matrix(rnorm(10000), 2500, 4)

# Have 3 different responses
y1 <- 0.9 * X[, 1] + 0.00 * X[, 2] + 1.2 * X[, 3] + 0.3 * X[, 4]
y2 <- 1.1 * X[, 1] + 0.02 * X[, 2] + 0.8 * X[, 3] + 0.4 * X[, 4]
y3 <- 0.9 * X[, 1] + 0.01 * X[, 2] + 1.3 * X[, 3] + 0.5 * X[, 4]

w = c(0.9, 0.00, 1.2, 0.3, 1.1, 0.02, 0.8, 0.4, 0.9, 0.01, 1.3, 0.5)
w = matrix(w, 4, 3)

Function for multi-task LASSO experiment

mt_LASSO <- function(Cov){
    set.seed(12345)
    # epsilon is the gaussian error. dimension: n * 3. 
    epsilon <- MASS::mvrnorm(n=2500, mu = rep(0, 3), Sigma = Cov)
    y <- cbind(y1, y2, y3)  
    y <- y + epsilon
    
    mt_fit = glmnet(X, y, family = "mgaussian")
    plot(mt_fit, type.coef = "2norm", label = TRUE)
    
    ## Use Cross Validation to select $\lambda$, and compare.
    cv.mt_fit = cv.glmnet(X, y, family = "mgaussian", intercept = FALSE)
    cv.fit1 = cv.glmnet(X, y[, 1], intercept = FALSE)
    cv.fit2 = cv.glmnet(X, y[, 2], intercept = FALSE)
    cv.fit3 = cv.glmnet(X, y[, 3], intercept = FALSE)
    #plot(cv.mt_fit)
    #plot(cv.fit1)
    #plot(cv.fit2)
    #plot(cv.fit3)
    
    # multi-task
    mt_beta = cbind(coef(cv.mt_fit)$y1, coef(cv.mt_fit)$y2, coef(cv.mt_fit)$y3)
    mt_beta = mt_beta[-1,]
    print(mt_beta)
    print( sum((mt_beta - w)^2) )
    
    # lasso for each response
    lasso_beta = cbind(coef(cv.fit1), coef(cv.fit2), coef(cv.fit3))
    lasso_beta = lasso_beta[-1,]
    print(lasso_beta)
    print( sum((lasso_beta - w)^2) )
}
  1. No correlation.
Cov1 = diag(3)
Cov1
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
mt_LASSO(Cov1)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8462746 0.9825700 0.8328112
## V2 .         .         .        
## V3 1.1802004 0.7286246 1.2499327
## V4 0.2874215 0.3621474 0.4446249
## [1] 0.03434121
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8071545 0.9576152 0.7964037
## V2 .         .         .        
## V3 1.1483408 0.6737016 1.2266681
## V4 0.2360065 0.3269934 0.4236760
## [1] 0.07937392
  1. Positive correlation 2.1 Identical
Cov2 = matrix(rep(1, 9), 3, 3) 
Cov2
##      [,1] [,2] [,3]
## [1,]    1    1    1
## [2,]    1    1    1
## [3,]    1    1    1
mt_LASSO(Cov2)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8160383 0.9943091 0.8163733
## V2 .         .         .        
## V3 1.0588833 0.6976364 1.1475391
## V4 0.2011229 0.2741342 0.3463963
## [1] 0.1285631
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8105887 1.0368798 0.8124457
## V2 .         .         .        
## V3 1.0689516 0.6949369 1.1700414
## V4 0.1741034 0.3007439 0.3756211
## [1] 0.1064173

2.2

Cov3 = matrix(rep(0.5, 9), 3, 3) + diag(0.5,3) 
Cov3
##      [,1] [,2] [,3]
## [1,]  1.0  0.5  0.5
## [2,]  0.5  1.0  0.5
## [3,]  0.5  0.5  1.0
mt_LASSO(Cov3)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8392249 1.0425384 0.8132784
## V2 .         .         .        
## V3 1.1182886 0.7195983 1.1565727
## V4 0.2403472 0.2846033 0.3831914
## [1] 0.07924777
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8256398 1.0572803 0.7947820
## V2 .         .         .        
## V3 1.1182657 0.6970610 1.1576019
## V4 0.2148866 0.2801171 0.3919066
## [1] 0.08977977

2.3

set.seed(12345)
Cov4 = diag(3) 
Cov4[lower.tri(Cov4, diag = FALSE)] = runif(3, min=0, max=1)
Cov4[upper.tri(Cov4, diag = FALSE)] = Cov4[lower.tri(Cov4, diag = FALSE)]
Cov4
##           [,1]      [,2]      [,3]
## [1,] 1.0000000 0.7209039 0.8757732
## [2,] 0.7209039 1.0000000 0.7609823
## [3,] 0.8757732 0.7609823 1.0000000
mt_LASSO(Cov4)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.7895786 1.0007739 0.7924897
## V2 .         .         .        
## V3 1.1181112 0.7614904 1.1882666
## V4 0.2546002 0.3116459 0.3918438
## [1] 0.07633566
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.7825212 1.0321000 0.7889642
## V2 .         .         .        
## V3 1.1278918 0.7494973 1.2092113
## V4 0.2354336 0.3232996 0.4183386
## [1] 0.06395371

2.4

set.seed(7749)
Cov5 = diag(3) 
Cov5[lower.tri(Cov5, diag = FALSE)] = runif(3, min=0, max=1)
Cov5[upper.tri(Cov5, diag = FALSE)] = Cov5[lower.tri(Cov5, diag = FALSE)]
print(Cov5)
##            [,1]      [,2]       [,3]
## [1,] 1.00000000 0.3981169 0.08865501
## [2,] 0.39811693 1.0000000 0.41728454
## [3,] 0.08865501 0.4172845 1.00000000
mt_LASSO(Cov5)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##           1         1         1
## V1 0.858229 1.0050388 0.7961722
## V2 .        .         .        
## V3 1.166639 0.7419670 1.2259252
## V4 0.252717 0.3279717 0.4454390
## [1] 0.04241121
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8318526 1.0160257 0.7600844
## V2 .         .         .        
## V3 1.1489868 0.7185787 1.2088129
## V4 0.2072310 0.3234016 0.4367745
## [1] 0.06778983
  1. Nested 3.1
Cov2n = Cov2
Cov2n[upper.tri(Cov2n, diag = FALSE)] = Cov2n[upper.tri(Cov2n, diag = FALSE)] * c(1/3, -1/3, 1/3)
Cov2n[lower.tri(Cov2n, diag = FALSE)] = Cov2n[upper.tri(Cov2n, diag = FALSE)]
print(Cov2n)
##            [,1]      [,2]       [,3]
## [1,]  1.0000000 0.3333333 -0.3333333
## [2,]  0.3333333 1.0000000  0.3333333
## [3,] -0.3333333 0.3333333  1.0000000
mt_LASSO(Cov2n)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##           1         1         1
## V1 0.794670 0.9915535 0.8525891
## V2 .        .         .        
## V3 1.117700 0.7416349 1.2660330
## V4 0.280915 0.3490821 0.4343120
## [1] 0.04420811
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.7701739 0.9955007 0.8299128
## V2 .         .         .        
## V3 1.1026309 0.7147418 1.2586116
## V4 0.2476114 0.3414947 0.4276283
## [1] 0.06305495

3.2

Cov3n = Cov3
Cov3n[upper.tri(Cov3n, diag = FALSE)] = Cov3n[upper.tri(Cov3n, diag = FALSE)] * c(1, -1, 1)
Cov3n[lower.tri(Cov3n, diag = FALSE)] = Cov3n[upper.tri(Cov3n, diag = FALSE)]
print(Cov3n)
##      [,1] [,2] [,3]
## [1,]  1.0  0.5 -0.5
## [2,]  0.5  1.0  0.5
## [3,] -0.5  0.5  1.0
mt_LASSO(Cov3n)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.7853107 0.9841137 0.8463054
## V2 .         .         .        
## V3 1.0933937 0.7529864 1.2440799
## V4 0.2665102 0.3542374 0.4161866
## [1] 0.05690908
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.7672388 1.0012575 0.8220369
## V2 .         .         .        
## V3 1.0855855 0.7387775 1.2359840
## V4 0.2379030 0.3623392 0.4077790
## [1] 0.06866989

3.3

Cov5n = Cov5
Cov5n[upper.tri(Cov5n, diag = FALSE)] = Cov5n[upper.tri(Cov5n, diag = FALSE)] * c(1, -1, 1)
Cov5n[lower.tri(Cov5n, diag = FALSE)] = Cov5n[upper.tri(Cov5n, diag = FALSE)]
print(Cov5n)
##             [,1]      [,2]        [,3]
## [1,]  1.00000000 0.3981169 -0.08865501
## [2,]  0.39811693 1.0000000  0.41728454
## [3,] -0.08865501 0.4172845  1.00000000
mt_LASSO(Cov5n)

## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8764248 1.0272527 0.8089080
## V2 .         .         .        
## V3 1.1388556 0.7040244 1.1929661
## V4 0.2278543 0.2965519 0.4199692
## [1] 0.06136338
## 4 x 3 sparse Matrix of class "dgCMatrix"
##            1         1         1
## V1 0.8734025 1.0275528 0.7636907
## V2 .         .         .        
## V3 1.1459483 0.6689585 1.1661989
## V4 0.2059279 0.2795383 0.4030489
## [1] 0.09579256