Решить огромную разреженную линейную систему, где A является результатом произведения крон - PullRequest
Новая
       39

Решить огромную разреженную линейную систему, где A является результатом произведения крон

2 голосов
/ 28 мая 2019

Как решить линейную систему (A ⊗ B + C ⊗ D) x = b в MATLAB без вычисления фактической матрицы, которая умножается на x (⊗ обозначает произведение Кронекера). Хотя A, B, C и D являются sparse матрицами, наивный подход, 

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

не помещается в памяти и вызывает сбой MATLAB для больших матриц (~ 1000 x 1000 элементов на матрицу).

Что можно с этим сделать?

1 Ответ

2 голосов
/ 28 мая 2019

Видя, что ваши матрицы обычно довольно разреженные, конечный результат тензорного произведения не должен занимать так много памяти.Это один из случаев, когда векторизация просто не может быть выполнена из-за огромных требований к памяти промежуточных вычислений, но это просто возможно при использовании циклов (и частичной векторизации).

Обратите внимание, что это " лучше, чем ничего, но не намного"тип решения.

Я буду использовать ndSparse представление , так как с ним легче работатьразреженные матрицы. 

% 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

Так что это всего лишь теоретическая демонстрация того, как в конечном итоге выполнитьвычисления, но, к сожалению, они очень неэффективны, так как они должны вызывать методы класса снова и снова, а также не гарантируют, что будет достаточно памяти для вычисления оператора \.

Одна из возможныхспособ улучшить это, вызывая matlab.internal.sparse.kronSparse некоторым блочным способом и сохранять промежуточные результаты в правильной позиции выходного массива, но это шоДля этого нужно немного подумать.

Кстати, я пытался использовать представление FEX, которое упомянул Андер (KronProd), но это не дает никакой выгоды, когда вам нужно вычислить kron(A,B) + kron(C,D) (хотя этоудивительно для kron(A,B)\b ситуаций).

...