Лучшая практика при работе с разреженными матрицами

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

Мой вопрос двоякий:

В ниже, A = full(S), где S - разреженная матрица.

Каков "правильный" способ доступа к элементу в разреженной матрице?

То есть, каков будет разреженный эквивалент var = A(row, col)?

Мой взгляд на эту тему: Ты бы не сделал ничего другого. var = S(row, col) так же эффективен, как и он. Мне объяснили это с объяснением:

Доступ к элементу в строке 2 и столбце 2, как вы сказали, S (2,2) то же, что и добавление нового элемента: var = S(2,2) = > A = full(S) = > var = A(2,2) = > S = sparse(A) => 4.

Может ли это утверждение действительно быть правильным?

Каков "правильный" способ добавления элементов в разреженную матрицу?

То есть, каков будет разреженный эквивалент A(row, col) = var? (Предполагая A(row, col) == 0 для начала)

Известно, что простое выполнение A(row, col) = var медленное для больших разреженных матриц. Из документация:

Если вы хотите изменить значение в этой матрице, вы можете искушаться использовать ту же индексацию:

B (3,1) = 42; % Этот код работает, однако он медленный.

Мой взгляд на эту тему: При работе с разреженными матрицами вы часто начинаете с векторов и используете их для создания матрицы следующим образом: S = sparse(i,j,s,m,n). Конечно, вы могли бы создать его следующим образом: S = sparse(A) или sprand(m,n,density) или что-то подобное.

Если вы начнете первый путь, вы просто выполните:

i = [i; new_i];
j = [j; new_j];
s = [s; new_s];
S = sparse(i,j,s,m,n);

Если вы начали с отсутствия векторов, вы сделали бы то же самое, но сначала используйте find:

[i, j, s] = find(S);
i = [i; new_i];
j = [j; new_j];
s = [s; new_s];
S = sparse(i,j,s,m,n);

Теперь вы, конечно, будете иметь векторы и можете повторно использовать их, если вы делаете эту операцию несколько раз. Однако было бы лучше добавить все новые элементы сразу, а не делать это в цикле, потому что растущие векторы медленные. В этом случае new_i, new_j и new_s будут векторами, соответствующими новым элементам.

Ответ 1

EDIT: ответ изменен в соответствии с предложениями Олега (см. комментарии).

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

Ssize = 10000;
NumIterations = 100;
NumInsertions = round(logspace(0, 4, 10));
NumInitialNZ = round(logspace(1, 4, 4));

NumTests = numel(NumInsertions) * numel(NumInitialNZ);
TimeDirect = zeros(numel(NumInsertions), numel(NumInitialNZ));
TimeIndices = zeros(numel(NumInsertions), 1);

