Spatial DP: Finding the Largest Square

In the past two weeks we’ve explored a couple different problems in dynamic programming. These were simpler 1-dimensional problems. But dynamic programming is often at its most powerful when you can work across multiple dimensions. In today’s problem, we’ll consider a problem that is actually a 2D spatial problem where we can use dynamic programming.

If you want to learn how to write dynamic programming solutions in Haskell from the ground up, take a look at our Solve.hs course. DP is one of several algorithmic approaches you’ll learn in Module 3!

The Problem

Today’s problem (Maximal Square) is fairly simple conceptually. We are given a grid of 1’s and 0’s like so:

10100
11111
00111
10101

We must return the size of the largest square in the grid composed entirely of 1’s. So in the example above, the answer would be 4. There are two 2x2 squares we can form, starting in the 2nd row, using either the 3rd or 4th column as the “top left” of the square.

We can do a couple small edits to change the answer here. For example, we can flip the second ‘0’ in the bottom row and we’ll get a 3x3 grid, allowing the answer 9:

10100
11111
00111
10111

We could instead flip the second ‘1’ in the third row, and now the answer is only 1, as there are no 2x2 squares remaining:

10100
11111
00101
10101

The Algorithm

To solve this, we can imagine a DP grid that has “layers” where each layer has the same dimensions as the original grid. Each layer has a number “k” associated with it. The index {row,column} at layer k tells us whether or not a square of size k exists in the original grid with size k x k, with the cell {row, column} as its top left cell.

To construct this grid, we would need a base case and a recursive case. The base case is to consider layer 1. This is identical to the original grid we receive. Any location with 1 in the original grid is the top left for a 1x1 square.

So how do we build the layer k+1? This requires one simple insight. Suppose we are dealing with a single index {r,c}. In order for this to be the top left of a square of size k+1, we just need to check that 4 cells begin squares of size k: {r,c}, {r+1,c}, {r,c+1},{r+1,c+1}.

So to form the next layer, we just loop through each index in the layer and fill it in with 1 if it meets that criterion. Once we reach a layer where each entry is 0, we are done. We should return the square of the last layer we found.

There are a few optimizations possible here. Thinking back to our first DP problem, we didn’t need to store the full DP array since each new step only depended on a couple prior values. This time, we don’t need a full grid with “k” layers. We could alternate with only two grids, saving new values from the prior grid, and then making our “new” grid the “old” grid for the next layer.

But even simpler than that, we can keep modifying a single grid in place. Each “new” value we calculate depends on numbers below and/or to its right. So as long as we loop through the grid from left to right and top to bottom, we are safe modifying its values in place. At least, that’s what we’ll do in Rust. In Haskell we could do this with the mutable array API, but we’ll stick with the more conventional, immutable, approach in this article. (You can learn more about Haskell’s mutable arrays in Solve.hs).

Rust Solution

Let’s start with the Rust solution, demonstrating the mutable array approach. We’ll start by defining a series of terms, like the dimensions of our input and our dp grid (which is initially a clone of the input). We’ll also define a boolean (found) to indicate if we’ve found at least a single 1 on the current layer. We’ll track level, the number of layers confirmed to have a 1.

pub fn maximal_square(matrix: Vec<Vec<char>>) -> i32 {
    let m = matrix.len();
    let n = matrix[0].len();
    let mut level = 0;
    let mut dp = matrix.clone();
    let mut found = true;
    ...
    return level * level;
}

Of course, our final answer is just the square of the final “level” we determine. But how do we find this? We’ll need an outer while loop that terminates once we hit a level that does not hold a 1. We reset found as false to start each loop, but at the end of the loop, we’ll increment the level if we have found something.

pub fn maximal_square(matrix: Vec<Vec<char>>) -> i32 {
    let m = matrix.len();
    let n = matrix[0].len();
    let mut level = 0;
    let mut dp = matrix.clone();
    let mut found = true;
    while (found) {
        found = false;
        ...
        if (found) {
            level += 1;
        }
    }
    return level * level;
}

