Optimal strategy for wandering robots

Question

If the agent is at the origin it always moves right when you choose it.

The goal is to get one agent to 10.

The cost is the sum of the square of the total number of moves (left or right) taken by each individual agent until the first agent gets to 10. You want to minimise the expected cost

Once you have fixed a strategy you should run your code ten thousand times (each time up to the point where the first agent gets to 10) and report the average cost.

Input and output

There is no input to this challenge. Your code need only output the mean cost over 10000 100,000 runs as specified. However please do include an explanation of your strategy for choosing which agent to move when.

In short, the cost of moving a particular agent effectively increases each time you move it. This is a result of the scoring system. As another example, if you move agent 1 twice for your first two moves and it ends up back at the origin, it is now cheaper to move any other agent.

If two answers have the same strategy then the first answer wins. If two different strategies give similar mean costs then they should be rerun with 100,000 a million (or even more) iterations to see which is better.

A very simple strategy gives a mean cost of 6369. So the first hurdle is to do better than that.

Winners

This question was first brilliantly solved optimally by Anders Kaseorg.

Answer 1

Ruby, modified minimum product, average cost ~4950

f=->a,b{(0..3).min_by{|q|(9-a[q])**3*b[q]}}

Try it online!

Always select the agent that minimizes (9-position)^3*(number of moves):

• On the first round, every agent makes a move.
• By using 9-a instead of 10-a, any agent in position 9 will try to move on the next iteration.
• Position seems to be more important than the cost projection, so after trying out different values, I found out that the third power is the best way to represent it.

The easiest strategy (round-robin) has a cost of about 7000.

If I let it run in a continuous loop, after a while the average stabilizes around 4950.

Answer 2

Rust, cost = 4928.296960314256

Computed exactly (modulo floating point imprecision) with a non-random algorithm, but also measured as \$4923.8 ± 2.6\$ (\$1σ\$ confidence) with \$10^7\$ random runs of the actual game as a check.

Strategy

Consider the following equivalent but wetter version of the game. There’s a tower of cells \$(p, m)\$ with eleven columns \$0 ≤ p ≤ 10\$ and infinitely many rows \$m ≥ 0\$ going downwards. We’ll call the first ten columns \$p = 0, \dotsc, 9\$ “live” and the last column \$p = 10\$ “dead”. Each live cell has a plugged drain in the floor, connected by equal-width pipes to its diagonal neighbors below. We start with \$1\$ liter of water in \$(0, 0)\$, with \$4\$ grains of sand floating in the water. Our goal is to control the plugs until one of the grains ends up in the dead column, while minimizing the expected sum of the squared row numbers of the grains at that time.

The reason to formulate the game like this is that our strategy does not need to look at the grains. Instead, we’re going to open the plugs in some fixed order, letting the grains float where they may, until 100% of the water flows into the dead column in the limit as all of the plugs are open. (We’re very good at opening plugs.)

I’ll explain in a moment how to decide which order we’ll use, but to motivate that, I first need to explain how to evaluate our expected score. Define “time \$t\$” to be the time when \$t\$ liters of water have flowed into the dead column. Then, because the grains of sand flow equivalently to water molecules, the distribution of time quadruples \$(t_1, \dotsc, t_4)\$ when the four grains enter the dead column is just the uniform distribution over \$[0, 1]^4\$. Therefore, the time \$t_\min = \min \{t_1, \dotsc, t_4\}\$ when the score should be evaluated has a beta distribution \$t_\min \sim \mathrm{Beta}(1, 4)\$, whose PDF is \$4(1 - t)^3\$.

Let \$q_{(p, m)}(t)\$ be the amount of water in cell \$(p, m)\$ at time \$t\$; let \$c_{(p, m)}(t) = m^2q_{(p, m)}(t)\$ be its weighted cost, and \$c_{\textrm{live}}(t)\$ (resp. \$c_{\textrm{dead}}(t)\$) the sum over all live (resp. dead) cells. Then we can write our expected score as

\begin{split} &\int_0^1 4(1 - t)^3 \left(3\frac{c_{\textrm{live}}(t)}{1 - t} + c'_{\textrm{dead}}(t)\right)\,dt \\ &= \int_0^1 12(1 - t)^2 c_{\textrm{live}}(t)\,dt + \int_0^1 4(1 - t)^3 c'_{\textrm{dead}}(t)\,dt \\ &= \Bigl[-4(1 - t)^3 c_{\textrm{live}}(t)\Bigr]_0^1 - \int_0^1 -4(1 - t)^3 c'_{\textrm{live}}(t)\,dt + \int_0^1 4(1 - t)^3 c'_{\textrm{dead}}(t)\,dt \\ &= \int_0^1 4(1 - t)^3 (c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t))\,dt. \end{split}

Therefore, we see that our goal at each step is to minimize the total growth rate of the weighted cost of all cells (\$c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t)\$), relative to the rate of flow into the dead column (i.e. growth of \$t\$ itself). \$4(1 - t)^3\$ is decreasing, so we prefer to have smaller ratios earlier. This makes our strategy obvious: at each step, we’ll open the plug that minimizes this ratio, wait until all water has flowed out of that cell, and leave it open while proceeding to the following cells.

Since \$c'_{\textrm{live}}(t) + c'_{\textrm{dead}}(t)\$ is piecewise constant, the integral evaluates to an (infinite) sum, one that converges much more quickly than the Monte Carlo simulation. That’s how I determined the precise answer above.

