将 Matlab 中矩阵的每一行与其余行进行比较

Question

我在 Matlab 中有一个矩阵 A，由零和维度为 BxM 的矩阵组成。

具体来说，A 包含M 空格中所有可能的 1 和 0 配置，同时考虑到顺序（因此，B=2^M）。

• 矩阵的构建使得对于任何i=1,...,N、A(i,:)>A(j,:) 或A(i,:)><A(j,:) 对于j=i+1,...,N。
• 写A(i,:)>=A(j,:) 意味着A(i,h)>=A(j,h) 代表h=1,...,M。
• 写成A(i,:)>A(j,:) 表示A(i,h)>=A(j,h) 对应h=1,...,M，至少有一个h 严格不等式。
• 写入A(i,:)><A(j,:) 意味着无法确定是A(i,:)>=A(j,:) 还是相反。

例如当M=3

  A=[1 1 1;
     1 1 0;
     1 0 1;
     1 0 0;
     0 1 1;
     0 1 0;
     0 0 1;
     0 0 0];

考虑

B=cell(2^M, 2^M);

对于A、A(i,:)>A(j,:) 的任意两个可比较行，我希望B{i,j} 包含所有行A(k,:)，这样A(j,:)<A(k,:)<A(i,:)。

在上面的示例中，所需的输出将是

B{1,4}=[1 1 0; 1 0 1];
B{1,6}=[1 1 0; 0 1 1];
B{1,7}=[1 0 1; 0 1 1];
B{1,8}=[1 1 0; 1 0 1; 1 0 0; 0 1 1; 0 1 0; 0 0 1];
B{2,8}=[1 0 0; 0 1 0];
B{3,8}=[1 0 0; 0 0 1];
B{5,8}=[0 0 1; 0 1 0];

这段代码做我想做的事

B=cell(2^M, 2^M); 
  for j=1:size(A,1)
      for h=1:size(A,1)
          if  sum(A(j,:)==A(h,:))~=M && sum(A(j,:)>=A(h,:))==M %if A(j,:)>A(h,:) according to the meaning indicated above
              B{j,h}=A(any(bsxfun(@ne, A,A(j,:)),2) & any(bsxfun(@ne, A,A(h,:)),2) &... 
              all((bsxfun(@le, A,A(j,:)) &  bsxfun(@ge, A,A(h,:))),2),:);
          end
      end
  end

但是，上面的代码是不可行的，因为在我的实际情况下M=20。您对是否可以加快速度有什么建议，如果可以，如何加快速度？

代码的一个主要问题是，对于M=20，它需要预先分配一个单元格(2^20)x(2^20)，这显然是做不到的。

另一方面，在双循环结束时，很多单元格是空的（因为许多行对不满足 if 条件），我真正需要的是跟踪 内容和仅非空坐标的坐标。因此，也许，一个“稀疏单元”可以帮助我。

任何建议都将不胜感激。

Answer 1

而不是预先计算和记忆 A 矩阵（“旧 Matlab 风格”），

编写一个“生成器函数”，为您提供任何行索引的 A 矩阵的行内容。通过这种方式，您可以将内存密集型问题转换为 CPU 密集型问题。除此之外，子结果彼此独立。

所以你需要的是

A_by_row=@(rowIdx)(....)

这样就不需要大量的 RAM，您可以将问题分发到 parfor、GPU 甚至多个节点，然后组合子范围的子结果。

试试这个：

ListIdx=0; 
for j=1:B 
 for h=1:B 
  if sum(A_by_row(j)==A_by_row(h))~=M, 
    ListIdx=ListIdx+1, 
    B_List(ListIdx).Coordinates=[j,h],     
    B_List(ListIdx).Result=**YourCodeThatMakesAnArbitraryLengthVec‌​tor**, 
  end, 
 end, 
end 

这样，在结构“列表”中，您将获得每个条目的坐标和答案向量。

祝你好运。

Answer 2

我们可以做更多的数学运算来预测测试的结果，然后保存预测，而不是暴力破解答案。

我们会做出以下假设：

1) 在数组 A 中，索引位置“n”存储“n”的二进制表示。例如：A(3) = [0,1,1]

2) 因为 A(j) 不能等于 A(h)，A(j) 的二进制表示也不能总是小于 A(h)，所以结果必须是下三角的。

3) 在绘制 M=5 的结果位置后，我们看到存在由合规条件创建的分形模式。可以使用卢卡斯对应定理“AND(NOT(j),h)=binomial_coefficient(j,h) mod(2)”重新创建此分形

4) 以下函数复制您的结果。

function out = user3285148()
    M = 5;
    maxm = 2^M-1;
    for j = 0:maxm
        for h = 0:maxm;
            if predict(j,h)
                B{j+1,h+1} = [binaryarray(j,M);binaryarray(h,M)];
            end
        end
    end
    out  = B;
end

function out = binaryarray(in, length)
    out = zeros(1,length);
    bit = 0;
    while in>0;
        out(end-bit)=mod(in,2);
        in = floor(in/2);
        bit =bit+1;
    end
