help-octave
[Top][All Lists]

## Re: Condensing indices in a matrix

 From: Kai Torben Ohlhus Subject: Re: Condensing indices in a matrix Date: Tue, 22 Dec 2020 14:49:52 +0900 User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:78.0) Gecko/20100101 Thunderbird/78.5.0

```On 12/21/20 4:22 AM, Brett Green wrote:
```
```
```
On Tue, Dec 8, 2020 at 1:19 AM Kai Torben Ohlhus <k.ohlhus@gmail.com <mailto:k.ohlhus@gmail.com>> 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

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

```