convlution matrix in four lines

This simple function generates the convolution matrix.

Description

% CONV_MTX Convolution matrix.
% If 'x' is a column vector of length 'nx', then conv_mtx(x,nh) gives a toeplitz
% matrix 'X' of size (nx+nh-1) times (nh). 'nh' is the length of the column vector
% 'h' with which 'x' is convolved. Thus finally 'X*h' is same as conv(x,h).
% If 'x' and 'h' are row vectors, then 'X' is of size (nh) times (nx+nh-1).
% Then 'h*X' is same as conv(x,h). This does not depend on the in built command
% 'TOEPLITZ'.

• MATLAB release: MATLAB 7.4 (R2007a)

[John D'Errico]

Without looking more deeply at this code, it produces a banded matrix. As such, it is best stored/returned in sparse form, especially important when the matrix may potentially be large. Since the sparse form of the matrix can be generated using a single call to sparse, I'd strongly recommend that alternative. The eventual matrix multiplications would of course be much more efficient too. The sparse code might look like this (please add the necessary additional code to make it friendly)

function M = convmtx(x,nh)
n = length(x);
M = sparse(bsxfun(@plus,(1:n)',0:(nh-1)), ...
    repmat(1:nh,n,1),repmat(x,1,nh),n+nh-1,nh);

Note that it could have been done easily without bsxfun, for users of older Matlab releases. See how it does for efficiency on a large matrix.

x = (1:5)';
timeit(@() conv_mtx(x,500))
ans =
     0.012447

timeit(@() convmtx(x,500))
ans =
    0.0029383

Thus for large convolutions, generating the full matrix is too much work. Plus, simply storing it will be impossible for larger convolutions unless you use a sparse form.

Next, look at the memory storage required.

M0 = conv_mtx(x,2500);
M1 = convmtx(x,2500);
whos M*
  Name Size Bytes Class Attributes
  M0 2504x2500 50080000 double
  M1 2504x2500 160004 double sparse

Had nh been much larger, the conv_mtx code would have overrun my available memory.

The final point to make is the time required for the final matrix*vector multiply when you might use the result.

timeit(@() M0*ones(2500,1))
ans =
      0.13361

timeit(@() M1*ones(2500,1))
ans =
   0.00088314

These are the reasons why we use sparse storage in Matlab for such banded matrices.