I want to compute all d-dimensional monomials of degree less than k and put them into a cell array - Pbase - in an ordered way (I use Matlab but this ploblem applies for other languages too). k and d are provided by the user but arbitrary.
What I have done so far,(the array called degreeindex show the location of monomials of a certain degree in Pbase):
n=nchoosek(d+k,k);%i know there are n over k possibilities for such monomials
Pbase=cell(1,n);
degreeindex=zeros(k+1,3);%showing [degree, start, end]
x=sym('x',[d,1]);
%the polynomials of degree 0
Pbase{1}=1;
degreeindex(1,1:3)=[0,1,1];%initialized for degree 0
%degree 1 monomials
degree=2;
degreeindex(degree,:)=[degree-1,degreeindex(degree-1,3)+1,degreeindex(degree-1,3)+1];
for i=1:d
Pbase{degreeindex(degree,2)+i-1}=x(i);
degreeindex(degree,3)=degreeindex(degree,3)+1;
end
%degree2 monomials
degree=3;
degreeindex(degree,:)=[degree-1,degreeindex(degree-1,3),degreeindex(degree-1,3)+1];
for i=1:d
for j=i:d
Pbase{degreeindex(degree,3)-1}=x(i).*x(j);
degreeindex(degree,3)=degreeindex(degree,3)+1;
end
end
%degree3 monomials
degree=4;
degreeindex(degree,:)=[degree-1,degreeindex(degree-1,3),degreeindex(degree-1,3)+1];
for i1=1:d
for i2=i1:d
for i3=i2:d
Pbase{degreeindex(degree,3)-1}=x(i1).*x(i2)*x(i3);
degreeindex(degree,3)=degreeindex(degree,3)+1;
end
end
end
...
The problem is that I don't find a way to do this for an arbitrary degree k. In the solution above I would have to include a new (and deeper) nested loop for every degree.
I know this seems like a faily trivial problem but I can't get my head around it. I appreciate all advice.
You are trying to do a type of partitioning, specifically you want to find all integer partitions of k into m or fewer parts (where order matters).
Using nchoosek
you can easily achive this (adapted from the Matlab File Exchange,here):
d = 3; % num dims
k = 4; % degree
x=sym('x',[d,1]);
m = nchoosek(k+d-1,d-1);
dividers = [zeros(m,1),nchoosek((1:(k+d-1))',d-1),ones(m,1)*(k+d)];
a = diff(dividers,1,2)-1;
PBase = cell(1, size(a,1));
for i = 1:size(a,1)
PBase{i} = prod(x.' .^ a(i,:));
end
Loop over whatever values of k
you need.