Recovery set of signals same sparsity

Hello. I shall plot the result of the following two actions:

1. Reconstruct a set (/series) of signals, e.g. 100 signals, all with same length n, and same sparsity k: what varies randomly is the position of non-zero values.
2. Compute probability of successful reconstruction of that set of signals, varying the number of observations (m), e.g. from 1 to 40% n.

While I coded/solved the two problems above with - 1 n-long - signal in input, I am experiencing issues dealing with a - set - of signals, I miss some Matlab/Simulink knowledge still. I am not sure that representing the set of signals in the shape of a matrix of n rows and 100 columns (i.e. each column is 1 signal per se) be the best solution, especially in terms of performances (e.g. application of BP/BPDN recovery).

When I use - one - signal in input, I generate it as:

support = randsample(n,k);
sparsesignal = zeros(n,1);
sparsesignal(support) = randn(k,1);

When I use - a set - of e.g. 100 signals in input, I generate it as:

support = randsample(n,k);
sparseset = zeros(n,100);
for j = 1:100
    sparseset(support,j) = randn(k,1);
end

Concerning action 1: I miss the logic of plotting the reconstruction of a - set of - signals, i.e.

• how would the plot look like?
• what algorithm would best suit the process, e.g. apply reconstruction on each matrix column separately, but then plot what ?

Concerning action 2: For each value m of observations, I compute the reconstruction error for each column of the signal matrix, then I average the sum error (e.g. MSE) over 100. The algorithm foresees tot iterations per each value of m (tot random observations per each m), e.g. tot=1000. I don't seem this solution to be efficient computationally. What to plot here is clear, but would you have any hint on how to improve the algorithm please?

Many thanks for your time and support.

Answer 1

[Stephen Jue]

If I understand correctly, you would like to reconstruct a set of signals with the same length “n” and sparsity “k”, but with the non-zero values randomly positioned. For each signal you are calculating the reconstruction error for the number of observations m = 1 to 40% of n.

You can make your code more efficient by replacing the “for” into a one-line vectorized statement, like so:

support = randsample(n,k);
sparseset = zeros(n,100);
sparseset(support,:) = randn(k,100);

Calculating the mean squared error can be vectorized as well, by using the “mean” function on the difference between your observed and predicted values. By default, the “mean” function acts on a matrix by calculating the mean of each column.

One thing I noticed in your code is that all of the signals have non-zero values in the same positions. Is this what you want? It seems by your problem statement that you want each signal to have nonzero values in random positions. To do this you’ll need to calculate “support” 100 times.

Also, depending on how sparse your signals are, your algorithm could be more space-efficient by using a sparse matrix rather than a normal matrix.

[VMat]

Thank you very much Stephen. You understood correctly - I should probably write that my question was about 'joint sparse recovery' of an ensemble of signals, and Multi Measurement Vectors (MMV).

Let's have: an input ensemble of q signals in the matrix X0, and the reconstructed ensemble in the matrix Xhat. Each column of those matrices represents 1 signal per se: same sparsity k for all signals/columns, but as you say, non-zeros shall appear randomly, not necessarily in the same indexes. That means applying q 'supports' of fixed k, please correct if wrong.

