快速找到全零答案的方法

Question

对于每个长度为 n+h-1 且值从 0 到 1 的数组，我想检查是否存在另一个长度为 n 且值从 -1,0,1 的非零数组，以便所有的 h内积为零。我的天真方法是

import numpy as np
import itertools
(n,h)= 4,3
for longtuple in itertools.product([0,1], repeat = n+h-1):
    bad = 0
    for v in itertools.product([-1,0,1], repeat = n):
        if not any(v):
            continue
        if (not np.correlate(v, longtuple, 'valid').any()):
            bad = 1
            break
    if (bad == 0):
        print "Good"
        print longtuple

如果我们设置 n = 19 和 h = 10 这是我想要测试的，这会非常慢。

我的目标是找到一个长度为n+h-1 的“好”数组。有没有 如何加快速度，使n = 19 和h = 10 可行？

当前的幼稚方法需要 2^(n+h-1)3^(n) 次迭代，每次迭代大约需要 n 次。那是 n = 19 和 h = 10 的 311,992,186,885,373,952 次迭代，这是不可能的。

注意 1 将 convolve 更改为 correlate 以便代码正确地考虑 v。

---

2015 年 7 月 10 日

问题仍然存在，没有足够快的解决方案来解决 n=19 和 h=10 的问题。

Answer 1

考虑以下“在中间相遇”的方法。

首先，将 leekaiinthesky 提供的矩阵公式中的情况重铸。

接下来，请注意，如果我们将问题更改为寻找h x n Hankel matrix H 的 0 和 1 使得 Hs1 永远不会等于 Hs2 对于 0 和 1 的两个不同的短向量。那是因为Hs1 = Hs2 暗示H(s1-s2)=0 暗示有一个向量v 由1、0 和-1 组成，即s1-s2，这样Hv = 0；相反，如果Hv = 0 在{-1,0,1}^n 中为v，那么我们可以在{0,1}^n 中找到s1 和s2 使得v = s1 - s2 和Hs1 = Hs2。

当n=19 时，{0,1}^n 中只有 524,288 个向量 s 可以尝试；散列结果Hs，如果相同的结果出现两次，则H 不好，然后尝试另一个H。就内存而言，这种方法是相当可行的。有2^(n+h-1)汉克尔矩阵H可以试试；当 n=19 和 h=10 是 268,435,456 个矩阵时。那是 2^38 测试，或 274,877,906,944，每个测试都有大约 nh 操作来将矩阵 H 和向量 s 相乘，大约 52 万亿次操作。这似乎可行，不是吗？

由于您现在只处理 0 和 1，而不是 -1，您还可以通过使用位操作（移位和计数 1）来加快处理速度。

更新

我用 C++ 实现了我的想法。我正在使用位运算来计算点积，将结果向量编码为长整数，并使用 unordered_set 检测重复项，当发现点积的重复向量时，提前退出给定的长向量。

几分钟后，对于 n=17 和 h=10，我获得了 00000000010010111000100100，对于 n=18 和 h=10，我在几分钟后获得了 000000111011110001001101011。我正要运行它 n=19 和 h=10。

#include <iostream>
#include <bitset>
#include <unordered_set>

/* Count the number of 1 bits in 32 bit int x in 21 instructions.
 * From /Hackers Delight/ by Henry S. Warren, Jr., 5-2
 */
int count1Bits(int x) {
  x = x - ((x >> 1) & 0x55555555);
  x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
  x = (x + (x >> 4)) & 0x0F0F0F0F;
  x = x + (x >> 8);
  x = x + (x >> 16);
  return x & 0x0000003F;
}

