Dynamic programming with large number of subproblems

Dynamic programming with large number of subproblems. So I'm trying to solve this problem from Interview Street:

Grid Walking (Score 50 points)
You are situated in an $N$-dimensional grid at position $(x_1,x_2,\dots,x_N)$. The dimensions of the grid are $(D_1,D_2,\dots,D_N$). In one step, you can walk one step ahead or behind in any one of the $N$ dimensions. (So there are always $2N$ possible different moves). In how many ways can you take $M$ steps such that you do not leave the grid at any point? You leave the grid if for any $x_i$, either $x_i \leq 0$ or $x_i > D_i$.

def number_of_ways(steps, starting_point):
global n, dimensions, mem
#print steps, starting_point
if (steps, tuple(starting_point)) in mem:
return mem[(steps, tuple(starting_point))]
val = 0
if steps == 0:
val = 1
else:
for i in range(0, n):
tuple_copy = starting_point[:]
tuple_copy[i] += 1
if tuple_copy[i] <= dimensions[i]:
val += number_of_ways(steps - 1, tuple_copy)
tuple_copy = starting_point[:]
tuple_copy[i] -= 1
if tuple_copy[i] > 0:
val += number_of_ways(steps - 1, tuple_copy)
mem[(steps, tuple(starting_point))] = val
return val


Big surprise: it fails for a large number of steps and/or dimensions due to a lack of memory.

So the next step is to improve my solution by using dynamic programming. But before starting, I'm seeing a major problem with the approach. The argument starting_point is an $n$-tuple, where $n$ is as large as $10$. So in fact, the function could be number_of_ways(steps, x1, x2, x3, ... x10) with $1 \leq x_i \leq 100$.

The dynamic programming problems I've seen in textbooks almost all have twp variables, so that only a two-dimensional matrix is needed. In this case, a ten-dimensional matrix would be needed. So $100^{10}$ cells in total.

With 2-D matrixes in dynamic programming, usually only the previous row of calculations is needed for the next calculation, hence reducing the spatial complexity from $mn$ to $\min(m,n)$. I'm not sure how I would do the same in this case. Visualizing a table isn't feasible, so the answer would have to come directly from the recursion above.

UPDATE

Using Peter Shor's suggestions, and making some minor corrections, notably the need to keep track of position in the $W(i, t_i)$ function, and rather than only splitting dimensions into two sets A and B, doing the splitting recursively, effectively using a divide-and-conquer method, until a base case is reached where only one dimension is in the set.

I came up with the following implementation, which passed all tests below the maximum execution time:

def ways(di, offset, steps):
global mem, dimensions
if steps in mem[di] and offset in mem[di][steps]:
return mem[di][steps][offset]
val = 0
if steps == 0:
val = 1
else:
if offset - 1 >= 1:
val += ways(di, offset - 1, steps - 1)
if offset + 1 <= dimensions[di]:
val += ways(di, offset + 1, steps - 1)
mem[di][steps][offset] = val
return val

def set_ways(left, right, steps):
# must create t1, t2, t3 .. ti for steps
global mem_set, mem, starting_point
#print left, right
#sleep(2)
if (left, right) in mem_set and steps in mem_set[(left, right)]:
return mem_set[(left, right)][steps]
if right - left == 1:
#print 'getting steps for', left, steps, starting_point[left]
#print 'got ', mem[left][steps][starting_point[left]], 'steps'
return mem[left][steps][starting_point[left]]
#return ways(left, starting_point[left], steps)
val = 0
split_point =  left + (right - left) / 2
for i in xrange(steps + 1):
t1 = i
t2 = steps - i
mix_factor = fact[steps] / (fact[t1] * fact[t2])
#print "mix_factor = %d, dimension: %d - %d steps, dimension %d - %d steps" % (mix_factor, left, t1, split_point, t2)
val += mix_factor * set_ways(left, split_point, t1) * set_ways(split_point, right, t2)
mem_set[(left, right)][steps] = val
return val

import sys
from time import sleep, time

fact = {}
fact[0] = 1
start = time()
accum = 1
for k in xrange(1, 300+1):
accum *= k
fact[k] = accum
#print 'fact_time', time() - start

