Fast/efficient way to decompose separable integer 2D filter coefficients

I would like to be able to quickly determine whether a given 2D kernel of integer coefficients is separable into two 1D kernels with integer coefficients. E.g.

 2   3   2
4   6   4
2   3   2


is separable into

 2   3   2


and

 1
2
1


The actual test for separability seems to be fairly straightforward using integer arithmetic, but the decomposition into 1D filters with integer coefficients is proving to be a more difficult problem. The difficulty seems to lie in the fact that ratios between rows or columns may be non-integer (rational fractions), e.g. in the above example we have ratios of 2, 1/2, 3/2 and 2/3.

I don't really want to use a heavy duty approach like SVD because (a) it's relatively computationally expensive for my needs and (b) it still doesn't necessarily help to determine integer coefficients.

Any ideas ?

FURTHER INFORMATION

Coefficients may be positive, negative or zero, and there may be pathological cases where the sum of either or both 1D vectors is zero, e.g.

-1   2  -1
0   0   0
1  -2   1


is separable into

 1  -2   1


and

-1
0
1

• I remember trying to figure this out way back in college. I almost succeeded, but I don't remember how. = ) I can't stop thinking about it now that you mentioned it! – Phonon Mar 28 '12 at 16:05
• @Phonon: heh - well keep thinking - I could use some inspiration on this one. ;-) – Paul R Mar 28 '12 at 16:20
• Is it possible to do the same thing but for double or float values ? – Diego Catalano Oct 6 '14 at 18:09
• @DiegoCatalano: see Denis's answer below, and the question he links to on math.stackexchange.com - I think that might work for the more general case of floating point coefficients. – Paul R Oct 6 '14 at 20:56
• @PaulR, How could one contact you on email? Thank You. – Royi Jul 19 '16 at 6:24

I have taken @Phonon's answer and modified it somewhat so that it uses the GCD approach on just the top row and left column, rather than on row/column sums. This seems to handle pathological cases a little better. It can still fail if the top row or left column are all zeroes, but these cases can be checked for prior to applying this method.

function [X, Y, valid] = separate(M)    % separate 2D kernel M into X and Y vectors
X = M(1, :);                          % init X = top row of M
Y = M(:, 1);                          % init Y = left column of M
nx = numel(X);                        % nx = no of columns in M
ny = numel(Y);                        % ny = no of rows in M
gx = X(1);                            % gx = GCD of top row
for i = 2:nx
gx = gcd(gx, X(i));
end
gy = Y(1);                            % gy = GCD of left column
for i = 2:ny
gy = gcd(gy, Y(i));
end
X = X / gx;                           % scale X by GCD of X
Y = Y / gy;                           % scale Y by GCD of Y
scale = M(1, 1) / (X(1) * Y(1));      % calculate scale factor
X = X * scale;                        % apply scale factor to X
valid = all(all((M == Y * X)));       % result valid if we get back our original M
end


Many thanks to @Phonon and @Jason R for the original ideas for this.

Got it! Posting MATLAB code, will post an explanation tonight or tomorrow

% Two original arrays
N = 3;
range = 800;
a = round( range*(rand(N,1)-0.5) )
b = round( range*(rand(1,N)-0.5) )

% Create a matrix;
M = a*b;
N = size(M,1);

% Sanity check
disp([num2str(rank(M)) ' <- this should be 1!']);

% Sum across rows and columns
Sa = M * ones(N,1);
Sb = ones(1,N) * M;

% Get rid of zeros
SSa = Sa( Sa~=0 );
SSb = Sb( Sb~=0 );

if isempty(SSa) | isempty(SSb)
break;
end

% Sizes of array without zeros
Na = numel(SSa);
Nb = numel(SSb);

% Find Greatest Common Divisor of Sa and Sb.
Ga = SSa(1);
Gb = SSb(1);

for l=2:Na
Ga = gcd(Ga,SSa(l));
end

for l=2:Nb
Gb = gcd(Gb,SSb(l));
end

%Divide by the greatest common divisor
Sa = Sa / Ga;
Sb = Sb / Gb;

%Scale one of the vectors
MM = Sa * Sb;
Sa = Sa * (MM(1) / M(1));

disp('Two arrays found:')
Sa
Sb
disp('Sa * Sb = ');
Sa*Sb
disp('Original = ');
M

