Когда мне следует использовать `sparse`?

Я просматривал sparse документацию Matlab, пытаясь найти, есть ли какие-либо рекомендации относительно того, когда это имеет смысл использовать разреженное представление, а не полное представление.

Например, у меня есть матрица data с примерно 30% ненулевых записей. Я могу проверить используемую память.

whos data
  Name             Size                 Bytes  Class     Attributes

  data      84143929x11            4394073488  double    sparse    

data = full(data);
whos data
  Name             Size                 Bytes  Class     Attributes

  data      84143929x11            7404665752  double              

Здесь я явно экономлю память, но будет ли это справедливо для любой матрицы с 30% ненулевых записей? А как насчет 50% ненулевых записей? Есть ли эмпирическое правило, при каком проценте мне следует перейти на полную матрицу?

А как насчет вычислений? Умножение матриц на разреженную матрицу, как правило, медленнее или быстрее? Операции с разреженной матрицей говорит, что

Вычислительная сложность разреженных операций пропорциональна nnz, количеству ненулевых элементов в матрице. Вычислительная сложность также линейно зависит от размера строки m и размера столбца n матрицы, но не зависит от произведения m * n, общего количества нулевых и ненулевых элементов.

Это трудно сравнить с полной матрицей, не зная более подробных сведений.

Библиотека разреженных матриц Scipy объясняет плюсы и минусы каждого разреженного формата. Например, для csc_matrix

Преимущества формата CSC

  • эффективные арифметические операции CSC + CSC, CSC * CSC и т. д.
  • эффективная нарезка столбцов
  • быстрые матричные векторные произведения (CSR, BSR могут быть быстрее)

Недостатки формата CSC

  • медленные операции нарезки строк (рассмотрите CSR)
  • изменения в структуре разреженности дороги (рассмотрите LIL или DOK)

Существует ли аналогичная информация о реализации sparse в Matlab? Если да, то где мне его найти?


person Cecilia    schedule 21.04.2015    source источник
comment
Дайте мне знать, если вы хотите, чтобы я подробно остановился на какой-либо части моего ответа   -  person krisdestruction    schedule 22.04.2015
comment
Я бы предложил награду за этот вопрос, но у меня очень низкая репутация. @LuisMendo, может быть, ты хочешь?   -  person krisdestruction    schedule 22.04.2015
comment
Из документации по разреженным матрицам выбранная библиография: Гилберт, Джон Р., Клив Молер и Роберт Шрайбер, Разреженные матрицы в MATLAB: проектирование и реализация, SIAM J. Matrix Anal. Appl., Vol. 13, No. 1, январь 1992 г., стр. 333-356. Вероятно, это хорошая отправная точка.   -  person Sam Roberts    schedule 22.04.2015


Ответы (3)


Многие операции с полными матрицами используют вызовы библиотеки BLAS / LAPACK, которые безумно оптимизированы и их трудно превзойти. На практике операции с разреженными матрицами будут превосходить операции с полными матрицами только в специализированных ситуациях, которые могут в достаточной степени использовать (i) разреженность и (ii) особую структуру матриц.

Просто случайное использование sparse, вероятно, ухудшит ваше положение. Пример: что быстрее, добавление полной матрицы 10000x10000 к полной матрице 10000x10000? Или добавить полную матрицу 10000x10000 к полностью разреженной (то есть все равно нулю) матрице 10000x10000? попытайся! В моей системе полный + полный быстрее!

Какие есть примеры ситуаций, когда разреженные ДАВЛЕНИЯ полны?

Пример 1: решение линейной системы A * x = b, где A - 5000x5000, но это блочно-диагональная матрица, построенная из 500 блоков 5x5. Код установки:

As = sparse(rand(5, 5));
for(i=1:999)
   As = blkdiag(As, sparse(rand(5,5))); 
end;                         %As is made up of 500 5x5 blocks along diagonal
Af = full(As); b = rand(5000, 1);

Затем вы можете проверить разницу в скорости:

As \ b % operation on sparse As takes .0012 seconds
Af \ b % solving with full Af takes about 2.3 seconds

В общем, линейная система с 5000 переменными несколько сложна, но 1000 отдельных линейных систем с 5 переменными - тривиально. Последнее, в основном, и решается в редком случае.

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

С другой стороны, случайное добавление разреженного кода без глубоких, внимательных размышлений почти наверняка замедлит ваш код.

person Matthew Gunn    schedule 07.11.2015

Я не эксперт в использовании матриц sparse, однако в Mathworks есть некоторая документация, относящаяся к работе и эффективности вычислений.

Описание их вычислительной сложности:

Вычислительная сложность разреженных операций пропорциональна nnz - количеству ненулевых элементов в матрице. Вычислительная сложность также линейно зависит от размера строки m и размера столбца n матрицы, но не зависит от произведения m * n, общего количества нулевых и ненулевых элементов.

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

Не утомляя вас алгоритмическими деталями, другой ответ предлагает вам не беспокоиться о разреженных для массива, состоящего только из 25% ненулевых. Они предлагают вам код для тестирования. Подробности смотрите в их сообщении.

A = sprand(2000,2000,0.25);
tic,B = A*A;toc
Elapsed time is 1.771668 seconds.

