When following a top-down approach we memoize answers and have explicit recursion. Usually this is done with a mutable look-up table. For instance the knapsack problem can be solved this way:

import scala.collection.mutable

def knapsack(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    val solutions = mutable.HashMap.empty[(Int,Int), Int]
    def f(i: Int, j: Int): Int = {
        if(!solutions.contains((i,j))) {
            val solution = 
                if( i == 0 || j == 0 ) {
                    0
                } else if( j - weight(i-1) >= 0 ) {
                    Math.max( 
                        f(i-1,j) , 
                        f(i-1, j - weight(i-1)) + value(i-1) 
                    )
                } else {
                    f(i-1, j)
                }
            solutions.put((i,j), solution)
        }
        solutions((i,j))
    }
    f(n, maxWeight)
}

What if we could do the same thing with a data structure that’s defined recursively, just as the solutions? Can this be done with a normal map? Well there’s the method tabulate that can be used to build a Vector using a function that describes every entry by it’s indexes. Let’s try something with it:

def knapsack2(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    lazy val solutions: Vector[Vector[Int]] = 
        Vector.tabulate(n + 1, maxWeight + 1) { (i,j) =>
            if( i == 0 || j == 0 ) {
                0
            } else if( j - weight(i-1) >= 0 ) {
                Math.max( 
                    solutions(i-1)(j) , 
                    solutions(i-1)(j - weight(i-1)) + value(i-1) 
                )
            } else {
                solutions(i-1)(j)
            }
        }
    solutions(n)(maxWeight)
}

Well, this even compiles but fails at runtime with a Stackoverflow error. This is because when we use the solutions matrix inside the tabulate function we demand for the solutions value to be already completely built. It’s kind of a chicken and egg problem. If we could relax this condition and make it possible to have some entries defined and others not…

Let’s define a LazyVector allowing precisely that!

// not threadsclass LazyVector[A](thunks: Vector[() => A]) { // (1)
    
    private val values: Array[Option[A]] = Array.fill(thunks.length)              (None) // (2)
    
    def apply(i: Int): A = { // (3)
        if(!values(i).isDefined) {
            values(i) = Some( thunks(i)() )
        }
        values(i).get
    }   
}

At (1) we declare that to build a LazyVector we need a Vector describing how to compute each entry. This is just a Vector of thunks (a thunk is a zero argument function). At (2) we see where we will save the computed values. Under the hood we will memoize the values using a mutable array. Yes, mutable! Don’t worry, our intentions are pure (get it?). At (3) we declare how we are going to retrieve values for this data structure: if it’s already computed just return it, otherwise compute it and save it.

Let’s also add some constructor functions in the companion object:

object LazyVector {
    
    def tabulate[A](n: Int)(f: Int => A): LazyVector[A] = {
        new LazyVector( Vector.tabulate(n)( i => () => f(i) ) )
    }

    def tabulate[A](m: Int, n: Int)(f: (Int,Int) => A): LazyVector[LazyVector[A]] = {
        tabulate(m)( i => tabulate(n)( j => f(i,j) ) )
    }

}

Let’s try the same thing we did with Vectors:

def knapsack(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    lazy val solutions: LazyVector[LazyVector[Int]] = 
        LazyVector.tabulate(n + 1, maxWeight + 1) { (i,j) =>
            if( i == 0 || j == 0 ) {
                0
            } else if( j - weight(i-1) >= 0 ) {
                Math.max( 
                    solutions(i-1)(j) , 
                    solutions(i-1)(j - weight(i-1)) + value(i-1) 
                )
            } else {
                solutions(i-1)(j)
            }
        }
    solutions(n)(maxWeight)
}

And this works! No stack overflow!

This has the same time complexity as we had before with the bottom-up approach. The space complexity, unfortunately, is O(nW) which is not the best we can do if we are only interested in computing the maximum total value.

There is an important advantage with this approach, though. The code is way more declarative. It doesn’t say how to iterate over the matrix of solutions in such a way that dependencies are pre-computed. It just says: here’s what the solution looks like, let the data structure’s laziness avoid re-computations.

Returning the concrete solution

If we also want to know which items conform the maximum total value then we have to record all the decisions we make, so the space complexity will be O(nW) anyway. This just implies returning more information when we are looking for the maximum:

