3
$\begingroup$

As a follow up to Generate the Matrix Form of 2D Convolution Kernel, could someone explain how to generate the matrix form of a 1D convolution kernel?

How different convolutions shapes are handled?
Are there cases where it is better to implement and apply 1D convolution using the matrix form?

$\endgroup$
5
$\begingroup$

The way to build the matrix is playing with indices of the signal data and the convolution kernel.

For example:

function [ mK ] = CreateConvMtx1D( vK, numElements, convShape )
% ----------------------------------------------------------------------------------------------- %
% [ mK ] = CreateConvMtx1D( vK, numElements, convShape )
% Generates a Convolution Matrix for 1D Kernel (The Vector vK) with
% support for different convolution shapes (Full / Same / Valid). The
% matrix is build such that for a signal 'vS' with 'numElements = size(vS
% ,1)' the following are equivalent: 'mK * vS' and conv(vS, vK,
% convShapeString);
% Input:
%   - vK                -   Input 1D Convolution Kernel.
%                           Structure: Vector.
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
%   - numElements       -   Number of Elements.
%                           Number of elements of the vector to be
%                           convolved with the matrix. Basically set the
%                           number of columns of the Convolution Matrix.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3, ...}.
%   - convShape         -   Convolution Shape.
%                           The shape of the convolution which the output
%                           convolution matrix should represent. The
%                           options should match MATLAB's conv2() function
%                           - Full / Same / Valid.
%                           Structure: Scalar.
%                           Type: 'Single' / 'Double'.
%                           Range: {1, 2, 3}.
% Output:
%   - mK                -   Convolution Matrix.
%                           The output convolution matrix. The product of
%                           'mK' and a vector 'vS' ('mK * vS') is the
%                           convolution between 'vK' and 'vS' with the
%                           corresponding convolution shape.
%                           Structure: Matrix (Sparse).
%                           Type: 'Single' / 'Double'.
%                           Range: (-inf, inf).
% References:
%   1.  MATLAB's 'convmtx()' - https://www.mathworks.com/help/signal/ref/convmtx.html.
% Remarks:
%   1.  The output matrix is sparse data type in order to make the
%       multiplication by vectors to more efficient.
%   2.  In case the same convolution is applied on many vectors, stacking
%       them into a matrix (Each signal as a vector) and applying
%       convolution on each column by matrix multiplication might be more
%       efficient than applying classic convolution per column.
% TODO:
%   1.  
%   Release Notes:
%   -   1.1.000     19/07/2021  Royi Avital
%       *   Updated to use modern MATLAB arguments validation.
%   -   1.0.000     20/01/2019  Royi Avital
%       *   First release version.
% ----------------------------------------------------------------------------------------------- %

arguments
    vK (:, :) {mustBeNumeric, mustBeVector}
    numElements (1, 1) {mustBeNumeric, mustBeReal, mustBePositive, mustBeInteger}
    convShape (1, 1) {mustBeNumeric, mustBeMember(convShape, [1, 2, 3])} = 1
end

CONVOLUTION_SHAPE_FULL         = 1;
CONVOLUTION_SHAPE_SAME         = 2;
CONVOLUTION_SHAPE_VALID        = 3;

kernelLength    = length(vK);

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        rowIdxFirst = 1;
        rowIdxLast  = numElements + kernelLength - 1;
        outputSize  = numElements + kernelLength - 1;
    case(CONVOLUTION_SHAPE_SAME)
        rowIdxFirst = 1 + floor(kernelLength / 2);
        rowIdxLast  = rowIdxFirst + numElements - 1;
        outputSize  = numElements;
    case(CONVOLUTION_SHAPE_VALID)
        rowIdxFirst = kernelLength;
        rowIdxLast  = (numElements + kernelLength - 1) - kernelLength + 1;
        outputSize  = numElements - kernelLength + 1;
end

mtxIdx = 0;

% The sparse matrix constructor ignores values of zero yet the Row / Column
% indices must be valid indices (Positive integers). Hence 'vI' and 'vJ'
% are initialized to 1 yet for invalid indices 'vV' will be 0 hence it has
% no effect.
vI = ones(numElements * kernelLength, 1);
vJ = ones(numElements * kernelLength, 1);
vV = zeros(numElements * kernelLength, 1);

for jj = 1:numElements
    for ii = 1:kernelLength
        if((ii + jj - 1 >= rowIdxFirst) && (ii + jj - 1 <= rowIdxLast))
            % Valid otuput matrix row index
            mtxIdx = mtxIdx + 1;
            vI(mtxIdx) = ii + jj - rowIdxFirst;
            vJ(mtxIdx) = jj;
            vV(mtxIdx) = vK(ii);
        end
    end
end

mK = sparse(vI, vJ, vV, outputSize, numElements);


end

Few remarks:

  1. When the length of the filter is smaller than the number of samples of the signal the matrix pattern is highly sparse. Hence I used sparse matrix in the code.
  2. The code supports 3 different shapes of the convolution: full, same and valid. I use the same interpretation as MATLAB's conv().
  3. The code has gone through a validation process and matches MATLAB's conv().

For a simple filter (Taken from Deblurring 1D data using Direct Inverse Filtering) vK = [1; 4; 6; 4; 1] / 16; and a signal of 50 samples we get:

enter image description here

As one can see, the different shapes are different on how they handle the edges.

In Deep Learning in some cases convolutions are implemented in a matrix form.
The reasoning is, if you apply the same kernel to multiple images you might gain performance:

$$ \boldsymbol{y}_{i} = \boldsymbol{k} \ast \boldsymbol{x}_{i} \Rightarrow Y = K X $$

Where $ \boldsymbol{y}_{i} $ is basically the i -th column of $ Y $ and the same for $ X $.

This form is also used often to solve optimization problems, yet there is a better way to implement this. But this is for another question.

The code is available at my StackExchange Signal Processing Q76344 GitHub Repository.

$\endgroup$
2
  • $\begingroup$ Great answer by the code. Could you elaborate on the optimization? How is it used? $\endgroup$
    – Thomas
    Jul 19 at 14:19
  • $\begingroup$ There is a lot to write on this so I think it deserves its own question. But it should be a focused one. $\endgroup$
    – Royi
    Jul 20 at 4:28
1
$\begingroup$

A convolution matrix is really just a diagonal band-structure matrix, where every row is all zeros but for the elements around the diagonal, which are identical (but shifted) for every row: the elements of the kernel.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.