As also described in Cormen, et al (2009) p. 65, in algorithm design, divide-and-conquer paradigm incorporates a recursive approach in which the main problem is:

Instead of constructing indefinite number of nested loops destroying the readability of the code and the performance of execution, the “recursive” way utilizes just one block of code which calls itself (hence the term “recursive”) for the smaller problem. The main point is to define a “stop” rule, so that the function does not sink into an infinite recursion depth. While nested loops modify the same object (or address space in the low level sense), recursion moves the “stack pointer”, so each recursion depth uses a different part of the stack (a copy of the objects will be created for each recursion). This illustrates a well-known trade-off in algorithm design: Memory versus performance; recursion enhances performance at the expense of using more memory.

Given the benefits of the recursive divide-and-conquer paradigm, however, the performance of the naive recursion may not still be as good as we may want, since there may be too many smaller problems to navigate through recursion. Sometimes, the same or a similar sub-problem may be visited multiple times and each sub-problem paves way to its sub-problems down to the stop rule – so a naive recursive approach solves the same sub-problems many and many times. The recipe to enhance performance is to “memoize” solved sub-problems: Simply to keep solutions in a look-up table, and first check the table if the value was memoized before commencing with the “divide-and-conquer”. Divide-and-conquer paradigm on steroids with memoization is known as “dynamic programming” or “dp” for short (Cormen, et al. 2009, p. 359).

Now we will illustrate an application of the “dynamic programming” approach using a question from the  PuzzlOR site – a site which publishes bi-monthly decision support puzzles for applied mathematicians. Since the answers to expired questions are disclosed by the website itself, there is no problem in describing a solution in detail. The question that appeared in October 2014 is called “Fighters”. The question involves a quadruel: A kind of duel fight that is among four participants instead of the classical two participant case.

The question goes as follows:

Four different fighters are having an all-out battle to determine who among them is the strongest.  The image above shows those four fighters: Allan, Barry, Charles, and Dan. 

 

Each fighter has varying attack and health abilities.  At the start of the battle, they have differing health points: Allan has 10, Barry has 12, Charles has 16, and Dan has 18.  Also, each fighter has differing attack points: Allan has 4, Barry has 3, Charles has 2, and Dan has 1.

 

The battle takes place over multiple rounds, each round consisting of a single attack.  In each round, one random attacker and one random defender are chosen.  When the attacker attacks a defender, the defender loses health points in the amount equivalent to the attacker’s attack points.  For example, if Allan is the attacker and Barry is the defender, Barry would lose 4 health points.

 

The fighters continue to randomly attack and defend in subsequent rounds until there is only one fighter left, who is then declared the winner.  A fighter is removed from the battle when his life points become zero (or less).

 

 

Question:  Which fighter is most likely to win the battle?

Starting with the initial health levels of 10, 12, 16, and 18, we can divide the main problem of finding the most likely winner into smaller ones, each of which is the ex-post health level configuration after a possible and equally likely attack between any “alive” fighters – “alive” meaning having a non-zero health level. We can divide the problem until we have only one alive fighter – which is the winner. We calculate the actual probability of winning here in this deepest recursion step at the stop rule (only one winner and the probability is zero for other fighters), so this is the “conquer” phase in which the smaller problem is solved.

Each recursion paves way into smaller problems – each one involving a different pair of fighters – number of which is the count of two-member permutations of alive fighters without replacement .  We calculate permutations instead of combinations since Fighter X attacking Fighter Y is different than the other way around – the attacker injures the attacked one,  so the attack is one-way and not symmetric. We calculate the permutations without replacement, since an attacker cannot attack herself/himself! So in the case of four alive fighters, there are 12 possible attacks. And each attack has a probability of 1/12. The winning probabilities of each fighter calculated from each of these 12 attacks are summed together in a bottom-up manner (from the deepest recursion level involving the last alive fighter, up to the initial question) and this step is the “combine” phase: In each step involving n number of ordered pairs of fighters, the 1/n probability is multiplied with the bottom-up probabilities to allow for carrying the joint probabilities up to the topmost level (the initial question).

We first start with the divide-and-conquer paradigm without memoization and see the performance:

## for non-repeated permutations
library(gtools) 

### non memoized version

# get the winning probabilities of each player
# win means only one player remains with non-zero health

puzzlor.fighters.nm <- function() { 
    
    health <- c(10, 12, 16, 18) # initial health points of players a:d
    attack <- c(4, 3, 2, 1) # attack points of players a:d

    # permutations of players. first column attackers, second column attacked.
    # self attack not allowed (no repetitions)
    at.combs <- gtools::permutations(4, 2, 1:4) 

    # get matrix of results through ply
    results <- t(apply(at.combs, 1, fight.rec.nm, health, attack, 1/12))

    # sum all probabilities of each player
    result <- colSums(results)

    # update player names
    names(result) <- c("Allan", "Barry", "Charles", "Dan")
 
    return(result)

}
    

# recursive function to recalculate the win probabilities after an attack,
# incorporating bottom-up cumulative joint probabilities