Code

Build and run with cargo run --release.

Cargo.toml

[package]
name = "wandering"
version = "0.1.0"
edition = "2021"

[dependencies]
rand = "0.8.5"

src/main.rs

use rand::prelude::*;
use rand::rngs::StdRng;
use std::collections::HashMap;

#[derive(Default, Debug)]
struct Rate {
    dead: f64,
    cost: f64,
    live: HashMap<(u32, u32), f64>,
}

fn add(map: &mut HashMap<(u32, u32), f64>, state: (u32, u32), value: f64) {
    if value != 0.0 {
        *map.entry(state).or_insert(0.0) += value;
    }
}

struct SortingHat {
    distribution: HashMap<(u32, u32), f64>,
    moves_frontier: [u32; 10],
    rates: HashMap<(u32, u32), Rate>,
    orders: HashMap<(u32, u32), usize>,
}

impl Default for SortingHat {
    fn default() -> SortingHat {
        let mut distribution = HashMap::new();
        distribution.insert((0, 0), 1.0);
        SortingHat {
            distribution,
            moves_frontier: [0_u32, 1, 2, 3, 4, 5, 6, 7, 8, 9],
            rates: HashMap::new(),
            orders: HashMap::new(),
        }
    }
}

impl SortingHat {
    fn get_order(&mut self, query_position: u32, query_moves: u32) -> usize {
        if let Some(&order) = self.orders.get(&(query_position, query_moves)) {
            return order;
        }

        loop {
            let (_, position, moves) = (0_u32..10)
                .zip(&self.moves_frontier)
                .filter(|&(position, &moves)| {
                    position + 1 == 10 || self.moves_frontier[position as usize + 1] > moves + 1
                })
                .map(|(position, &moves)| {
                    let mut cost = 0.0;
                    let mut dead = 0.0;
                    for new_position in [position.abs_diff(1), position + 1] {
                        cost += moves as f64 + 1. / 2.;
                        if new_position == 10 {
                            dead += 1. / 2.;
                        } else if let Some(new_rate) = self.rates.get(&(new_position, moves + 1)) {
                            dead += new_rate.dead / 2.;
                            cost += new_rate.cost / 2.;
                        }
                    }
                    (cost / dead, position, moves)
                })
                .min_by(|a, b| a.partial_cmp(&b).unwrap())
                .unwrap();

            self.orders.insert((position, moves), self.orders.len());

            let mut rate = Rate::default();
            for new_position in [position.abs_diff(1), position + 1] {
                rate.cost += moves as f64 + 1. / 2.;
                if new_position == 10 {
                    rate.dead += 1. / 2.;
                } else if let Some(new_rate) = self.rates.get(&(new_position, moves + 1)) {
                    rate.dead += new_rate.dead / 2.;
                    rate.cost += new_rate.cost / 2.;
                    for (&newer_state, live_rate) in &new_rate.live {
                        add(&mut rate.live, newer_state, live_rate / 2.);
                    }
                } else {
                    add(&mut rate.live, (new_position, moves + 1), 1. / 2.);
                }
            }

            if let Some(mass) = self.distribution.remove(&(position, moves)) {
                for (&new_state, live_rate) in &rate.live {
                    add(&mut self.distribution, new_state, live_rate * mass);
                }
            }

            for old_rate in self.rates.values_mut() {
                if let Some(live_rate) = old_rate.live.remove(&(position, moves)) {
                    old_rate.dead += rate.dead * live_rate;
                    old_rate.cost += rate.cost * live_rate;
                    for (&new_state, new_live_rate) in &rate.live {
                        *old_rate.live.entry(new_state).or_insert(0.0) += new_live_rate * live_rate;
                    }
                }
            }

            self.rates.remove(&(position + 1, moves + 1));
            self.rates.insert((position, moves), rate);
            self.moves_frontier[position as usize] += 2;

            if (position, moves) == (query_position, query_moves) {
                return self.orders.len() - 1;
            }
        }
    }

    fn strategy(&mut self, positions: &[u32; 4], num_moves: &[u32; 4]) -> usize {
        (0..4)
            .min_by_key(|&agent| self.get_order(positions[agent], num_moves[agent]))
            .unwrap()
    }
}

fn test(mut rng: impl Rng, mut strategy: impl FnMut(&[u32; 4], &[u32; 4]) -> usize, runs: u32) {
    let mut count = 0;
    let mut mean = 0.0;
    let mut m2 = 0.0;
    let mut mean_stdev = f64::NAN;

    for _ in 0..runs {
        let mut positions = [0; 4];
        let mut num_moves = [0; 4];
        loop {
            let agent = strategy(&positions, &num_moves);
            if positions[agent] == 0 || rng.gen() {
                positions[agent] += 1;
            } else {
                positions[agent] -= 1;
            }
            num_moves[agent] += 1;
            if positions[agent] == 10 {
                break;
            }
        }
        let sample: f64 = num_moves.into_iter().map(|m| (m as f64).powi(2)).sum();

        count += 1;
        let delta = sample - mean;
        mean += delta / count as f64;
        m2 += delta * delta * (1. - 1. / count as f64);
        mean_stdev = (m2 / (count as f64 * (count - 1) as f64)).sqrt();
        if count % 100 == 0 {
            eprint!("\r\x1b[K{count}: {mean} ± {mean_stdev}");
        }
    }
    eprint!("\r\x1b[K");
    println!("{count}: {mean} ± {mean_stdev}");
}

