How to find all multi-indices satisfying given contraints?

Questions
, ,
1

Hey there,

First, in case someone isn’t familiar with multi-indices, they are essentially lists of non-negative integers, e.g. \alpha = (\alpha_1,\ldots,\alpha_n). Further, denote \left| \alpha \right|_1 = \sum_{i=1}^n \left|\alpha_i\right|.

What I would like is some efficient way of finding all multi-indices for a given n satisfying

a \leq \left| \alpha \right|_1 \leq b,

with the additional constraint that \alpha_i > 0 for each 1 \leq i \leq n. The integers a and b are of course given.

More specifically, I need to sum over all such multi-indices (passing each as an argument to a function), but that’s simple enough once I have a list of the multi-indices.

-Jonatan

2

Have you tried Subsets (if you want the indices in order) or Tuples (for the more general case)? You may then have to add a shift.

The most common case is to pick k indices from n in ascending order. Say 3 indices from 5:

Subsets[Range[5], {3}]

Giving:

{{1, 2, 3}, {1, 2, 4}, {1, 2, 5}, {1, 3, 4}, {1, 3, 5}, {1, 4, 5}, {2, 3, 4}, {2, 3, 5}, {2, 4, 5}, {3, 4, 5}}
3

You could just have a1 go from 1 to b-n+1, a2 from 1 to (b-n+1)-a1, a3 from 1 to (b-n+1)-(a1+a2), etc.

One can use Reduce to find all actual multiindex sets, but that 's overkill. Likewise, could look at exponents in the expansion of (1+a1+a2+...+an)^(b-a)*(a1+...+an)^a. Again, overkill.

4

Just table it out

With[{b= ...., n = ...},
allowedInd = Table[{a1,a2,a3,...}, {a1,1,b-n+1},{a2, 1, b-n+2-a1},{a3, 1, b-n+3-a1-a2},...]
]

and use it afterwards. Be careful to formulate correct Table-iterators.

Another simple way is to include If-statements into the generic summand term, giving 0 if the multi-index is not allowed to contribute.

5

Thank you, Daniel and Udo, that seems to work. I do have a few issues with this solution though, I hope there’s some easy fix for these.

Firstly, the constraint a \leq \left|\alpha\right|_1 is missing from these solutions. As you say, this is fixable with a simple If-statement, and that’s fine, however I can’t see a way of doing this without leaving nulls/zeros/empty sets in the output. Secondly, I think I failed to mention this earlier but ideally I’d like the output to be as a list of lists, such as the output from David’s suggested method.

I’m sure both of these problems are easily fixable, I’m just new to Mathematica so I don’t yet know how to do it easily. My first thought was Flatten, but then I would need to know for sure how many levels the output list has; also, that still leaves the task of weeding out the nulls/whatnots.

Finally, there’s the issue that the suggested method only works if n is fixed – I would prefer to implement a function which simply takes a, b and n as input. Again, this is mostly a case of my question being unclear, sorry for that.

Thank you for the help so far. =)

-Jonatan

P.S. As a side note, in case someone else ends up needing this for something: the appropriate upper limits for the iterators are b-n+1, b-n+2-a1, b-n+3-a1-a2, and so forth. Otherwise, if we have b = 5 and n = 3, say, then the last limit becomes at most 5 - 3 + 1 - 1 - 1 = 1, when it should be 3 in the case where \alpha_1 = \alpha_2 = 1.

6

Lower conditions as well as the whole job can be done by brute force up to a certain n

In[41]:= With[{b = 7, n = 5, a = 1},
 If[a < 0 || n < 0 || b < 1 || b - n < a, Print["You will get the empty set!"]];
 If[a < n, Print["Your conditions run empty on lower boundary!"]];
 Union[Flatten[
    DeleteMissing[
     Table[If[a <= Plus @@ {a1, a2, a3, a4, a5} <= b, 
       X[a1, a2, a3, a4, a5], Missing[]],
      {a1, 1, b - 1}, {a2, 1, b - 1}, {a3, 1, b - 1}, {a4, 1, 
       b - 1}, {a5, 1, b - 1}], n]]] /. X -> List
 ]

During evaluation of In[41]:= Your conditions run empty on lower boundary!

