Быстрое умножение разреженных матриц
для класса я должен написать свой собственный решатель линейных уравнений для разреженных матриц.
Я свободен использовать любой тип структуры данных для разреженных матриц, и мне нужно реализовать несколько решений, в том числе сопряженный градиент.
Мне было интересно, есть ли известный способ хранения разреженных матриц, чтобы умножение с вектором было относительно быстрым.
Прямо сейчас мои разреженные матрицы в основном реализованы std::map< std::pair<int, int>, double> который хранит данные, если таковые имеются. Это преобразует умножение матрицы из вектора в сложность O (n²) в O (n²log (n)), так как я должен выполнить поиск для каждого элемента матрицы.
Я посмотрел в формате матрицы Yale Sparse, и кажется, что поиск элемента также в O (log (n)), поэтому я не уверен, что это будет намного быстрее.
Для справки у меня есть матрица 800×800, которая заполнена 5000 записей. Требуется примерно до 450 секунд, чтобы решить такую систему с помощью метода сопряженных градиентов.
Как вы думаете, возможно ли сделать это намного быстрее с другой структурой данных?
Спасибо!
Решение
Наиболее распространенные варианты CSC или CSR хранилище. Они оба эффективны для умножения матрицы на вектор. Также очень легко кодировать эти процедуры умножения, если вам придется делать это самостоятельно.
Тем не менее, хранилище Yale также дает очень эффективное умножение матрицы на вектор. Если вы выполняете поиск элементов матрицы, вы неправильно поняли, как использовать формат. Я предлагаю вам изучить некоторые стандартные разреженные библиотеки, чтобы узнать, как реализовано умножение матрицы на вектор.
Даже с вашим текущим хранилищем вы можете выполнять матричное умножение в сложности O (n). Все редкие алгоритмы умножения матрицы на вектор, которые я когда-либо видел, сводятся к одним и тем же шагам. Например, рассмотрим y = Ax.
- Обнулить вектор результата, у.
- Инициализируйте итератор для ненулевых элементов матрицы A.
- Получите следующий ненулевой элемент матрицы, скажем, A [i, j]. Обратите внимание, что шаблон i, j не следует регулярному шаблону. Это просто отражает порядок, в котором ненулевые элементы A хранятся.
- y [i] + = A [i, j] * x [j]
- Если элементов A больше, переходите к 3.
Я подозреваю, что вы пишете классический двойной код для умножения с плотным циклом:
for (i=0; i<N; i++)
for (j=0; j<N; j++)
y[i] += A[i,j]*x[j]
и это то, что ведет вас к поискам.
Но я не предлагаю вам придерживаться std::map место хранения. Это не будет супер эффективным. Я бы рекомендовал CSC в основном потому, что он наиболее широко используется.
Другие решения
Других решений пока нет …