How to find the indices of non-zero elements in a matrix

Hello,

I have an Na-by-Nt matrix which is sparse, i.e.: most elements are zeros. I want to find the indices of the non-zeros elements in the form of (i,j) where i is the row and j is the column.

I have then these two Na-by-1 vector a and Nt-by-1 vector tau. i corresponds to the ith element in a and j the jth element in tau, and I want to find them as well.

How can I do the above in an efficient way without the need for for loops?

Thanks

Answer 1

[Cedric Wannaz]

If it is for storing only non-zero elements, MATLAB supports SPARSE matrices. If you need to get row/column indices of non-zero elements of A:

 [rId, cId] = find( A ) ;

and if you need values as well:

 [rId, cId, val] = find( A ) ;

If you just need values, you can get them simply with

 vals = A(A ~= 0) ;

[S. David]

Thanks, but in this way I have repetitive indices, because the same row/column may contain more than one non-zero element. I was looking for something like this:

1- Find the indices of non-zero elements as (i,j).

2- b=[b a(i)].

3- t=[t tau(j)].

But I am trying to avoid using for loops since Na*Nt is relatively large.

[Cedric Wannaz]

You will end up roughly re-implementing the internal implementation of sparse matrices if you do this, but for what purpose?

[S. David]

OK, what do you suggest? Now I have rId and cId. How can I extract the corresponding elements from a and tau without using for loops?

[Cedric Wannaz]

For data storage/manipulation of matrices which are sparse enough, I suggest that you use MATLAB built-in sparse matrices, and that you avoid building your own, less efficient and less versatile data structure. Most functions from the MATLAB base package do support sparse matrices and will run algorithms optimized for them when they are given sparse matrices as input.

I recommend that you build sparse matrices using calls of the form

 S = sparse(i,j,s,m,n) ;

and that you avoid building them dense first and then converting them to sparse with

 S = sparse( A )

At this point, I don't understand why you need indices if you want to get all non-zero values, as you can get them with

 val = A(A ~= 0) ;

If you have a vector of row indices and a vector of column indices (which are subsets of the sets of all indices), you can get corresponding values by converting the "subscript" indices to linear indices using SUB2IND, and indexing the array with the linear index.

[S. David]

Please have a look here

http://www.mathworks.com/matlabcentral/answers/138011-sparse-recovery-problem-solution

What I am doing is that after finding x which is Na*Nt-by-1 vector I form the Na-by-Nt matrix

xMtx=transpose(reshape(x,Nt,Na));

From this matrix I want to find Np=Na unique triplets. This implies that I should have one and only one non-zero element in each row of xMtx, and these non-zero elements occur at different column indices.

[Matt J]

This implies that I should have one and only one non-zero element in each row of xMtx, and these non-zero elements occur at different column indices.

But you haven't explained what changes to xMtx you are free to make in order to satisfy this property. Do you want to just discard non-zero elements, replacing them with zeros, until xMtx has the pattern of non-zeros that you want?

Incidentally, the property you describe is impossible if Na>Nt.

[S. David]

Na is always less than Nt.

This is the problem, I was expecting that I will have Np non-zero elements in the solution. Now, I don't know how to handle it. Is there any sparse algorithm that force the solution to have exactly Np non-zero elements?

[Cedric Wannaz]

Could you build a small numeric example?

[S. David]

Here is my MATLAB code for A and v in the sparse problem Ax=v:

clear all;
clc;

SNRdB=0:2:12;
SNR=10.^(SNRdB./10);

N=512;
f0=30*10^3;
BW=4*10^3;
df=BW/N;
T=1/df;
Ts=T/N;%Sample time
nsamp=100;%Oversampling
Ta=Ts/nsamp;%Sample time after oversampling

tau=[0 5 7 20].*10^-3;%Path Delays
Np=length(tau);
a=[0 0.001 0.0012 0.0024];
h=[1 0.8 0.7 0.5];
h=h./sqrt(h*h');

CP_ms=30*10^-3;%CP in ms
CP=round(CP_ms/Ta);%CP in samples 
Tg=CP*Ta;

t =( -Tg : Ta : (T+Tg-Ta) )';

fr=(f0+(0:N-1)./T)'*ones(1,N);%Subcarrier indices increase with rows
fc=transpose((f0+(0:N-1)./T)'*ones(1,N));%Subcarrier indices increase with column

H=zeros(N,Np);

for pp=1:Np
    H(:,pp)=H(:,pp)+h(pp).*exp( -1i * 2 * pi * (f0 + (0:N-1)/T)' * tau(pp) );
end

Phi=zeros(N);
for pp=1:Np
    Hc_Col_Incr=transpose(H(:,pp)*ones(1,N));
    Theta=(fc.*kron(1+a(pp),ones(N))-fr);
    Phi=Phi+Hc_Col_Incr.*exp(1i.*pi.*T.*Theta).*sinc(T.*Theta);
end

lambda=2;
Ntau=lambda*N*Tg/T;
NtauL=1:Ntau;
Tau_Tentative_Set=NtauL.*(T/(lambda*N));

aMax=max(a);%5*10^-4;
Doppler_Tentative_Set=linspace(0,aMax,Np);%-aMax:1*10^-4:aMax;
% Na=length(Doppler_Tentative_Set);

s=2.*(rand(N,1)>0.5)-1;

v=Phi*s;

for pp=1:Np
    for tt=1:Ntau
        theta_p=fc.*kron(1+Doppler_Tentative_Set(pp),ones(N))-fr;
        Gamma_p=exp(1i*pi*T.*theta_p).*sinc(T.*theta_p);
        Lambda_p=diag(exp(-1i*2*pi.*(f0+(0:N-1)/T).*Tau_Tentative_Set(tt)));
        if pp==1 && tt==1
            A=Gamma_p*Lambda_p*s;
        else
            A=[A Gamma_p*Lambda_p*s];
        end
    end
end

[xVec,inform] = as_bpdn(A,v);

The function as_bpdn is attached.

Thanks

[Cedric Wannaz]

I meant a simple/minimal numeric example showing which triplets you need to ID in which data structure, because I am not sure that anybody will have time to understand your code/algorithm.

If it is working but too slow, the first thing that comes to my mind is to replace the last double FOR loop with

 A = complex( zeros( numel(s), Np * Ntau )) ;
 colId = 0 ;
 for pp=1:Np
     theta_p=fc.*kron(1+Doppler_Tentative_Set(pp),ones(N))-fr;
     Gamma_p=exp(1i*pi*T.*theta_p).*sinc(T.*theta_p);
     for tt=1:Ntau
         Lambda_p_diag = exp(-1i*2*pi.*(f0+(0:N-1)/T).*Tau_Tentative_Set(tt));
         colId = colId + 1 ;
         A(:,colId) = bsxfun( @mtimes, Gamma_p, Lambda_p_diag ) * s ;
     end
 end

[S. David]

It is difficult to give a non real example. However, I think it is not necessary to understand the code/algorithm. At the end, the code will give A, v, and the sparse solution xVec. That is all is needed I guess.