num_tests = int(data.pop(0))
for ignore in xrange(0, num_tests):
n_and_steps = data.pop(0)
n, steps = map(lambda x: int(x), n_and_steps.split())
starting_point = map(lambda x: int(x), data.pop(0).split())
dimensions = map(lambda x: int(x), data.pop(0).split())
mem = {}
for di in xrange(n):
mem[di] = {}
for i in xrange(steps + 1):
mem[di][i] = {}
ways(di, starting_point[di], i)
start = time()
#print 'mem vector is done'
mem_set = {}
for i in xrange(n + 1):
for j in xrange(n + 1):
mem_set[(i, j)] = {}
#print time() - start

• "it fails for a large number of steps and/or dimensions" -- what does "fail" mean here?
– Raphael
Oct 8, 2012 at 7:54
• Welcome! I edited your question to a) use proper Markdown and LaTeX formatting (please to so yourself in the future) and b) remove superfluous gutter. We don't care for blurbs of C-code; please restrict yourself to ideas, that is pseudo code of the central things.
– Raphael
Oct 8, 2012 at 8:01
• Fails means it exhausts all available system memory by filling up the mem[] dictionary. And thank you for cleaning up my answer. Not too familiar with LaTeX but will make an effort next time. Oct 8, 2012 at 8:06
• You can find help on Markdown next to the editor box; see here for a primer on LaTeX.
– Raphael
Oct 8, 2012 at 8:14

The different dimensions are independent. What you can do is compute, for each dimension j, how many different walks there are in just that dimension which take $t$ steps. Let us call that number $W(j,t)$. From your question, you already know how to compute these numbers with dynamic programming.

Now, it's easy to count the number of walks that take $t_i$ steps in dimension $i$. You have $N \choose t_1,t_2, \ldots, t_M$ ways of interspersing dimensions so that the total number of steps taken in dimension $i$ is $t_i$, and for each of these ways you have $\Pi_1^N W(i,t_i)$ walks. Sum over these to get $$\sum_{t_1+t_2+\ldots+t_N=M}{M \choose t_1,t_2,\ldots,t_N}\ \Pi_{i=1}^N W(i,t_i).$$ Now, the memory is under control, since you only need to remember the values $W(j,t)$. The time grows superpolynomially for large $N$, but most computers have a lot more time than memory.

You can do even better. Recursively divide the dimensions into two subsets, $A$ and $B$, and compute how many walks there are using just the dimensions in subset $A$, and just those in $B$. Call these numbers $W_A(t)$ and $W_B(t)$, respectively. You get the total number of walks by

$$\sum_{t_1 + t_2 = M} {M \choose t_1} \, W_A(t_1) W_B(t_2).$$

