• 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)