Как лучше всего создать разреженный массив на C++?
https://stackoverflow.com/questions/4306
-
08-06-2019 - |
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianВопрос
Я работаю над проектом, который требует манипулирования огромными матрицами, в частности пирамидального суммирования для расчета копулы.
Короче говоря, мне нужно отслеживать относительно небольшое количество значений (обычно значение 1, а в редких случаях больше 1) в море нулей в матрице (многомерном массиве).
Разреженный массив позволяет пользователю хранить небольшое количество значений и считать все неопределенные записи заданными значениями.Поскольку физически невозможно хранить все значения в памяти, мне нужно сохранить только несколько ненулевых элементов.Это может быть несколько миллионов записей.
Скорость является огромным приоритетом, и мне также хотелось бы динамически выбирать количество переменных в классе во время выполнения.
В настоящее время я работаю над системой, которая использует двоичное дерево поиска (b-дерево) для хранения записей.Кто-нибудь знает лучшую систему?
Решение
Для 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;
Другие советы
Общий совет: метод, использующий строки в качестве индексов, на самом деле является очень медленный.Гораздо более эффективным, но в остальном эквивалентным решением было бы использование векторов/массивов.Совершенно нет необходимости записывать индексы в строку.
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.
Boost имеет шаблонную реализацию BLAS, называемую uBLAS, которая содержит разреженную матрицу.
http://www.boost.org/doc/libs/1_36_0/libs/numeric/ublas/doc/index.htm
Небольшая деталь в сравнении индексов.Вам нужно выполнить лексикографическое сравнение, иначе:
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
Хэш-таблицы имеют быструю вставку и поиск.Вы можете написать простую хэш-функцию, поскольку знаете, что в качестве ключей будете иметь дело только с парами целых чисел.
Эйген — это библиотека линейной алгебры C++, имеющая выполнение разреженной матрицы.Он даже поддерживает матричные операции и решатели (LU-факторизация и т. д.), оптимизированные для разреженных матриц.
Полный список решений можно найти в Википедии.Для удобства я процитировал соответствующие разделы следующим образом.
https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_.28DOK.29
Словарь ключей (ДОК)
DOK состоит из словаря, который отображает (ряд, столбец) -PAIRS до значения элементов.Элементы, которые отсутствуют в словаре, воспринимаются как ноль.Формат хорош для постепенного строительства разреженной матрицы в случайном порядке, но плохо для итерации над ненулевыми значениями в лексикографическом порядке.Один обычно строит матрицу в этом формате, а затем преобразует в другой более эффективный формат для обработки. [1
Список списков (LIL)
LIL хранит один список на строку, причем каждая запись содержит индекс столбца и значение.Как правило, эти записи соблюдаются с помощью индекса столбца для более быстрого поиска.Это еще один формат хорошего для инкрементной матричной конструкции. [2
Координационный список (COO)
COO хранит список кортежей (строка, столбец, значение).В идеале записи сортируются (по индексу строк, а затем столбцом индекса) для улучшения случайного времени доступа.Это еще один формат, который хорош для инкрементной матричной конструкции. [3
Сжатая разреженная строка (формат CSR, CRS или Yale)
Формат сжатой разреженной строки (CSR) или формат сжатой строки (CRS) представляет матрицу M на три (одномерные) массивы, которые соответственно содержат ненулевые значения, экстенты строк и индексы столбцов.Это похоже на COO, но сжимает индексы строк, отсюда и название.Этот формат обеспечивает быстрый доступ к ряду и умножению матрицы-вектора (MX).
Лучший способ реализовать разреженные матрицы — не реализовывать их, по крайней мере, не самостоятельно.Я бы посоветовал BLAS (который, как мне кажется, является частью LAPACK), который может обрабатывать действительно огромные матрицы.
Поскольку только значения с [a][b][c]...[w][x][y][z] имеют значение, мы сохраняем только сам индекс, а не значение 1, которое есть почти везде - всегда то же самое + нет возможности хешировать.Отмечая, что проклятие размерности присутствует, предложите использовать какой-нибудь установленный инструмент NIST или Boost, по крайней мере, прочитайте источники, чтобы избежать ненужной ошибки.
Если в работе необходимо отразить распределения временных зависимостей и параметрические тенденции неизвестных наборов данных, то карта или B-дерево с однозначным корнем, вероятно, нецелесообразно.Мы можем хранить только сами индексы, хешированные, если порядок (чувствительность к представлению) может подчиняться сокращению временной области во время выполнения, для всех 1 значений.Поскольку ненулевых значений, отличных от единицы, мало, очевидным кандидатом для них является любая структура данных, которую вы можете легко найти и понять.Если набор данных действительно имеет размер огромной вселенной, я предлагаю какое-то скользящее окно, которое самостоятельно управляет файлом/диском/постоянным вводом-выводом, перемещая части данных в область видимости по мере необходимости.(написание кода, который вы можете понять). Если вы берете на себя обязательство предоставить реальное решение рабочей группе, невыполнение этого требования оставляет вас во власти операционных систем потребительского уровня, единственная цель которых - отобрать у вас обед.
Вот относительно простая реализация, которая должна обеспечивать разумный быстрый поиск (с использованием хеш-таблицы), а также быструю итерацию по ненулевым элементам в строке/столбце.
// 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 если вам нужны разумные эффективные «вставки» (замена нуля на ненулевое).
Я бы предложил сделать что-то вроде:
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.
