Improving Speed and Reducing Memory Consumption with Creation of 2D Sparse Convolution Matrix

In a previous question of mine, Creating Convolution Matrix of 2D Kernel for Different Shapes of Convolution, among answers of the great Matt I came up with the following code:

function [ mK ] = CreateConvMtx2DSparse( mH, numRows, numCols, convShape )

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

numColsKernel   = size(mH, 2);
numBlockMtx     = numColsKernel;

cBlockMtx = cell(numBlockMtx, 1);

for ii = 1:numBlockMtx
    cBlockMtx{ii} = CreateConvMtxSparse(mH(:, ii), numRows, convShape);
end

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        % For convolution shape - 'full' the Doubly Block Toeplitz Matrix
        % has the first column as its main diagonal.
        diagIdx     = 0;
        numRowsKron = numCols + numColsKernel - 1;
    case(CONVOLUTION_SHAPE_SAME)
        % For convolution shape - 'same' the Doubly Block Toeplitz Matrix
        % has the first column shifted by the kernel horizontal radius.
        diagIdx     = floor(numColsKernel / 2);
        numRowsKron = numCols;
    case(CONVOLUTION_SHAPE_VALID)
        % For convolution shape - 'valid' the Doubly Block Toeplitz Matrix
        % has the first column shifted by the kernel horizontal length.
        diagIdx     = numColsKernel - 1;
        numRowsKron = numCols - numColsKernel + 1;
end

vI = ones(min(numRowsKron, numCols), 1);
mK = kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{1});
for ii = 2:numBlockMtx
    diagIdx = diagIdx - 1;
    mK = mK + kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{ii});
end


end

The code is running pretty fast.

But I think there might be ways to improve it farther, specifically in the following lines:

vI = ones(min(numRowsKron, numCols), 1);
mK = kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{1});
for ii = 2:numBlockMtx
    diagIdx = diagIdx - 1;
    mK = mK + kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{ii});
end

This code snippet basically creates a block matrix form the sparse matrices in the cell array cBlockMtx. Where the diagonal of the first element in cBlockMtx is defiend by diagIdx.

Is there a more efficient way to generate this matrix? The main issue is that each time generating kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{ii}) requires a lot of memory.

[Matt J]

I believe you already found a more efficient alternative here?

[Royi Avital]

Well, I tried this:

function [ mK ] = CreateConvMtx2DI( mH, numRows, numCols, convShape )

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

numColsKernel   = size(mH, 2);
numBlockMtx     = numColsKernel;

cBlockMtx = cell(numBlockMtx + 1, 1);

for ii = 1:numBlockMtx
    cBlockMtx{ii} = CreateConvMtx1D(mH(:, ii), numRows, convShape);
end

cBlockMtx{numBlockMtx + 1} = sparse([], [], [], size(cBlockMtx{numBlockMtx}, 1), size(cBlockMtx{numBlockMtx}, 2), 0);

switch(convShape)
    case(CONVOLUTION_SHAPE_FULL)
        % For convolution shape - 'full' the Doubly Block Toeplitz Matrix
        % has the first column as its main diagonal.
        diagIdx     = 0;
        numRowsKron = numCols + numColsKernel - 1;
    case(CONVOLUTION_SHAPE_SAME)
        % For convolution shape - 'same' the Doubly Block Toeplitz Matrix
        % has the first column shifted by the kernel horizontal radius.
        diagIdx     = floor(numColsKernel / 2);
        numRowsKron = numCols;
    case(CONVOLUTION_SHAPE_VALID)
        % For convolution shape - 'valid' the Doubly Block Toeplitz Matrix
        % has the first column shifted by the kernel horizontal length.
        diagIdx     = numColsKernel - 1;
        numRowsKron = numCols - numColsKernel + 1;
end

vC = (numBlockMtx + 1) * ones(numRowsKron, 1);
vR = (numBlockMtx + 1) * ones(numCols, 1);

vC(1:(numBlockMtx - diagIdx)) = (diagIdx + 1):numBlockMtx;
vR(1:(diagIdx + 1)) = (diagIdx + 1):-1:1;

mK = cell2mat(cBlockMtx(toeplitz(vC, vR)));


end

Namely I replaced:

vI = ones(min(numRowsKron, numCols), 1);
mK = kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{1});
for ii = 2:numBlockMtx
    diagIdx = diagIdx - 1;
    mK = mK + kron(spdiags(vI, diagIdx, numRowsKron, numCols), cBlockMtx{ii});
end

With:

vC = (numBlockMtx + 1) * ones(numRowsKron, 1);
vR = (numBlockMtx + 1) * ones(numCols, 1);

vC(1:(numBlockMtx - diagIdx)) = (diagIdx + 1):numBlockMtx;
vR(1:(diagIdx + 1)) = (diagIdx + 1):-1:1;

mK = cell2mat(cBlockMtx(toeplitz(vC, vR)));

It turned out to be slower (x2 slower). So I will be happy to have another idea...