• Hi Peter. Alright, that was the missing insight. Now I only have one doubt left. The outer sum iterates over all possible combinations of t1,t2,...tn that sum to M. Unfortunately, the number of such combinations is C(M+1, N-1), which can be as high as C(300+1, 10-9). Very large number... :( Oct 8, 2012 at 15:44
• @Alexandre: My second algorithm (starting with "You can do even better") doesn't have that problem. I left the first algorithm in my answer because it's the first one I came up with, and because I think it's much easier to explain the second algorithm as a variant of the first than just giving it without any motivation. Oct 8, 2012 at 16:00
• I implemented the second algorithm. It's faster, but still too low for the largest bounds. The problem with the first one was iterating over all possibilites of t1, t2, t3, ... tn that summed to M. The second algorithm only iterates over solutions to t1 + t2 = M. But then the same must be done for Wa(t1), iterating over solutions to t1' + t2' = t1. And so forth recursively. Here's the implementation in case you're insterested: pastebin.com/e1BLG7Gk . And in the second algorithm, there multinomial should be M over t1,t2 no? Oct 9, 2012 at 10:36
• Nevermind! Solved it! Had to use memoization in the set_ways function as well. Here's the final solution, which is blazing fast! pastebin.com/GnkjjpBN Thank you for your insight Peter. You made both key observations: problem independence and divide-and-conquer. I recommend people look at my solution because there are some things which are not in the answer above, such as the W(i,ti) function needing a third argument, which is the position. That has to be calculated for values combinations of i, ti, and position. If you can, also add t2 the the multinomial in your second algorithm. Oct 9, 2012 at 12:33

Let's extract a formula for $\newcommand{\now}{\operatorname{now}}\now(s,x_1,\dots,x_n)$ from your code (for an inner cell, that is ignoring border cases):

\qquad \begin{align} \now(s,x_1,\dots,x_n) = \phantom{+} &\sum_{i=0}^n \now(s-1,x_1,\dots,x_{i-1},x_i + 1, x_{i+1},\dots,x_n) \\ + &\sum_{i=0}^n\now(s-1,x_1,\dots,x_{i-1},x_i - 1, x_{i+1},\dots,x_n) \end{align}

Here are some ideas.

• We see that once you have computed all values for $s=k$, you can drop all computed values for $s<k$.
• For a fixed $s$, you should compute the table entries in lexicographic order (just because it's simple). Then, note that every cell only needs such cells inside a "radius of one", that is no coordinate can be farther away than one. Therefore, once your iteration hits $x_1=i$, you can drop all values for $x_1\leq i-2$. If that is not enough, do the same for $x_2$ -- for fixed $x_1=i$, drop values with $x_1 = i$ and $x_2 \leq j-2$ when $x_2=j$ is reached -- and so on.
• Note that "so there are always $2N$ possible different moves" holds only in the middle of the grid, that is if $x_i - M > 0$ and $x_i + M < D_i$ for all $i$. But that also means that the answer is easy in the middle: it's just $(2N)^M$. If you had a working dynamic programming recurrence, that alone would allow you shave away most of the table (if $M\ll N$).
• Another thing to note is that you don't have to compute the whole table; most values will be filled with $0$ anyway (if $M\ll N$). You can restrict yourself to the (hyper)cube of edge-length $2M$ around $x$ (note that it will be dented because of paths leaving the grid).

That should be sufficient to keep memory usage quite low.

• Hi Raphael, let's say our goal is now(3, 3, 3, 3), on a grid 5x5x5. Using dynamic programming and using lex order as you suggested, we would calculate now(0, 0, 0, 0), then (0, 0, 0, 1), ... now(0, 5, 5, 5). At what point we could discard now(0, 0, 0, 0) (which is more than a radius of one away from (5, 5, 5), since we will now need it to calculate now(1, 0, 0, 0), now(1, 0, 0, 1), etc? You mentioned M << N a couple of times, but the bounds are 1 <= M <= 300, and 1 <= N <= 10. So, at the extremes, it doesn't seem that 1 << 300. Oct 8, 2012 at 9:51
• 1) What is unclear in my second bullet? As soon as you calculate $(2,0,0,0)$, you can discard $(0,\*,\*,\*)$. That is not the earliest point where you can discard $(0,0,0,0)$, though; the last cell you need it for is $(1,0,0,0)$. 2) I am not too much concerned about your specific values for $M$ and $N$, to be honest. I'd rather look at the general problem. If you don't have $M \ll N$, the last two bullets won't help you much. $M=1$ and $N=10$ should suffice to notice the effect, though, and neither strategy ever hurts.
– Raphael
Oct 8, 2012 at 10:13
• The 1) bullet I understand. That reduces spatial complexity from M*D^N to D^N, but D^N is still too much. I'm not quite seeing how the 2) bullet works though. Can you use the example in my comment to illustrate it? Oct 8, 2012 at 10:19
• @Alexandre I did in my earlier comment. If I read $D$ as meaning $\max_{i=1,\dots,N} D_i$, then applying the second bullet once reduces space complexity to $D^{N-1}$, the second time to $D^{N-2}$ and so on. (More precisely, it goes from $\prod_{i=1}^N D_i$ to $\prod_{i=2}^N D_i$ and so on.)
– Raphael
Oct 8, 2012 at 10:24
• Not quite grasping how to do it... Let's say I did understand, and that I reduced spatial complexity to D. Fundamentally, won't M*D^N subproblems still need to be solved? Isn't an additional property needed to make the problem polynomial? Oct 8, 2012 at 10:33