fn main() {
    let rng = StdRng::seed_from_u64(0);
    let mut sorting_hat = SortingHat::default();
    test(
        rng,
        |positions, num_moves| sorting_hat.strategy(positions, num_moves),
        10000000,
    )
}

Answer 3

Python 3 (PyPy), 4957.6 +/- 2.6 (10M trials)

def expco(positions, num_moves,depth=0) -> float:
    if depth==0:
        return min((10**2 - positions[i] ** 2)
            * (5 * 10**2 - positions[i] ** 2 + 6 * num_moves[i] - 2)/3 for i in range(4))
    return (10 not in positions)and min(
       expco([abs(p-(i==j)) for j,p in enumerate(positions)],
             [m+(i==j) for j,m in enumerate(num_moves)],depth-1)/2 +
       expco([p+(i==j) for j,p in enumerate(positions)],
             [m+(i==j) for j,m in enumerate(num_moves)],depth-1)/2 +
       num_moves[i]*2+1 for i in range(4))
def strategy(positions, num_moves, depth=1) -> int:
    return min(range(4), key = lambda i:
       expco([abs(p-(i==j)) for j,p in enumerate(positions)],
             [m+(i==j) for j,m in enumerate(num_moves)],depth)/2 +
       expco([p+(i==j) for j,p in enumerate(positions)],
             [m+(i==j) for j,m in enumerate(num_moves)],depth)/2 +
       num_moves[i]*2+1)

Try it online!

This is a simple brute force search on top of @Anders Kaseorg's formula (and code).

The score given was obtained with search depth 1, i.e. for each agent compute the expected cost if we can switch agents once immediately after the first move and take the best.

Apologies for the messy code.

Answer 4

Python 3, cost ≈ 5034.3 ± 2.7

\$1σ\$ confidence, evaluated with \$10^7\$ runs.

def strategy(positions, num_moves):
    return min(
        range(4),
        key=lambda i: (10**2 - positions[i] ** 2)
        * (5 * 10**2 - positions[i] ** 2 + 6 * num_moves[i] - 2),
    )

Try it online! (\$10^5\$ runs)

How it works

If we were to pick one agent at position \$p\$ that’s previously moved \$m\$ times, and use it for the rest of the game, the expected additional cost would be exactly

$$f(p, m) = \frac13(10^2 - p^2)(5⋅10^2 - p^2 + 6m - 2).$$

(You can verify that

$$m^2 + f(p, m) = (m + 1)^2 + \frac12(f(\lvert p - 1\rvert, m + 1) + f(p + 1, m + 1))$$

and \$f(10, m) = 0\$.)

At each step, we pick the agent that minimizes this quantity. (We don’t restrict ourselves to using it for the rest of the game.)

Answer 5

5050? ± 26.64

def strategy(positions, num_moves):
  E = [100, 99, 96, 91, 84, 75, 64, 51, 36, 19]
  _, n = min(((E[p] + m) ** 2 - m ** 2, i) for p, m, i in zip(positions, num_moves, range(4)))
  return n

Try it online!

Use minimal \$\left(E(P_i)+M_i\right)^2-M_i^2\$, where \$E(X)\$ is a pre-calculated array $$\left[E(0),\dots,E(9)\right]=\left[100, 99, 96, 91, 84, 75, 64, 51, 36, 19\right]$$

generated by following codes

syms E0 E1 E2 E3 E4 E5 E6 E7 E8 E9 EA
solve([E0-E1-1, E1-(E0+E2)/2-1, E2-(E1+E3)/2-1, E3-(E2+E4)/2-1, E4-(E3+E5)/2-1, E5-(E4+E6)/2-1, E6-(E5+E7)/2-1, E7-(E6+E8)/2-1, E8-(E7+E9)/2-1, E9-(E8+EA)/2-1, EA])

Try it online!

The driver code is used from other posts.

Answer 6

Python 3, cost ≈ 4928.4 ± 1

I believe this is essentially optimal. Thank you to @usul for the inspiration.

import functools

f = lambda p, m: (1. / 3) * (100 - p * p) * (500 - p * p + 6 * m - 2)

@functools.lru_cache(None)
def gittins_value(p, m, l):
    if m >= 400:
        return f(p, m)
    if p == 0:
        v_next = gittins_value(p + 1, m + 1, l)
    else:
        v_next = 0.5 * gittins_value(p - 1, m + 1, l)
        if p != 9:
            v_next += 0.5 * gittins_value(p + 1, m + 1, l)
    return min((2 * m + 1) + v_next, l)

@functools.lru_cache(None)
def gittins_index(p, m):
    low = 2 * m + 1
    high = f(p, m)
    while True:
        mid = 0.5 * (low + high)
        v = gittins_value(p, m, mid)
        if v < mid:
            high = mid
        else:
            low = mid
        if abs(high - low) <= 1e-6:
            break
    return 0.5 * (low + high)

def strategy(p, m):
    indices = [gittins_index(p_i, m_i) for p_i, m_i in zip(p, m)]
    return indices.index(min(indices))

Answer 7

R / Rcpp: 5006 ± 7.31, or 4932 ± 3.63

Score computed by taking the average of 20 million simulated games. Bounds are a 95% bootstrap confidence interval with 10,000 iterations.

