Kronecker Tensor Product for very large matrices

I need to build a 2D kronecker tensor product of with a size of approximately 10^6 x 10^6. The problem is that I get the following message:

%Out of memory. Type HELP MEMORY for your options.

Is there any way to get around this memory limitation? My code is posted below:

nx = 1000;ny = 1000;                                                    %Number of points along x and y
dx = 1/(nx-1);dy = 1/(ny-1);                                            %Spacing between each point
diag_block = eye(ny-1) * (-2/dx^2-2/dy^2);                              %Identity Matrix of ny-1 points                           
diag_block = diag_block + diag(ones(ny-2,1)/dy^2,1);                    %
diag_block = sparse(diag_block) + sparse(diag(ones(ny-2,1)/dy^2,-1));   %
Matrix = kron(sparse(eye(nx-1)),sparse(diag_block));                    %Using sparse() made this line work
Matrix = Matrix + diag(ones((nx-2) * (ny-1),1)/dx^2, ny-1);             %FAILS HERE!!!
Matrix = Matrix + diag(ones((nx-2) * (ny-1),1)/dx^2, -(ny-1));          %AND HERE!!!

Answer 1

[Matt J]

You need the matrices generated by diag() in sparse form, too. Use spdiags() instead, e.g.,

nx = 3;ny =3;     dx=1;                                               %Number of points along x and y


Qfull=diag(ones((nx-2) * (ny-1),1)/dx^2, ny-1)

Qfull = 4×4

0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0

Bin=[0;0;ones((nx-2) * (ny-1),1)/dx^2];
Q=spdiags(Bin,ny-1,4,4)

Q = 
   (1,3)        1
   (2,4)        1

full(Q)

ans = 4×4

0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0

[William Kett]

I apologize for taking so long to get back to you! A lot of things came up since your reply. I went ahead and tried implementing spdiags() and it did drastically reduce the size of matrix that I am trying to add onto the original matrix. The only issue is that it seems to revert back to its full size when I try to add said sparse matrix to the original matrix. The full code will be posted below with comments to highlight where changes were made that implement your code.

nx = 1000;ny = 1000;
dx = 1/(nx-1);dy = 1/(ny-1);
diag_block = eye(ny-1) * (-2/dx^2-2/dy^2);
diag_block = diag_block + diag(ones(ny-2,1)/dy^2,1);
diag_block = diag_block + diag(ones(ny-2,1)/dy^2,-1);
%% Issues occur below
Matrix = kron(sparse(eye(nx-1)),sparse(diag_block));
Matrix = Matrix + spdiags(ones((nx-2) * (ny-1),1)/dx^2, ny-1); %Implemented spdiags here!!! Also breaks here!!!
Matrix = Matrix + spdiags(ones((nx-2) * (ny-1),1)/dx^2, -(ny-1)); %Implemented spdiags here!!! Also breaks here!!!

I went ahead and broke up the 2 lines where the code implemented spdiags as well as using whos to show just how much spdiags() cuts down on the size of the matrices below:

%% First line
Matrix = kron(sparse(eye(nx-1)),sparse(diag_block)); % The Matrix we are trying to add onto.
whos Matrix %yields a size of roughly 56 Megabytes

%% Second line broken up
% Original line: Matrix = Matrix + spdiag(ones((nx-2) * (ny-1),1)/dx^2, ny-1);
A = spdiags(ones((nx-2) * (ny-1),1)/dx^2, ny-1); %Works
whos A %Yields a size of only 8 bytes
Matrix = Matrix + A;
%Error using  + 
%Requested 998001x998001 (7420.8GB) array exceeds maximum array size preference. Creation of arrays greater
%than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or
%preference panel for more information.

%% Third line broken up
B = spdiags(ones((nx-2) * (ny-1),1)/dx^2, -(ny-1)); %Works
whos B %Yields a size of only 8 bytes
Matrix = Matrix + B;
%Error using  + 
%Requested 998001x998001 (7420.8GB) array exceeds maximum array size preference. Creation of arrays greater
%than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or
%preference panel for more information.

So, we should be taking a 56 Megabyte matrix and adding an 8 Byte array onto it, but it seems as though the "+" operation reverts it back to its massive 7420 Gigabyte size. Is there a way to add the matrices together so that they retain their sparsity?

[Matt J]

Your A and B matrices are full scalars, not sparse matrices:

>> A

A =

     0

>> B

B =

      998001

[William Kett]

Ah! So when I create A and B using spdiags(), I am instead getting a scalar instead of getting a matrix and by adding a scalar to the Matrix, I guess it fills in the empty space and creates a non-sparse matrix? I did not change my arguments or the order in which I passed them in from when I changed from diag() to spdiags().

Could this be why?

Changing nx = ny = 4 yields the following Matrix:

Matrix =

-36 9 0 9 0 0 0 0 0

9 -36 9 0 9 0 0 0 0

0 9 -36 0 0 9 0 0 0

9 0 0 -36 9 0 9 0 0

0 9 0 9 -36 9 0 9 0

0 0 9 0 9 -36 0 0 9

0 0 0 9 0 0 -36 9 0

0 0 0 0 9 0 9 -36 9

0 0 0 0 0 9 0 9 -36

and converting that to sparse by running sparse(Matrix) yields:

ans =

(1,1) -36

(2,1) 9

(4,1) 9

(1,2) 9

(2,2) -36

(3,2) 9

(5,2) 9

(2,3) 9

(3,3) -36

(6,3) 9

(1,4) 9

(4,4) -36

(5,4) 9

(7,4) 9

(2,5) 9

(4,5) 9

(5,5) -36

(6,5) 9

(8,5) 9

(3,6) 9

(5,6) 9

(6,6) -36

(9,6) 9

(4,7) 9

(7,7) -36

(8,7) 9

(5,8) 9

(7,8) 9

(8,8) -36

(9,8) 9

(6,9) 9

(8,9) 9

(9,9) -36

So, I know that I am supposed to be getting a non-scalar result for my A and B, so it's hopefully the case that I just input my arguments incorrectly.

[Matt J]

Could this be why?

Yes, you should review the documentation for spdiags() and/or re-examine the example I gave. you It's syntax is very different from that of diag().