Что мы можем сделать для реализации SparseMatrix?

Question

Я использовал матрицу Eigen и библиотеку SparseVector. У меня были проблемы с производительностью, и все, что мне нужно, это продукт Sparse vector dot. Поэтому я развернул свою собственную реализацию SparseMatrix, надеясь, что это будет как-то немного быстрее:

Немного примера кода:

#include <map>
using namespace std ;

struct SparseMatrix
{
  map<int, Vector> vals ;

  Vector dot( SparseMatrix& o )
  {
    SparseMatrix *LS, *RS ;

    // iterate over the smaller of the 2 
    if( vals.size() < o.vals.size() )
    {
      // walk vals
      LS = this ;
      RS = &o ;
    }
    else
    {
      LS = &o ;
      RS = this ;
    }

    // walk LS
    Vector sum = 0 ;

    for( map<int,Vector>::iterator LSIter = LS->vals.begin() ; LSIter != LS->vals.end() ; ++LSIter )
    {
      const int& key = LSIter->first ;

      // use the key, see if RS has a similar entry.
      map<int,Vector>::iterator RSIter = RS->vals.find( key );
      if( RSIter != RS->vals.end() )
        sum += RSIter->second * LSIter->second ;
    }

    return sum ;
  }
} ;

Итак, скалярное произведение двух векторов, скажем, имеет такие записи, как:

+---------------+
| vec 1         |
| index   value |
|   2      18   |
|   7       4   |
|  18      33   |
+---------------+

+---------------+
| vec 2         |
| index   value |
|   2       1   |
|  15      87   |
|  21      92   |
+---------------+

Точечный продукт тогда равен 18.

Итак, как вы можете видеть, я использовал std::map для поиска элементов, чтобы увидеть, находится ли элемент из одного вектора в другом векторе.

Поскольку я использую только целочисленную индексацию и 1d-массивы, есть ли способ ускорить поиск? Мое умножение разреженных матриц все еще является узким местом (производительность моего кода лишь незначительно выше, чем у Eigen)

Answer 1

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

Если, конечно, вы не сильно меняете данные.

т.е. в псевдокоде:

for i = 0, j = 0; i < vec1.size && j < vec2.size:
    if vec1[i].first == vec2[j].first:
        result += vec1[i++].second * vec2[j++].second
    elif vec1[i].first < vec2[j].first:
        i++
    else
        j++