Thanks to Anders Kaseorg for pointing out a typo in my code that caused the solution to accidentally be optimal. Reader, may all your future bugs behave thusly. The following description explains the 5006 solution with a note for how to get the 4932 solution.

Strategy:

Let's consider a simplified version of the environment in order to construct a risk heuristic. We'd like to be able to treat n consecutive movements of an agent as independent, but that won't be true if the agent ever reaches position 0 or 10 during those movements.

Suppose instead that each of the number lines is infinite in both directions, and that the probability of moving right is always 0.5 (even at position 0). In this simplified case, if we move a single agent n consecutive times, we may treat the resulting k number of movements to the right as binomially distributed.

At each step, let n = 10 - x, where x is the agent's position on the number line. Since the agent's movements are independent in the simplified environment, we can calculate the probability that the agent has moved backwards given k moves to the right as

$$Prob = 1 - \sum_{k=\lceil \frac{n}{2} \rceil}^n \dbinom n k \cdot 0.5^k \cdot (1-0.5)^{n-k}$$

Note: If you incorrectly move the one-minus inside the sum, the solution becomes near optimal for reasons that are unclear.

Let the cost function be the increase in cost accrued by moving the agent n times:

$$ Cost = (x + n)^2 - (x)^2 $$

Then define the risk of moving an agent as

$$ Risk = Prob \cdot Cost $$

and at each step move the agent with the lowest risk.

Intuition:

The idea is to calculate the risk that over the minimum number of steps an agent would need to reach the goal, the agent ends up moving backwards instead of forwards, i.e. We try to maximize the chance that the chosen agent will move closer to the goal. Minimizing this risk heuristic appears to perform fairly well even though the method does not take into account the effect of the order of the movements, nor the effect of reaching the zero position.

#include <Rcpp.h>
using namespace Rcpp;

//* Initialize Environment
// [[Rcpp::export]]
NumericMatrix initEnv() {
  return(NumericMatrix (4, 2));
}


//* Update Environment
// [[Rcpp::export]]
NumericMatrix updateEnv(NumericMatrix x, int a) {

  //Update move counter
  x(a, 1) = x(a, 1) + 1;
  
  //Update position
  if(x(a, 0) == 0){
    x(a, 0) = 1;
  } else {
    RNGScope scope;
    NumericVector r = Rcpp::runif(1, 0.0, 1.0);
    if(r[0] <= 0.5){
      x(a, 0) = x(a, 0) - 1;
    } else {
      x(a, 0) = x(a, 0) + 1;
    }
  }
  
  //Return
  return(x);
  
}


//* Calculate Loss
// [[Rcpp::export]]
double envCost(NumericMatrix x) {
  double lambda = 0;
  for(int i = 0; i < 4; i++){
    lambda = lambda + pow(x(i, 1), 2);
  }
  return(lambda);
}


//* Calculate Risk
// [[Rcpp::export]]
double agentRisk(NumericMatrix x, int a) {
  
  //Parameters
  int pos = x(a, 0);
  int n = 10 - pos;
  
  //Calculate risk correctly
  double prob = 0;
  double cost = pow(x(a, 1) + n, 2) - pow(x(a, 1), 2);
  for(int k = ceil(n / 2); k <= n; k++){
    prob = prob + (Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k));
  }
  prob = 1 - prob;

  //Calculate risk with typo
  //double prob = 0;
  //double cost = pow(x(a, 1) + n, 2) - pow(x(a, 1), 2);
  //for(int k = ceil(n / 2); k <= n; k++){
  //  prob = prob + (1 - Rf_choose(n, k) * pow(0.5, k) * pow(0.5, n - k));
  //}
  
  //Return
  return(prob * cost);
  
}


//* Simulate Game
// [[Rcpp::export]]
List simGame() {
  
  //Generate environment
  NumericMatrix en = initEnv();
  
  //Complete game
  bool done = false;
  while(!done){
    
    //Parameters
    int a = 0;
    double lambda = agentRisk(en, 0);
    
    //Calculate risks and choose action
    for(int i = 1; i < 4; i++){
      double lambda_new = agentRisk(en, i);
      if(lambda_new < lambda){
        a = i;
        lambda = lambda_new;
      }
    }
    
    //Update game
    en = updateEnv(en, a);
    
    //Check victory conditions
    for(int i = 0; i < 4; i++){
      if(en(i, 0) >= 10){
        done = true;
        break;
      }
    }
    
  }
  
  //Return
  return(
    List::create(
      Named("state") = en,
      Named("cost") = envCost(en)
    )
  );
  
}

Answer 8

Rust, cost = ~5032 (1M iterations)

const depth: usize = 1200;
const iterations: usize = 1_000_000;

lazy_static! {
    static ref EXPECTED_SCORES: [[usize; 10]; depth] = {
        let mut expected_scores = [[1_000_000usize; 10]; depth];

        for i in (0..depth-1).rev() {
            for j in 0..10 {
                if j == 0 {
                    expected_scores[i][j] = expected_scores[i + 1][j + 1].saturating_add(2 * i + 1);
                } else if j == 9 {
                    expected_scores[i][j] =
                        (expected_scores[i + 1][j - 1] / 2).saturating_add(2 * i + 1);
                } else {
                    expected_scores[i][j] = (expected_scores[i + 1][j + 1]/ 2)
                    .saturating_add(expected_scores[i + 1][j - 1] / 2)
                    .saturating_add(2 * i + 1);
                }
            }
        }

        expected_scores
    };
}

