A fast cell array generation

Dear all,

in our optimization procedure, we need to create a cell array 'v_indices' from a vector 'v' and another cell array of indices 'indices'. Here is my code:

repets=4000;       %repetitions
n=1000;           %vector length

%defines a cell of indices
indices = mat2cell([0:n-1; 1:n; 2:n+1]',ones(n,1),3);
indices{1}=[1 2];  %correction of the first entry
indices{n}=[n-1 n]; %correction of the last entry

tic
for r=1:repets
    %defines a random vector
    v = randi(n*10,n,1);
    v_indices=indices;          %preallocation
    for i=1:numel(indices)   %loop over cell (inefficient?)
        v_indices{i}=v(indices{i}); 
    end
end
toc

The parameter 'repets' can be set higher in tests. Perhaps a for loop generating 'v_indices' is not very efficient there. Can you suggest some improvement, if possible?

Thank you, Jan Valdman

[Bruno Luong]

CELL and INDEXING is inherently slow (not the for-loop which is fast). You won't have a faster code structured like this.

If you explain what you would do with v_indices mighte be there is a fast way to carry out the calculation without breaking up sliding 3-window in cell (such as conv, im2col, blockproc, etc...).

[Jan Valdman]

Thank you for your comment. My aim is to compute an approximate gradient vector of an energy functional with variable 'v'. Since the gradient can be computed locally (due to the specific problem structure and underlyig finite element discretization), there is a cell array structure behind. The attached code shows the implementation for 1D case of the p-Laplacian operator. In fact, the cell array is not needed in 1D, but we would like to keep the cell array for an extension to 2D, where it can not be avoided at all. Therefore, this is probably a right benchmark for it. Here is the code:

function test_energy_gradient
repets=2000;       %repetitions
n=1000;            %vector length
p=2;               %power>1
epsilon=1e-5;      %central difference tolerance

h=1/numel(n);            %mesh size for energy computation

%defines a constant cell of indices
indices = mat2cell([0:n-1; 1:n; 2:n+1]',ones(n,1),3);
indices{1}=[1 2];  %correction of the first entry
indices{n}=[n-1 n]; %correction of the last entry

tic
v_indices_p=indices;          %preallocation
v_indices_m=indices;          %preallocation
for r=1:repets  
    v = randi(n*10,n,1);      %defines a random vector
   
    for i=1:numel(indices)   %loop over cell (inefficient?)
        v_loc_p=v(indices{i}); 
        v_loc_m=v(indices{i}); 
        if i==1
           v_loc_p(1)=v_loc_p(1)+epsilon;
           v_loc_m(1)=v_loc_m(1)-epsilon;
        else
           v_loc_p(2)=v_loc_p(2)+epsilon;
           v_loc_m(1)=v_loc_m(1)-epsilon;
        end      
        v_indices_p{i}=v(indices{i}); 
        v_indices_m{i}=v(indices{i}); 
    end
    
    e_local_p=local_energies(v_indices_p);  %energy plus
    e_local_m=local_energies(v_indices_m);  %energy minus 
    
    De=(e_local_p-e_local_m)/epsilon/2;       %output gradient vector 
end
toc

    function e_local = local_energies(v_indices)  
        e_local = zeros(numel(v_indices),1);  
                
        for j=1:numel(v_indices)          %loop over gradient components   
            v_loc_p = v_indices{j};       %local values of y values 
            Dv_local = diff(v_loc_p)./h;  %local gradients
            e_local(j) = (1/p)*h*sum(abs(Dv_local).^p);           
        end           
    end
end

It takes about 10 seconds on my desktop and I am wondering, whether there is a chance for some speedup. The key feature of our method (research article) is an efficient local evaluation of the gradient vector, which then becomes a part of an optimization routine.

Best, Jan Valdman

[Bruno Luong]

I'll leave out the detail of plus/minus potential and let you sort out, but instead of cell you can form an n x 3 array and work on this array. (NOTE: Your original code looks odd to me if not buggy).

The idea is this code:

repets=2000;       %repetitions
n=1000;            %vector length
p=2;
h=1/numel(n); % to me it's not correct
for r=1:repets  
    v = randi(n*10,n,1);      %defines a random vector
    
    vpad = v([1 1:end end]);
    v_loc = [vpad(1:end-2) vpad(2:end-1) vpad(3:end)];
    Dv_local = diff(v_loc,1,2) / h;
    e_local = (1/p)*h*sum(Dv_local.^p,2); % h would cancel out, p norm require power of (1/p), not devided by p

end

[Bruno Luong]

"we would like to keep the cell array for an extension to 2D, where it can not be avoided at all."

Usually the right data structure for (mesh) connectivity topology if double array vertex/face. One does not need CELL.

Avoid calculation with CELL when you can if speed matter for you.

[Jan Valdman]

Dear Bruno, thank you for your tips, I am thinking about them.

In 2D, we end up with a local vector with has a variable size (depending on number of elements a node is shared by), so it is for instance someting between 3 and 'na', where 'na' is the number of adjacent nodes, for simplicity let us take na=20 for a nice geometry. In 1D, it is, as you are suggesting, bounded by 3 (it is in fact 2 for boundary nodes and 3 for internal nodes). So we have to figure out how to address and process these local contributions efficiently. I am thinking of precomputing the value 'na' (it follows from the mesh topology) and perhaps putting zero value at positions not used. Let us see.

[Bruno Luong]

And when you connect all the adjadcent nodes do you get a constant-vertex polygonals such as triangles, rectangles, pentagonals, heaxagonals ? If it the case the best is vertexe/face number then you can "loop" on face to compute the local enegy operator.

If not there is another structure more generic is graph/digraph. Not sure about the speed though.

Under the hood of graph there is a sparse adjacent matrix. It can be of class logical if you are only interested in adjacency relation ship and you can work with such data structure as well. For example you can imgaine you build a sparse matrix such that when you multiply by the potential it returns the gradient field. Then you can compute the local p-norm, etc...

Bottom line is use two/three bigs arrays to store the topology. Avoid cell.

Answer 1

[Walter Roberson]

You hae v_indices=indices inside your for r loop, so nothing assigned to v_indices is kept for the next iteration of r. You do not assign to indices within the loop, so running multiple times is not making iterative changes to indices.

Therefore, your final result in v_indices is going to be the same as if you had only done one (the last) iteration of the for r loop, so you might as well not have a for r loop there.

[Jan Valdman]

Thank you, you are right: the cell array 'v_indices' is recomputed inside the loop. The code is only a simplified benchmark, there is an extra operation missing on 'v_indices' and needed - it is an energy computation expained in my answer to Bruno above.