Каков наилучший способ создания разреженного массива в C ++?

Для C ++ карта работает хорошо. Несколько миллионов объектов не будут проблемой. 10 миллионов элементов заняли около 4,4 секунды и около 57 мегабайт на моем компьютере.

Мое тестовое приложение выглядит следующим образом:

#include <stdio.h>
#include <stdlib.h>
#include <map>

class triple {
public:
    int x;
    int y;
    int z;
    bool operator<(const triple &other) const {
        if (x < other.x) return true;
        if (other.x < x) return false;
        if (y < other.y) return true;
        if (other.y < y) return false;
        return z < other.z;
    }
};

int main(int, char**)
{
    std::map<triple,int> data;
    triple point;
    int i;

    for (i = 0; i < 10000000; ++i) {
        point.x = rand();
        point.y = rand();
        point.z = rand();
        //printf("%d %d %d %d\n", i, point.x, point.y, point.z);
        data[point] = i;
    }
    return 0;
}

Теперь, чтобы динамически выбрать количество переменных, самое простое решение - представить индекс в виде строки , а затем использовать строку в качестве ключа для карты. Например, элемент, расположенный в [23] [55], может быть представлен строкой «23,55». Мы также можем расширить это решение для более высоких измерений; например, для трех измерений произвольный индекс будет выглядеть как «34,45,56». Простая реализация этого метода заключается в следующем:

std::map data<string,int> data;
char ix[100];

sprintf(ix, "%d,%d", x, y); // 2 vars
data[ix] = i;

sprintf(ix, "%d,%d,%d", x, y, z); // 3 vars
data[ix] = i;

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

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

Маленькая деталь в сравнении индекса. Вам нужно сделать лексикографическое сравнение, иначе:

a= (1, 2, 1); b= (2, 1, 2);
(a<b) == (b<a) is true, but b!=a

Изменить: Таким образом, сравнение должно быть:

return lhs.x<rhs.x
    ? true 
    : lhs.x==rhs.x 
        ? lhs.y<rhs.y 
            ? true 
            : lhs.y==rhs.y
                ? lhs.z<rhs.z
                : false
        : false

Boost имеет шаблонную реализацию BLAS, которая называется uBLAS и содержит разреженную матрицу.

https://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm

Как общий совет, метод, использующий строки в качестве индексов, на самом деле очень медленный. Гораздо более эффективным, но в остальном эквивалентным решением было бы использование векторов / массивов. Абсолютно нет необходимости записывать индексы в строку.

typedef vector<size_t> index_t;

struct index_cmp_t : binary_function<index_t, index_t, bool> {
    bool operator ()(index_t const& a, index_t const& b) const {
        for (index_t::size_type i = 0; i < a.size(); ++i)
            if (a[i] != b[i])
                return a[i] < b[i];
        return false;
    }
};

map<index_t, int, index_cmp_t> data;
index_t i(dims);
i[0] = 1;
i[1] = 2;
// … etc.
data[i] = 42;

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

Поскольку только значения с [a] [b] [c] ... [w] [x] [y] [z] имеют значение, мы храним только сам индекс, а не значение 1, которое почти везде - всегда тот же + нет способа хэшировать. Отмечая, что проклятие размерности присутствует, предложите пойти с каким-то установленным инструментом NIST или Boost, хотя бы прочитайте источники для этого, чтобы обойти ненужную ошибку.

Если в работе необходимо зафиксировать распределения временной зависимости и параметрические тенденции неизвестных наборов данных, то карта или B-дерево с однозначным корнем, вероятно, не практичны. Мы можем хранить только сами индексы, хэшированные, если порядок (чувствительность к представлению) может подчиняться сокращению временной области во время выполнения для всех 1 значений. Поскольку ненулевые значения, отличные от единицы, немногочисленны, очевидным кандидатом для них является любая структура данных, которую вы можете легко найти и понять. Если набор данных действительно огромного размера, то я предлагаю какое-то скользящее окно, которое самостоятельно управляет файлом / диском / persistent-io, перемещая части данных в область по мере необходимости. (Написание кода, который вы можете понять) Если вы обязуетесь предоставить реальное решение для рабочей группы,

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

// Copyright 2014 Leo Osvald
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#ifndef UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_
#define UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_

#include <algorithm>
#include <limits>
#include <map>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

// A simple time-efficient implementation of an immutable sparse matrix
// Provides efficient iteration of non-zero elements by rows/cols,
// e.g. to iterate over a range [row_from, row_to) x [col_from, col_to):
//   for (int row = row_from; row < row_to; ++row) {
//     for (auto col_range = sm.nonzero_col_range(row, col_from, col_to);
//          col_range.first != col_range.second; ++col_range.first) {
//       int col = *col_range.first;
//       // use sm(row, col)
//       ...
//     }
template<typename T = double, class Coord = int>
class SparseMatrix {
  struct PointHasher;
  typedef std::map< Coord, std::vector<Coord> > NonZeroList;
  typedef std::pair<Coord, Coord> Point;

 public:
  typedef T ValueType;
  typedef Coord CoordType;
  typedef typename NonZeroList::mapped_type::const_iterator CoordIter;
  typedef std::pair<CoordIter, CoordIter> CoordIterRange;

  SparseMatrix() = default;