def max[A](a: (Int, A), b: (Int, A)): (Int, A) =
    if(a._1 > b._1)
        a
    else
        b

def knapsackWithItems(
    maxWeight: Int, value: Vector[Int], weight: Vector[Int]): List[Int] = {
    val n = value.length
    lazy val solutions: LazyVector[LazyVector[(Int, List[Int])]] = 
        LazyVector.tabulate(n + 1, maxWeight + 1) { (i,j) =>
            if( i == 0 || j == 0 ) {
                (0, List.empty)
            } else if( j - weight(i-1) >= 0 ) {
                val (including, itemsIds) = solutions(i-1)(j - weight(i-1))
                max( 
                    solutions(i-1)(j) ,
                    (including + value(i-1), (i-1) :: itemsIds )
                )
            } else {
                solutions(i-1)(j)
            }
        }
    val (_,items) = solutions(n)(maxWeight)
    items
}

Similarly we can solve Project Euler’s 81-th problem like this:

sealed trait Instruction
case object GoRight extends Instruction
case object GoDown extends Instruction

def minPath(values: Vector[Vector[Int]]): List[Instruction] = {
    val m = values.length
    val n = values(0).length
    lazy val solutions: LazyVector[LazyVector[(Int, List[Instruction])]] = 
        LazyVector.tabulate(m, n) { 
            case (0,0) => (values(0)(0), List.empty)
            case (0,j) => 
                val (leftVal, leftInsts) = solutions(0)(j-1)
                (values(0)(j) + leftVal, leftInsts ++ List(GoRight))
            case (i,0) => 
                val (upVal, upInsts) = solutions(i-1)(0)
                (values(i)(0) + upVal, upInsts ++ List(GoDown))
            case (i,j) =>
                val (leftVal, leftInsts) = solutions(i)(j-1)
                val (upVal  , upInsts  ) = solutions(i-1)(j)
                if(leftVal < upVal) {
                    (leftVal + values(i)(j), leftInsts ++ List(GoRight))
                } else {
                    (upVal + values(i)(j)  , upInsts ++ List(GoDown))                    
                }
        }
    val (_,insts) = solutions(m-1)(n-1)
    insts
}

Conclusions

This top-down approach has the disadvantage that it may consume more memory than we need (if we are just interested in the maximum or minimum value). But for the kind of problems where we need to have all the decisions then it’s worth it anyway. It produces very clear and declarative code. The only disadvantage I can think of is that it may produce stack-overflows for large inputs.

— Miguel Vilá

For instance in this project Euler problem the solutions depend on each other differently. Once again, let’s start with a problem definition that allows us to break the problem into subproblems:

Let’s consider some position (i,j). To get to that position we only have two possibilities: we either get to that position by coming from the upper row and going down from there or we come from the previous column and go right from there. These two possibilities are mutually exclusive thus we can compute the minimum between them to compute the optimal solution for (i,j) :

Here v(i,j) is the value stored at row i and column j.

And, once again, let’s draw a diagram showing the dependencies:

There is an important difference with respect to the base cases here. The minimal path for each value in the first row is simply the path that starts at (0,0) and goes right from there (this is the formula R(i,j), which is the accumulative row sum). Similarly, the minimal path for each value in the first column is the part that starts at and goes down from there (this is the formula C(i,j) , which is the accumulative column sum).

Now, unlike the 0/1 knapsack problem, the solution for (i,j) not only depends on a value at the previous row. It also depends on the previous value in the same row.

Let’s start with the base cases: how to compute the accumulative sum of an array of numbers? fold can be used to combine all the values using a reducer function. But the intermediate results of that function aren’t kept, which is what we want: we don’t want to have the total sum. We want to get the accumulative sum. scan does this. It has the same parameters as fold but instead of returning just one value it’s going to return a vector of all the intermediate values:

def accumulativeSum(nums: Vector[Int]): Vector[Int] =
    nums.scanLeft(0)(_ + _).drop(1)

The drop(1) is so we don’t include the entry with 0 at the beginning.

So the first row and the first column can be computed like this:

