MATLAB: Finding all possible polynomial combinations of n variables

combinationdiophantineMATLABmeshgridndgridpolynomial

I am trying to create additional polynomial features (with polynomial degree up to p where p can be anywhere between 1 to 10) to my current dataset of n variables/features.
For example, I have three variables/features: x1, x2, and x3
If I want to expand the dataset to include polynomial degree up to 2,
then I would have: x1, x1^2, x2, x1x2, x2^2, x3, x1x3, x2x3, x3^2
Once I have all possible combination of exponents the variables can take, I can easily do element-wise exponentiation.
Solution I came up with for the toy example above is the following:
% 3 variables
% 2 deg max polynomial
[var1_exp, var2_exp, var3_exp] = ndgrid(0:2);
% all exponent combinations
exp_comb = [var1_exp(:), var2_exp(:), var3_exp(:)];
% logic to remove combinations that aren't valid
exp_comb = exp_comb(sum(exp_comb,2)<=2 & sum(exp_comb,2)>0, :);
The solution indeed matches what I had above:
>> exp_comb
exp_comb =
     1     0     0
     2     0     0
     0     1     0
     1     1     0
     0     2     0
     0     0     1
     1     0     1
     0     1     1
     0     0     2
However, I am not sure how to generalize this easily to large n cases.
Any help would be greatly appreciated.
UPDATES:
Following solutions have been posted so far with its own pros and cons
Solution 1 – works but doesn't scale well to both n and p in terms of memory (n=30 and p=2 requires 1534009.5GB memory)
n = 3; % Number of variables/features

p = 2; % Maximum polynomial degree

vars = cell(n,1);
[vars{1:n}] = ndgrid(0:p);
exp_comb = reshape(cat(n+1, vars{:}), [], n);
exp_comb = exp_comb(sum(exp_comb,2)<=p & sum(exp_comb,2)>0, :);
Solution 2 – works without demanding tremendous memory, but doesn't generalize to p other than 2
n = 3; % Number of variables/features
p = 2; % Maximum polynomial degree
combs=nchoosek(1:n,p);
m=size(combs,1);
S=sparse([1:m;1:m],combs.',1,m,n);
exp_comb = [S;speye(n);2*speye(n)];

Best Answer

This should be better:
n=30;
combs=nchoosek(1:n,2);
m=size(combs,1);
S=sparse([1:m;1:m],combs.',1,m,n);
exp_comb = [S;speye(n);2*speye(n)];
Related Question