Create 3D matrix with 1s in certain position and all other elements 0s (NO FOR LOOP)

Hello,

I'm trying to compute a 3D matrix with all 0s except for columns of 1s defined by random position vector. I want this solution WITHOUT FOR LOOPS or any loops.

So I have a 0s matrix of dimension 3x5x5

and a random position vector of scalars 1x5 where I put the index of columns that has to be 1s.

So for example if this ranom vector is equal to

% The first number 5 is only for semplicity, it depends on the size of my system that
% can reach huge dimension (above 10e3)
site = randi(5,1,5);

I want to have a 3D matrix (called A) in this form

A(:,:,1) =
1    0    0    0    0
1    0    0    0    0
1    0    0    0    0
A(:,:,2) = 
0    1    0    0    0
0    1    0    0    0
0    1    0    0    0
A(:,:,3) = 
0    0    1    0    0
0    0    1    0    0
0    0    1    0    0
A(:,:,4) = 
0    0    0    1    0
0    0    0    1    0
0    0    0    1    0
A(:,:,5) = 
0    0    0    0    1
0    0    0    0    1
0    0    0    0    1

[Jan]

Why doyou want to avoid loops? They might be the most efficient solution.

[Pietro Marabotti]

Because I have to do 10^8 iterations and for loops are very very slow.

[Jan]

Please be a specific as possible: Which dimension is 1e8? Perhaps it is only a problem of your code, so please post the slow loops.

Does the contents of site concerns the 2nd or the 3rd or both dimensions? [1,2,3,4,5] might be a misleading example.

Are you sure, that you really need such a huge array, which contain the same information as your vector site?

[Pietro Marabotti]

The number of iterations gives the 3rd dimension. I have modified the answer and write site in the proper way. I need this such big array beacuase I'm computing a simulation and I need sucha a huge number of iterations. This is only a part of a code that I cannot unfortunately publish.

[Jan]

@Pietro: Now the inputs are less clear than before. What is the meaning of the elements of site?

Is the size of the problem 1e8 or 1e3 - and what does this number define? One of the dimensions of A or the 2nd and 3rd one? Then A(3, 1e8, 1e8) will exhaust your RAM massively and the loop is not the problem: 80'000 TerraByte of RAM are less than you find in modern computers.

n = 1e3;
site = randi(n, 1, n);
tic;
A = zeros(3,n,n);
for k = 1:n
    A(:,site(k),site(k)) = 1;
end
toc

This takes about 0.01 seconds. With n=1e4 it needs 0.06 seconds. So what exactly is your problem? You explain, that "the

The need of huge arrays, which contain a few 1s only at well defined positions is a secure mark of a design flaw of the program. It would be much easier to use site directly. Even the need for a huge number of iterations does not explain, why you need to create this almost empty input matrix statically.

Answer 1

[Luna]

Hi Pietro

A = zeros(3,5,5);
site = [1,2,3,4,5]';
for i = 1:numel(site)
A(:,site(i),site(i)) = ones(1,3)';
end

This might solve your problem but you need to be sure that if any element of site vector is bigger than any dimensions of A, it will change A matrix sizes. For example site cannot be [1,2,3,4,30].

I don't understand why you don't want for loop.

[Pietro Marabotti]

Because I have to do 10^8 iterations and for loops are very very slow.

[Jan]

@Luna: ones(3,1) is smarter than ones(1,3)', because it avoids the transposition. But here you can define the right hand side by a scalar also, which is even better:

A(:,site(i),site(i)) = 1;

[Pietro Marabotti]

Yes, I see, but I have to repeat this code many many times so the foor loop take me much more time than operations between matrices.

[Luna]

Thanks Jan! I am a big fan of you! :)