Generating all possible pairs of polynomial interaction combinations?

I have a matrix X where column one represents a time series, lets call it “ a ”, and column two represents another time series, that we can call “ b ”.

What I want to produce is all possible combinations of each time series multiplied with polynomial versions of the other series (“the other series” will be several other series in the future), after specifying a maximum polynomial order. I would also like to receive the polynomial versions of each variable separately (multiplied with nothing) if possible.

So say I want the polynomial order to be no higher than 3 and I only have the two variables a and b, what I would like to produce is the following (provided I didn’t miss any combination when doing this series manually):

a^1
a^2
a^3
b^1
b^2
b^3
a^1*b^1
a^2*b^1
a^3*b^1
a^1*b^2
a^1*b^3

Preferably I would like to be able to do this even when the polynomial orders aren’t integers, and when I have several more variables than a and b. So if I were to add a variable c I would also like all possible combinations of that with a and b, but I would still like the combinations expressed as pairs, or rather combinations of max length 2, so a^1*b^1*c^1 for example is too long and therefore not a valid combination.

Most combination generating approaches I can think of will probably also generate combinations such as: a^1*a^2, but since that’s the same as just a^3 I don’t want it, I could of course generate it and then remove it by a function like “unique” etc. but if possible I would prefer to not generate those combinations to begin with in the interest of saving memory.

Does anyone have any ideas for solving this?

Thanks in advance!

Answer 1

[John D'Errico]

Hint: All you need to do is generate the set of exponents. For example, here column 1 of uv might represent the exponent of a, column 2 the exponent of b.

[u,v] = meshgrid(0:3);
uv = [u(:),v(:)]
uv =
   0     0
   0     1
   0     2
   0     3
   1     0
   1     1
   1     2
   1     3
   2     0
   2     1
   2     2
   2     3
   3     0
   3     1
   3     2
   3     3

If you don't wish a constant term, then delete the first row. You state that the polynomial order should be no greater than 3, but you included a term in there for a*b^3, which would have a total order of 4. So I'm not sure if you were careful in your example.

But if you do not wish to have terms in the model with total order greater than 3, that is trivially done:

uv(sum(uv,2)>3,:) = []
uv =
   0     0
   0     1
   0     2
   0     3
   1     0
   1     1
   1     2
   2     0
   2     1
   3     0

Dropping the constant term, we get:

uv(sum(uv,2) == 0,:) = []
uv =
   0     1
   0     2
   0     3
   1     0
   1     1
   1     2
   2     0
   2     1
   3     0

I'm not sure what you mean by non-integer orders of the polynomials, but the above example does not require the exponents to be integers. As well, I don't see any issue with the fact that you must delete some terms. The memory involved is trivial. Finally, actual computation of the polynomial is as trivial as:

syms a b
sum(a.^uv(:,1).*b.^uv(:,2))
ans =
a^3 + a^2*b + a^2 + a*b^2 + a*b + a + b^3 + b^2 + b

Higher orders?

help ndgrid
help meshgrid

[Peta]

Thank you for the advice, that did absolutely give me a push in the right direction. But I still find one thing problematic with this approach: when I run my script the number of time series ( a,b,c…) can vary from time to time. As I understand it you have to specify manually how many outputs you want from the ndgrid command by filling in additional variables.

So if I run the script one time I might need:

[X1,X2] = ndgrid(0:3)

But if I run it a second time I could instead need:

[X1,X2,X3] = ndgrid(0:3)

And so on. It could be that I need up to 50 or so variables, it depends on the width of my x matrix. Is there any way to write that so it can dynamically adjust itself and create all combinations given the dimensions of x?