Now the core of the “layer” loop is to loop through each cell, left to right and top to bottom.

pub fn maximal_square(matrix: Vec<Vec<char>>) -> i32 {
    ...
    while (found) {
        found = false;
        for i in 0..m {
            for j in 0..n {
                ...
            }
        }
        if (found) {
            level += 1;
        }
    }
    return level * level;
}

So what happens inside the loop? When we hit a 0 cell, we don’t need to do anything. It always remains a 0 and we haven’t “found” anything. But interesting things happen if we hit a 1.

First, we note that found is now true - this layer is not empty. We have found a k x k square. But second, we should now reset this cell as 0 if it does not make a square of size k+1. We need to first check the dimensions to make sure we don’t go out of bounds, but then also check the 3 spaces, to the right, below, and diagonally away from us. If any of these are 0, we reset this cell as 0.

pub fn maximal_square(matrix: Vec<Vec<char>>) -> i32 {
    ...
    while (found) {
        found = false;
        for i in 0..m {
            for j in 0..n {
                if (dp[i][j] == '1') {
                    found = true;
                    if (i + 1 >= m || 
                        j + 1 >= n || 
                        dp[i][j+1] == '0' ||
                        dp[i+1][j] == '0' ||
                        dp[i+1][j+1] == '0') {
                        dp[i][j] = '0';
                    }
                }
            }
        }
        if (found) {
            level += 1;
        }
    }
    return level * level;
}

And just by filling in this logic, our function is suddenly done! Our inner loop is complete, and our outer loop will break once we find no more increasingly large squares. Here is the full Rust solution:

pub fn maximal_square(matrix: Vec<Vec<char>>) -> i32 {
    let m = matrix.len();
    let n = matrix[0].len();
    let mut level = 0;
    let mut dp = matrix.clone();
    let mut found = true;
    while (found) {
        found = false;
        for i in 0..m {
            for j in 0..n {
                if (dp[i][j] == '1') {
                    found = true;
                    if (i + 1 >= m || 
                        j + 1 >= n || 
                        dp[i][j+1] == '0' ||
                        dp[i+1][j] == '0' ||
                        dp[i+1][j+1] == '0') {
                        dp[i][j] = '0';
                    }
                }
            }
        }
        if (found) {
            level += 1;
        }
    }
    return level * level;
}

Haskell Solution

Now let’s write this in Haskell. We’ll start with a few definitions, including a type alias for our DP map. We’ll take an Array as the problem input, but we want a HashMap for our stateful version since we can “mutate” a HashMap efficiently:

type SquareMap = HM.HashMap (Int, Int) Bool

maximalSquare :: A.Array (Int, Int) Bool -> Int
maximalSquare grid = ...
  where
    ((minRow,minCol), (maxRow, maxCol)) = A.bounds grid
    initialMap = HM.fromList [vs | vs <- A.assocs grid]

    ...

Now we’ll define two loop functions - one for the inner loop, one for the outer loop. The “state” for the inner loop is our current level number, as well as the map of the previous layer. The inner loop (coordLoop) should return us an updated map, as well as the found bool value telling us if we’ve found at least a single 1 in the prior layer.

maximalSquare :: A.Array (Int, Int) Bool -> Int
maximalSquare grid = ...
  where
    ((minRow,minCol), (maxRow, maxCol)) = A.bounds grid
    initialMap = HM.fromList [vs | vs <- A.assocs grid]

    coordLoop :: (Bool, SquareMap) -> (Int, Int) -> (Bool, SquareMap)
    coordLoop (found, mp) coord@(r, c) = ...

    loop :: Int -> HM.HashMap (Int, Int) Bool -> Int
    loop level mp = ...

    ...

Notice that coordLoop has the argument pattern for foldl, rather than foldr. We want to loop through our coordinates in the proper order, from left to right and top down. If we use a right fold over the indices of the grid, it will go in reverse order.

Let’s start by filling in the inner loop. The first thing to do is determine if the found value needs to change. This is the case if we discover a True value at this index:

