在MATLAB中找到矩阵的最大公约数

Question

我正在寻找一种方法将某个矩阵元素与其最低公约数相除。

例如，我有向量

[0,0,0; 2,4,2;-2,0,8]

我知道最小公约数是2，所以除法后的矩阵就是

[0,0,0;1,2,1;-1,0,4]

可以计算这个的内置方法是什么？

提前致谢

附言我个人不喜欢使用循环进行这种计算，似乎有内置的计算可以执行矩阵元素除法。

Answer 1

既然你不喜欢循环，那么递归函数呢？

iif = @(varargin) varargin{2 * find([varargin{1:2:end}], 1, 'first')}();
gcdrec=@(v,gcdr) iif(length(v)==1,v, ...
                     v(1)==1,1, ...
                     length(v)==2,@()gcd(v(1),v(2)), ...
                     true,@()gcdr([gcd(v(1),v(2)),v(3:end)],gcdr));
mygcd=@(v)gcdrec(v(:)',gcdrec);

A=[0,0,0; 2,4,2;-2,0,8];
divisor=mygcd(A);
A=A/divisor;

第一个函数iif 将定义一个内联条件构造。这允许定义递归函数gcdrec，以找到数组的最大公约数。这个iif 是这样工作的：它测试第一个参数是否为true，如果是，则返回第二个参数。否则，它会测试第三个参数，如果 那是 true，则返回第四个参数，依此类推。您需要使用@() 保护递归函数，有时还需要保护其中出现的其他数量，否则可能会出错。

使用iif 递归函数gcdrec 的工作方式如下：

• 如果输入向量是标量，则返回它
• 否则，如果向量的第一个分量为 1，则无法恢复，因此返回 1（允许快速返回大型矩阵）
• 否则，如果输入向量的长度为 2，则通过 gcd 返回最大公约数
• 否则，它会使用一个缩短的向量来调用自身，其中前两个元素被替换为它们的最大公约数。

函数mygcd只是一个方便的前端。

应该很快，我猜只有递归深度可能是非常大的问题的问题。我做了一个快速的时间检查来与@Adriaan 的循环版本进行比较，使用A=randi(100,N,N)-50，N=100、N=1000 和N=5000 和tic/toc。

1. N=100:
   • 循环 0.008 秒
   • 递归 0.002 秒
2. N=1000:
   • 循环 0.46 秒
   • 递归 0.04 秒
3. N=5000:
   • 循环 11.8 秒
   • 递归 0.6 秒

更新：有趣的是，我没有超出递归限制（默认为 500）的唯一原因是我的数据没有公约数。设置一个随机矩阵并将其加倍将导致已经达到N=100 的递归限制。因此，对于大型矩阵，这是行不通的。再说一次，对于小矩阵@Adriaan 的解决方案非常好。

我还尝试在每个递归步骤中将其重写为输入向量的一半：这确实解决了递归限制问题，但它非常慢（N=100 2 秒，261 秒N=1000)。某处可能有一个中间地带，矩阵大小很大（ish）并且运行时还不错，但我还没有找到它。

Answer 2

 A = [0,0,0; 2,4,2;-2,0,8];
 B = 1;
 kk = max(abs(A(:))); % start at the end
 while B~=0 && kk>=0
     tmp = mod(A,kk);
     B = sum(tmp(:));
     kk = kk - 1;
 end
 kk = kk+1;

这可能不是最快的方法，但现在可以了。我在这里所做的是初始化一些计数器B，以在获取mod 后存储矩阵中所有元素的总和。 kk 只是一个遍历整数的计数器。 mod(A,kk) 计算 A 中每个元素除法后的模数。因此，如果您的所有元素都可以被 2 整除，那么它将为每个元素返回 0。 sum(tmp(:)) 然后从模矩阵中生成一列，将其相加以获得一些数字。当且仅当该数字为 0 时，存在公约数，因此 A 中的所有元素都可以被 kk 整除。一旦发生这种情况，您的循环就会停止，您的公约数就是kk 中的数字。由于kk 每次计数都会减少一个值，因此实际上是一个值太低，因此添加了一个。

注意：我刚刚编辑了循环以运行向后，因为您正在寻找最大的 cd，而不是最小的 cd。如果你有一个像[4,8;16,8] 这样的矩阵，它将停在2，而不是4。对此深表歉意，这现在有效，尽管这里的其他两种解决方案都快得多。

最后，除法矩阵可以这样完成：

divided_matrix = A/kk;

Answer 3

同意，我也不喜欢循环！让我们杀了他们-

unqA = unique(abs(A(A~=0))).';             %//'
col_extent = [2:max(unqA)]';               %//'
B = repmat(col_extent,1,numel(unqA));
B(bsxfun(@gt,col_extent,unqA)) = 0;
divisor = find(all(bsxfun(@times,bsxfun(@rem,unqA,B)==0,B),2),1,'first');
if isempty(divisor)
    out = A;
else
    out = A/divisor;
end

样本运行

案例#1：

A =
     0     0     0
     2     4     2
    -2     0     8
divisor =
     2
out =
     0     0     0
     1     2     1
    -1     0     4

案例#2：

A =
     0     3     0
     5     7     6
    -5     0    21
divisor =
     1
out =
     0     3     0
     5     7     6
    -5     0    21

Answer 4

这是另一种方法。让A 成为您的输入数组。

1. 获取A 的非零值并取其绝对值。调用结果向量B。
2. 测试从1 到max(B) 的每个数字，看看它是否除以B 的所有条目（即除法的余数是否为零）。
3. 取最大的数。

代码：

A = [0,0,0; 2,4,2;-2,0,8];                  %// data
B = nonzeros(abs(A));                       %// step 1
t = all(bsxfun(@mod, B, 1:max(B))==0, 1);   %// step 2
result = find(t, 1, 'last');                %// step 3