def minPath(values: Vector[Vector[Int]]): Int = {
    val m = values.length
    val n = values(0).length
    val firstRow = accumulativeSum(values(0))
    val firstColumn = accumulativeSum( (0 to m-1).map(values(_)(0)) )
    ???
}

Let’s use the same strategy as before. We’ll compute every new row from the bottom up, starting with the first row as the initial state:

def minPath(values: Vector[Vector[Int]]): Int = {
    val m = values.length
    val n = values(0).length
    val firstRow = accumulativeSum(values(0))
    val firstColumn = accumulativeSum( (0 to m-1).map(values(_)(0)) )
    (1 to m-1).foldLeft(firstRow) { (upperRow, i) =>
        ??? // How to compute the new row?
    }.last
}

Now to compute say M(i,1) we’ll need to use two things: M(i — 1, 1)

which is in the upper row and M(i, 0) which is a value in the very first column, at the left of M(i,1). Then to compute M(i,2) we will need to use M(i — 1,2) and M(i,1) which is the value that we previously computed. The function fold can help us with using the last value to compute the next. But we also need to return all the list of intermediate values, not just the last one, so scan is the one we want. Putting it all together:

def minPath(values: Vector[Vector[Int]]): Int = {
    val m = values.length
    val n = values(0).length
    val firstRow = accumulativeSum(values(0))
    val firstColumn = accumulativeSum( (0 to m-1).map(values(_)(0)) )
    (1 to m-1).foldLeft(firstRow) { (upperRow, i) =>
        val leftmostSolution = firstColumn(i)
        (1 to n-1).scanLeft(leftmostSolution) { (leftSolution,j) =>
            val upperSolution = upperRow(j)
            Math.min(leftSolution, upperSolution) + values(i)(j)
        }.toVector
    }.last
}

That’s it! No mutation whatsoever!

Conclusions

Even for a class of algorithms that looks inherently imperative there are ways to avoid mutation and to end up with the very same algorithm using immutable data structures and higher-order functions like fold or scan. You could just translate the mutations into State monad values but that’s just too complex and unnecessary. There are simpler ways to functionally express the same thing.

Is this worth it? Maybe not if you are used to the imperative approach, which already yields a purely functional solution up to observable mutation.

In the next post we will see how to do the same with a top-down approach that’s both functional and expressive.

— Miguel Vilá

Translating every mutation to a state modification

In functional programming one way to handle state is the State monad. Then we could mindlessly replace the modifications by creating State values that represent those modifications. Using the State monad implementation in Cats and an immutable data structure like Vector we can reimplement the O(nW) space version like this:

First some imports:

import cats._, cats.instances.all._, cats.syntax.traverse._, cats.syntax.foldable._
import cats.data.State

This function will prove useful when updating the state:

def setSolution(i: Int, j: Int, newVal: Int)
               (solutions: Vector[Vector[Int]]): Vector[Vector[Int]] =
    solutions.updated(i, solutions(i).updated(j, newVal))

and the whole thing:

def knapsack(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    val initialState: Vector[Vector[Int]] = Vector.fill(n+1, maxWeight + 1)( 0 )
    val st: State[Vector[Vector[Int]], Unit] = ( 1 to n ).toList.traverseU_ { i =>
        ( 1 to maxWeight ).toList.traverseU_ { j =>
            for {
                solutions <- State.get[Vector[Vector[Int]]]
                newVal = if( j - weight(i-1) >= 0 ) {
                    Math.max(
                        solutions(i-1)(j) , 
                        solutions(i-1)(j - weight(i-1)) + value(i-1) 
                    )
                } else {
                    solutions(i-1)(j)
                }
                _ <- State.modify(setSolution(i,j,newVal))
            } yield ()
        }
    }
    val solution = st.runS(initialState).value
    solution(n)(maxWeight)
}

Not too interesting. This even looks too much like the original code. In a similar way we can do the same for the O(W)space version:

def knapsack(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    val initialState: Vector[Int] = Vector.fill(maxWeight + 1)( 0 )
    val st: State[Vector[Int], Unit] = ( 1 to n ).toList.traverseU_ { i =>
        for {
            solutions <- State.get[Vector[Int]]
            newSolutions = 0 +: ( 1 to maxWeight ).map { j: Int =>
                if( j - weight(i-1) >= 0 ) {
                    Math.max(
                        solutions(j) , 
                        solutions(j - weight(i-1)) + value(i-1) 
                    )
                } else {
                    solutions(j)
                }
            }.toVector
            _ <- State.set(newSolutions)
        } yield ()
    }
    val solution = st.runS(initialState).value
    solution(maxWeight)
}