Out[41]= {{1, 1, 1, 1, 1}, {1, 1, 1, 1, 2}, {1, 1, 1, 1, 3}, {1, 1, 1, 2, 1}, {1, 1, 1, 2, 2}, {1, 1, 1, 3, 1},
          {1, 1, 2, 1, 1}, {1, 1, 2, 1, 2}, {1, 1, 2, 2, 1}, {1, 1, 3, 1, 1}, {1, 2, 1, 1, 1}, {1, 2, 1, 1, 2}, 
          {1, 2, 1, 2, 1}, {1, 2, 2, 1, 1}, {1, 3, 1, 1, 1}, {2, 1, 1, 1, 1}, {2, 1, 1, 1, 2}, {2, 1, 1, 2, 1}, 
          {2, 1, 2, 1, 1}, {2, 2, 1, 1, 1}, {3, 1, 1, 1, 1}}

if you want multi-indices disregarding the order, it is simply

In[42]:= With[{b = 7, n = 5, a = 1},
 If[a < 0 || n < 0 || b < 1 || b - n < a,  Print["You will get the empty set!"]];
 If[a < n, Print["Your conditions run empty on lower boundary!"]];
 Union[Sort /@ 
    Flatten[DeleteMissing[
      Table[If[a <= Plus @@ {a1, a2, a3, a4, a5} <= b, 
        X[a1, a2, a3, a4, a5], Missing[]],
       {a1, 1, b - 1}, {a2, 1, b - 1}, {a3, 1, b - 1}, {a4, 1, 
        b - 1}, {a5, 1, b - 1}], n]]] /. X -> List
 ]

During evaluation of In[42]:= Your conditions run empty on lower boundary!

Out[42]= {{1, 1, 1, 1, 1}, {1, 1, 1, 1, 2}, {1, 1, 1, 1, 3}, {1, 1, 1, 2, 2}}

this last thing brings you to an efficient way to do it: Start with the trivial allowed lower bound multiindex {1,1,1,..., Max[1,a-n+1]}, raise the last position until the limit b is reached ({1,1,1,...,b-n+1}) now distribute the IntegerPartitions/@ Range[Max[1,a-n+1],b-n+1-Max[1,a-n+1]] over the trivial allowed lower bound multiindex, with other words

