A graded lexicographic index, part 2

The previous post on this topic covered basic code for producing a graded lexicographic index (a/k/a the state index) and its inverse. This post will cover some more advanced techniques. The first one of these is simple: as de Boor pointed out (and it took me quite a while to catch on) the state index I(\alpha) of \alpha \in X_{B,b} can be given by a simple formula. Let PA denote the matrix given by the following MATLAB code:

mb=max(b+1,B+1);  P1=pascal(mb);  P=P1(1:B-1,:);  PA=[zeros(1,mb);P];

For example, with B = 4, b = 8 this is

PA = \left( \begin{smallmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 \\ 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 \\ 1 & 3 & 6 & 10 & 15 & 21 & 28 & 36 & 45 \end{smallmatrix} \right).

(Identifying the elements of the Pascal matrix with binomial coefficients would introduce some minor technical complications w/r/t the top row, and so we avoid this.) Now let

\alpha \equiv \sum_{k=1}^K \alpha_{j_k} \delta_{j_k}, \ \ \alpha_{j_k} > 0 \ \ \forall k.

Then

I(\alpha) = 1 + \sum_{k=1}^K \sum_{m=1}^{\alpha_{j_k}} PA \left( B + 1 - j_k, y_{km}(\alpha)\right)

where

y_{km}(\alpha) = b + 2 - m - \sum_{\ell = 0}^{k-1} \alpha_{j_\ell}

and an empty sum is evaluated as zero. This is embodied in the following MATLAB code:

J = find(alpha);
K = length(J);

for k = 1:K
   temp1 = sum(alpha(J(1:(k-1))));
   for m = 1:alpha(J(k))
      stateindex = stateindex + PA(B+1-J(k),b+2-temp1-m);
   end
end

stateindex = stateindex + 1;

This is hard to beat in general. However, for some applications it is sufficient to compute the state index at a neighboring site. And this can actually be done much more effectively in practice.

A deceptively simple algorithm always produces a superset of the neighborhood of a state in bucketspace. Write

\tau_k (\alpha) := 1+ \sum_{j=1}^k \alpha_{B-j+1}

and for a finite subset Y of \mathbb{R}, let \mu_k(Y) and \mu^k(Y) respectively denote the k least and greatest elements of Y. Now set I_1 := \{I(\alpha)\} and given I_k for k > 0 form

I_{k+1} := I_{k+1}^- \cup \{I_k\} \cup I_{k+1}^+

where

I_{k+1}^- := \{ \left( \mu_{k-1}(I_k) \cup I_1 \right) - PA(k, \tau_k) \}

and

I_{k+1}^+ := \{ \left( \mu^{k-1}(I_k) \cup I_1 \right) - PA(k, \tau_k - 1) \}.

Then the state indices of the nearest neighbors of \alpha are contained in the set I_B. It is not hard to see that this algorithm requires \sim 3B^2 simple operations (e.g., lookup, addition/subtraction, min/max) if we store the (fairly small) Pascal matrix. On “interior” states (i.e., states with no zero entries), this algorithm produces the exact neighborhood. However, on “boundary” states (i.e., states with at least one zero entry), some spurious indices are returned. It is not obvious (especially given its structure) how the algorithm might be refined to deal with this problem directly.

However, the spurious indices can themselves be checked efficiently. Given B and b, the following MATLAB code illustrates how to invert the state index:

PA2 = flipud(fliplr(PA));
alpha = zeros(1,B);
ind0 = 0; % init
for j = 1:B
   temp1 = cumsum(PA2(j,(ind0+1):(mb-1)));
   ind1 = max(find(temp1<s));
   if ind1
      temp2 = temp1(ind1);
      s = s - temp2;
      alpha(j) = ind1;
      ind0 = ind0 + ind1;
   end
end

The inversion essentially takes O(b) operations. With this tool in hand, we can trivially establish which state indices are spurious and which are not.

In the next post on this topic, we’ll produce a technique which is actually much faster than any of these for updating state indices in both MATLAB and Java, even after accounting for considerable care in the optimization (such as accounting for the order in which bitshifts are executed).

This entry was posted on Monday, October 19th, 2009 at 00:01 and is filed under Equilibrium Networks, Mathematics, Nonrandom bits. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

One Response to “A graded lexicographic index, part 2”

  1. A graded lexicographic index, part 3 « Equilibrium Networks Says:
    26 October 2009 at 16:35 | Reply

    [...] graded lexicographic index, part 3 In an earlier post we covered some more sophisticated ways of implementing a graded lexicographic index, a/k/a the [...]

Leave a Reply