# Fast / Efficient Way to Decompose a Separable integer 2D Filter Coefficients without the SVD

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! Mar 28, 2012 at 16:05
• @Phonon: heh - well keep thinking - I could use some inspiration on this one. ;-) Mar 28, 2012 at 16:20
• Is it possible to do the same thing but for double or float values ? Oct 6, 2014 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. Oct 6, 2014 at 20:56
• Seems someone calls this the «inverse outer product»: Apr 7, 2020 at 9:18

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. Mar 29, 2012 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. Mar 29, 2012 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 ! Mar 29, 2012 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. Mar 28, 2012 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? Mar 28, 2012 at 17:59
• @JimClay: Yes, that's effectively what you're doing at the end, if you have that functionality available. Mar 28, 2012 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. Apr 21, 2013 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? Apr 21, 2013 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. Apr 22, 2013 at 18:39