Af = full(A);
tic,B = Af*Af;toc
Elapsed time is 0.499045 seconds.
person krisdestruction    schedule 21.04.2015
comment
Я нашел такую ​​же сложную документацию. Но как вы можете сравнить это со сложностью работы с полной матрицей? Например, предположим, что полное матричное умножение для матрицы размера nxm, а матрица mxp равна O (nmp). Если в документации для разреженных матриц указано, что сложность линейна с m и n, могу ли я перейти к выводу, что умножение равно O (m + n) или O (nnz (m) + nnz (n))? - person Cecilia; 22.04.2015
comment
Что касается использования памяти, я мог бы расширить этот фрагмент кода, чтобы он соответствовал функции с изменяющимся процентом ненулевых значений. Это может быть полезно, но похоже, что эти характеристики также должны быть присущи реализации, и мы должны быть в состоянии доказать, что матричная реализация будет использовать определенный объем памяти с учетом некоторого количества ненулевых записей. - person Cecilia; 22.04.2015
comment
По сути. Я бы хотел изрядно скучных алгоритмических деталей. :) - person Cecilia; 22.04.2015
comment
@Cecilia Я ценю тех, кому нравятся детали, но, к сожалению, я этого не знаю (как я уже сказал, я не эксперт). Теперь появляется раздел «Перестановки и перестановки», в котором описывается, как он переставляет умножение. Может быть, кто-то может быть более полезным? : / - person krisdestruction; 22.04.2015

Если у вас есть матрица фиксированного размера, лучший способ получить надежный ответ - это просто метод проб и ошибок. Однако, если вы не знаете размеров своих матриц / векторов, то практические правила таковы:

Ваши разреженные векторы должны иметь эффективное постоянное количество ненулевых записей

что для матриц будет означать

Ваша N x N разреженная матрица должна иметь <= c * N ненулевые элементы, где c - это константа, "намного меньшая", чем N.

Дадим псевдотеоретическое объяснение этому правилу. Мы рассмотрим довольно простую задачу построения скалярного (или точечного) произведения двух векторов с двузначными координатами. Теперь, если у вас есть два плотных вектора одинаковой длины N, ваш код будет выглядеть так:

//define vectors vector, wector as double arrays of length N 
double sum = 0;
for (int i = 0; i < N; i++)
{
    sum += vector[i] * wector[i];
}

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

UPD: Фактически, в for цикле вы рискуете выбрать неправильную ветвь только один раз, в конце вашего цикла, поскольку по определению выбираемая ветвь по умолчанию будет входить в цикл. Это составляет не более 1 перезапуска конвейера на операцию скалярного произведения.

Давайте теперь посмотрим, как разреженные векторы реализуются в BLAS. Здесь каждый вектор кодируется двумя массивами: одним из значений и одним из соответствующих индексов, что-то вроде

1.7    -0.8    3.6
171     83     215

(плюс одно целое число, указывающее предполагаемую длину N). В документации BLAS указано, что порядок индексов здесь не играет роли, поэтому данные

-0.8    3.6    1.7
 83     215    171

кодирует тот же вектор. Это замечание дает достаточно информации для восстановления алгоритма скалярного произведения. Учитывая два разреженных вектора, закодированных данными int[] indices, double[] values и int[] jndices, double[] walues, один будет вычислять их скалярное произведение в строках этого «кода»:

double sum = 0;
for (int i = 0; i < indices.length; i++)
{
    for (int j = 0; j < jndices.length; j++)
    {
        if(indices[i] == jndices[j])
        {
            sum += values[indices[i]] * walues[jndices[j]];
        }
    }
}

что дает нам общее количество indices.length * jndices.length * 2 + indices.length условных переходов. Это означает, что для того, чтобы справиться с плотным алгоритмом, ваши векторы должны иметь не более sqrt(N) ненулевых записей. Дело в том, что зависимость от N уже нелинейна, поэтому нет смысла спрашивать, нужно ли вам заполнение 1%, 10% или 25%. 10% идеально подходит для векторов длины 10, все еще вроде нормально для длины 50 и уже полное разорение для длины 100.

UPD. В этом фрагменте кода у вас есть ветвь if, и вероятность выбрать неверный путь составляет 50%. Таким образом, скалярное произведение двух разреженных векторов будет примерно в 0,5–1 раза больше среднего числа ненулевых записей на один разреженный вектор), в зависимости от того, насколько разрежены ваши векторы. Числа должны быть скорректированы: в операторе if без else самая короткая инструкция будет выполняться первой, то есть «ничего не делать», но все же.

Обратите внимание, что наиболее эффективная операция - это скалярное произведение разреженного и плотного векторов. Учитывая разреженный вектор indices и values и плотный вектор dense, ваш код будет выглядеть так:

double sum = 0;
for (int i = 0; i < indices.length; i++)
{
    sum += values[indices[i]] * dense[indices[i]];
}

то есть у вас будет indices.length условных ветвей, и это хорошо.

UPD. Еще раз, я уверен, что у вас будет не более одного перезапуска конвейера за операцию. Также обратите внимание, что afaik в современных многоядерных процессорах обе альтернативы выполняются параллельно на двух разных ядрах, так что в альтернативных ветвях вам нужно только дождаться завершения самого длинного из них.

Теперь при умножении матрицы на вектор вы в основном берете скалярные произведения векторов #rows. Умножение матрицы на сумму матрицы при взятии # ((ненулевых) столбцов во второй матрице) матрицы на векторные умножения. Приглашаем разобраться в сложности самостоятельно.

И вот здесь начинается вся черная магия глубокая теория хранения различных матриц. Вы можете сохранить свою разреженную матрицу как плотный массив разреженных строк, как разреженный массив плотных строк или разреженный массив разреженных строк. То же самое и со столбцами. Все забавные сокращения от Scipy, процитированные в вопросе, имеют к этому отношение.

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

person Kolya Ivankov    schedule 07.12.2017