Take the 2-minute tour ×

Say I have a matrix with 3 columns: c1, c2,c3, and I want to create a new matrix in which each column is any possible product of two of the columns of this matrix.

So, if I had a matrix with d columns, I would like to create a new matrix with d+d(d-1)/2+d columns. For example, consider the matrix with 3 columns c1, c2,c3. The matrix that I would like to create should have the columns c1, c2,c3, c1xc2, c2xc3,c1xc3, c1^2, c2^2 and c3^2.

Is there any efficient way to do this?

share|improve this question
look into nchoosek –  slayton Oct 9 '12 at 0:36

3 Answers 3

up vote 4 down vote accepted

I'm embarrassed to post this - I'm sure there must be a simpler way (there is a MUCH simpler way - see my december update at the bottom of the answer), but this will do the job:

A = [1 2 3; 4 5 6];
n =  size(A, 2);

B = A(:, reshape(ones(n, 1) * (1:n), 1, n^2)) .* repmat(A, 1, n);
Soln = [A, B(:, logical(reshape(tril(toeplitz(ones(n, 1))), 1, n^2)'))];

The calculation is not efficient, since in the B step I actually calculate double the number of combinations that I need (ie I get c1.*c1, c1.*c2, c1.*c3, c2.*c1, c2.*c2, c2.*c3, c3.*c1, c3.*c2, c3.*c3), and then in the second step I pull out only the columns that I need (eg I get rid of c3.*c1 as I've already got c1.*c3 and so on).

UPDATE: Was just out driving and a much better method occurred to me. You just need to construct two index vectors of the form: I1 = [1 1 1 2 2 3] and I2 = [1 2 3 2 3 3], then (A(:, I1) .* A(:, I2)) will get you all the column products you are after. I'm away from my computer at the moment, but will come back later and work out a general way to construct the index vectors. I think it can be fairly easily accompished using the tril(toeplitz) construction. Cheers. Will update in a few hours.

UPDATE: Rody's second solution (+1) is exactly what I had in mind with my previous update so I won't bother repeating what he has done there now. Yoda's is quite neat too actually, so another +1.

DECEMBER UPDATE: Funnily enough, after working on it here, I had to revisit this problem for my own research (coding up White's test for heteroscedasticity). I'm actually favoring a new approach now, recommended (somewhat cryptically) by @slayton in the comments. Specifically, using nchoosek. My new solution looks like this:

T = 20; K = 4;
X = randi(100, T, K);
Index = nchoosek((1:K), 2);
XAll = [X, X(:, Index(:, 1)) .* X(:, Index(:, 2)), X.^2];

nchoosek yields exactly the indices we need to construct the cross-products quickly and easily!

share|improve this answer
I just completed my answer and only then saw your edit; we're clearly on the same page here :) +1 for you. I didn't know the toeplitz and hankel commands, those will definitely come in useful in the future. –  Rody Oldenhuis Oct 9 '12 at 6:25
Great Guys. Your help made my day. –  FabianG Oct 11 '12 at 5:04

The following is somewhat decent:

B = arrayfun(@(x) circshift(A, [0 -x]), 0:size(A,2)-1, 'UniformOutput', false);
B = cat(2, ones(size(A)), B{:});

C = repmat(A, 1, size(A,2)+1) .* B

This will however result in the matrix

c1 c2 c3 (c1.*c1) (c2.*c2) (c3.*c3) (c1.*c2) (c2.*c3) (c3.*c1) (c1.*c3) (c2.*c1) (c3.*c2)

of which the sequence is different than what you asked, and the products are not unique. If you only want all unique products, use this:

[sA1, sA2] = size(A);

aa = repmat(1:sA2, sA2,1);

C = [A, A(:,nonzeros(triu(aa))).*A(:,nonzeros(triu(aa.')))]  

which results in

c1 c2 c3 (c1.*c1) (c2.*c1) (c2.*c2) (c3.*c1) (c3.*c2) (c3.*c3)

which is a different sequence than you asked for, but it contains only the unique products.

Does the sequence matter for your purpose?

share|improve this answer
Your second solution +1 is exactly what I had in mind, but I wasn't confident enough that I could write it down without testing it in Matlab first :-) –  Colin T Bowers Oct 9 '12 at 10:18

Here's another alternative. First, define a function that returns all possible pairs (as per your requirements in the question) for a given number of columns:


The function should be fairly self-explanatory. Try looking at the output for a few small values of n and see for yourself. With this in place, you can proceed as follows:

s = RandStream('twister','Seed',1); %for reproducibility 
x = rand(s, 4, 3) %your matrix

%    0.4170    0.1468    0.3968
%    0.7203    0.0923    0.5388
%    0.0001    0.1863    0.4192
%    0.3023    0.3456    0.6852

o = cellfun(@(c)prod(x(:,c),2),cols(size(x,2)),'UniformOutput',0);
out = cat(2,o{:})

%    0.4170    0.1468    0.3968    0.1739    0.0215    0.1574    0.0612    0.1655    0.0582
%    0.7203    0.0923    0.5388    0.5189    0.0085    0.2903    0.0665    0.3881    0.0498
%    0.0001    0.1863    0.4192    0.0000    0.0347    0.1757    0.0000    0.0000    0.0781
%    0.3023    0.3456    0.6852    0.0914    0.1194    0.4695    0.1045    0.2072    0.2368
share|improve this answer

Your Answer


