Coding multiple linear regression in R when sample size is smaller than variable size

Question

Here is my multiple linear regression code where $Xtrain$ and $Xtest$ are matrices containing respectively the train and test samples in rows and explanatory variables in columns, and $ytrain$ and $ytest$ are the output column vectors containing the response variable values respectively for the train and test samples.

ginv() is the generalized inverse function since normal inverse (function solve()) would not work here since the number of variables (1000) is larger than the number of samples (30).

$pred$ is the prediction of the response for the test samples using the fitted model.

library(MASS)
w = ginv(t(Xtrain) %*% Xtrain) %*% t(Xtrain) %*% ytrain
pred = Xtest %*% w

Interestingly, the prediction has 100% accuracy when I compare with the true $ytest$. The responses are correct in all decimal points!

How come this can happen? I know this is impossible with high-dimensional data, and also, when I use glmnet() (with no regularization) instead of my own code, it gives 22% accuracy. So, do you have any idea about what is wrong with above code? Is there something wrong with ginv() (this is the first time I use it)? Does glmnet() function use generalized inverse as well when the sample size is smaller than the variable size?

Answer 1

What do you expect with 1000 variables and 30 observations? There is nothing wrong with your code, with fewer obs than variables you can always find a perfect fit with least squares (the solution is not unique, but using ginv you try to get the solution with smallest norm.)

glmnet does not do this, it does not compute matrix inverses (generalized or not), it uses regularization, you could look through the tag regularization. So it does not find a least-squares solution, in stead it uses cross-validation as a help to choose the best regularization parameter (you also get a solution with regularization parameter close to zero, but this solution does not need to equal that you find with ginv.)

Your best bet is to go with glmnet and use regularization, or better, get more data.