求解巨大的稀疏线性系统，其中 A 是 kron 乘积的结果

Question

如何在MATLAB中求解线性系统(A ⊗ B + C ⊗ D) x = b而不计算乘以x的实际矩阵（⊗表示kronecker乘积）。即使 A、B、C 和 D 是 sparse 矩阵，天真的方法

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

不适合内存并导致 MATLAB 因大型矩阵（每个矩阵约 1000 x 1000 个元素）而崩溃。

对此可以做些什么？

Answer 1

看到你的矩阵通常很稀疏，张量积的最终结果应该不会占用那么多内存。这是由于中间计算的巨大内存需求而根本无法完成矢量化的情况之一，但使用循环（和部分矢量化）可能是可能的。

请注意，这是一种“总比没有好，但好不了多少”的解决方案。

我将使用 ndSparse submission，因为它可以更轻松地处理稀疏矩阵。

% 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

从这里可以根据可用内存继续：

% 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

所以这只是一个如何最终执行计算的理论演示，但不幸的是它非常低效，因为它必须一遍又一遍地调用class方法再次，它也不能保证有足够的内存来评估 \ 运算符。

一种可能的改进方法是以块方式调用 matlab.internal.sparse.kronSparse 并将中间结果存储在输出数组的正确位置，但这需要仔细考虑。

顺便说一句，我尝试使用 Ander 提到的 FEX 提交 (KronProd)，但是当你需要计算 kron(A,B) + kron(C,D) 时它没有任何好处（尽管它在 kron(A,B)\b 情况下很了不起） .

求解巨大的稀疏线性系统，其中 A 是 kron 乘积的结果

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

matlab

lazy-evaluation

sparse-matrix

equation-solving

tensor