maximalSquare :: A.Array (Int, Int) Bool -> Int
maximalSquare grid = ...
  where
    ((minRow,minCol), (maxRow, maxCol)) = A.bounds grid
    initialMap = HM.fromList [vs | vs <- A.assocs grid]

    coordLoop :: (Bool, SquareMap) -> (Int, Int) -> (Bool, SquareMap)
    coordLoop (found, mp) coord@(r, c) =
      let found' = found || mp HM.! coord
          ...
      in  (found’, ...)

Now we need the 5 conditions that tell us if this cell should get cleared. Calculate all these, and insert False at the cell if any of them match. Otherwise, keep the map as is!

maximalSquare :: A.Array (Int, Int) Bool -> Int
maximalSquare grid = ...
  where
    ((minRow,minCol), (maxRow, maxCol)) = A.bounds grid
    initialMap = HM.fromList [vs | vs <- A.assocs grid]

    coordLoop :: (Bool, SquareMap) -> (Int, Int) -> (Bool, SquareMap)
    coordLoop  (found, mp) coord@(r, c) =
      let found' = found || mp HM.! coord
          tooRight = c >= maxCol
          tooLow = r >= maxRow
          toRight = mp HM.! (r, c + 1)
          under = mp HM.! (r + 1, c)
          diag = mp HM.! (r + 1, c + 1)
          failNext = tooLow || tooRight || not toRight || not under || not diag
          mp' = if failNext then HM.insert coord False mp else mp
      in  (found', mp')

    ...

Now for the outer loop, we use foldl to go through our coordinates using the coordLoop. If we’ve found at least 1 square at this size, then we recurse with the new map and an incremented size. Otherwise we return the square of the current level. Then we just need to call this loop with initial values:

```haskell
type SquareMap = HM.HashMap (Int, Int) Bool

maximalSquare :: A.Array (Int, Int) Bool -> Int
maximalSquare grid = loop 0 initialMap
  where
    ((minRow,minCol), (maxRow, maxCol)) = A.bounds grid
    initialMap = HM.fromList [vs | vs <- A.assocs grid]

    coordLoop :: (Bool, SquareMap) -> (Int, Int) -> (Bool, SquareMap)
    coordLoop  (found, mp) coord@(r, c) = ...

    loop :: Int -> HM.HashMap (Int, Int) Bool -> Int
    loop level mp =
      let (found, mp') = foldl coordLoop (False, mp) (A.indices grid)
      in  if found then loop (level + 1) mp' else (level * level)

This completes our Haskell solution!

type SquareMap = HM.HashMap (Int, Int) Bool

maximalSquare :: A.Array (Int, Int) Bool -> Int
maximalSquare grid = loop 0 initialMap
  where
    ((minRow,minCol), (maxRow, maxCol)) = A.bounds grid
    initialMap = HM.fromList [vs | vs <- A.assocs grid]

    coordLoop :: (Bool, SquareMap) -> (Int, Int) -> (Bool, SquareMap)
    coordLoop  (found, mp) coord@(r, c) =
      let found' = found || mp HM.! coord
          tooRight = c >= maxCol
          tooLow = r >= maxRow
          toRight = mp HM.! (r, c + 1)
          under = mp HM.! (r + 1, c)
          diag = mp HM.! (r + 1, c + 1)
          failNext = tooLow || tooRight || not toRight || not under || not diag
          mp' = if failNext then HM.insert coord False mp else mp
      in  (found', mp')

    loop :: Int -> HM.HashMap (Int, Int) Bool -> Int
    loop level mp =
      let (found, mp') = foldl coordLoop (False, mp) (A.indices grid)
      in  if found then loop (level + 1) mp' else (level * level)

Conclusion

Next week we’ll look at one more multi-dimensional DP problem where the dimensions aren’t quite as obvious in this spatial way. The best way to understand DP is to learn related concepts from scratch, including your basic use-it-or-lose-it problems and memoization. You’ll study all these concepts and learn Haskell implementation tricks in Module 3 of Solve.hs. Enroll in the course now!

Next
Next

Making Change: Array-Based DP