Have we learned anything from this? Not much, really. And this is boring precisely because the state monad can be used to implement the very same imperative logic. This solves the problem in a functional way but it’s kind of an over-engineered approach. Is there a simpler, more elegant and functional way to do this?

Just use fold!

In the imperative version at each iteration we are computing a new row based on the last one. And we are interested in the last row because it holds the final answer. The first row is full of zeros, that’s our “base case”. And a perfect function for building incremental state from an initial state is `fold`! `fold` has two parameters: the initial state and a function that computes a new state based on the previous one and on the current element of the structure we are folding over. This fits us perfectly because to compute each new row we must have two ingredients: the previous row and which row number is that we are processing. We could start by folding over a list of the row numbers starting with a row full of zeros:

val firstRow = Vector.fill(maxWeight + 1)( 0 )
val lastRow = (1 to n).foldLeft(firstRow) { (upperRow, i) =>
    ??? // How to compute each new row?
}
val answer = lastRow.last

To compute the new row we must do the same that we were doing before. Well, kind of: instead of modifying the input we are going to return a new row. To do this we can compute each value based on the column number, the row number and the values of the previous row:

val firstRow = Vector.fill(maxWeight + 1)( 0 )
val lastRow = (1 to n).foldLeft(firstRow) { (upperRow, i) =>
    0 +: (1 to maxWeight).map { j => 
        if( j - weight(i-1) >= 0 ) {
            Math.max( 
                upperRow(j) ,
                upperRow(j - weight(i-1)) + value(i-1)
            )
        } else {
            upperRow(j)
        }
    }.toVector
}
val answer = lastRow.last

And that’s pretty much it! There’s no mutation and we accomplish the same thing with very similar code! This even has O(W) space complexity.

What we have seen so far

In this post we saw two functional solutions to a dynamic programming problem. The first approach is very equivalent to the imperative solution and the second one uses a higher order function. In the next post we will see what happens when the dependencies between the solutions have another shape.

— Miguel Vilá

Overlapping subproblems are subproblems that may have the same dependencies. DP helps us here because we want to compute these common dependencies just once. There are two usual approaches:

• The top-down approach starts trying to solve the bigger problem and uses memoization.
• The bottom-up approach starts with the subproblems to build up the solution to the bigger problem and uses some data structure to save the subproblems solutions.

If you want to learn more about dynamic programming you can look into some of the links at the end of this post.

In my opinion the best way to understand Dynamic Programming is by looking at multiple examples. The first usual example is the Fibonacci sequence but I think that’s a bad instance of the general idea because nobody would call Dynamic Programming the “technique” of storing the two last solutions in two variables. But in more complex examples you will be able to see the “pattern”.

This post is not so much an introduction to DP. As you will see it’s more about how to implement DP algorithms in a functional way when taking the bottom-up approach.

The imperative perspective

Dynamic programming problems, when solved in a bottom-up fashion, usually rely in some mutable data structure which holds the solutions for smaller subproblems. For example let’s explore the 0/1 knapsack problem. This and this are good videos about it and describe a solution. Here’s my summary of the problem:

Given n items indexed 1, 2, ⋯, n with values Vi and weights Wi find a subset of them with maximum total value that doesn’t exceed the input weight W. Each item may be included just once.

For example, for this list of items with their respective value:

and maximum target weight W=9.

You can select the items 3, 4 and 5 for a total weight of 9 and a total value of 15. This is good but you could do even better by selecting items 1, 2, 4 and 5. In this case we would get a total value of 16 which is in fact the maximum you can get without exceeding the maximum weight of 9.

To solve this problem in the general case we can first break it by considering subsets. Let’s start with a definition that allows us to express subsets of the problem:

For each item we have two choices: either we include it and gain it’s value but we have reduced our available weight. Or, on the other hand, we ignore it and never use it. These two possibilities totally exclude each other so if we compute the maximum between them for every item (and for every possible maximum weight ≤ W) then we can get the overall maximum. And these two possibilities correspond to two subproblem definitions.

