Re: Condensing indices in a matrix

On 12/21/20 4:22 AM, Brett Green wrote:
> 
> On Tue, Dec 8, 2020 at 1:19 AM Kai Torben Ohlhus <[email protected] 
> <mailto:[email protected]>> 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