Bounty Ended with 50 reputation awarded by amoeba
5 Substantial change reflecting comment.
source | link

Why is the minimum-norm OLS solution so (comparably) good in this case? I think it is related to the following behavior that I found very counter-intuitive, but on reflection makes a lot of sense.

I don't know why thisThis happens for a trivial reason. There could$y$ can be some geometric argument showing that very high dimensional OLS solution surfacesexpressed exactly as a linear combination of columns of $X$. $\hat{\beta}$ is the minimum-norm vector of coefficients. As more columns are likelyadded the norm of $\hat{\beta}$ must decrease or remain constant, because a possible linear combination is to have some point closekeep the previous coefficients the same and set the new coefficients to 0$0$.

Why is the minimum-norm OLS solution so (comparably) good in this case? I think it is related to the following behavior that I found very counter-intuitive.

I don't know why this happens. There could be some geometric argument showing that very high dimensional OLS solution surfaces are likely to have some point close to 0.

Why is the minimum-norm OLS solution so (comparably) good in this case? I think it is related to the following behavior that I found very counter-intuitive, but on reflection makes a lot of sense.

This happens for a trivial reason. $y$ can be expressed exactly as a linear combination of columns of $X$. $\hat{\beta}$ is the minimum-norm vector of coefficients. As more columns are added the norm of $\hat{\beta}$ must decrease or remain constant, because a possible linear combination is to keep the previous coefficients the same and set the new coefficients to $0$.

4 added 34 characters in body
source | link
library(glmnet)
set.seed(1846)
noise <- 10
N <- 80
num.vars <- 100
target <- runif(N,-1,1)
training.data <- matrix(nrow = N, ncol = num.vars)
for(i in 1:num.vars){
  training.data[,i] <- target + rnorm(N,0,noise)
}
plot(cv.glmnet(training.data, target, alpha = 0,
               lambda = exp(seq(-10, 10, by = 0.1))))
library(glmnet)
set.seed(1846)
noise <- 10
N <- 80
num.vars <- 100
target <- runif(N,-1,1)
training.data <- matrix(nrow = N, ncol = num.vars)
for(i in 1:num.vars){
  training.data[,i] <- target + rnorm(N,0,noise)
}
plot(cv.glmnet(training.data, target, alpha = 0,
               lambda = exp(seq(-10, 10, by = 0.1))))
max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)

  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
}


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")

abline(v =80)
max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)

  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
}


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")

abline(v = 80)
library(glmnet)
set.seed(1846)
noise <- 10
N <- 80
num.vars <- 100
target <- runif(N,-1,1)
training.data <- matrix(nrow = N, ncol = num.vars)
for(i in 1:num.vars){
  training.data[,i] <- target + rnorm(N,0,noise)
}
plot(cv.glmnet(training.data, target, alpha = 0,
               lambda = exp(seq(-10, 10, by = 0.1))))
max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)

  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
}


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")

abline(v =80)
library(glmnet)
set.seed(1846)
noise <- 10
N <- 80
num.vars <- 100
target <- runif(N,-1,1)
training.data <- matrix(nrow = N, ncol = num.vars)
for(i in 1:num.vars){
  training.data[,i] <- target + rnorm(N,0,noise)
}
plot(cv.glmnet(training.data, target, alpha = 0,
               lambda = exp(seq(-10, 10, by = 0.1))))
max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)

  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
}


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")

abline(v = 80)
3 Replaced glmnet edit with explicit solution
source | link
max.beta.random <- function(num.vars, s){
  # s between 1 and 201, with 1 corresponding to log(lambda) = 10, 201 log(lambda) = -10
  num.vars <- round(num.vars)
  s <- round(s)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars) 

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
 
  model.cvudv <- cv.glmnetsvd(training.data, target, alpha = 0,)
                    
  U <- udv$u
  S <- diag(udv$d)
 lambda =V exp(seq(<-10, 10, by = 0.1)))udv$v

  beta.valhat <- V %*% solve(model.cv$glmnet.fit$beta[,s]S) %*% t(U) %*% target

  max(abs(beta.valhat))
} 


curve(Vectorize(max.beta.random)(x, 201), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients") 

abline(v =80)

Plot of max magnitude of coefficients as number of predictors increasesPlot of max magnitude of coefficients as number of predictors increases

I don't know ifwhy this is a computational artifact, either caused by glmnet or because my "small lambda" should decrease as p increaseshappens. But if not, then thereThere could be some geometric argument showing that very high dimensional OLS solution surfaces are likely to have some point close to 0.

max.beta.random <- function(num.vars, s){
  # s between 1 and 201, with 1 corresponding to log(lambda) = 10, 201 log(lambda) = -10
  num.vars <- round(num.vars)
  s <- round(s)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars)
  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
 
  model.cv <- cv.glmnet(training.data, target, alpha = 0,
                        lambda = exp(seq(-10, 10, by = 0.1)))

  beta.val <- (model.cv$glmnet.fit$beta[,s])

  max(abs(beta.val))
}

curve(Vectorize(max.beta.random)(x, 201), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients")
abline(v =80)

Plot of max magnitude of coefficients as number of predictors increases

I don't know if this is a computational artifact, either caused by glmnet or because my "small lambda" should decrease as p increases. But if not, then there could be some geometric argument showing that very high dimensional OLS solution surfaces are likely to have some point close to 0.

max.beta.random <- function(num.vars){
  num.vars <- round(num.vars)
  set.seed(1846)
  noise <- 10
  N <- 80
  target <- runif(N,-1,1)
  training.data <- matrix(nrow = N, ncol = num.vars) 

  for(i in 1:num.vars){
    training.data[,i] <- rnorm(N,0,noise)
  }
  udv <- svd(training.data)
 
  U <- udv$u
  S <- diag(udv$d)
  V <- udv$v

  beta.hat <- V %*% solve(S) %*% t(U) %*% target

  max(abs(beta.hat))
} 


curve(Vectorize(max.beta.random)(x), from = 10, to = 1000, n = 50,
      xlab = "Number of Predictors", y = "Max Magnitude of Coefficients") 

abline(v =80)

Plot of max magnitude of coefficients as number of predictors increases

I don't know why this happens. There could be some geometric argument showing that very high dimensional OLS solution surfaces are likely to have some point close to 0.

2 Expanded answer and struck through incorrect reasoning
source | link
1
source | link