C++でスパース配列を作成する最も良い方法は何ですか？

Question

私は巨大な行列の操作、特に、コピュラ計算のためのピラミッドな加算を必要とするプロジェクトに取り組んでいます。

要するに、行列（多次元配列）内のゼロの海で比較的少数の値（通常は1の値、まれには1以上の値）を追跡する必要があります。

スパース配列を使用すると、少数の値を格納し、すべての未定義レコードをプリセット値とみなすことができます。物理的にはすべての値をメモリに格納するわけではないので、ゼロ以外の要素はほとんど保存する必要はありません。これは数百万エントリになる可能性があります。

速度は非常に重要です。実行時にクラス内の変数の数を動的に選択することもできます。

私は現在、エントリを格納するためにバイナリ検索ツリー（bツリー）を使用するシステムで動作します。誰かがより良いシステムを知っていますか？

Answer 1

C++では、マップがうまく機能します。数百万のオブジェクトが問題になることはありません。私のコンピュータには約10秒間で約4.4秒かかっていた。次のように

私のテストアプリケーションは次のとおりです。

#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"文字列によって表すことができる。より高次元のためにこのソリューションを拡張することもできます。例えば3次元の場合、任意のインデックスは「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; 

Answer 2

ハッシュテーブルには高速挿入とルックアップがあります。整数のペアだけをキーとして扱うことがわかっているので、単純なハッシュ関数を書くことができます。

Answer 3

スパース行列を実装する最善の方法は、それを実装しないことです - 少なくともあなた自身ではありません。私はBLAS（私はLAPACKの一部であると思う）に、本当に巨大な行列を扱えることを提案します。

Answer 4

インデックス比較の細部。

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 

Answer 5

ブーストは、BLASのテンプレート実装では、スパース行列が含まれているuBLASを呼びかけています。ただ、アドバイスとして

http://www.boost.org/doc/libs/1_36_0/libs/numeric/ublas/doc/index.htm

Answer 6

：インデックスとして文字列を使用する方法は、実際には非常に遅いです。はるかに効率的ですが、等価なソリューションはベクトル/配列を使用することです。インデックスを文字列に書く必要は全くありません。

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を使用するため、平衡二分探索木の面で実装の実際には非常に効率的ではありません。この場合のデータ構造のパフォーマンスは、（ランダム化された）ハッシュテーブルのほうがはるかに優れています。

Answer 7

[a] [b] [c] ... [w] [x] [y] [z]の値は重要なので、値1だけではなく、いつでもどこでも - それをハッシュする方法は全くありません。次元の呪いが存在することを銘記すれば、NISTやBoostの確立されたツールを使うことをお勧めします。少なくとも、不必要な大失敗を回避するための情報源を読んでください。

未知のデータセットの時間依存性分布とパラメトリック傾向を取得する必要がある場合、単値ルートを持つマップまたはBツリーはおそらく実用的ではありません。実行時に1つの値すべてに対して、順序付け（提示の感性）が実行時の時間領域の縮小に従属する場合、ハッシュしたインデックスだけを格納できます。 1以外のゼロ以外の値は少ないので、それらの明白な候補は、簡単に見つけて理解できるデータ構造です。データセットが本当に巨大なユニバースサイズの場合は、ファイル/ディスク/永続性を自分で管理し、必要に応じてデータの一部を範囲に移動する、ある種のスライディングウィンドウを提案します。 （あなたが理解できる記述コード）実際の解決策をワーキンググループに提供しようと約束している場合、そうしなければ、昼食を取ることを唯一の目標とする消費者グレードのオペレーティングシステムの慈悲が残されます。

Answer 8

これは、妥当な高速検索（ハッシュテーブルを使用）と、行/列の非ゼロ要素に対する高速反復を提供する比較的簡単な実装です。簡単にするために

// 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に変更してください。

Answer 9

Eigenは、疎行列のimplementationを有するC++線形代数ライブラリである。これは、疎行列に最適化された行列演算とソルバ（LU分解など）をサポートします。

Answer 10

私が何かやってお勧めします：スパースデータを保つのを助けるために

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

Answer 11

解決策の完全なリストは、ウィキペディアで見つけることができます。便宜上、関連セクションを以下のように引用しました。

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

キーの辞典（DOK）

DOK要素の 値（行、列）-pairsをマッピングする辞書から成ります。辞書 から欠落している要素はゼロと見なされます。この形式は漸進的に漸進的に疎な行列を構築するのに適していますが、字句的な順序でゼロ以外の値を反復するのには貧弱です（ ）。 1つの典型的には、 は、このフォーマットのマトリックスを構築し、さらに処理するために効率的なフォーマットに変換する。リスト（LIL）

LIL店の[1]

リストカラム 指標と値を含む各エントリに行ごとにリスト。通常、これらのエントリは、検索の高速化のために、 列インデックスでソートされています。これは、 インクリメンタルマトリックス構造に適した別の形式です。[2]

座標リスト（COO）

COOストア（行、列、値）タプルのリスト。理想的には、エントリ はランダムアクセスを改善するために（行インデックス、列インデックスによって）並べ替えられます 回。これは増分マトリックス の構成に適した別のフォーマットです。[3]

圧縮スパース行（CSR、CRSまたはエール・フォーマット）

圧縮スパース行（CSR）または圧縮行ストレージ（CRS）形式 三（ワンにより行列Mを表します にはそれぞれ非ゼロ値、行のエクステント、および列 インデックスが含まれています。これはCOOに似ていますが、行インデックスを圧縮するため、 という名前になります。このフォーマットは、高速の行アクセスと行列ベクトル の乗算（Mx）を可能にします。