  // Reads a matrix stored in MatrixMarket-like format, i.e.:
  // <num_rows> <num_cols> <num_entries>
  // <row_1> <col_1> <val_1>
  // ...
  // Note: the header (lines starting with '%' are ignored).
  template<class InputStream, size_t max_line_length = 1024>
  void Init(InputStream& is) {
    rows_.clear(), cols_.clear();
    values_.clear();

    // skip the header (lines beginning with '%', if any)
    decltype(is.tellg()) offset = 0;
    for (char buf[max_line_length + 1];
         is.getline(buf, sizeof(buf)) && buf[0] == '%'; )
      offset = is.tellg();
    is.seekg(offset);

    size_t n;
    is >> row_count_ >> col_count_ >> n;
    values_.reserve(n);
    while (n--) {
      Coord row, col;
      typename std::remove_cv<T>::type val;
      is >> row >> col >> val;
      values_[Point(--row, --col)] = val;
      rows_[col].push_back(row);
      cols_[row].push_back(col);
    }
    SortAndShrink(rows_);
    SortAndShrink(cols_);
  }

  const T& operator()(const Coord& row, const Coord& col) const {
    static const T kZero = T();
    auto it = values_.find(Point(row, col));
    if (it != values_.end())
      return it->second;
    return kZero;
  }

  CoordIterRange
  nonzero_col_range(Coord row, Coord col_from, Coord col_to) const {
    CoordIterRange r;
    GetRange(cols_, row, col_from, col_to, &r);
    return r;
  }

  CoordIterRange
  nonzero_row_range(Coord col, Coord row_from, Coord row_to) const {
    CoordIterRange r;
    GetRange(rows_, col, row_from, row_to, &r);
    return r;
  }

  Coord row_count() const { return row_count_; }
  Coord col_count() const { return col_count_; }
  size_t nonzero_count() const { return values_.size(); }
  size_t element_count() const { return size_t(row_count_) * col_count_; }

 private:
  typedef std::unordered_map<Point,
                             typename std::remove_cv<T>::type,
                             PointHasher> ValueMap;

  struct PointHasher {
    size_t operator()(const Point& p) const {
      return p.first << (std::numeric_limits<Coord>::digits >> 1) ^ p.second;
    }
  };

  static void SortAndShrink(NonZeroList& list) {
    for (auto& it : list) {
      auto& indices = it.second;
      indices.shrink_to_fit();
      std::sort(indices.begin(), indices.end());
    }

    // insert a sentinel vector to handle the case of all zeroes
    if (list.empty())
      list.emplace(Coord(), std::vector<Coord>(Coord()));
  }

  static void GetRange(const NonZeroList& list, Coord i, Coord from, Coord to,
                       CoordIterRange* r) {
    auto lr = list.equal_range(i);
    if (lr.first == lr.second) {
      r->first = r->second = list.begin()->second.end();
      return;
    }

    auto begin = lr.first->second.begin(), end = lr.first->second.end();
    r->first = lower_bound(begin, end, from);
    r->second = lower_bound(r->first, end, to);
  }

  ValueMap values_;
  NonZeroList rows_, cols_;
  Coord row_count_, col_count_;
};

#endif  /* UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_ */

Для простоты это так immutable, но вы можете сделать его изменчивым; не забудьте изменить std::vectorна, std::setесли вы хотите разумные эффективные «вставки» (изменение нуля на ненулевой).

Eigen - это библиотека линейной алгебры C ++, в которой реализована разреженная матрица. Он даже поддерживает матричные операции и решатели (факторизация LU и т. Д.), Которые оптимизированы для разреженных матриц.

Я бы предложил сделать что-то вроде:

typedef std::tuple<int, int, int> coord_t;
typedef boost::hash<coord_t> coord_hash_t;
typedef std::unordered_map<coord_hash_t, int, c_hash_t> sparse_array_t;

sparse_array_t the_data;
the_data[ { x, y, z } ] = 1; /* list-initialization is cool */

for( const auto& element : the_data ) {
    int xx, yy, zz, val;
    std::tie( std::tie( xx, yy, zz ), val ) = element;
    /* ... */
}

Чтобы сохранить разреженность данных, вы можете написать подкласс unorderd_map, чьи итераторы автоматически пропускают (и стирают) любые элементы со значением 0.

Полный список решений можно найти в википедии. Для удобства я процитировал соответствующие разделы следующим образом.

https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_.28DOK.29

Словарь ключей (ДОК)

DOK состоит из словаря, который сопоставляет пары (строка, столбец) со значением элементов. Элементы, отсутствующие в словаре, принимаются равными нулю. Формат подходит для постепенного построения разреженной матрицы в случайном порядке, но не подходит для итерации по ненулевым значениям в лексикографическом порядке. Один обычно создает матрицу в этом формате, а затем преобразует в другой, более эффективный формат для обработки. [1]

Список списков (LIL)

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

Список координат (COO)

COO хранит список (строка, столбец, значение) кортежей. В идеале, записи сортируются (по индексу строки, затем по индексу столбца), чтобы улучшить время произвольного доступа. Это еще один формат, который подходит для построения инкрементных матриц. [3]

Сжатый разреженный ряд (CSR, CRS или формат Yale)

Формат сжатой разреженной строки (CSR) или хранилища сжатой строки (CRS) представляет собой матрицу M из трех (одномерных) массивов, которые соответственно содержат ненулевые значения, экстенты строк и индексы столбцов. Он похож на COO, но сжимает индексы строк, отсюда и название. Этот формат обеспечивает быстрый доступ к строкам и умножение матрицы на вектор (Mx).