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

Question

我有一个尺寸为 BxM 的由零和一组成的 Matlab 矩阵 A。

具体来说，A 包含了所有在考虑顺序的情况下 M 个位置上可能的 0 和 1 的排列方式（因此，B=2^M）。

以下是该矩阵的构建规则：

- 对于任意 i=1,...,N，都有 A(i,:)>A(j,:) 或 A(i,:)><A(j,:) 对于 j=i+1,...,N。 - 表示 A(i,:)>=A(j,:) 意味着对于所有的 h=1,...,M，A(i,h)>=A(j,h)。 - 表示 A(i,:)>A(j,:) 意味着对于所有的 h=1,...,M，至少存在一个 h 满足 A(i,h)>A(j,h)。 - 表示 A(i,:)><A(j,:) 意味着无法确定 A(i,:)>=A(j,:) 还是 A(j,:)>=A(i,:)。

例如，当 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 条件未满足许多行对），我真正需要的是仅跟踪非空单元格的内容和坐标。因此，也许，“稀疏单元格”可以帮助我解决这个问题。

非常感谢任何建议。

- TEX

---

一个值没有给出的单元格为空。您不需要预分配 (2^20)x(2^20) 的单元格，只需预分配一个准备填充的结构即可。一个单元格已经是一个“稀疏”的单元格。您只是有太多的数据。请尝试大致估算您需要填充其中的(2^20)x(2^20)值的数量。 - Ander Biguri

@AnderBiguri 我明白了，但是预测满足条件的坐标并不容易。我不知道应该如何做到这一点。 - TEX

进行近似数学计算。如果一个大小为(2^20)x(2^20)的矩阵只有1%是稀疏的，并且在每个非零值中存储一个double（而不是单元格），则需要80Gb的RAM。请注意，您的建议似乎需要比这更多的值。 - Ander Biguri

你的代码没有产生你所期望的输出。 你的代码产生了不正确的 B{1,6}, B{1,7}，B{1,8}，B{3,8}，B{5,8}。 - Sardar Usama

现在应该是正确的。谢谢。每个单元格内行的顺序无关紧要。 - TEX

Answer 1

不要像"旧版Matlab"那样预先计算和记忆A矩阵，而是编写一个"生成器函数"，为任何行索引提供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 

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

- George Rey

---

你真的需要一个稀疏的B矩阵吗？还是只需要一个包含坐标和该行答案列表的列表就足够了？如果是后者，我可以帮助你。 - George Rey

user3285148，你是认真的吗？ - George Rey

你说“我在Matlab中有一个由0和1组成的矩阵A，其维度为BxM”，使得size(A,1) = B。使用常量比调用函数size更快。 - George Rey

关于双重循环，这不是三分钟就能解决的问题。你需要将方程式实际转换成不同的形式。我以前用Mathematica中的符号操作做过这样的事情。如果没有符号转换，可以尝试使用parallel-for和/或在C中编写内部循环。普通的Matlab在处理布尔值时特别低效（它们被存储为double(!)），类型转换通常会使情况更糟。 - George Rey

@excaza: 是的，但只有在你明确告诉它的情况下才会这样。如果您执行 q=1，则 class(q)=='double'；必须执行 q=true 才能获得 logical 类型。 - George Rey

Answer 2

不必采用暴力方法得出答案，我们可以进行更多的数学计算，来预测测试结果，然后只需保存预测结果即可。

我们将做出以下假设：

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

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

3）在M=5的情况下，绘制结果的位置后，我们发现遵从条件创建了分形模式。这个分形模式可以使用Lucas Correspondence Theorem重新创建，“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个单元格的信息的情况下按需获取数据。请告知我是否符合您的要求。

- Alea Kootz

---

你正在有效地为空间B重新创建左下方的Sierpinski筛子。 - Alea Kootz

你的函数在M=3时创建了19个结果，而在M=3时只有7个有效对。 - Hadi

