找出2d矩阵是否是另一个2d矩阵的子集

Question

最近我参加了一个黑客马拉松,我开始了解一个试图在2d矩阵中找到网格形式的问题.一个模式可能是U,H和T,并将用3*3矩阵表示如果我想呈现H和U.

+--+--+--+            +--+--+--+
|1 |0 |1 |            |1 |0 |1 |
+--+--+--+            +--+--+--+
|1 |1 |1 |  --> H     |1 |0 |1 |    -> U
+--+--+--+            +--+--+--+
|1 |0 |1 |            |1 |1 |1 |
+--+--+--+            +--+--+--+

现在我需要搜索到10*10 matrix containing 0s and 1s.Closest和唯一的解决方案我可以得到它的强力算法O(n ^ 4).在MATLAB和R这样的语言中有非常微妙的方法来做到这一点但不是在C,C++中.我尝试了很多在谷歌和SO上搜索这个解决方案.但我最接近的是这个SO POST,它讨论了实现Rabin-Karp字符串搜索算法.但是没有伪代码或任何帖子解释这个.可以任何人帮忙或者提供任何链接,pdf或一些逻辑来简化这个？

编辑

作为Eugene Sh.评论说如果N是大矩阵(NxN)的大小而k是小矩阵(kxk),则buteforce算法应该取O((Nk)^ 2).由于k是固定的,它减少到O(N ^ 2).是绝对正确的.但如果N和K很大,有没有任何一般化的方法？

Answer 1

好吧,这是2D Rabin-Karp方法.

对于以下讨论,假设我们想要在(n,n)矩阵内找到(m,m)子矩阵.(这个概念也适用于矩形矩阵,但我的指数用完了.)

我们的想法是,对于每个可能的子矩阵,我们计算一个哈希值.只有当哈希匹配我们想要找到的矩阵的哈希时,我们才会按元素进行比较.

为了提高效率,我们必须避免每次重新计算子矩阵的整个哈希值.因为今晚我睡不着觉,我唯一可以弄清楚如何做到这一点的哈希函数就是相应子矩阵中1的总和.我将它作为练习留给比我更聪明的人来找出更好的滚动哈希函数.

现在,如果我们刚刚检查了从(i,j)到(i + m - 1,j + m - 1)的子矩阵,并且知道它内部有x 1,我们可以计算子中的1的数量.矩阵一到右 - 即从(i,j + 1)到(i + m - 1,j + m) - 通过从(i,j)到(i)减去子矢量中的1的数量+ m - 1,j)并在子矢量中添加1的数量(i,j + m)到(i + m - 1,j + m).

如果我们点击大矩阵的右边距,我们将窗口向下移动一个然后再向左移动边距然后再向下移动一个然后再向右移动,依此类推.

注意,这需要O(m)操作,而不是每个候选者的O(m ²).如果我们为每对索引执行此操作,我们得到O(m n ²)的工作.因此,通过巧妙地将潜在子矩阵的大小的窗口移动通过大矩阵,我们可以将工作量减少m倍.也就是说,如果我们没有得到太多的哈希冲突.

这是一张图片:

当我们将当前窗口向右移动时,我们在左侧的红色列向量中减去1的数量,并在右侧的绿色列向量中添加1的数量以获得新的1的数量.窗口.

我使用伟大的Eigen C++模板库实现了这个想法的快速演示.该示例还使用了Boost中的一些内容,但仅用于参数解析和输出格式化,因此如果您没有Boost但想要尝试代码,则可以轻松摆脱它.索引摆弄有点乱,但我会在此处不做进一步解释.上述散文应充分涵盖.

#include <cassert>
#include <cstddef>
#include <cstdlib>
#include <iostream>
#include <random>
#include <type_traits>
#include <utility>

#include <boost/format.hpp>
#include <boost/lexical_cast.hpp>

#include <Eigen/Dense>

#define PROGRAM "submatrix"
#define SEED_CSTDLIB_RAND 1

using BitMatrix = Eigen::Matrix<bool, Eigen::Dynamic, Eigen::Dynamic>;
using Index1D = BitMatrix::Index;
using Index2D = std::pair<Index1D, Index1D>;

std::ostream&
operator<<(std::ostream& out, const Index2D& idx)
{
  out << "(" << idx.first << ", " << idx.second << ")";
  return out;
}

