Finding the greatest common divisor of a matrix in MATLAB

Question

Im looking for a way to divide a certain matrix elements with its lowest common divisor.

for example, I have vectors

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

I can tell the lowest common divisor is 2, so the matrix after the division will be

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

What is the built in method that can compute this?

Thanks in advance

p.s. I personally do not like using loops for this computation, it seems like there is built in computation that can perform matrix element division.

Answer 1

Since you don't like loops, how about recursive functions?

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;

The first function iif will define an inline conditional construct. This allows to define a recursive function, gcdrec, to find the greatest common divisor of your array. This iif works like this: it tests whether the first argument is true, if it is, then it returns the second argument. Otherwise it tests the third argument, and if that's true, then it returns the fourth, and so on. You need to protect recursive functions and sometimes other quantities appearing inside it with @(), otherwise you can get errors.

Using iif the recursive function gcdrec works like this:

• if the input vector is a scalar, it returns it
• else if the first component of the vector is 1, there's no chance to recover, so it returns 1 (allows quick return for large matrices)
• else if the input vector is of length 2, it returns the greatest common divisor via gcd
• else it calls itself with a shortened vector, in which the first two elements are substituted with their greatest common divisor.

The function mygcd is just a front-end for convenience.

Should be pretty fast, and I guess only the recursion depth could be a problem for very large problems. I did a quick timing check to compare with the looping version of @Adriaan, using A=randi(100,N,N)-50, with N=100, N=1000 and N=5000 and tic/toc.

1. N=100:
   • looping 0.008 seconds
   • recursive 0.002 seconds
2. N=1000:
   • looping 0.46 seconds
   • recursive 0.04 seconds
3. N=5000:
   • looping 11.8 seconds
   • recursive 0.6 seconds

Update: interesting thing is that the only reason that I didn't trip the recursion limit (which is by default 500) is that my data didn't have a common divisor. Setting a random matrix and doubling it will lead to hitting the recursion limit already for N=100. So for large matrices this won't work. Then again, for small matrices @Adriaan's solution is perfectly fine.

I also tried to rewrite it to half the input vector in each recursive step: this indeed solves the recursion limit problem, but it is very slow (2 seconds for N=100, 261 seconds for N=1000). There might be a middle ground somewhere, where the matrix size is large(ish) and the runtime's not that bad, but I haven't found it yet.

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;

This is probably not the fastest way, but it will do for now. What I did here is initialise some counter, B, to store the sum of all elements in your matrix after taking the mod. the kk is just a counter which runs through integers. mod(A,kk) computes the modulus after division for each element in A. Thus, if all your elements are wholly divisible by 2, it will return a 0 for each element. sum(tmp(:)) then makes a single column out of the modulo-matrix, which is summed to obtain some number. If and only if that number is 0 there is a common divisor, since then all elements in A are wholly divisible by kk. As soon as that happens your loop stops and your common divisor is the number in kk. Since kk is decreased every count it is actually one value too low, thus one is added.

Note: I just edited the loop to run backwards since you are looking for the Greatest cd, not the Smallest cd. If you'd have a matrix like [4,8;16,8] it would stop at 2, not 4. Apologies for that, this works now, though both other solutions here are much faster.

Finally, dividing matrices can be done like this:

divided_matrix = A/kk;

Answer 3

Agreed, I don't like the loops either! Let's kill them -

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

Sample runs

Case #1:

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

Case #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

Here's another approach. Let A be your input array.

1. Get nonzero values of A and take their absolute value. Call the resulting vector B.
2. Test each number from 1 to max(B), and see if it divides all entries of B (that is, if the remainder of the division is zero).
3. Take the largest such number.

Code:

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