Condensing indices in a matrix

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Condensing indices in a matrix

BGreen
Is it possible to construct a matrix of the following form without for loops?

Start with a multidimensional array
M( T1 , n1 , m1 , T2 , n2 , m2 )
where we have the same size in the T1 and T2 dimensions, the n1 and n2 dimensions, and the m1 and m2 dimensions.

Then, for each n and m, set
M( 1 , n , m , 2 , n , m ) = a;
M( 2 , n , m , 1 , n+1 , m ) = b;
M( 2 , n , m , 1 , n , m+1 ) = c;
for some numbers a, b and c. I can do this with e.g.
M( 1 , : , : , 2 , : , : ) = a;
so this is not a problem.

Then condense the multidimensional array to a 2D array
M( u(T1,n1,m1) , u(T2,n2,m2) )
where u is some function that condenses the indices from tuples of 3 to single indices. The function u can be anything that condenses our indices this way, though it'd be especially nice if I could define it myself.

I am familiar with using reshape to condense all indices. Is there also a way to condense only certain indices? Or, is there another way to construct this matrix besides reshaping a multidimensional array?

- Brett Green


Reply | Threaded
Open this post in threaded view
|

Re: Condensing indices in a matrix

siko1056
On 12/8/20 6:10 AM, Brett Green wrote:

> Is it possible to construct a matrix of the following form without for
> loops?
>
> Start with a multidimensional array
> M( T1 , n1 , m1 , T2 , n2 , m2 )
> where we have the same size in the T1 and T2 dimensions, the n1 and n2
> dimensions, and the m1 and m2 dimensions.
>
> Then, for each n and m, set
> M( 1 , n , m , 2 , n , m ) = a;
> M( 2 , n , m , 1 , n+1 , m ) = b;
> M( 2 , n , m , 1 , n , m+1 ) = c;
> for some numbers a, b and c. I can do this with e.g.
> M( 1 , : , : , 2 , : , : ) = a;
> so this is not a problem.
>
> Then condense the multidimensional array to a 2D array
> M( u(T1,n1,m1) , u(T2,n2,m2) )
> where u is some function that condenses the indices from tuples of 3 to
> single indices. The function u can be anything that condenses our
> indices this way, though it'd be especially nice if I could define it
> myself.
>
> I am familiar with using reshape to condense all indices. Is there also
> a way to condense only certain indices? Or, is there another way to
> construct this matrix besides reshaping a multidimensional array?
>
> - Brett Green
>


Was it possible to send some minimal working example how array M is
created?  For u(x,y,z), maybe `sub2ind` [1] can be helpful.

HTH,
Kai

[1] https://octave.org/doc/v6.1.0/Advanced-Indexing.html



Reply | Threaded
Open this post in threaded view
|

Re: Condensing indices in a matrix

BGreen

On Tue, Dec 8, 2020 at 1:19 AM Kai Torben Ohlhus <[hidden email]> wrote:

Was it possible to send some minimal working example how array M is
created?  For u(x,y,z), maybe `sub2ind` [1] can be helpful.

HTH,
Kai

[1] https://octave.org/doc/v6.1.0/Advanced-Indexing.html


Thank you and please excuse the slow response. Here is an MWE:

ind = @(T,n,m) 2*N*(m-1) + 2*(n-1) + T + 1;
M = zeros(2*N^2, 2*N^2);

for n=1:N
for m=1:N
M( ind(0,n,m) , ind(1,n,m) ) = 1;
M( ind(1,n,m) , ind(0,n,m) ) = 1;

M( ind(0,n,m) , ind(1,n+1,m) ) = 1;
M( ind(1,n+1,m) , ind(0,n,m) ) = 1;

M( ind(0,n,m) , ind(1,n,m+1) ) = 1;
M( ind(1,n,m+1) , ind(0,n,m) ) = 1;
end
end

While it's nice that sub2ind exists, I prefer to do the indexing myself so that it is explicitly readable from the code, and so that I can control it directly. The key problem for me is the assignment, which I would like to do without for loops if possible.


Reply | Threaded
Open this post in threaded view
|

Re: Condensing indices in a matrix

siko1056
On 12/21/20 4:22 AM, Brett Green wrote:

