Monday, June 30, 2014

ICPC Problem H: Pachinko

Problem A has been covered quite well elsewhere, so I won't discuss it. That leaves only problem H. I started by reading the Google translation of a Polish writeup by gawry. The automatic translation wasn't very good, but it gave me one or two ideas I borrowed. I don't know how my solution compares in length of code, but I consider it much simpler conceptually.

This is a fairly typical Markov process, where a system is in some state, and each timestep it randomly selects one state as a function of the current state. One variation is that the process stops once the ball reaches a target, whereas Markov processes don't terminate. I was initially going to model that as the ball always moving from the target to itself, but things would have become slightly complicated.

Gawry has a nice way of making this explicitly a matrix problem. Set up the matrix M as for a Markov process i.e., \(M_{i,j}\) is the probability of a transition from state j to state i. However, for a target state j, we set \(M_{i,j}=0\) for all i. Now if \(b\) is our initial probability vector (equal probability for each empty spot in the first row), then \(M^t b\) represents the probability of the ball being in each position (and the game not having previously finished) after \(t\) timesteps. We can then say that the expected amount of time the ball spends in each position is given by \(\sum_{t=0}^{\infty} M^t b\). The sum of the elements in this vector is the expected length of a game and we're told that it is less than \(10^9\), so we don't need to worry about convergence. However, that doesn't mean that the matrix series itself converges: Gawry points out that if there are unreachable parts of the game with no targets, then the series won't converge. We fix that by doing an initial flood-fill to find all reachable cells and only use those in the matrix. Gawry then shows that under the right conditions, the series converges to \((I - M)^{-1} b\).

This is where my solution differs. Gawry dismisses Gaussian elimination, because the matrix can be up to 200,000 square. However, this misses the fact that it is banded: by numbering cells left-to-right then top-to-bottom, we ensure that every non-zero entry in the matrix is at most W cells away from the main diagonal. Gaussian elimination (without pivoting) preserves this property. We can exploit this both to store the matrix compactly, and to perform Gaussian elimination in \(O(W^3H)\) time.

One concern is the "without pivoting" caveat. I was slightly surprised that my first submission passed. I think it is possible to prove correctness, however. Gaussian elimination without pivoting is known (and easily provable) to work on strictly column diagonally dominant matrices. In our case the diagonal dominance is weak: columns corresponding to empty cells have a sum of zero, those corresponding to targets have a 1 on the diagonal and zeros elsewhere. However, the matrix is also irreducible, which I think is enough to guarantee that there won't be any division by zero.

EDIT: actually it's not guaranteed to be irreducible, because the probabilities can be zero and hence it's possible to get from A to B without being able to get from B to A. But I suspect that it's enough that one can reach a target from every state.

4 comments:

mahbub said...

Wondering if you can slightly describe gawry's method, since it is a bit difficult to understand from the translated text.

Bruce Merry said...

@mahbub: I didn't follow it fully either, but Petr has a description on his blog.

Godmar Back said...

Can you give a hint as to what representation you're using for M and b?

I tried Gaussian with a COO-style representation for both, along with column indices, but can't get it under 1:22m for test2.in. Though (I-M) stays reasonable sparse, the 'b' side of the augmented matrix quickly fills in.

Also, do you exploit the fact that you need at most the first 20 values of (I-M)^-1 * b?

I also tried LU factorization of the entire matrix (which has roughly 2M nonzeros for both L/U), which would require roughly 18 GFlop to solve for ~8,000 rhs - so no go.

Godmar Back said...

Oh, wait a minute. I think I see the issue. I was - naively - trying to compute the absorption probabilities for all initial states to all absorbing states (of which there are over 8,000 for the largest test case), then obtain the aggregate absorption probabilities by summing over the start states.

After reading SnapDragon's post, I realize that you instead compute the expected number of times the chain is in a state directly via a linear system of equations with only one right-hand side, rather than 20x8000 individual probabilities. That should fit the time envelope.