Is there a way to recursively iterate through all possible sub-matrices of a matrix while preventing some sub-matrices from being visited?

Question

My task is to find all sub-matrices inside a matrix such that each sub-matrix counted satisfies a certain condition and also is not a part of another sub-matrix that works.

My first thought was to write a recursive procedure so that we could simply return from the current sub-matrix whenever we find that it works (to prevent any sub-matrices of that sub-matrix from being tested). Here is my code that attempts to do this:

void find(int xmin, int xmax, int ymin, int ymax){
    if(xmin > xmax || ymin > ymax){return;}
    else if(works(xmin, xmax, ymin, ymax)){++ANS; return;}

    find(xmin + 1, xmax, ymin, ymax);
    find(xmin, xmax - 1, ymin, ymax);
    find(xmin, xmax, ymin + 1, ymax);
    find(xmin, xmax, ymin, ymax - 1);
}

The problem with my current method seems to be the fact that it allows sub-matrices to be visited more than once, meaning that the return statements are ineffective and don't actually prevent sub-matrices of working sub-matrices from being counted, because they are visited from other matrices. I think I have the right idea with writing a recursive procedure, though. Can someone please point me in the right direction?

Answer 1

Obviously, you need a way to check if a submatrix was evaluated before or is contained in a larger solution. Also, you need to account that after a solution is found, a larger solution may be found later which covers the currently found solution.

One way of doing this is to utilize a structure called R*-tree, which allows to query spatial (or geometric) data efficiently. For this purpose, you could use R-tree implementation from boost. By using a box (rectangle) type to represent a submatrix, you can then use the R-tree with queries boost::geometry::index::contains (to find previously found solutions which include the submatrix considered) and boost::geometry::index::within (to find previously found solutions which are contained within a newly found solution).

Here is a working example in C++11, which is based on your idea:

#include <vector>
#include <numeric>
#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/register/point.hpp>
#include <boost/geometry/geometries/register/box.hpp>
#include <boost/geometry/index/rtree.hpp>

namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;

struct Element
{
    int x, y;
};

struct Box
{
    Element lt, rb;
};

BOOST_GEOMETRY_REGISTER_POINT_2D(Element, long, cs::cartesian, x, y)
BOOST_GEOMETRY_REGISTER_BOX(Box, Element, lt, rb)

template<class M>
bool works(M&& matrix, Box box) {
    // Dummy function to check if sum of a submatrix is divisible by 7
    long s = 0;
    for (int y=box.lt.y; y<box.rb.y; y++)
        for (int x=box.lt.x; x<box.rb.x; x++)
            s += matrix[y][x];
    return s % 7 == 0;
}

template <class T, class M>
void find(T& tree, M&& matrix, Box box, T& result) {
    if (box.lt.x >= box.rb.x || box.lt.y >= box.rb.y) return;
    for (auto it = tree.qbegin(bgi::contains(box)); it != tree.qend(); ++it) return;
    if (works(matrix, box)) {
        // Found a new solution

        // Remove any working submatrices which are within the new solution
        std::vector<Box> temp;
        for (auto it = result.qbegin(bgi::within(box)); it != result.qend(); ++it)
            temp.push_back(*it); 
        result.remove(std::begin(temp), std::end(temp));

        // Remember the new solution
        result.insert(box);
        tree.insert(box);
        return;
    }

    // Recursion
    find(tree, matrix, Box{{box.lt.x+1,box.lt.y},{box.rb.x,box.rb.y}}, result);
    find(tree, matrix, Box{{box.lt.x,box.lt.y+1},{box.rb.x,box.rb.y}}, result);
    find(tree, matrix, Box{{box.lt.x,box.lt.y},{box.rb.x-1,box.rb.y}}, result);
    find(tree, matrix, Box{{box.lt.x,box.lt.y},{box.rb.x,box.rb.y-1}}, result);

    tree.insert(box);
}

template <class T>
void show(const T& vec) {
    for (auto box : vec) {
        std::cout << "Start at (" << box.lt.x << ", " << box.lt.y << "), width="
                  << box.rb.x-box.lt.x << ", height=" << box.rb.y-box.lt.y << "\n";
    }
}

int main()
{
    // Initialize R-tree
    const size_t rtree_max_size = 20000;
    bgi::rtree<Box, bgi::rstar<rtree_max_size> > tree, result;

    // Initialize sample matrix
    const int width = 4;
    const int height = 3;
    int matrix[height][width];
    std::iota((int*)matrix, (int*)matrix + height*width, 1);

    // Calculate result
    find(tree, matrix, Box{{0,0},{width,height}}, result);

    // Output
    std::cout << "Resulting submatrices:\n";
    show(result);

    return 0;
}

In this example the following matrix is considered:

1  2  3  4
5  6  7  8
9 10 11 12

And the program will find all submatrices for which the sum of their elements is divisible by 7. Output:

Resulting submatrices:
Start at (0, 2), width=4, height=1
Start at (1, 0), width=3, height=3
Start at (0, 0), width=2, height=2
Start at (0, 1), width=1, height=2

What I liked about your recursive approach is that it works very fast even for large matrices of 1000×1000 elements.