%% Single insertion operation (non-incremental)
% Method A: Direct insertion
for iInitialNZ = 1:numel(NumInitialNZ)
    disp(['Running with initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]);

    for iInsertions = 1:numel(NumInsertions)
        tSum = 0;
        for jj = 1:NumIterations
            S = spalloc(Ssize, Ssize, NumInitialNZ(iInitialNZ));
            r = randi(Ssize, NumInsertions(iInsertions), 1);
            c = randi(Ssize, NumInsertions(iInsertions), 1);

            tic
            S(r,c) = 1;
            tSum = tSum + toc;
        end

        disp([num2str(NumInsertions(iInsertions)) ' direct insertions: ' num2str(tSum) ' seconds']);
        TimeDirect(iInsertions, iInitialNZ) = tSum;
    end
end

% Method B: Rebuilding from index vectors
for iInsertions = 1:numel(NumInsertions)
    tSum = 0;
    for jj = 1:NumIterations
        i = []; j = []; s = [];
        r = randi(Ssize, NumInsertions(iInsertions), 1);
        c = randi(Ssize, NumInsertions(iInsertions), 1);
        s_ones = ones(NumInsertions(iInsertions), 1);

        tic
        i_new = [i; r];
        j_new = [j; c];
        s_new = [s; s_ones];
        S = sparse(i_new, j_new ,s_new , Ssize, Ssize);
        tSum = tSum + toc;
    end

    disp([num2str(NumInsertions(iInsertions)) ' indexed insertions: ' num2str(tSum) ' seconds']);
    TimeIndices(iInsertions) = tSum;
end

SingleOperation.TimeDirect = TimeDirect;
SingleOperation.TimeIndices = TimeIndices;

%% Incremental insertion
for iInitialNZ = 1:numel(NumInitialNZ)
    disp(['Running with initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]);

    % Method A: Direct insertion
    for iInsertions = 1:numel(NumInsertions)
        tSum = 0;
        for jj = 1:NumIterations
            S = spalloc(Ssize, Ssize, NumInitialNZ(iInitialNZ));
            r = randi(Ssize, NumInsertions(iInsertions), 1);
            c = randi(Ssize, NumInsertions(iInsertions), 1);

            tic
            for ii = 1:NumInsertions(iInsertions)
                S(r(ii),c(ii)) = 1;
            end
            tSum = tSum + toc;
        end

        disp([num2str(NumInsertions(iInsertions)) ' direct insertions: ' num2str(tSum) ' seconds']);
        TimeDirect(iInsertions, iInitialNZ) = tSum;
    end
end

% Method B: Rebuilding from index vectors
for iInsertions = 1:numel(NumInsertions)
    tSum = 0;
    for jj = 1:NumIterations
        i = []; j = []; s = [];
        r = randi(Ssize, NumInsertions(iInsertions), 1);
        c = randi(Ssize, NumInsertions(iInsertions), 1);

        tic
        for ii = 1:NumInsertions(iInsertions)
            i = [i; r(ii)];
            j = [j; c(ii)];
            s = [s; 1];
            S = sparse(i, j ,s , Ssize, Ssize);
        end
        tSum = tSum + toc;
    end

    disp([num2str(NumInsertions(iInsertions)) ' indexed insertions: ' num2str(tSum) ' seconds']);
    TimeIndices(iInsertions) = tSum;
end

IncremenalInsertion.TimeDirect = TimeDirect;
IncremenalInsertion.TimeIndices = TimeIndices;

%% Plot results
% Single insertion
figure;
loglog(NumInsertions, SingleOperation.TimeIndices);
cellLegend = {'Using index vectors'};
hold all;
for iInitialNZ = 1:numel(NumInitialNZ)
    loglog(NumInsertions, SingleOperation.TimeDirect(:, iInitialNZ));
    cellLegend = [cellLegend; {['Direct insertion, initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]}];
end
hold off;
title('Benchmark for single insertion operation');
xlabel('Number of insertions'); ylabel('Runtime for 100 operations [sec]');
legend(cellLegend, 'Location', 'NorthWest');
grid on;

% Incremental insertions
figure;
loglog(NumInsertions, IncremenalInsertion.TimeIndices);
cellLegend = {'Using index vectors'};
hold all;
for iInitialNZ = 1:numel(NumInitialNZ)
    loglog(NumInsertions, IncremenalInsertion.TimeDirect(:, iInitialNZ));
    cellLegend = [cellLegend; {['Direct insertion, initial nzmax = ' num2str(NumInitialNZ(iInitialNZ))]}];
end
hold off;
title('Benchmark for incremental insertions');
xlabel('Number of insertions'); ylabel('Runtime for 100 operations [sec]');
legend(cellLegend, 'Location', 'NorthWest');
grid on;

Я запустил это в MATLAB R2012a. Результаты для отдельных операций вставки суммируются на этом графике:

Это показывает, что использование прямой вставки намного медленнее, чем использование индексных векторов, если выполняется только одна операция. Рост в случае использования индексных векторов может быть либо из-за роста самих векторов, либо из более длинной разреженной матрицы, я не уверен, что. Начальный nzmax, используемый для построения матриц, по-видимому, не влияет на их рост.

Результаты для инкрементных вставок суммируются на этом графике:

Здесь мы видим противоположную тенденцию: использование индексных векторов происходит медленнее из-за накладных расходов их постепенного роста и восстановления разреженной матрицы на каждом шаге. Способ понять это - посмотреть на первую точку на предыдущем графике: для вставки одного элемента более эффективно использовать прямую вставку, а не перестраивать с использованием векторов индекса. В инкрементальном случае эта единственная вставка выполняется повторно, и поэтому становится жизнеспособным использовать прямую вставку, а не индексные векторы, против предложения MATLAB.

Это понимание также предполагает, что если бы мы постепенно добавляли, скажем, 100 элементов за раз, то эффективным выбором было бы использовать индексные векторы, а не прямую вставку, так как первый график показывает, что этот метод будет быстрее для вставки этот размер. Между этими двумя режимами есть область, где вы, вероятно, должны поэкспериментировать, чтобы определить, какой метод более эффективен, хотя, вероятно, результаты покажут, что разница между методами там невелика.

Нижняя строка: какой метод я должен использовать?

Мой вывод состоит в том, что это зависит от характера ваших предполагаемых операций вставки.

Если вы собираетесь вставлять элементы по одному, используйте прямую вставку.
Если вы намереваетесь вставлять большое ( > 10) количество элементов за раз, перестройте матрицу из индексных векторов.

Ответ 2

MATLAB хранит разреженные матрицы в сжатом формате столбца. Это означает, что когда вы выполняете операции типа A (2,2) (чтобы получить элемент в строке 2, столбец 2), MATLAB сначала получает доступ ко второму столбцу, а затем находит элемент в строке 2 (индексы строк в каждом столбце сохраняются в порядке возрастания). Вы можете думать об этом как:

 A2 = A(:,2);
 A2(2)

Если вы получаете доступ только к одному элементу разреженной матрицы, выполняющему var = S(r,c), это хорошо. Но если вы перебираете элементы разреженной матрицы, вы, вероятно, захотите получить доступ к одному столбцу за раз, а затем перебирать ненулевые индексы строк через [i,~,x]=find(S(:,c)). Или используйте что-то вроде spfun.

Вам следует избегать построения плотной матрицы A, а затем делать S = sparse(A), так как эти операции просто выжимают нули. Вместо этого, как вы заметили, гораздо эффективнее построить разреженную матрицу с нуля, используя триплетную форму и вызов sparse(i,j,x,m,n). MATLAB имеет приятную страницу, в которой описывается, как эффективно создавать разреженные матрицы.

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