Solve huge sparse linear system, where A is a result of a kron product

Question

How can I solve the linear system (A ⊗ B + C ⊗ D) x = b in MATLAB without calculating the actual matrix that multiplies x (⊗ denotes the kronecker product). Even though A,B,C and D are sparse matrices, the naive approach,

x = (kron(A,B) + kron(C,D))\b

does not fit in memory and causes MATLAB to crash for large Matrices (~1000 x 1000 elements per matrix).

What can be done about this?

Do you know of any algorithms out there to solve this sort of problems, and are looking for their implementation in MATLAB? Or are you asking how, in general, one can solve in MATLAB problems that don't fit in memory? — Dev-iL, May 28 '19 at 10:48
How to solve problems that don't fit in memory. Or how the create a "matrix object" that stores the result of (A x B + C x D) without explicity calculating it in the first place which can be used with the backslash operator, if this is possible — Julian H, May 28 '19 at 10:55
[_Tall arrays_](https://www.mathworks.com/help/matlab/import_export/tall-arrays.html) provide a general solution for such problems. However, if this equations represents some known problem, perhaps there are some efficient implementations out there, which take `A, B, C, D` and `b` and return `x` w/o solving the full system. You might want to mention which field this equation appears in.... — Dev-iL, May 28 '19 at 11:15
The equations arise the space time finite element method where B and D are mass and stiff matrices and and A and C are matrices for the discretication in time. The latter two are sparse diagonal matrices with only one or two secondary diagonals. One approche to solve such equations is presenten in [link](https://www.igpm.rwth-aachen.de/Download/reports/pdf/IGPM309.pdf) but I'm searching for an alternative — Julian H, May 28 '19 at 11:21
According to Mathworks, they do not have a sparse `kron`. https://uk.mathworks.com/matlabcentral/answers/111918-is-it-possible-to-improve-performance-of-the-kron-operator-by-making-it-a-built-in-function-rather-t . In the same thread there is a link to a user-implemented sparse kron, I suggest using that one — Ander Biguri, May 28 '19 at 12:11

score 2 · Accepted Answer · answered May 28 '19 at 14:51

Seeing how your matrices are generally quite sparse, the end result of the tensor product shouldn't take that much memory. This is one of the cases when vectorization simply cannot be done due to the huge memory requirements of the intermediate computations, but it just might be possible using loops (and partial vectorization).

Note that this is a "better than nothing, but not by much" type of a solution.

I'll be using the ndSparse submission, as it makes it easier to work with sparse matrices.

% Create some matrices
[A,B] = deal(sparse(1000,1000));
A(randperm(1000^2, 10000)) = randn(1E4, 1);
B(randperm(1000^2, 10000)) = randn(1E4, 1);
A = ndSparse(A); B = ndSparse(B);

% Kronecker tensor product, copied from kron.m
[ma,na] = size(A);
[mb,nb] = size(B);
A = reshape(A,[1 ma 1 na]);
B = reshape(B,[mb 1 nb 1]);
% K = reshape(A.*B,[ma*mb na*nb]);                    % This fails, too much memory.
K = ndSparse(sparse(mb*ma*nb*na,1),[mb, ma, nb, na]); % This works

From here one can proceed depending on the available memory:

% If there's plenty of memory (2D x 1D -> 3D):
for ind1 = 1:mb
  K(ind1,:,:,:) = bsxfun(@times, A, B(ind1, :, :, :));
end

% If there's less memory (1D x 1D -> 2D):
for ind1 = 1:mb
  for ind2 = 1:ma
    K(ind1,ind2,:,:) = bsxfun(@times, A(:, ind2, :, :), B(ind1, :, :, :));
  end
end

% If there's even less memory (1D x 0D -> 1D):
for ind1 = 1:mb
  for ind2 = 1:ma
    for ind3 = 1:nb
      K(ind1,ind2,ind3,:) = bsxfun(@times, A(:, ind2, :, :), B(ind1, :, ind3, :));
    end
  end
end

% If there's absolutely no memory (0D x 0D  -> 0D):
for ind1 = 1:mb
  for ind2 = 1:ma
    for ind3 = 1:nb
      for ind4 = 1:na
        K(ind1,ind2,ind3,ind4) = A(:, ind2, :, ind4) .* B(ind1, :, ind3, :);
      end
    end
  end
end

K = sparse(reshape(K,[ma*mb na*nb])); % Final touch

So this is just a theoretical demonstration of how to eventually perform the computation, but unfortunately it is very inefficient since it has to call class methods over and over again, and also it doesn't guarantee that there will be enough memory to evaluate the \ operator.

One possible way to improve this is by calling matlab.internal.sparse.kronSparse in some blockwise manner and store the intermediate results in the correct position of the output array, but this should require some careful thought.

BTW, I tried using the FEX submission that Ander mentioned (KronProd) but it provides no benefit when you need to compute kron(A,B) + kron(C,D) (although it's amazing for kron(A,B)\b situations).

Solve huge sparse linear system, where A is a result of a kron product

1 Answers1