Recovery set of signals same sparsity

6 views (last 30 days)
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.

Accepted Answer

Stephen Jue
Stephen Jue on 1 Sep 2016
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.
  1 Comment
VMat
VMat on 5 Sep 2016
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!

Sign in to comment.

More Answers (2)

Stephen Jue
Stephen Jue on 7 Sep 2016
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%.
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.

VMat
VMat on 1 Sep 2016
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!

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!