Main Content

Matching Pursuit

Sensing Dictionary Creation and Visualization

This example shows how to create and visualize a sensing dictionary consisting of a Haar wavelet at level 2 and the discrete cosine transform-II (DCT-II) basis.

Create the sensing dictionary. The size of the dictionary matrix is 100-by-200. For each basis type, the size of the associated submatrix is 100-by-100. The matrix columns are the basis elements.

xdict = sensingDictionary(Type={'dwt','dct'},Level=[2 0])
xdict = 
  sensingDictionary with properties:

                Type: {'dwt'  'dct'}
                Name: {'db1'  ''}
               Level: [2 0]
    CustomDictionary: []
                Size: [100 200]

Visualize Haar

Extract the 100-by-100 submatrix associated with the Haar wavelet basis.

xmatHaar = subdict(xdict,1:100,1:100);

Display the matrix.

imagesc(xmatHaar)
colormap(gray)
title("Haar Basis")

Figure contains an axes object. The axes object with title Haar Basis contains an object of type image.

Choose and plot two basis elements.

belement = [10 40];
figure
subplot(1,2,1)
plot(xmatHaar(:,belement(1)))
title("Haar basis element: " + num2str(belement(1)))
ylim([-0.5 0.5])
subplot(1,2,2)
plot(xmatHaar(:,belement(2)))
title("Haar basis element: " + num2str(belement(2)))
ylim([-0.5 0.5])

Figure contains 2 axes objects. Axes object 1 with title Haar basis element: 10 contains an object of type line. Axes object 2 with title Haar basis element: 40 contains an object of type line.

To plot all the wavelet basis elements, set showAll to true.

showAll = false;
if showAll
    figure
    for k=1:100
        plot(xmatHaar(:,k))
        title("Haar basis element: " + num2str(k))
        ylim([-0.5 0.5])
        pause(0.1)
    end
end

Visualize DCT-II

Extract the 100-by-100 submatrix associated with the DCT-II basis.

xmatDCT = subdict(xdict,1:100,101:200);

Display the matrix.

figure
imagesc(xmatDCT)
colormap(gray)
title("DCT-II Basis")

Figure contains an axes object. The axes object with title DCT-II Basis contains an object of type image.

Choose and plot four basis elements.

belement = [10 40 60 90];
for k=1:4
    subplot(2,2,k)
    plot(xmatDCT(:,belement(k)))
    title("DCT-II basis element: " + num2str(belement(k)))
    ylim([-0.5 0.5])
end

Figure contains 4 axes objects. Axes object 1 with title DCT-II basis element: 10 contains an object of type line. Axes object 2 with title DCT-II basis element: 40 contains an object of type line. Axes object 3 with title DCT-II basis element: 60 contains an object of type line. Axes object 4 with title DCT-II basis element: 90 contains an object of type line.

To plot all the DCT-II basis elements, set showAll to true.

showAll = false;
if showAll
    figure
    for k=1:100
        plot(xmatDCT(:,k))
        title("DCT-II basis element: " + num2str(k))
        ylim([-0.5 0.5])
        pause(0.1)
    end
end

Orthogonal Matching Pursuit on 1-D Signal

This example shows how to perform orthogonal matching pursuit on a 1-D input signal that contains a cusp.

Load the cuspamax signal. Construct a dictionary consisting of Daubechies extremal phase wavelets at level 2 and the DCT-II basis.

load cuspamax
lsig = length(cuspamax);
A = sensingDictionary(Size=lsig, ...
    Type={'dwt','dct'}, ...
    Name={'db4'},Level=[2 0]);

Use orthogonal matching pursuit to obtain an approximation of the signal in the sensing dictionary with 100 iterations.

[Xr,YI,I,R] = matchingPursuit(A,cuspamax, ...
    maxIterations=100, ...
    Algorithm="OMP",maxerr={"L2",1});

Extract the vectors from the dictionary that correspond to the approximation. Multiply with the associated coefficients and confirm the result is equal to the approximation.

sMat = subdict(A,1:1024,I);
x = sMat*Xr(I,:);
max(abs(x-YI))
ans = 0

Plot the signal, approximation, and residual.

subplot(2,1,1)
plot(cuspamax)
hold on
plot(YI)
hold off
legend("Signal","Approx.")
title("Signal and Approximation")
axis tight
subplot(2,1,2)
plot(R)
title("Residual")
axis tight

Figure contains 2 axes objects. Axes object 1 with title Signal and Approximation contains 2 objects of type line. These objects represent Signal, Approx.. Axes object 2 with title Residual contains an object of type line.

Plot the proportion of retained signal energy for each iteration in the matching pursuit result.

cuspEnergy = cuspamax*cuspamax';
leni = length(I);
recEnergy = zeros(1,leni);
basisElts = subdict(A,1:1024,I);
for k=1:leni
    rec = basisElts(:,1:k)*Xr(I(1:k),:);
    recEnergy(k) = norm(rec)^2/cuspEnergy;
end
figure
plot(recEnergy)
title("Quality / Iteration")
xlabel("Iteration")
ylabel("Quality")

Figure contains an axes object. The axes object with title Quality / Iteration contains an object of type line.

Electricity Consumption Analysis Using Matching Pursuit

This example shows how to compare matching pursuit with a nonlinear approximation in the discrete Fourier transform basis. The data is electricity consumption data collected over a 24-hour period. The example demonstrates that by selecting vectors from a dictionary, matching pursuit is often able to approximate a signal more efficiently with M vectors than any single basis.