int main () {
  const int n = 19;
  const int h = 10;
  std::unordered_set<long> dotProductSet;

  // look at all 2^(n+h-1) possibilities for longVec
  // upper n bits cannot be all 0 so we can start with 1 in pos h
  for (int longVec = (1 << (h-1)); longVec < (1 << (n+h-1)); ++longVec) {

    dotProductSet.clear();
    bool good = true;

    // now look at all n digit non-zero shortVecs
    for (int shortVec = 1; shortVec < (1 << n); ++shortVec) {

      // longVec dot products with shifted shortVecs generates h results
      // each between 0 and n inclusive, can encode as h digit number in
      // base n+1, up to (n+1)^h = 20^10 approx 13 digits, need long
      long dotProduct = 0;

      // encode h dot products of shifted shortVec with longVec
      // as base n+1 integer
      for(int startShort = 0; startShort < h; ++startShort) {
        int shortVecShifted = shortVec << startShort;
        dotProduct *= n+1;
        dotProduct += count1Bits(longVec & shortVecShifted);
      }

      auto ret = dotProductSet.insert(dotProduct);
      if (!ret.second) {
        good = false;
        break;
      }
    }

    if (good) {
      std::cout << std::bitset<(n+h-1)>(longVec) << std::endl;
      break;
    }
  }

  return 0;
}

第二次更新

n=19 和 h=10 的程序在我的笔记本电脑的后台运行了两周。最后，它只是退出而没有打印任何结果。除非程序中出现某种错误，否则看起来没有具有您想要的属性的长向量。我建议寻找为什么没有这么长的向量的理论原因。也许某种计数论点会起作用。

Answer 2

这只是部分答案，因为这似乎仍然太慢而无法检查n=19, h=10 的情况（在这种情况下甚至可能不存在“好的”向量）。

这是 @vib 描述的检查算法的实现，并在 Mathematica 中对所有 2^(n+h-1) 向量使用随机抽样。

TestFullNullSpace[A_, ns_] := Module[{dim, n},
  {dim, n} = Dimensions[ns];
  For[i = 1, i < 3^dim, i++,
   testvec = Mod[IntegerDigits[i, 3, dim].ns, 3] /. {2 -> -1};
   If[Norm[A.testvec] == 0,
    Return[False]];
   ];
  Return[True];
]

n = 17;
h = 10;

Do[
 v = Table[RandomChoice[{0, 1}], {n + h - 1}];
 A = Table[v[[i ;; i + n - 1]], {i, 1, h}];
 ns = NullSpace[A, Modulus -> 3] /. {2 -> -1};
 If[TestFullNullSpace[A, ns],
  Print[v]];,
 {1000}]

经过几秒钟的计算，上述运行的示例输出：

{0,0,1,1,0,0,0,0,0,0,1,0,1,1,1,0,1,0,1,1,0,0,0,1,1,0}
{1,1,0,1,0,0,0,1,1,0,1,1,1,1,1,0,1,0,1,0,0,1,1,0,0,0}
{1,0,1,1,1,1,1,0,0,0,1,1,0,1,0,1,1,0,0,1,1,0,1,1,1,0}
{0,0,0,0,1,0,1,1,1,0,1,1,0,0,1,1,1,1,0,1,0,1,0,0,1,1}

所以从 1000 个检查向量中，有 4 个是“好”的（除非我有错误）。不幸的是，对于n=18，我运行了几分钟，但仍然没有找到“好的”向量。我不知道它们是不存在还是非常罕见。

Answer 3

可能有更快的方法*

您要查找的内容与kernel or null space of a matrix 的概念有关。

特别是，对于每个 n+h-1 “longtuple”并给定 n，构造一个 h×n 矩阵，其行是 longtuple 的 n 个子元组。换句话说，如果你的 longtuple 是 [0,0,0,1,0,0] 并且 n = 3，那么你的矩阵是：

[[0 0 0]
 [0 0 1]
 [0 1 0]
 [1 0 0]]

将此矩阵称为A。您正在寻找一个矢量 x 使得 Ax = 0，其中 0 是一个全为 0 的矢量.如果存在这样的x（即不是本身全为0）并且可以缩放到只包含{-1,0,1}，那么你想抛出A 并继续下一个 longtuple。