#[derive(Copy, Clone, Eq, PartialEq, Debug)]
struct Agent {
    position: u8,
    actions: u32,
}

fn strategy(agents: [Agent; 4]) -> usize {
    return (0..4)
        .min_by_key(|i| EXPECTED_SCORES[(agents[*i].actions as usize).min(depth - 1)][agents[*i].position as usize])
        .unwrap();
}

Basically pre-computes a table of expected scores then uses that.

Playground link Need to reduce iterations if you don't want it to time out.

Answer 9

Python 3 (PyPy), 5027.1 ± 8.3

def strategy(positions, num_moves) -> int:
    return min(
        range(4),
        key=lambda i: (10 - positions[i] + 2*num_moves[i]) * (10 - positions[i]),
    )

Try it online!

Copies the function/testing template from Anders Keorg's answer.

Chooses the agent with the minimal best-case cost \$(10-\mathit{position}+\mathit{moves})^2-\mathit{moves}^2=(10-\mathit{position}+2\mathit{moves})(10-\mathit{position})\$.

Answer 10

C++: ~4953 (100 million iterations)

#include <bits/stdc++.h>

using namespace std;

mt19937 rng(chrono::steady_clock::now().time_since_epoch().count());
int rand(int a, int b) { return uniform_int_distribution<int>(a, b)(rng);}

int main() {
    const int NRIT = 100'000'000;
    auto go = [&] {
        int l[4] = { 10,10,10,10 }, t[4] = {}, sol = 0;
        auto cost = [&](int nr) {return l[nr] * (4 * t[nr] * l[nr] + 6); };
        while (1) {
            pair<int, int> best = { INT_MAX, -1 };
            for (int i = 0; i < 4; ++i)
                best = min(best, { cost(i), i });
            int opt = best.second;
            sol += 2 * t[opt]++ + 1;
            if (l[opt] == 10)
                --l[opt];
            else {
                if (rand(0,1))
                    --l[opt];
                else
                    ++l[opt];
            }
            if (!l[opt])
                return sol;
        }
    };
    long long total = 0;
    for (int i = 0; i < NRIT; ++i)
        total += go();
    cout << (double)(total) / NRIT;

}

I misread the question at first and thought it had a 1/2 chance to move backwards or stay still. This chooses the agent with the least expected cost to finish, according to my calculations and the (wrong) assumption it can't go backwards.

\$ \text{cost} = (10-p)(4 \cdot (10-p) \cdot m+6)\$

I want to add 10 000 iterations is definitely not enough, I could cherrypick 10k iterations to get 4790, which is far from the real average. Standardised testing also encounters the issue of different number generators, which can actually make a small difference (it did from my testing).

Answer 11

Python: ~4973

TL;DR By Gittins theory, there is an optimal solution of the following form, but it has some complicated replacement for approx_grades and this is a heuristic.

import random,numpy

def approx_grade(x,n):
  if x <= 0:
    return 2*n+1 + approx_grade(1, n+1)
  else:
    return 2**(10-x) * (9 - x) + (2**(11-x) -2)*(2*n + 1)

def sim(num_players):
  xs = [0]*num_players
  ns = [0]*num_players
  grades = [approx_grade(0,0)]*num_players
  while True:
    i = numpy.argmin(grades)  # STRATEGY HERE
    ns[i] += 1
    if xs[i] == 0:
      xs[i] = 1
    elif random.random() <= 0.5:
      xs[i] += 1
      if xs[i] == 10:
        return sum([n*n for n in ns])
    else:
      xs[i] -= 1
    grades[i] = approx_grade(xs[i], ns[i])

T = 100000
scores = [sim(4) for t in range(T)]
print(numpy.mean(scores))

This is a special case of the Gittins index theorem, and there is a very nice exposition of essentially this game in the 2003 paper On Playing Golf with Two Balls[1]

[1] https://epubs.siam.org/doi/abs/10.1137/S0895480102408341

Notice we can rephrase the cost as: when we move an agent that has taken n steps so far, we pay 2n+1 immediately, and the total cost at the end is the sum of all payments made along the way. Note 2n+1 = (n+1)^2 - n^2.

The key point is that for any given agent i and state of the agent (i.e. location x[i] and number of steps so far n[i]), there is an index or "grade" we can give the agent that perfectly captures how desirable it is to choose that agent next. The optimal strategy is to pick the agent with the best grade. One agent's grade does not depend on the others' states!

The grade is the following number. Imagine we had a fifth agent with the ability to magically teleport to the end for a cost g. And imagine we are playing the game with only agent i and the magic agent. If g is very small, we will just choose the magic agent. But if g is large, we will take a risk and play agent i, see what happens, and then decide whether to use the magic agent or play agent i again.

There is some magic number g(x[i],n[i]) where we are indifferent between the magic agent and i. That number is the grade.

---

Now, in theory the grade is polynomial-time computable, but it's pretty complex. So I used the following approximation: suppose the optimal strategy in the magic game is to keep trying agent i as long as they move right, but as soon as they move left, pay g to stop. Then the index for step 9 would satisfy

$$ g(9,n) = 2n + 1 + \frac{1}{2}g(9,n) . $$

The left side is if we pay g(9,n) for the magic agent to halt immediately, and the right side is if we pay 2n+1 to advance agent i: with probability half we terminate and pay nothing, and with probability half we go left, at which point we pay g(9,n) to stop.

