2008-08-07 2 views
43

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

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

Редкий массив позволяет пользователю хранить небольшое количество значений и считать все неопределенные записи заданными значениями. Поскольку физически невозможно сохранить все значения в памяти, мне нужно сохранить только несколько ненулевых элементов. Это может быть несколько миллионов записей.

Скорость является огромным приоритетом, и я также хотел бы динамически выбирать количество переменных в классе во время выполнения.

В настоящее время я работаю над системой, использующей двоичное дерево поиска (b-tree) для хранения записей. Кто-нибудь знает о лучшей системе?

ответ

25

Для 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; 
+1

Что можно сказать о том, как получить диапазон элементов от этого или проверить, полностью ли диапазон в массиве? – 2012-03-20 11:18:10

+2

Выполнение оператора <неверно. Рассмотрим Triple {1,2,3} и Triple {3,2,1}, и не будет меньше, чем другой. Правильная реализация проверяет x, затем y, затем z последовательно, а не сразу. – Whanhee 2014-06-27 15:14:45

+1

Поскольку это не было исправлено в течение длительного времени, я взял на себя смелость заменить его правильной реализацией. – celtschk 2015-05-06 07:34:07

3

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

1

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

4

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

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

Edit: Так что сравнение, вероятно, следует:

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 
16

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

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 это на самом деле не очень эффективным, поскольку реализации с точки зрения сбалансированного двоичного дерева поиска.В этом случае гораздо более эффективными структурами данных будет хеш-таблица (рандомизированная).

0

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

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

0

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

// 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, если вы хотите разумные эффективные «вставки» (изменение нуля на ненулевое).

3

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

0

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

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.

2

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

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

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

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

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

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

координат лист (СОО)

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

Сжатые редкие строки (КСО, СКИ или формат Йельского)

Сжатых редкие строки (КСО) или сжатое хранение строк (СКИ) Формат представляет собой матрицу М три (одно- мерных) массивов, то соответственно содержат ненулевые значения, экстенты строк и столбцы . Он похож на COO, но сжимает индексы строк, поэтому имя. Этот формат обеспечивает быстрый доступ к строке и матричный вектор умножения (Mx).

Смежные вопросы