end

function out = predict(j,h)
    out =0;
    if h<j
        out = mod(nchoosek(j,h),2);
    end
end

建议您以紧凑的形式存储 B，或使用预测在存储 (2^20)^2 个信息单元的 leu 中按需提取数据。

如果这符合您的要求，请告诉我。

Answer 3

我会把你的问题分成 3 部分。

1. 找到 B 所需的尺寸
2. 在 A 中找到所有可比较的行对，其中至少有一行（B 的左侧部分）
3. 在每个可比对之间查找 A 的所有行（B 的右侧部分）

我将使用的几个概念

• 汉明权重：连续 1 的个数
• 汉明距离：行间不同值的数量

第 1 部分

我们知道行 i 和 j 是可比较的，并且当且仅当存在 i>j i 中的 0 也是 j 中的 0，j 中的任何 1 也是 i 中的 1。这意味着 i 的汉明权重总是会大于 j

的汉明权重

我们还知道，如果 i 和 j 之间的汉明距离为 1，则其间的行集将为空。

因此，对于任何 i，我们可以通过将其至少两个 1 翻转为 0 来找到其可行的 j。

由于我们还不关心顺序，任何具有相同汉明权重的 i 将具有相同数量的可行 j。如果 i 的汉明权重是 1 或 0，则没有可行的 j，因为没有两个 1 可以翻转。

现在，我们可以找到给定汉明权重 i w 的可行 j 的数量，它们在给定的汉明距离 d ，使用函数nchoosek(w,d)。我们知道有用 i 的汉明权重是从 2 到 M，而有用的汉明距离是从 2 到每个 w。

因此我们可以生成大小为MxM 的矩阵n_j，其中行表示汉明权重，列表示汉明距离。

n_j = zeros(M);

for w=2:20
    for d=2:w
        n_j(w,d) = nchoosek(w,d);
    end
end 

如果我们对每一行求和，我们将得到具有给定汉明权重的每个 i 的可用 j 数量。

我们还可以找到在 A 中有多少具有给定有用汉明权重的 i，因为它是完整的并且具有所有可能的组合。我们可以再次使用nchoosek

n_i = zeros(M,1);

for w=2:M
    n_i(w) = nchoosek(M,w);
end

现在我们可以通过将每个汉明权重上可能的 i 的数量乘以可能的 jij /em> 在每个汉明重量上。

sum(n_i .* sum(n_j,2))

在 M=20 的情况下是巨大的，它是 3.4753 e09。这只是对的数量，我们仍然需要找到每对之间的行数。它仍然比 2^20 x 2^20 有所改进，但我现在会考虑我的选择。

继续，我们可以找到对之间的行数作为每对汉明距离的函数。基本上它的2^d -2 它来自这样一个事实，即汉明距离表示您可以翻转以获得新的有效行的值的数量，减去两个原始行 i 和 j em>。

而且由于我们还通过汉明权重将原始 n_j 分开，所以我们可以相乘。

n_r = 2.^(1:M) - 2;

sum(sum((n_i*n_r).*n_j))

对于 M=20 的情况是 1.0925 e12。这意味着 B 的长度必须为 3.4753 e09 并保持 1.0925 e12 值，这是我不会尝试的顺序

---

第 2 部分

现在生成所有对ij。我相信生成一个我们输入 i 并返回有效 j 的函数会更容易和更快。

你说你已经有了A，所以你可以直接索引它，但我会使用de2bi和bi2de函数将索引i关联到每一行A_row_i

A_row_i = de2bi(2^M - i, M)

i = 2^M - bi2de(A_row_i)

首先要计算行i的汉明权重，并用它来获得我们期望的j的数量。

w = sum(A_row_i);

j_s = n_j(w,:);

j_s 将是一个大小为 M 的向量，它告诉我们在给定的汉明距离下我们可以得到多少个 j。

可以修改code from this answer 以创建大小为sum(j_s)xw 的矩阵flips，并对A_row_i 中的1 进行所有必要的翻转以获得A_row_j 的矩阵

flips=[];

