Optimal way to creating new features from training set (in R) My training dataset has 1500 features - all numeric. (But only about 200 data points). I want to create additional features and then use the exhaustive list for feature selection.
For creating additional features, I considered generating all 2-way interactions between the features. I ran the following code in R - but it is taking a lot of time. For 1500 features, selecting all 2-way products will result in about 1.2 million features. After about 2 hours of running, the result had only about 300k features generated. Is there a more efficient way to generate this? Could this be done faster?
Here is the relevant code:
for(i in 1:ncol(train)){
    for(j in i:ncol(train)){
        two_way_features <- cbind(two_way_features,train[,i]*train[,j])
}}

EDIT:
Apologies - I didn't add in the context. We are running an online campaign - sent to 200 companies. 6 measures were observed post the campaign. I have about 1500 features (including both about the campaign and about the companies). The aim is to learn from this. (And unfortunately, in B2B space, there's only so much we can do to get data. While I can try to get more features, can't expand to more than 200 customers/companies). 
Learning from this data, we would like to predict the 6 measures for another set of 100 companies (on whom we plan to run the campaign post thanksgiving).  
The entire dataset is numeric(both the measured target variables and the feature vectors). All of them are non-negative. No categorical variables too. 
For now, glmnet is what I am doing. After playing with glmnet for a while with the given feature vectors, the cross-validation score hasn't been improving. I am trying to investigate if combination of feature vectors provide a better CV score. Given that there's really limited data points, I haven't explored outside glmnet. (strictly LASSO). 
 A: First of all, I don't see how a reliable feature selection among > 1e6 features can take place on the basis of only 200 cases.  
Just to be sure we're talking about the same thing: modeling 200 cases with 1500 features did not lead to sufficient degrees of freedom in your model!?
Feature selection is a really difficult task, which usually leads to massive multiple comparison situations. You may get away with so few cases in extremely benign regression problems, but I don't see any chance for exhaustive search feature selection for classification (not even with proper scoring rules). But I'd be happy to learn the opposite :-)  

Now about the more programming-related part of the question:


*

*repeatedly using cbind is for sure going to be very slow, because a lot of copying will be going on.

*one way to speed up things would be to set up the matrix (or data.frame) with the correct number of columns, and then just fill the columns

*for loops are also slow in R because there's a whole lot of overhead to make sure you're not messing around with the loop variable.


Here are 2 possibilities to speed up things.
Idea 1: more efficient looping
You want to generate 1.2e6 features, but only for 200 cases. So if you have to loop, loop over the cases rather than over the features.  
Feature generation for one row can be done very efficiently by calculating the outer product of the row with itself, and then just taking the lower (or upper) triangle.
feature.gen <- function (x){
  x <- outer (x, x)
  x [lower.tri (x)]
}

outer can also help you generating appropriate column names:
colnames.gen <- function (X){
  x <- outer (colnames (X), colnames (X), paste, sep = ".")
  x [lower.tri (x)]
}

Now you have 2 possibilities:


*

*apply:
xnew <- apply (X, 1, feature.gen)


*for-loop:
xnew <- matrix (NA_real_, nrow = (ncol (X)^2 - ncol (X)) / 2, ncol = nrow (X))
for (i in 1 : nr) 
    xnew [,i] <- feature.gen (X [i,])

I tested both with a data matrix of size 200 x 1000 (memory!), and 


*

*runtimes were basically the same (15 s) but

*the 2nd approach needs considerably less memory.

*colnames.gen needed 11 s for a matrix with 1500 columns.


Note that for both approaches you need to cbind your original data matrix to the transpose of xnew to have the original and interactions in the new data matrix.
Idea 2: Use a kernel
If you don't care whether the square terms are included or not, maybe your modeling is available in a kernel version. In that case, you could use a polynomial kernel of degree 2.
A: Apart from what cbeleites told you, remember that $log( x \cdot y ) = \log( x ) + \log( y )$, and that addition is much more computationally effective than multiplication. So instead of multiplying the columns, logarithmize them first, add them, and exp them.
But hey, if you have addition, you might as well leave that to your ML algorithm: just do 
new_matrix <- cbind( old_matrix, log( old_matrix )

In my field (biology) what sometimes really helps is to create new features based on relative proportions (so $\frac{a}{b}$ rather than $a\cdot b$); often not the feature itself, but its proportion relative to another feature is informative, as in "relative level of amino acid X to general level of amino acids".
However, this will also be taken care of when you add the logarithmized features, since $\frac{a}{b}=\log(a)-\log(b)$