The recursive formula that express is this:

According to this formula we can see that the problem (i,j) depends on two subproblems: (i — 1, j - Wi) and (i — 1, j)

Here’s a diagram that describes the subproblem dependencies:

The blue square represents a subproblem we are trying to solve. It depends on two other subproblems (the ones pointed by the arrows). And the final answer is the red square at the lower-right corner (the value

). It represents the solution when we consider all the elements and with the problem’s input weight. We can also see the “base” cases here: when considering no item the maximum value we can produce is and the same happens when the maximum weight is 0 (this assumes all Wi> 0).

Drawing this kind of diagrams can be useful because it helps us to determine the shape of the iterations. In this case, an element of a row depends on the values of the previous row. Thus, we can build up solutions going down one row at a time.

Let’s use Scala to solve this. A traditional imperative solution can be:

def knapsack(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    val solutions: Array[Array[Int]] = Array.fill(n+1, maxWeight + 1)( 0 )
    (1 to n) foreach { i =>
        (1 to maxWeight) foreach { j =>
            solutions(i)(j) = if( j - weight(i-1) >= 0 ) {
                Math.max( 
                    solutions(i-1)(j) , 
                    solutions(i-1)(j - weight(i-1)) + value(i-1) 
                )
            } else {
                solutions(i-1)(j)
            }
        }
    } 
    solutions(n)(maxWeight)
}

This iterates one row at a time and computes the solution given the recursive formula. The time complexity of this is O(nW) where n is the number of items to consider and W is the maximum weight. And the space complexity is also O(nW). The previous solution can be improved given that we only need to remember the last row of solutions for each iteration:

def knapsack(maxWeight: Int, value: Vector[Int], weight: Vector[Int]): Int = {
    val n = value.length
    var solutions: Array[Int] = Array.fill(maxWeight + 1)( 0 )
    (1 to n) foreach { i =>
        val newSolutions = Array.fill(maxWeight + 1)( 0 )
        (1 to maxWeight) foreach { j =>
            newSolutions(j) = if( j - weight(i-1) >= 0 ) {
                Math.max( 
                    solutions(j) , 
                    solutions(j - weight(i-1)) + value(i-1) 
                )
            } else {
                solutions(j)
            }
        }
        solutions = newSolutions
    } 
    solutions(maxWeight)
}

The space complexity of this will be O(W).

How to do the same thing functionally?

How could we do the same thing functionally? One of the main objectives in functional programming is avoiding side-effects. The most prominent side-effect that we have here is the mutation of a data structure, and this can be avoided with immutable data structures.

But before doing that, we should ask ourselves: is it worth it? Should someone calling the function knapsack care if it is implemented using mutable variables? My opinion is that no: it doesn’t matter if knapsack is implemented with or without mutability, even from a functional programming point of view. Let me explain.

This is because the implementation of knapsack that I listed before is already “purely functional” with respect to observable mutations! Will it return the same results when called with the same arguments whether it is implemented with mutable or immutable data structures? Yes! Is it implemented with any side-effect? Well, it has a mutable variable but does that mutable variable escape the function definition? No. Is it a global variable whose mutation can be observed from the outside? No. Thus, with respect to “observable mutations”, the function knapsack is referentially transparent.

Having said that, the motivation behind implementing knapsack functionally is just an exercise on how we can model this problem in a functional fashion. The idea is just to stretch our “functional” muscles: to see how far the paradigm can take us in an apparently mutable domain. But my opinion is that, even when following the functional programming ideas, mutation is not bad as long as it’s just an implementation detail and not something that can be observed by the bigger system.

What we have seen so far

In this post we saw how to solve a common dynamic programming problem using an imperative approach. This involved mutation and a for loop. In the next post we will see an initial approach on how to do the same thing functionally.

Additional resources

If you want to know more about dynamic programming here are some good resources:

• The Wikipedia page on DP is actually very comprehensive.
• The second part of Stanford’s Algorithm’s course has a section on DP.
• This chapter of Dasgupta, Papadimitriou & Vazirani’s “Algorithms”.

— Miguel Vilá