Быстрое умножение разреженных матриц

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

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

Прямо сейчас мои разреженные матрицы в основном реализованы как завернутый std::map< std::pair<int, int>, double>, в котором хранятся данные, если таковые имеются. Это преобразует умножение матрицы с вектором на сложность O (n²) в O (n²log (n)), поскольку мне нужно выполнить поиск для каждого элемента матрицы. Я изучил формат матрицы Yale Sparse, и кажется, что получение элемента также находится в O (log (n)), поэтому я не уверен, будет ли это намного быстрее.

Для справки у меня есть матрица 800x800, заполненная 5000 записями. Решение такой системы с помощью метода сопряженных градиентов занимает примерно 450 секунд.

Как вы думаете, можно ли сделать это намного быстрее с другой структурой данных?

Благодарность!


c++ matrix sparse-matrix matrix-multiplication linear-algebra
person lezebulon    schedule 28.03.2013    source источник
comment
Сначала прочтите википедию. en.wikipedia.org/wiki/Sparse_matrix содержит хороший список распространенных методов хранения, которые даст вам эффективные операции.   -  person Daniel Williams    schedule 29.03.2013
comment
@Song Wang: цель класса в основном состоит в том, чтобы развернуть наш собственный решатель метода конечных элементов   -  person lezebulon    schedule 29.03.2013

Ответы (1)


arrow_upward
25
arrow_downward

Наиболее распространенный выбор - CSC или хранилище CSR. Оба они эффективны для умножения матрицы на вектор. Также очень легко запрограммировать эти процедуры умножения, если вам придется делать это самому.

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

Даже с вашим текущим хранилищем вы можете выполнять матричное умножение со сложностью O (n). Все алгоритмы умножения разреженных матриц на вектор, которые я когда-либо видел, сводятся к одним и тем же шагам. Например, рассмотрим y = Ax.

  1. Обнулить результирующий вектор y.
  2. Инициализируйте итератор для ненулевых элементов матрицы A.
  3. Получите следующий ненулевой элемент матрицы, скажем, A [i, j]. Обратите внимание, что образец i, j не соответствует обычному образцу. Он просто отражает порядок, в котором хранятся ненулевые элементы A.
  4. y[i] += A[i,j]*x[j]
  5. Если элементов A больше, перейдите к 3.

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

for (i=0; i<N; i++)
    for (j=0; j<N; j++)
        y[i] += A[i,j]*x[j]

и это то, что побуждает вас выполнять поиск.

Но я не предлагаю вам использовать std::map хранилище. Это не будет суперэффективным. Я бы порекомендовал CSC главным образом потому, что он наиболее широко используется.

person David Heffernan    schedule 28.03.2013
comment
Если вы выполняете поиск элементов матрицы, значит, вы неправильно поняли, как использовать формат, да, вы совершенно правы, это была моя первоначальная проблема. Например, для метода CSR у меня будет такая же сложность поиска, но на самом деле мне не нужно выполнять поиск, чтобы выполнить умножение матрицы на вектор, просто чтобы один раз пройти по массиву - person lezebulon; 29.03.2013
comment
Что касается вашего редактирования: вы правы, это тоже сработает, но только потому, что мои баллы уже отсортированы по строкам, верно? - person lezebulon; 29.03.2013
comment
Просто пройти по массиву один раз. Это именно то, что нужно. Вы получили это сейчас! Это требует некоторого сдвига в уме, но если вы так думаете, вы идете прямо домой! - person David Heffernan; 29.03.2013
comment
FWIW 800x800 с nz = 5000 должен быть почти мгновенным. Вы должны сделать около 1000 линейных решений для такой матрицы за секунду. - person David Heffernan; 29.03.2013
comment
Спасибо за информацию! Я быстро реализовал умножение в O (n), о котором вы упомянули, и общее время решения теперь составляет ~ 3 секунды, с выполнением 3000 умножений матрицы на вектор. Поэтому я думаю, что следующим шагом будет улучшение моего решателя, но это уже другой вопрос! - person lezebulon; 29.03.2013
comment
Это действительно зависит от формы ваших матриц. Мой код FE работает с длинными балками, и я думаю, что проблемы действительно легко решить. Я лично использую CSparse Тима Дэвиса, и это превосходно. Его книга-компаньон также превосходна. - person David Heffernan; 29.03.2013