• Thanks - this is great - I was lying awake last night thinking about factorising the coefficients etc but just using the GCD like this is a lot more simple and elegant. Unfortunately there is still one wrinkle to iron out - it needs to work with both positive and negative coefficients and this can lead to degenerate cases, e.g. A=[-2 1 0 -1 2]; B=[2 -3 6 0 -1]; M=A'*B;. The problem here is that sum(A) = 0 so Sb = [0 0 0 0 0]. I'm going to try modifying your algorithm so that it uses the sum of absolute values of the coefficients and see if that helps. Thanks again for your help. – Paul R Mar 29 '12 at 9:13
• OK - it looks like you can still get the GCDs and do the scaling by using abs(M), i.e. Sa=abs(M)*ones(N,1); Sb=ones(1,N)*abs(M); and then continue as above, but I can't yet see how to restore the signs to Sa, Sb at the end. I've added a pathological example that illustrates the problem in the original question above. – Paul R Mar 29 '12 at 9:59
• I think I have a working solution now - I've posted it as a separate answer, but the credit goes to you for the underlying idea. Thanks again ! – Paul R Mar 29 '12 at 14:50

Maybe I'm trivializing the problem, but it seems like you could:

• Break the $N$-by-$M$ matrix $\mathbf{A}$ into rows $\mathbf{a_i}$, $i = 0, 1, \ldots , N-1$.
• For each row index $j > 0$:

• Elementwise divide $\mathbf{a_j}$ by $\mathbf{a_0}$ to yield the ratio of each element in row $j$ to its counterpart in row zero $\mathbf{r_j}$. This would need to be done using floating-point or some other fractional arithmetic.
• Inspect the elements in $\mathbf{r_j}$ to determine if they are constant. You'll need to allow for some tolerance in this comparison to allow for floating-point rounding.
• If the values in $\mathbf{r_j}$ are constant, then row $\mathbf{a_j}$ is a scalar multiple of row $\mathbf{a_0}$. Store the ratio of row $j$ to row $0$ in a list $\mathbf{x}$.
• If the ratio vector is not constant, then row $\mathbf{a_j}$ is a scalar multiple of row $\mathbf{a_0}$, so you will not be able to decompose the matrix in this way. Stop checking.
• If all of the rows were deemed to be constant multiples of row 0 in the tests above, then take the list of row ratio scalars $\mathbf{x}$ and normalize it by the smallest ratio:

$$x_{k,norm} = \frac{x_k}{\min_{i=0}^{N-1} x_i}$$

• After doing so, the normalized list of ratios $\mathbf{x_{norm}}$ will contain a value of 1 for the row that has the smallest norm. You want to see if you can scale this list in some way to yield a list that contains all integer coefficients. A brute-force approach could just do a linear search, that is, calculate: $$\mathbf{x_{scaled}} = K\mathbf{x_{norm}}, K = 1, 2, \ldots, M$$ For each value of $K$, after calculating the scaled list of ratios, you would need to check each to see if they are integer-valued (again, with some tolerance). Set $M$ equal to the largest denominator that you're willing to look for in the inter-row ratios.

Not the most elegant method, and it's likely that there's a better way, but it should work, is pretty simple to implement, and should be relatively speedy for modestly-sized matrices.

• Thanks - I think I was probably heading in something like this direction before I got bogged down in the details. It's not 100% clear to me that you'll always arrive at a solution using this method, but anyway, I should probably code this up and try it with a few examples. I have a hunch that it may need to be applied both row-wise and column-wise to see which yields the "best" solution. Thanks for taking the time to spell out the details - I'll get busy with it and let you know how it works out. – Paul R Mar 28 '12 at 17:16
• Couldn't you find the greatest common divisor of the first elements of the rows, and use that to determine your basis vector? – Jim Clay Mar 28 '12 at 17:59
• @JimClay: Yes, that's effectively what you're doing at the end, if you have that functionality available. – Jason R Mar 28 '12 at 18:07

Another way is to find a separable approximation $x \otimes y \otimes z$ to your kernel $A$, and see how close it is. Two ways of doing this, i.e. to minimize $|A - x \otimes y \otimes z|$:
1) brute-force optimize $x\ y\ z$; this takes time ~ the sum, not the product, of their lengths
2) for fixed $y$ and $z$, the optimal $x$ is just a projection, so optimize $x\ y\ z\ x\ y\ z\ ...$ in turn.

(From approximate-a-convolution-as-a-sum-of-separable-convolutions on math.stackexchange.)

• Try not to answer questions with unexplained links. It is better to explain the necessary details in your answer and include the link only for reference; that way if the link breaks the essential details of the answer are still there. – Sam Maloney Apr 21 '13 at 0:08
• @SamMaloney: I see no reason why that's necessary. The link explains everything in detail. It will still pop up in Q&A search. So why not? – Naresh Apr 21 '13 at 2:46
• @Naresh I only mention it because one of the goals of the stack exchange sites is to build up a repository of answered questions for future reference. So while I understand that this particular link is to another SE site and should be fairly safe, it is a general best practice not to count on links still working several years from now. Giving a general outline of these "two easy methods' in the answer would ensure that information is retained even if something happens to the linked question. As I said though, this was more of a general comment on best practices regarding links in answers. – Sam Maloney Apr 22 '13 at 18:39