This solves to g(9,n) = 4n+2.

Similarly, at location 8, we should be indifferent between paying g(8,n) to stop, or paying 2n+1 to get a half chance to go left -- when we pay g(8,n) to stop -- and a half chance to go right -- when we repeat the game at step 9.

$$ g(8,n) = 2n+1 + \frac{1}{2} g(8,n) + \frac{1}{2}\left(2(n+1)+1 + \frac{1}{2}g(8,n)\right) . $$

If we repeat this logic, we get an approximation to the grade which is

$$ g(x,n) = 2^{10-x}(9-x) + (2^{11-x}-2)(2n+1) . $$

After dealing with the special case of x=0, I get this code. Apparently it's still pretty suboptimal as another answer gets ~4927.

---

Note: I used multiprocessing to run 8 million total reps in 36min on a laptop, and got an average of just under 4973. I also reimplemented in C++ with probably hundreds of millions of reps and it seemed to converge to about 4974.5. But the quality of the RNG and precision arithmetic matters! C's built in rand() gave basically nonsense.

Answer 12

R / Rcpp: 4937 ± 4.22

Average of 15 million simulated games. Bounds are a 95% bootstrap confidence interval with 15,000 iterations.

As a supplement to the wonderful Gittins index answers, I wanted to show an empirical method that produces the same result. This solution predicts which robot to move using a neural net that's been trained by a genetic algorithm.

Strategy:

In order to optimally move our robots across the number lines, we need to determine a function that maps each robot's state information to a value that represents our decision to move it. This is called an action-value function. For simplicity, let's choose to always move the robot with the highest action-value.

At this point we could use the Gittins index as our action-value, but let's empirically estimate an action-value function instead.

Estimating the Action-Value Function

I decided to use a neural net (NN) to approximate the action-value function. The NN is fully connected with two hidden layers, each having ten neurons. Weights are initialized using the He method, and biases are initialized uniformly on \$[-1,1]\$. The activation functions are all rectified linear units (ReLU), except the activation function for the output, which is a ReLU clamped to \$[0, 1]\$.

Neural nets often have better performance when the inputs are symmetrically distributed and on the interval \$[-1, 1]\$, but in this problem the input state information of a robot is its current position on the interval \$[0 .. 9]\$ and the number of times it has been moved on the interval \$[0 .. \infty)\$.

I transformed the position of the robot to the \$[-1, 1]\$ interval by shifting and scaling.

To transform the total number of moves I first simulated 10 million games in which I always moved the same robot. This gave a distribution for the total number of moves, which I used to estimate an upper bound. The distribution of the total number of moves is heavily right skewed, so I applied a log transform to reduce skew, and then normalized the result using the log of my estimated upper bound in order to arrive at the \$[-1, 1]\$ interval for total number of moves. Note: Because the total number of moves can be zero, I added one to the total number of moves before taking the log.

Training the Neural Net

In reinforcement learning there are a number of well-established methods for training neural nets that predict action-values, but I have decided to completely ignore them because I've been having a lot of fun with genetic algorithms recently. There's no fun like gradient-free fun.

The genetic algorithm that I used updated NN weights and biases, but not the NN architecture itself. The initial population was 225 chromosomes, and fitness for each chromosome was calculated by evaluating its corresponding action-value function over 7500 games.

It is often desirable for neural net weights and biases to be relatively small and for the model's learned representation to be distributed across the neurons. To achieve this, I added batch normalization to the hidden layers of the neural nets. Regularizing the models in this way greatly improved performance and helped to ensure that the mutation function continued to be effective late in training by keeping the weights and biases close to the mutation function's generating distributions.

The batch normalization's \$\gamma\$ and \$\beta\$ parameters were also updated by the genetic algorithm.

Learned Action-Values:

This graph shows the final neural net's predicted action-values for robots in the most commonly traversed subset of the state space. Note: Because these estimated action-values are not with reference to any rewards, they only have meaning with respect to each other. I have arbitrarily scaled them to be between zero and one.

Code:

R Code to Generate Policy

# Parameters
iters <- 100
popNum <- 225
layers <- c(10, 10)
neval <- 7500
mutRate <- 0.001

# Evolve
population <- generatePopulation(popNum, layers, neval)
out <- rep(NA, iters)
for(i in 1:iters){
  
  # Evaluate population fitness
  fitness = evaluatePopulation(population, neval)
  out[i] = mean(fitness)

  # Keep track of best policy so far
  minFit = fitness[length(fitness)]
  if(min(fitness) <= minFit){
    minFit = min(fitness)
    minPolicy = population[[which.min(fitness)]]
  }
  
  # Selection
  cutoff = quantile(fitness, 0.5)
  population = population[fitness <= cutoff]
  
  # Crossover
  offspringNum = (popNum - length(population))
  offspring = lapply(1:offspringNum, function(i){
    inds = sample(1:length(population), 2, FALSE)
    crossover(population[[ inds[1] ]], population[[ inds[2] ]])
  })
  population = c(population, offspring)
  
  # Mutate
  population = lapply(1:length(population), function(i){
    mutate(population[[i]], rate = mutRate)
  })
  
  # Elitism
  eliteInds = c(round(0.99 * popNum):popNum)
  population[eliteInds] = lapply(eliteInds, function(i){minPolicy})
  
}

R Code Helper Functions

#' Evaluate Population
#' @export
evaluatePopulation <- function(population, neval){
  sapply(population, function(x){
    return(mean(evaluatePolicy(x, neval)))
  })
}

