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

Question

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

Answer 1

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.

Answer 2

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

Answer 3

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.

Answer 4

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.)