Now, what would be the best matlab coding to compute the MSE/RMSE between X0 and Xhat ? Considering only 1 signal in the ensemble, cases are that: if the 2 corresponding columns match (given a tolerance) then 'error' is 0 and probability of exact recovery 1; if they don't match, I'd need to know 'how much they don't match', i.e. the error (and so the probability of recovery < 1). This is the 'error' for 1 single signal: to compute the error at matrix level I'd average all the single errors on q. Current pseudo code (doesn't work) for this logic would be:

   % outer for cycle sets flag = 0;
      for i=1:iter
          equi = all(abs(Xhat-X0) < 1.0e-6);
          if (equi == 1)
              flag = flag+1;
          else
              % suppose to have q signals in the input ensemble, i.e. vector equi has q elements, whose value can be 1 or 0
              flag = flag+(sum(equi)/q);
          end 
        end
      % j is index of outer for cycle, indicates measurements
    % vector res (to plot) shall tell the probability of successful recovery or the MSE/RMSE for each number of measurements
      res(j) = flag/iter;
Don't know how to use the 'mean' function here, I found 'immse' function but results are not that good.

Finally, I am not clear on your statement: "depending on how sparse your signals are, your algorithm could be more space-efficient by using a sparse matrix rather than a normal matrix". Seems relevant, could you please re-explain that? maybe with a code example? Sparsity k shall range around 3% / 8% of n (length of each signal in the ensemble).

I would most appreciate any additional help: many thanks!

Answer 2

[Stephen Jue]

Your code looks good, but be careful with the construction of your signals. Since you do want each signal to have random numbers in random positions, you will want to calculate "support" within the for-loop, like so:

sparseset = zeros(n,100);
for j = 1:100
    support = randsample(n,k);
    sparseset(support,j) = randn(k,1);
end

The "immse" function should work fine.

To elaborate further on the "sparse matrix" that I mentioned:

MATLAB has a data structure called the "sparse matrix", which only stores nonzero values of the matrix and their positions. This data structure makes sense when dealing with sparse data like yours, since it only ranges from about 3-8%.

This should get you started: http://www.mathworks.com/help/matlab/ref/sparse.html

To build the matrix that you want, you can use code similar to this:

n = 100; % length of each signal
k = 3; % sparsity (number of nonzero values per signal)
numSignals = 1000; % number of signals

row_indices = [];
column_indices = [];
random_numbers = [];

for j = 1:numSignals
    row_indices = [row_indices; randperm(n, k)']; % Random row indices (1 to n inclusive)
    column_indices = [column_indices; j*ones(k, 1)]; % Column indices
    random_numbers = [random_numbers; rand(k, 1)]; % Random numbers at each point
end
signalMatrix = sparse(row_indices, column_indices, random_numbers, n, numSignals);

"signalMatrix" should contain 1000 signals of length 100 and sparsity of 3. This matrix takes up ~56 kB on my machine, versus a 1000-by-100 normal matrix, which takes up ~800 kB.

I hope that helps you.

Answer 3

[VMat]

Thank you very much Stephen. You understood correctly - I should probably write that my question was about 'joint sparse recovery' of an ensemble of signals, and Multi Measurement Vectors (MMV).

Let's have: an input ensemble of q signals in the matrix X0, and the reconstructed ensemble in the matrix Xhat. Each column of those matrices represents 1 signal per se: same sparsity k for all signals/columns, but as you say, non-zeros shall appear randomly, not necessarily in the same indexes. That means applying q 'supports' of fixed k, please correct if wrong.

Now, what would be the best matlab coding to compute the MSE/RMSE between X0 and Xhat ? Considering only 1 signal in the ensemble, cases are that: if the 2 corresponding columns match (given a tolerance) then 'error' is 0 and probability of exact recovery 1; if they don't match, I'd need to know 'how much they don't match', i.e. the error (and so the probability of recovery < 1). This is the 'error' for 1 single signal: to compute the error at matrix level I'd average all the single errors on q. Current pseudo code (doesn't work) for this logic would be:

   % outer for cycle sets flag = 0;

      for i=1:iter
          equi = all(abs(Xhat-X0) < 1.0e-6);
          if (equi == 1)
              flag = flag+1;
          else
              % suppose to have q signals in the input ensemble, i.e. vector equi has q elements, whose value can be 1 or 0
              flag = flag+(sum(equi)/q);
          end 
        end
      % j is index of outer for cycle, indicates measurements
    % vector res (to plot) shall tell the probability of successful recovery or the MSE/RMSE for each number of measurements
      res(j) = flag/iter;

Don't know how to use the 'mean' function here, I found 'immse' function but results are not that good.

Finally, I am not clear on your statement: "depending on how sparse your signals are, your algorithm could be more space-efficient by using a sparse matrix rather than a normal matrix". Seems relevant, could you please re-explain that? maybe with a code example? Sparsity k shall range around 3% / 8% of n (length of each signal in the ensemble).

I would most appreciate any additional help: many thanks!