|
|
|
|
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?
|
|
|
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
|
|
|
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.
|
|
|
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
|
|
|
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.
|
|
Locked