Matching Pursuit Using DCT, Fourier, and Wavelet Dictionaries

Load the dataset and plot the data. The dataset contains 35 days of electric consumption. Choose day 32 for further analysis. The data is centered and scaled, so the actual units of usage are not relevant.

load elec35_nor
x = signals(32,:);

plot(x)
xlabel('Minutes')
ylabel('Usage')

Figure contains an axes object. The axes object contains an object of type line.

The electricity consumption data contains smooth oscillations punctuated by abrupt increases and decreases in usage.

Zoom in on a time interval from 500 minutes to 1200 minutes.

xlim([500 1200])

Figure contains an axes object. The axes object contains an object of type line.

You can see the abrupt changes in the slowly-varying signal at approximately 650, 760, and 1120 minutes. In many real-world signals like these data, the interesting and important information is contained in the transients. It is important to model these transient phenomena.

Construct a signal approximation using 50 vectors chosen from a dictionary with orthogonal matching pursuit (OMP). The dictionary consists of the Haar wavelet at level 2, the discrete cosine transform (DCT) basis, a Fourier basis, and the Kronecker delta basis. Then, use OMP to find the best 50-term greedy approximation of the electric consumption data. Plot the result.

Dict = sensingDictionary('Size',length(x), ...
    'Type',{'dwt','dct','fourier','eye'}, ...
    'Name',{'db1'}, ...
    'Level',[2]);
[Xr,YI,I,R] = matchingPursuit(Dict,x, ...
    maxIterations=50, ...
    Algorithm="OMP",maxerr={"L2",1});

plot(x)
hold on
plot(real(YI))
hold off
xlabel('Minutes')
ylabel('Usage')
legend('Original Signal','OMP')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Original Signal, OMP.

You can see that with 50 coefficients, orthogonal matching pursuit approximates both the smoothly-oscillating part of the signal and the abrupt changes in electricity usage.

Determine how many vectors the OMP algorithm selected from each of the subdictionaries.

for k=1:4
    basezind{k} = (k-1)*1440+(1:1440);
end
dictvectors = cellfun(@(x)intersect(I,x),basezind,'UniformOutput',false);
for k=1:4
    fprintf("Subdictionary %d: %d vectors\n",k,numel(dictvectors{k}));
end
Subdictionary 1: 0 vectors
Subdictionary 2: 40 vectors
Subdictionary 3: 2 vectors
Subdictionary 4: 8 vectors

The majority of the vectors come from the DCT and Fourier basis. Given the overall slowly-varying nature of the electricity consumption data, this is expected behavior. The additional vectors capture the abrupt signal changes.

Matching Pursuit Using DCT vs. Full Dictionary

Repeat the OMP with only the DCT subdictionary. Set the OMP to select the 50 best vectors from the DCT dictionary. Construct the dictionary and perform the OMP. Compare the OMP with the DCT dictionary to the OMP with the additional subdictionaries. Notice that adding the Kronecker delta subdictionary shows the abrupt changes in electricity usage more accurately. The advantage of including the Kronecker delta basis is especially clear especially in approximating the upward and downward spikes in usage at approximately 650 and 1120 minutes.

Dict2 = sensingDictionary('Size',length(x),'Type',{'dct'});
[Xr2,YI2,I2,R2] = matchingPursuit(Dict2,x, ...
    maxIterations=50, ...
    Algorithm="OMP",maxerr={"L2",1});

plot(x)
hold on
plot(real(YI2),'linewidth',2)
hold off
title('DCT Dictionary')
xlabel('Minutes')
ylabel('Usage')
xlim([500 1200])

Figure contains an axes object. The axes object with title DCT Dictionary contains 2 objects of type line.

figure
plot(x)
hold on
plot(real(YI),'linewidth',2)
hold off
xlim([500 1200])
title('Full Dictionary')
xlabel('Minutes')
ylabel('Usage')

Figure contains an axes object. The axes object with title Full Dictionary contains 2 objects of type line.

Obtain the best 50-term nonlinear approximation of the signal in the discrete Fourier basis. Obtain the DFT of the data, sort the DFT coefficients, and select the 50 largest coefficients. The DFT of a real-valued signal is conjugate symmetric, so only consider frequencies from 0 (DC) to the Nyquist (1/2 cycles/minute).

xdft = fft(x);
[~,I] = sort(xdft(1:length(x)/2+1),'descend');
ind = I(1:50);

Examine the vector ind. None of the indices correspond to 0 or the Nyquist. Add the corresponding complex conjugate to obtain the nonlinear approximation in the DFT basis. Plot the approximation and the original signal.

indconj = length(xdft)-ind+2;
ind = [ind indconj];
xdftapp = zeros(size(xdft));
xdftapp(ind) = xdft(ind);
xrec = ifft(xdftapp);

plot(x)
hold on
plot(xrec,'LineWidth',2)
hold off
xlabel('Minutes')
ylabel('Usage')
legend('Original Signal','Nonlinear DFT Approximation')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Original Signal, Nonlinear DFT Approximation.

Similar to the DCT dictionary, the nonlinear DFT approximation performs well at matching the smooth oscillations in electricity consumption data. However, the nonlinear DFT approximation does not approximate the abrupt changes accurately. Zoom in on the interval of the data containing the abrupt changes in consumption.

plot(x)
hold on
plot(xrec,'LineWidth',2)
hold off
xlabel('Minutes')
ylabel('Usage')
legend('Original Signal','Nonlinear DFT Approximation')
xlim([500 1200])

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Original Signal, Nonlinear DFT Approximation.