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.