#' Generate Random Population
#' @export
generatePopulation <- function(size, layers, neval, maxCost = Inf){
  
  # Population
  population = list()
  
  # Increase yield
  while(TRUE){
    
    # Generate seed population
    out = lapply(1:size, function(i){
      generateChromosome(layers)
    })
    
    # Calculate fitness
    fitness = evaluatePopulation(out, neval)
    fitness[fitness > maxCost] = NA
    inds = order(fitness)
    out = out[inds]
    fitness = fitness[inds]
    
    # Seed population
    population = c(population, out[!is.na(fitness)])
    
    # Report
    if(length(population) >= size){
      population = population[1:size]
      break
    } else {
      message(paste0("Found: ", length(population)))
    }
    
  }
  
  # Return
  return(population)
  
}

Rcpp Code for Genetic Algorithm

#include <Rcpp.h>
using namespace Rcpp;

//* Initialize Environment
// [[Rcpp::export]]
NumericMatrix initEnv() {
  return(NumericMatrix (4, 2));
}

//* Update Environment
// [[Rcpp::export]]
NumericMatrix updateEnv(
    NumericMatrix env,
    int a
) {
  
  //Increase movement count
  env(a, 1) += 1;
  
  //Move robot
  if(env(a, 0) == 0){
    env(a, 0) = 1;
  } else {
    NumericVector r = Rcpp::runif(1, 0.0, 1.0);
    if(r[0] >= 0.5){
      env(a, 0) += 1;
    } else {
      env(a, 0) -= 1;
    }
  }
  
  //Return
  return(env);
  
}

//* Calculate Total Cost
// [[Rcpp::export]]
double envCost(NumericMatrix env) {
  double lambda = 0;
  for(int i = 0; i < env.nrow(); i++){
    lambda += pow(env(i, 1), 2);
  }
  return(lambda);
}

//* Generate Chromosome
// [[Rcpp::export]]
List generateChromosome(
  IntegerVector layers
) {
  
  //Create container for the NN matrices
  List chrom;
  
  //Create weight and bias matrices for each layer
  for(int i = 0; i < layers.length() + 1; i++){
    
    //Determine matrix dimensions
    int inputNum;
    int outputNum;
    if(i == 0){
      inputNum = 2;
      outputNum = layers[i];
    } else if(i == layers.length()){
      inputNum = layers[i - 1];
      outputNum = 1;
    } else {
      inputNum = layers[i - 1];
      outputNum = layers[i];
    }
    
    //Initialize weights and biases
    NumericVector wtList = Rcpp::rnorm(outputNum * inputNum, 0.0, sqrt(0.5 * inputNum));
    NumericVector biList = Rcpp::runif(outputNum, -1.0, 1.0);
    NumericMatrix wt (outputNum, inputNum, wtList.begin());
    NumericMatrix bi (outputNum, 1, biList.begin());
    
    //Record
    List chLayer = List::create(
      Named("wt") = wt,
      Named("bi") = bi
    );
    chrom.push_back(chLayer);
    
  }
  
  //Initialize batch normalization matrix
  NumericMatrix batch (layers.length(), 2);
  NumericVector batchVals = Rcpp::rnorm(batch.size(), 0.0, 1.0);
  for(int i = 0; i < batch.size(); i++){
    batch[i] = batchVals[i];
  }
  
  //Return
  List out = List::create(
    Named("layers") = layers,
    Named("chrom") = chrom,
    Named("batch") = batch
  );
  return(out);
  
}

//* Crossover
// [[Rcpp::export]]
List crossover(
    List x,
    List y,
    double rate = -1,
    double scale = -1
) {
  
  //Set parameters
  NumericVector probs = Rcpp::runif(2, 0.0, 1.0);
  if(rate >= 0){
    probs[0] = rate;
  }
  if(scale >= 0){
    probs[1] = scale;
  }
  
  //Crossover the NN
  List out = clone(x);
  
  //Access the components
  List chrom = out["chrom"];
  List xChrom = x["chrom"];
  List yChrom = y["chrom"];
  
  //Crossover weights and biases
  for(int i = 0; i < chrom.length(); i++){
    
    //Reference the layer
    List chLayer = chrom[i];
    List xLayer = xChrom[i];
    List yLayer = yChrom[i];
    
    //Access the matrices
    NumericMatrix wt = chLayer["wt"];
    NumericMatrix bi = chLayer["bi"];
    NumericMatrix xWt = xLayer["wt"];
    NumericMatrix xBi = xLayer["bi"];
    NumericMatrix yWt = yLayer["wt"];
    NumericMatrix yBi = yLayer["bi"];
    
    //Crossover
    LogicalVector doCross_wt = (Rcpp::runif(wt.size(), 0.0, 1.0) < probs[0]);
    LogicalVector doCross_bi = (Rcpp::runif(bi.size(), 0.0, 1.0) < probs[0]);
    for(int j = 0; j < wt.size(); j++){
      if(doCross_wt[j]){
        wt[j] = probs[1] * (xWt[j] - yWt[j]) + xWt[j];
      }
    }
    for(int j = 0; j < bi.size(); j++){
      if(doCross_bi[j]){
        bi[j] = probs[1] * (xBi[j] - yBi[j]) + xBi[j];
      }
    }
    
  }
  
  //Crossover gamma and beta
  NumericMatrix batch = out["batch"];
  NumericMatrix xBatch = x["batch"];
  NumericMatrix yBatch = y["batch"];
  LogicalVector doCross_yb = (Rcpp::runif(batch.size(), 0.0, 1.0) < probs[0]);
  for(int i = 0; i < batch.size(); i++){
    if(doCross_yb[i]){
      batch[i] = probs[1] * (xBatch[i] - yBatch[i]) + xBatch[i];
    }
  }
  
  //Return
  return(out);
  
}