In[51]:= With[{b = 7, n = 5, a = 1},
 If[a < 0 || n < 0 || b < 1 || b - n < a, Print["You will get the empty set!"]];
 If[a < n, Print["Your conditions run empty on lower boundary!"]];
 Join[{{1, 1, 1, 1, Max[1, a - n + 1]}},
  {1, 1, 1, 1, Max[1, a - n + 1]} + PadLeft[#, n] & /@ 
   Flatten[IntegerPartitions /@ 
     Range[Max[1, a - n + 1], b - n + 1 - Max[1, a - n + 1]], 1]]
 ]

During evaluation of In[51]:= Your conditions run empty on lower boundary!

Out[51]= {{1, 1, 1, 1, 1}, {1, 1, 1, 1, 2}, {1, 1, 1, 1, 3}, {1, 1, 1, 2, 2}}

if they are needed with respect to order, the result must be shuffled back to un-order.

Check for errors please, and find possibly other, even more efficient ways to do it!

7

This might be what you have in mind.

indicesFromLowToHigh[l_, h_, n_, var_] := Module[
  {vars = Array[var, n]},
  Table[{var[j], 1, h - n + j - Sum[var[k], {k, j - 1}]}, {j, n}]
  ]

Example:

indicesFromLowToHigh[2, 7, 3, x]

Out[361]= {{x[1], 1, 5}, {x[2], 1, 6 - x[1]}, {x[3], 1, 
  7 - x[1] - x[2]}}

For actual use as a sequence of iterators one would have something like Sequence@@indicesFromLowToHigh[...].

8

Thank you, that got me just what I needed! There was some bug fixing needed, but that was quite trivial once I had the right idea from your code. This is what I finally ended up with:

With[{b=10,n=5,a=8},
    If[a < 0 || n < 0 || b < 1 || b < a, Print["You will get the empty set!"]];
    If[a < n, Print["Your conditions run empty on lower boundary!"]];
    Flatten[Permutations /@ (ConstantArray[1,n] + PadLeft[#,n]& /@ 
        Flatten[IntegerPartitions[#,n]& /@ 
            Range[Max[0, a-n], b-n], 1]), 1]
]

Basically, your code worked fine when a < n, so all I had to do was fix the Range-call and replace the “initial vector” with all ones, and everything worked perfectly. As a nice side-effect, replacing the initial vector with a call to ConstantArray removed the only place in the function with a hard-coded value for n, allowing it to work for any n given as input.

Thanks again, your help was invaluable. =)

-Jonatan

P.S. Since the solution ended up taking quite a different form, Daniel’s solution wasn’t needed in the end, but thank you for that as well nonetheless, it was nice to learn how I might go about doing it if I find myself needing it later on.

P.P.S. I tested the efficiency, and with (a,b,n) = (10,20,5), your “brute force” solution took 15.7 seconds to run, the code above took 0.02 seconds. Very nice improvement, if I had gone about solving this myself I doubt I would’ve found anything better than the brute force solution with my limited knowledge.

9

This disrespects the lower bound if a > n

In[9]:= (* Lichtblau *)
Remove[indicesFromLowToHigh]
indicesFromLowToHigh[l_, h_, n_, var_] := 
 Module[{vars = Array[var, n]}, Table[{var[j], 1, h - n + j - Sum[var[k], {k, j - 1}]}, {j, n}]]

In[12]:= indicesFromLowToHigh[8, 10, 5, x]
Out[12]= {{x[1], 1, 6}, {x[2], 1, 7 - x[1]}, {x[3], 1, 8 - x[1] - x[2]}, 
          {x[4], 1, 9 - x[1] - x[2] - x[3]}, {x[5], 1, 10 - x[1] - x[2] - x[3] - x[4]}}

In[18]:= Select[
 Flatten[Table[X[x1, x2, x3, x4, x5], {x1, 1, 6}, {x2, 1, 7 - x1}, {x3, 1, 8 - x1 - x2}, {x4, 1, 9 - x1 - x2 - x3},
      {x5, 1, 10 - x1 - x2 - x3 - x4}]] /.  X -> List, (Plus @@ # < 8 || Plus @@ # > 10) &]

Out[18]= {{1, 1, 1, 1, 1}, {1, 1, 1, 1, 2}, {1, 1, 1, 1, 3}, {1, 1, 1, 2, 1}, {1, 1, 1, 2, 2}, {1, 1, 1, 3, 1}, 
          {1, 1, 2, 1, 1}, {1, 1, 2, 1, 2}, {1, 1, 2, 2, 1}, {1, 1, 3, 1, 1}, {1, 2, 1, 1, 1}, {1, 2, 1, 1, 2}, 
          {1, 2, 1, 2, 1}, {1, 2, 2, 1, 1}, {1, 3, 1, 1, 1}, {2, 1, 1, 1, 1}, {2, 1, 1, 1, 2}, {2, 1, 1, 2, 1},
          {2, 1, 2, 1, 1}, {2, 2, 1, 1, 1}, {3, 1, 1, 1, 1}}

but if they become deleted afterwards, it does fit

In[19]:= (* Lehtonen *)
lehto = With[{b = 10, n = 5, a = 8},
   If[a < 0 || n < 0 || b < 1 || b < a, Print["You will get the empty set!"]];
   If[a < n, Print["Your conditions run empty on lower boundary!"]];
   Flatten[Permutations /@ (ConstantArray[1, n] + PadLeft[#, n] & /@ 
       Flatten[IntegerPartitions[#, n] & /@ Range[Max[0, a - n], b - n], 1]), 1]
   ];

In[26]:= lichtb = 
  DeleteCases[
    Flatten[Table[
      If[8 <= Plus @@ {x1, x2, x3, x4, x5} <= 10, X[x1, x2, x3, x4, x5], Y[x1, x2, x3, x4, x5]], 
      {x1, 1, 6}, {x2, 1, 7 - x1}, {x3, 1, 8 - x1 - x2}, {x4, 1, 9 - x1 - x2 - x3}, 
      {x5, 1, 10 - x1 - x2 - x3 - x4}]], _Y] /. X -> List;

In[28]:= Sort[lehto] == Sort[lichtb]
Out[28]= True

so one could call this subtle brute force.

10

Sorry if I come back to this, but the Parma Polyhedron Library might be of interest in cases where IntegerPartitions[] might take too big memory requirements (usual around b - n = 80 in the language of the current problem on a 8 GB RAM computer). I tried to fit a Henrich package Multidimensional Convex Polyhedra and Searching for Magic Squares to this multi-index problem, but all I fiddled out was the null space, possibly because I’m unable to define an equation in this case, only inequalities are available.

11

An interesting addition, but in my work computing the multi-indices will never actually be the bottleneck, since both the computation time and memory usage for the remainder of my program is exponential in b.