@darenshan，如果你按照它们的 A(j) 索引将这 19 对数据分组，你会得到 7 个集合。 - Alea Kootz

谢谢。为什么当M=5时，我得到的单元格B是32x31而不是32x32？此外，当M=20时，它需要无限的时间。 - TEX

你能否写出只需要给定 j 和 h 就可以得到 B{j,h} 的代码行。例如，像我的例子中一样，B{1,6}=[1 1 0; 0 1 1]。 - TEX

Answer 3

我将把你的问题分成三个部分。

1. 找到B所需的大小 2. 找到A中至少有一行在中间的可比较行对（B的左半部分） 3. 找到每个可比较对之间的A行（B的右半部分）

我会使用以下几个概念：

- 汉明重量：一行中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可以翻转。

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

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

n_j = zeros(M);

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

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

由于A是完整的并具有可能的每个组合，因此我们还可以找到在A中具有给定有用海明重量的 i 的数量。我们可以再次使用nchoosek。

n_i = zeros(M,1);

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

现在，我们可以通过将每个汉明重量上可能的i的数量乘以每个汉明重量上可能的j的数量来计算ij总对数。

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

在M等于20的情况下，它非常庞大，为3.4753亿。这仅仅是成对数量，我们还需要找到每一对之间行数的数量。虽然这比2^20 x 2^20有所改进，但是此时我会考虑我的选择。

继续查找，我们可以将每对之间的行数作为每对汉明距离的函数来查找。基本上是2^d -2，这是因为汉明距离告诉你可以翻转多少个值以获得一个新的有效行，减去两个原始的行i和j。

由于我们按汉明重量保留了原始的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个值，这是一个非常大的数字。

---

第二部分 现在要生成所有的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。

可以修改这个答案中的代码，创建一个大小为 sum(j_s) x w 的矩阵 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

值得注意的是，此时的翻转仅依赖于，而不依赖于，因此只有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，这需要一段时间。

---

第三部分

如果我们有有效的 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，所需的内存量是不实际的。

- Noel Segura Meraz

---

谢谢。假设我只想获取B{1,3}，我该怎么做？ - TEX

假设我有索引 i 和 j（我不知道行是否可比较），如何获取它们之间的行索引？ - TEX

要获取B{1,3}中的行，只需将i=1和j=3设置为代码中第3部分并运行，您将在k中获得这些行。（您需要使用第1部分中定义的n_r，但它只是n_r = 2.^(1:M) - 2;。） - Noel Segura Meraz

我的示例中，它比循环内的行慢。 - TEX

哦，我明白你的意思了。你仍然可以使用我的答案的第二部分来减少双重1:2^M循环，变成单个1:2^M循环。这与你的代码结合起来一定会大大减少时间。 - Noel Segura Meraz

Answer 4

简单的数学问题！

A：

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,:)

B：

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 

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 集所需的时间是巨大的（即使我们神奇地知道在进行任何计算之前哪些元素是非空的）。

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值；（如果想要像我用通配符样式那样存储B值，则至少需要20~60 GB的内存，而如果想要写出所有B值，所需内存则更多）

最后要提到的是，非空B的确切数量可以使用一种简单的算法来计算！

- Hadi

---

@user3285148 是的，但在进行任何计算之前，你不可能知道 B{i,j} 是否为空。我的回答只是说，你可以在短时间内计算出任何 B{i,j}（即使是 M20），那么为什么要计算并存储所有的 B 值呢？ - Hadi

关于可行性，这是可行的，但需要大约121天的计算来计算和存储所有非空B值 :) - Hadi

谢谢。让j=4，i=5。您能在回答末尾总结一下获取B{4,5}的整个代码（不包括解释）吗？ - TEX

好的，非常感谢。 - TEX

我不确定我理解了。像我的问题一样说M=3，假设我想得到B{1,6}=[1 1 0; 0 1 1];（参见我的示例）。我该怎么做？ - TEX