//* Mutate
// [[Rcpp::export]]
List mutate(
    List x,
    double rate = -1
) {
  
  //Set parameters
  NumericVector probs = Rcpp::runif(1, 0.0, 1.0);
  if(rate >= 0){
    probs[0] = rate;
  }
  
  //Access components
  List out = clone(x);
  List chrom = out["chrom"];
  
  //Mutate weights and biases
  for(int i = 0; i < chrom.length(); i++){
    
    //Access the weights and biases
    List chLayer = chrom[i];
    NumericMatrix wt = chLayer["wt"];
    NumericMatrix bi = chLayer["bi"];
    
    //Mutate
    int numInputs = wt.ncol();
    LogicalVector doMut_wt = (Rcpp::runif(wt.size(), 0.0, 1.0) <= probs[0]);
    LogicalVector doMut_bi = (Rcpp::runif(bi.size(), 0.0, 1.0) <= probs[0]);
    for(int j = 0; j < wt.size(); j++){
      if(doMut_wt[j]){
        wt[j] = Rcpp::rnorm(1, 0.0, sqrt(0.5 * numInputs))[0];
      }
    }
    for(int j = 0; j < bi.size(); j++){
      if(doMut_bi[j]){
        bi[j] = Rcpp::runif(1, -1.0, 1.0)[0];
      }
    }
    
  }
  
  //Mutate gamma and beta
  NumericMatrix batch = out["batch"];
  LogicalVector doMut_yb = (Rcpp::runif(batch.size(), 0.0, 1.0) <= probs[0]);
  for(int i = 0; i < batch.size(); i++){
    if(doMut_yb[i]){
      batch[i] = Rcpp::rnorm(1, 0.0, 1.0)[0];
    }
  }
  
  //Return
  return(out);
  
}

//* Neural Net Prediction
//* 
//* Environment:
//* env[0] = Position
//* env[1] = Total moves
//* 
// [[Rcpp::export]]
NumericVector predictNN(
    NumericVector env,
    List chromosome
) {
  
  //Initialize the output
  NumericVector pred (1);
  
  //Process environment
  //Normalizing to be in [-1, 1]
  NumericVector vIn (2);
  vIn[0] = (env[0] / 5 - 1);
  vIn[1] = 2 * (log(env[1] + 1) / log(2000)) - 1;
  
  //Compute NN outputs
  List chrom = chromosome["chrom"];
  NumericMatrix batch = chromosome["batch"];
  for(int i = 0; i < chrom.length(); i++){
    
    //Access the layer
    List chLayer = chrom[i];
    NumericMatrix wt = chLayer["wt"];
    NumericMatrix bi = chLayer["bi"];
    
    //Calculate activation
    NumericVector vOut (wt.nrow());
    for(int j = 0; j < wt.nrow(); j++){
      NumericVector wtVec = wt(j, _);
      NumericVector biVec = bi(j, _);
      NumericVector act = {sum(wtVec * vIn) + biVec};
      vOut[j] = act[0];
    }
    
    //Apply batch normalization
    if(vOut.size() > 1){
      vOut = (vOut - mean(vOut)) / sd(vOut);
      vOut = batch(i, 0) * vOut + batch(i, 1);
    }
    
    //Apply ReLU
    for(int j = 0; j < vOut.length(); j++){
      NumericVector relu = {0.0, vOut[j]};
      vOut[j] = max(relu);
    }
    
    //Copy
    vIn = clone(vOut);
    
  }
  
  //Record
  pred[0] = vIn[0];
  
  //Return
  return(Rcpp::clamp(0.0, pred, 1.0));
  
}

//* NN Prediction to Action
// [[Rcpp::export]]
int predToResponse(NumericVector pred){
  
  //Calculate response
  double maxPred = max(pred);
  IntegerVector response;
  for(int i = 0; i < 4; i++){
    if(pred[i] >= maxPred){
      response.push_back(i);
    }
  }
  int action = Rcpp::sample(response, 1)[0];
  
  //Return
  return(action);
  
}

//* Evaluate Chromosome
// [[Rcpp::export]]
NumericVector evaluatePolicy(
    List chromosome,
    int iters
) {
  
  //Initialize parameters
  NumericVector out (iters);
  NumericVector pred (4);
  
  //Run episodes
  for(int i = 0; i < iters; i++){
    
    //Episode parameters
    NumericMatrix env = initEnv();
    
    //Complete episode
    bool done = false;
    while(!done){
      
      //Update state
      for(int j = 0; j < 4; j++){
        NumericVector state = env(j, _);
        pred[j] = predictNN(state, chromosome)[0];
      }
      int a = predToResponse(pred);
      env = updateEnv(env, a);
      
      //Check victory condition
      for(int j = 0; j < 4; j++){
        if(env(j, 0) >= 10){
          out[i] = envCost(env);
          done = true;
          break;
        }
      }
      
    }
    
  }
  
  //Return
  return(out);
  
}