This is an idea to take advantage of matrix multiplication for speed. The order of the columns is a bit different than what you wanted, but hopefully this is still useful
import pandas as pd
import itertools
import scipy as scp
# Create a large matrix, A, to multiply df1 against
# Df1*A = Y
#
# Y needs to have (df1.rows) rows by (df1.cols*df2.cols*df2.rows) columns
#
# So A will have the shape of (df1.cols) rows by (df1.cols*df2.cols*df2.rows) columns
#
# A is just diagonal matrices hstacked together, where the diagonal element is the value in Df2
#
# Using scp.sparse matrices for A to avoid too much memory use
np.random.seed(1)
df1 = pd.DataFrame({
'col1':[1,2,3,4,5,6,7,8,9],
'col2':[11,12,13,14,15,16,17,18,19],
'col3':[21,22,23,24,25,26,27,28,29],
})
data = {}
n_df2_rows = 3 #control how many rows for Df2
n_df2_cols = 3 #control how many cols for Df2
for i in range(n_df2_rows):
data['C'+str(i)] = np.random.randint(1,10,n_df2_cols)
df2 = pd.DataFrame(
data,
index=['R{}'.format(i) for i in range(n_df2_cols)],
)
r1,c1 = df1.shape
r2,c2 = df2.shape
A = scp.sparse.hstack([scp.sparse.eye(c1,dtype=np.int8)*v for v in df2.values.flatten('F')])
Y = pd.DataFrame(df1.values @ A) #@ is sparse matrix multiplication
Y.columns = ['{}_{}_{}'.format(c,v,i) for i,v,c in itertools.product(df2.columns,df2.index,df1.columns)]
Y
Df2 using np.random.seed(1)

Output

Even when Df2 is 100 by 100 (and the output has 30,000 columns) this code runs in less than a minute for me