for d=2:M
    c = nchoosek(1:w,d);
    out = ones(j_s(d),w);
    out(sub2ind([j_s(d),w],(1:j_s(d))'*ones(1,d),c))=0;
    flips = [flips;out];
end

此时值得注意的是flips只依赖于w，而不依赖于i，所以只有M-2个不同的。从理论上讲，您可以预先计算它们以加快代码速度。

现在我们有了所需的 1 翻转，我们只需要将它们应用到 A_row_i 的副本并恢复 j 的值

A_row_j = repmat(A_row_i,[sum(j_s),1]);
A_row_j(A_row_j==1) = flips;

j = 2^M - bi2de(A_row_j);

现在j 包含相应i 的所有有效j。

您现在可以循环从 1 到 2^M 的所有值并继续构建 B，这将是对双循环的改进，但对于 M=20，它仍然需要一段时间。

---

第 3 部分

如果我们有有效的 i 和 j，那么在行之间获取索引 k 并不难。首先我们计算A_row_i和A_row_j之间的汉明距离d，并用它来创建一个类似于flips的矩阵，这将帮助我们构建A_row_k的所有值并恢复索引k

要翻转的值是A_row_i 和A_row_j 不同的值，我们可以通过bitxor 轻松找到这些值

A_row_i = de2bi(2^M - i, M);
A_row_j = de2bi(2^M - j, M);

d = pdist([A_row_i ; A_row_j],'cityblock');

flips = de2bi(1:2^d-2);

A_row_k = repmat(A_row_i,[n_r(d),1]);

flipable = repmat(bitxor(A_row_i,A_row_j),[n_r(d),1]);

A_row_k(flipable==1) = flips;

k = 2^M - bi2de(A_row_k);

现在您可以循环 B，或者每次为每个 i 获取 j 时运行此代码，无论哪种方式更适合您。

作为结论，我真的建议将其保留为生成器，并在需要值时调用它们，因为对于 M=20，所需的内存量是不切实际的。

Answer 4

简单的数学题！

答：

Function out=A(row , number_of_bits)
        out=~bitget(row-1,number_of_bits:-1:1);
end

A(3,3)
ans =
     1     0     1 %which is identical with your A(3,:)

乙：

function [out]=B(n1,n2,number_of_bits)
   c1=A(n1,number_of_bits);
   c2=A(n2,number_of_bits);
   out=zeros(1,number_of_bits);
   out((c1==0)&(c2==1))=-1;
   if sum(out)>=0
      out=zeros(1,number_of_bits)+2;
      out(c1==0)=0;
      out(c2==1)=1; 
      out=int8(out);
   else
      out=[];
   end

%example:
B(1,4,3)
ans =
     1     2     2 %number 2 serves as wildcard ("*") and means that element could be either 0 or 1, you should discard A(1,:) and A(4,:) from the results so the actual size_of_B is (2^ number_of_wildcards)-2

示例：

B(1,1000000,20) 
ans =
  Columns 1 through 10
     2     2     2     2     1     2     1     1     1     1
  Columns 11 through 20
     2     1     1     1     2     2     2     2     2     2
  %which means that B(1,1000000) has (2^12)-2=4094 rows 

B(1,2,20)
ans =
  Columns 1 through 10
     1     1     1     1     1     1     1     1     1     1
  Columns 11 through 20
     1     1     1     1     1     1     1     1     1     2
  %which means that B(1,2) has (2^1)-2=0 rows 

size_of_B：

function out=size_of_B(b)
     if numel(b)==0
        out=0;
     else
        b(~(b==2))=0;
        b(~(b==0))=1;
        out=(2^sum(b))-2;
     end
%example:
size_of_B(B(1,1000000,20))
ans =
        4094

---

有 (2^19)*((2^20)-1)~= 5e11 对 i, j 进行迭代 如果你有一个什么都不做的for循环：1100秒！所以计算所有非空B的时间是巨大的（即使我们神奇地知道B的哪个元素在任何计算之前都是非空的）

tic; 
 for i=1:1:10^9 
 end 
 toc
Elapsed time is 2.218067 seconds.
tic; 
 for i=1:1:10^10
 end 
 toc
Elapsed time is 22.165445 seconds.

---

%example:
 tic;
 M=10;
 n=1;
 for i=1:2^M
     for j=1:2^M
         if i<j
             c=B(i,j,M);
             if ~(size_of_B(c)==0)
                 d(n,:)=c;
                 n=n+1;
             end
         end
     end
 end
 toc
 [s,~]=size(d);
 s
Elapsed time is 10.095085 seconds. 
s =
          52905 %for M=10 there is 52905 pair of i,j with non-empty-B values;

对于 M=20 大约需要 10*4^10~=10M 秒 ~= 121 天！ （如果内存允许）

---

%script for plotting number of non-empty-B versus M
  d=[];
 for M=1:10
     tic;
     n=1;
     for i=1:2^M
         for j=1:2^M
             if i<j
                 c=B(i,j,M);
                 if ~(size_of_B(c)==0)
                     d(n,:)=c;
                     n=n+1;
                 end
             end
         end
     end
     toc
     [s,~]=size(d);
     d=[];
     h(M)=s;
 end
Elapsed time is 0.005418 seconds.
Elapsed time is 0.002238 seconds.
Elapsed time is 0.000933 seconds.
Elapsed time is 0.003058 seconds.
Elapsed time is 0.011650 seconds.
Elapsed time is 0.044217 seconds.
Elapsed time is 0.170286 seconds.
Elapsed time is 0.696161 seconds.
Elapsed time is 3.024212 seconds.
Elapsed time is 15.836280 seconds.
>> plot(h,'-o')

从图中估计，我们可以发现对于 M=20，大约有 ~3e9 个非空 B 值； （并且至少需要大约 20~60 GB 的内存来存储 B 值，例如我使用的通配符样式，如果您想编写没有通配符样式的所有 B 则更多）

最后我可以提一下，非空 B 的确切数量可以用一个简单的算法计算出来！