我不太确定计算内核的（理论上最有效的）计算复杂度是多少，但它似乎在 O(h+n)^3 左右，无论如何都是比 O(3^n) 好很多。有关如何计算内核的一些示例，请参阅上面的 Wikipedia 链接或 Python (NumPy, SciPy), finding the null space of a matrix。

无论如何，一旦你确定了内核，你就必须做一些额外的工作来确定是否有任何形式为 {-1, 0, 1}^n 的向量存在于其中，但我不认为那是计算负担很大。

*NB：在 cmets 中，@vib 指出这实际上可能是一个很大的计算负担。我不确定确定这些向量是否与内核相交的最佳算法是什么。也许它不能在多项式时间内解决，在这种情况下，这个答案并不能加速原始问题！

示例代码

根据您在 cmets 中给出的示例，从上面链接的其他 Stack Overflow 问题改编代码：

import scipy
from scipy import linalg, matrix
def null(A, eps=1e-15):
    u, s, vh = scipy.linalg.svd(A)
    null_mask = (s <= eps)
    null_space = scipy.compress(null_mask, vh, axis=0)
    return scipy.transpose(null_space)

A = matrix([[0,0,0,1],[0,0,1,0],[0,1,0,0],[0,0,0,0]])
print null(A)

#> [[-1.]
#>  [ 0.]
#>  [ 0.]
#>  [ 0.]]

代码给出了一个 n 元组的示例（实际上，与您给出的示例相同），它使 [0, 0, 0, 1, 0, 0] 无效，因为它是一个“好”的长元组。如果代码返回[]，那么大概没有这样的n 元组，并且longtuple 是“好”的。 （如果代码确实返回了一些东西，你仍然需要检查 {-1, 0, 1} 部分。）

进一步的想法

这样的x是否存在，暂时忽略{-1, 0, 1}约束，相当于nullityA的/strong>（内核维度）大于0。这相当于询问的rank是否A 等于 n。因此，如果您找到了一些巧妙处理 {-1, 0, 1} 约束的方法并将其分解为只需要计算 A 的等级，我相信这可以更快地完成.

顺便说一句，你（或给你这个问题的人）似乎很可能已经知道这一切......否则你为什么要称 longtuple 的长度为“n+h-1”，如果你还没有从高度 h 的矩阵开始...！

Answer 4

这是一种将其减少到O(n*h*3^(n/2 + 1)) 的方法。这可扩展性很差，但对于您的用例来说已经足够了。

遍历向量前半部分的所有可能性。创建一个字典字典 ... 数组字典，其键是每个移位内积的值，其最终值是产生该序列的向量的前半部分的数组。

现在遍历向量后半部分的所有可能性。当你计算它的每个内积时，遍历嵌套字典，看看是否有相应的前半部分对内积的贡献仍然取消。如果你一直遍历到最后，那么你可以把你找到的前半部分和你也找到的后半部分放在一起，你就有答案了。

不要忘记忽略全为 0 的答案！

Answer 5

下面是一个将复杂度从 3^n 降低到 3^{n-h} 的算法。

让 v_1、v_2、..、v_h 成为您需要与之正交的向量。

考虑向量空间 (Z/3Z)^n。令 v'_1, .., v'_h 为 v_1, .., v_h 在这个空间中的自然包含。

现在令 w 为系数在 {-1,0,1} 中的向量，令 w' 为 (Z/3Z)^n 的向量，通过自然地将 w 视为 (Z/3Z) 的向量而获得^n。那么 w 与 v_1, .., v_h (in R) 的标量积为零的必要条件是 w' 与 v'_1, .., v 的标量积为零 (in (Z/3Z)^n ) '_h.

现在您可以很容易地确定与 v'_1、..、v'_h 的标量积为零的 w'。它们将形成一个大小为 3^{n-h} 的空间。然后，您需要检查它们中的每一个，关联的 w 是否实际上与所有 v_i 正交。