>
> On Tue, Dec 8, 2020 at 1:19 AM Kai Torben Ohlhus <[hidden email]
> <mailto:[hidden email]>> wrote:
>
>
>     Was it possible to send some minimal working example how array M is
>     created?  For u(x,y,z), maybe `sub2ind` [1] can be helpful.
>
>     HTH,
>     Kai
>
>     [1] https://octave.org/doc/v6.1.0/Advanced-Indexing.html
>     <https://octave.org/doc/v6.1.0/Advanced-Indexing.html>
>
>
> Thank you and please excuse the slow response. Here is an MWE:
>
> ind = @(T,n,m) 2*N*(m-1) + 2*(n-1) + T + 1;
> M = zeros(2*N^2, 2*N^2);
>
> for n=1:N
> for m=1:N
> M( ind(0,n,m) , ind(1,n,m) ) = 1;
> M( ind(1,n,m) , ind(0,n,m) ) = 1;
>
> M( ind(0,n,m) , ind(1,n+1,m) ) = 1;
> M( ind(1,n+1,m) , ind(0,n,m) ) = 1;
>
> M( ind(0,n,m) , ind(1,n,m+1) ) = 1;
> M( ind(1,n,m+1) , ind(0,n,m) ) = 1;
> end
> end
>
> While it's nice that sub2ind exists, I prefer to do the indexing myself
> so that it is explicitly readable from the code, and so that I can
> control it directly. The key problem for me is the assignment, which I
> would like to do without for loops if possible.


It has been a longer project, but I hope it was worth the outcome.
Looking at your matrices M with the "spy" command, you can clearly see a
sparse band structure and create them more efficiently with my function
"createMfast" below.

Your function needs for dimension N = 50 (5100 x 5100 dense matrix)
about 4 seconds and about 200 MB on my machine.  From N=100 the
computational cost become too huge.

The function "createMfast" uses sparse banded and boolean matrices (as
your values only deal with {0,1} and for N=500 (501000 x 501000 sparse
banded matrix) only 0.1 seconds creation time are required and it uses
only 17 MB of main memory.  Thus you can grow larger if you want.

The following experiment shows a comparison up to dimension 50 and a
small timing benchmark.

````````````````````````````
clear all; clc;

function M = createM (N)
   ind = @(T,n,m) 2*N*(m-1) + 2*(n-1) + T + 1;
   M = zeros(2*N^2, 2*N^2);

   for n=1:N
     for m=1:N
       M( ind(0,n,m) , ind(1,n,m) ) = 1;
       M( ind(1,n,m) , ind(0,n,m) ) = 1;

       M( ind(0,n,m) , ind(1,n+1,m) ) = 1;
       M( ind(1,n+1,m) , ind(0,n,m) ) = 1;

       M( ind(0,n,m) , ind(1,n,m+1) ) = 1;
       M( ind(1,n,m+1) , ind(0,n,m) ) = 1;
     end
   end
endfunction

function M = createMfast (N)
   d = false (2*N^2, 1);
   d(1:2:end) = true;
   dd = false (2*(N^2 + N), 6);
   dd(1:length(d),1:3) = [d d d];
   dd(2:length(d)+1,4) = d;
   dd(4:length(d)+3,5) = d;
   dd(end-length(d)+2:end,6) = d(1:end-1);
   M = spdiags (dd, [-2*N-1 -3 -1 1 3 2*N+1], 2*(N^2 + N), 2*(N^2 + N));
endfunction


NN = [1:10, 15, 20, 30, 40, 50, 100, 200, 500];
T1 = [];
T2 = [];
for N = NN
   if (N <= 50)
     tic ();
     M = createM (N);
     T1(end+1) = toc ();
   else
     T1(end+1) = inf;
   end

   ## Display the matrix structure.
   #if (N == 10)
   #  figure ();
   #  spy (M);
   #end
   ## Get diagonals for inspection.
   #[a, b] = spdiags (M)

   tic ();
   MM = createMfast (N);
   T2(end+1) = toc ();

   if (N <= 50)
     Mdiff = MM - sparse (M);
     if (any (any (Mdiff)))
       figure ();
       spy (Mdiff);
     end
   end
end

figure ();
plot (NN, T1, 'b', NN, T2, 'r')
````````````````````````````

HTH,
Kai


Reply | Threaded
Open this post in threaded view
|

Re: Condensing indices in a matrix

BGreen
On Tue, Dec 22, 2020 at 12:49 AM Kai Torben Ohlhus <[hidden email]> wrote:

It has been a longer project, but I hope it was worth the outcome.
Looking at your matrices M with the "spy" command, you can clearly see a
sparse band structure and create them more efficiently with my function
"createMfast" below.

Your function needs for dimension N = 50 (5100 x 5100 dense matrix)
about 4 seconds and about 200 MB on my machine.  From N=100 the
computational cost become too huge.

The function "createMfast" uses sparse banded and boolean matrices (as
your values only deal with {0,1} and for N=500 (501000 x 501000 sparse
banded matrix) only 0.1 seconds creation time are required and it uses
only 17 MB of main memory.  Thus you can grow larger if you want.

The following experiment shows a comparison up to dimension 50 and a
small timing benchmark..

HTH,
Kai

Thank you very much! spdiags is the key I'd been missing. spy will also be very useful as I learn to use spdiags and troubleshoot.
... and to say that the runtime comparison was impressive is an understatement.