BitMatrix
get_random_bit_matrix(const Index1D rows, const Index1D cols)
{
  auto matrix = BitMatrix {rows, cols};
  matrix.setRandom();
  return matrix;
}

Index2D
findSubMatrix(const BitMatrix& haystack,
              const BitMatrix& needle,
              Index1D *const collisions_ptr = nullptr) noexcept
{
  static_assert(std::is_signed<Index1D>::value, "unsigned index type");
  const auto end = Index2D {haystack.rows(), haystack.cols()};
  const auto hr = haystack.rows();
  const auto hc = haystack.cols();
  const auto nr = needle.rows();
  const auto nc = needle.cols();
  if (nr > hr || nr > hc)
    return end;
  const auto target = needle.count();
  auto current = haystack.block(0, 0, nr - 1, nc).count();
  auto j = Index1D {0};
  for (auto i = Index1D {0}; i <= hr - nr; ++i)
    {
      if (j == 0)  // at left margin
        current += haystack.block(i + nr - 1, 0, 1, nc).count();
      else if (j == hc - nc)  // at right margin
        current += haystack.block(i + nr - 1, hc - nc, 1, nc).count();
      else
        assert(!"this should never happen");
      while (true)
        {
          if (i % 2 == 0)  // moving right
            {
              if (j > 0)
                current += haystack.block(i, j + nc - 1, nr, 1).count();
            }
          else  // moving left
            {
              if (j < hc - nc)
                current += haystack.block(i, j, nr, 1).count();
            }
          assert(haystack.block(i, j, nr, nc).count() == current);
          if (current == target)
            {
              // TODO: There must be a better way than using cwiseEqual().
              if (haystack.block(i, j, nr, nc).cwiseEqual(needle).all())
                return Index2D {i, j};
              else if (collisions_ptr)
                *collisions_ptr += 1;
            }
          if (i % 2 == 0)  // moving right
            {
              if (j < hc - nc)
                {
                  current -= haystack.block(i, j, nr, 1).count();
                  ++j;
                }
              else break;
            }
          else  // moving left
            {
              if (j > 0)
                {
                  current -= haystack.block(i, j + nc - 1, nr, 1).count();
                  --j;
                }
              else break;
            }
        }
      if (i % 2 == 0)  // at right margin
        current -= haystack.block(i, hc - nc, 1, nc).count();
      else  // at left margin
        current -= haystack.block(i, 0, 1, nc).count();
    }
  return end;
}

int
main(int argc, char * * argv)
{
  if (SEED_CSTDLIB_RAND)
    {
      std::random_device rnddev {};
      srand(rnddev());
    }
  if (argc != 5)
    {
      std::cerr << "usage: " << PROGRAM
                << " ROWS_HAYSTACK COLUMNS_HAYSTACK"
                << " ROWS_NEEDLE COLUMNS_NEEDLE"
                << std::endl;
      return EXIT_FAILURE;
    }
  auto hr = boost::lexical_cast<Index1D>(argv[1]);
  auto hc = boost::lexical_cast<Index1D>(argv[2]);
  auto nr = boost::lexical_cast<Index1D>(argv[3]);
  auto nc = boost::lexical_cast<Index1D>(argv[4]);
  const auto haystack = get_random_bit_matrix(hr, hc);
  const auto needle = get_random_bit_matrix(nr, nc);
  auto collisions = Index1D {};
  const auto idx = findSubMatrix(haystack, needle, &collisions);
  const auto end = Index2D {haystack.rows(), haystack.cols()};
  std::cout << "This is the haystack:\n\n" << haystack << "\n\n";
  std::cout << "This is the needle:\n\n" << needle << "\n\n";
  if (idx != end)
    std::cout << "Found as sub-matrix at " << idx << ".\n";
  else
    std::cout << "Not found as sub-matrix.\n";
  std::cout << boost::format("There were %d (%.2f %%) hash collisions.\n")
    % collisions
    % (100.0 * collisions / ((hr - nr) * (hc - nc)));
  return (idx != end) ? EXIT_SUCCESS : EXIT_FAILURE;
}

在编译和运行时,请将上述内容视为伪代码.我几乎没有尝试过优化它.这只是我自己的概念验证.