- 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.

*(This is what IMHO keeps the question appropriate for cross validated)*