Fastest way to create matrix of interactions many times

I have this code, where I have 8 columns in my matrix of regressors X (36 unique interactions):

interactions=zeros(size(X,1),36);
iii=0;
for kk=1:8
    for jj=kk:8
        iii=iii+1;
        interactions(:,iii)=X(:,kk).*X(:,jj);
    end
end
X_all=[ones(size(X,1),1),X,interactions];

I have to run this 500,000 times for different Xs.

The code above is much faster when doing it many times than running the matlab function below many times (due to overhead I think):

 X_all=x2fx(X,'quadratic');

Is there some way to make my code faster? Perhaps avoid the loop somehow?

The profiler shows that most of the time is spent on the multiplication part, followed by the last line.