fight.rec.nm <- function(at, health, attack, prob) { 

    # at: vector of attacker and attacked no's
    # at1: attacker no
    # at2: attacked no
    # health: vector for current health levels of all players
    # attack: vector for attack points of all players
    # prob: carried joint probability

    at1 <- at[1] # attacker no
    at2 <- at[2] # attacked no

    # update the health of at2 after attack
    health[at2] <- max(health[at2] - attack[at1], 0)

    players <- which(health > 0) # index of alive players

    # if only one player survives, so is the winner
    if (length(players) == 1) { 
        result <- rep(0,4) # vector for results
        result[players] <- prob # update the probability of the winner
        return(result)

    } else {
 
        # permutations of indices of attackers and attackeds.
        # self attack not allowed (no repetitions in permutations)
        at.combs <- gtools::permutations(length(players), 2, players)
 
        probs <- nrow(at.combs) # total count of attack permutations
        prob.new <- 1 / probs # probability of each attack permutation

        # collect matrix of results recursing the function
        results <- t(apply(at.combs, 1, fight.rec.nm, health, attack, prob.new)) 
        result <- colSums(results) # sum the probabilities for each player
        result <- prob * result # multiply with the joint probability
        return(result)

    }

}

The key points here in the above code are:

We construct two functions: The first one (“puzzlor.fighters.nm”) (line 9) is a wrapper to initiate necessary objects and to start the initial problem calling the recursive function. The second function (“fight.rec.nm”) (line 35) is the recursive function which does all the steps described in detail above.

The performance of this non-memoized version is … well it is still running as at the time I am writing these sentences.

The reason for the performance bottleneck is the fact that too many sub-problems are generated until only one winner remains and the number of “tree branches” grows geometrically with each recursion (multiplied by the number of possible ordered attack pairs in each step). And those sub-problems are revisited many times since the fights can yield the same health level configuration at different instances . We can enhance the performance without redefining the problem but by augmenting it with an additional data structure:

## for non-repeated permutations
library(gtools) 

# get the winning probabilities of each player
# win means only one player remains with non-zero health

puzzlor.fighters.5D <- function() {     

    health <- c(10, 12, 16, 18) # initial health points of players a:d
    attack <- c(4, 3, 2, 1) # attack points of players a:d

    ## global 5D object for memoization
    result.ar <<- array(dim = c(health + 1, 4), dimnames = NULL)

    # permutations of players. first column attackers, second column attacked.
    # self attack not allowed (no repetitions).
    at.combs <- gtools::permutations(4, 2, 1:4)
 
    # get matrix of results through ply
    results <- t(apply(at.combs, 1, fight.rec.5D, health, attack, 1/12)) 
    result <- colSums(results) # sum all probabilities of each player
    names(result) <- c("Allan", "Barry", "Charles", "Dan") # update player names
    return(result)

}
    

# recursive function to recalculate the win probabilities after an attack,
# incorporating bottom-up cumulative joint probabilities
fight.rec.5D <- function(at, health, attack, prob) { 

    # at: vector of attacker and attacked no's
    # at1: attacker no
    # at2: attacked no
    # health: vector for current health levels of all players
    # attack: vector for attack points of all players
    # prob: carried joint probability

    at1 <- at[1] # attacker no
    at2 <- at[2] # attacked no

    # update the health of at2 after attack
    health[at2] <- max(health[at2] - attack[at1], 0)

    # create a matrix for indexing the 5D array 
    mat.index <- cbind(matrix(rep(health + 1, 4), nrow = 4, byrow = T), 1:4)

    mat.row <- "["(result.ar, mat.index) # get the memoized health configuration

    # if the probabilities for the health conf. are already memoized
    if (!is.na(mat.row[1])) {

        # get the result and multiply with the prob
        result <- prob * mat.row 
        return(result)

    } else {

        players <- which(health > 0) # index of alive players

        # if only one player survives, so is the winner
        if (length(players) == 1) { # 

            result <- rep(0,4) # vector for results
            result[players] <- prob # update the probability of the winner
            return(result)

        } else {

            # permutations of indices of attackers and attackeds.
            # self attack not allowed (no repetitions in permutations)
            at.combs <- gtools::permutations(length(players), 2, players) 
            probs <- nrow(at.combs) # total count of attack permutations
            prob.new <- 1 / probs # probability of each attack permutation

            # collect matrix of results recursing the function
            results <- t(apply(at.combs, 1, fight.rec.5D, health, attack, prob.new))
            result <- colSums(results) # sum the probabilities for each player

            # memoize the health configuration
            "["(result.ar, mat.index) <<- result 

            # multiply with the joint probability
            result <- prob * result 
            return(result)

        }

    }

}

The main points of the modified code are as follows:

> microbenchmark(puzzlor.fighters.5D(), times = 1)
Unit: seconds
                  expr      min       lq     mean   median       uq      max
 puzzlor.fighters.5D() 7.983507 7.983507 7.983507 7.983507 7.983507 7.983507
 neval
     1

The code executes in less than 8 seconds. You can check the correct solutions from the PuzzlOR web page.

A downside to the dp approach is that, it cannot be utilized with the “mcmapply” and similar R functions of the “parallel” package for multicore parallelism, since these functions create copies of the global objects as well as the local objects for each thread. So all values cannot be memoized into the same object. However the gains from the memoization outweigh the disadvantage of using a single core.

You can view the whole code (with an additional option with a 4D array and a matrix, which slightly underperforms the 5D case) from my Github page.

Notes: