c++ - 在 C++ 中创建稀疏数组的最佳方法是什么？

Question

我正在研究一个需要处理巨大矩阵的项目，特别是用于 copula 计算的金字塔求和。

简而言之，我需要在矩阵（多维数组）中的大量零中跟踪相对少量的值（通常为 1，在极少数情况下超过 1）。

稀疏数组允许用户存储少量值，并将所有未定义的记录假定为预设值。由于物理上不可能将所有值存储在内存中，因此我只需要存储少数非零元素。这可能是几百万个条目。

速度是重中之重，我还想在运行时动态选择类中变量的数量。

我目前在一个使用二叉搜索树（b-tree）来存储条目的系统上工作。有人知道更好的系统吗？

score 29 · Accepted Answer

对于 C++，映射效果很好。几百万个对象不会是一个问题。1000 万个项目在我的计算机上花费了大约 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;
}

现在要动态选择变量的数量，最简单的解决方案是将index 表示为 string，然后使用 string 作为 map 的键。例如，位于 [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;

score 23 · Accepted Answer

接受的答案建议使用字符串来表示多维索引。

但是，为此构造字符串是不必要的浪费。如果在编译时不知道大小（因此std::tuple不起作用），则std::vector可以很好地用作索引，无论是哈希映射还是有序树。对于std::map，这几乎是微不足道的：

#include <vector>
#include <map>

using index_type = std::vector<int>;

template <typename T>
using sparse_array = std::map<index_type, T>;

对于std::unordered_map（或类似的基于哈希表的字典），它的工作量稍大一些，因为std::vector不专门std::hash：

#include <vector>
#include <unordered_map>
#include <numeric>

using index_type = std::vector<int>;

struct index_hash {
    std::size_t operator()(index_type const& i) const noexcept {
        // Like boost::hash_combine; there might be some caveats, see
        // <https://stackoverflow.com/a/50978188/1968>
        auto const hash_combine = [](auto seed, auto x) {
            return std::hash<int>()(x) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
        };
        return std::accumulate(i.begin() + 1, i.end(), i[0], hash_combine);
    }
};

template <typename T>
using sparse_array = std::unordered_map<index_type, T, index_hash>;

无论哪种方式，用法都是相同的：

int main() {
    using i = index_type;

    auto x = sparse_array<int>();
    x[i{1, 2, 3}] = 42;
    x[i{4, 3, 2}] = 23;

    std::cout << x[i{1, 2, 3}] + x[i{4, 3, 2}] << '\n'; // 65
}

score 13 · Accepted Answer

Boost 有一个名为 uBLAS 的 BLAS 模板化实现，其中包含一个稀疏矩阵。

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

score 5 · Accepted Answer

Eigen是一个 C++ 线性代数库，具有稀疏矩阵的实现。它甚至支持针对稀疏矩阵进行优化的矩阵运算和求解器（LU 分解等）。

score 4 · Accepted Answer

指数比较中的小细节。您需要进行字典比较，否则：

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

score 4 · Accepted Answer

完整的解决方案列表可以在维基百科中找到。为方便起见，我将相关部分引述如下。

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

密钥字典 (DOK)

DOK 包含一个将（行、列）对映射到元素值的字典。字典中缺少的元素被视为零。该格式适用于以随机顺序增量构造稀疏矩阵，但不适用于按字典顺序迭代非零值。通常以这种格式构造一个矩阵，然后转换为另一种更有效的格式进行处理。 [1]

列表列表 (LIL)

LIL 每行存储一个列表，每个条目包含列索引和值。通常，这些条目按列索引排序，以便更快地查找。这是另一种适用于增量矩阵构造的格式。 [2]

坐标列表（COO）

COO 存储（行、列、值）元组的列表。理想情况下，对条目进行排序（按行索引，然后按列索引）以缩短随机访问时间。这是另一种有利于增量矩阵构造的格式。 [3]

压缩稀疏行（CSR、CRS 或耶鲁格式）

压缩稀疏行 (CSR) 或压缩行存储 (CRS) 格式表示由三个（一维）数组组成的矩阵 M，它们分别包含非零值、行的范围和列索引。它类似于 COO，但压缩了行索引，因此得名。这种格式允许快速行访问和矩阵向量乘法 (Mx)。

score 3 · Accepted Answer

哈希表具有快速插入和查找功能。您可以编写一个简单的散列函数，因为您知道您将只处理整数对作为键。

score 1 · Accepted Answer

实现稀疏矩阵的最好方法是不要实现它们——至少不是你自己。我建议 BLAS（我认为它是 LAPACK 的一部分）可以处理非常大的矩阵。

score 0 · Accepted Answer

因为只有带有 [a][b][c]...[w][x][y][z] 的值是重要的，所以我们只存储索引本身，而不是几乎无处不在的值 1 - 总是相同+无法对其进行哈希处理。注意到存在维度的诅咒，建议使用一些已建立的工具 NIST 或 Boost，至少阅读源代码以避免不必要的错误。

如果工作需要捕获未知数据集的时间依赖性分布和参数趋势，那么具有单值根的 Map 或 B-Tree 可能不实用。对于所有 1 值，我们只能存储索引本身，如果排序（表示的敏感性）可以从属于运行时时域的缩减，则进行散列处理。由于除 1 之外的非零值很少，因此显而易见的候选者是您可以轻松找到并理解的任何数据结构。如果数据集真的是巨大的宇宙大小，我建议使用某种滑动窗口来自己管理文件/磁盘/持久性io，根据需要将部分数据移动到范围内。（编写您可以理解的代码）如果您承诺为工作组提供实际解决方案，

score 0 · Accepted Answer

这是一个相对简单的实现，它应该提供合理的快速查找（使用哈希表）以及对行/列中非零元素的快速迭代。

// 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

score 0 · Accepted Answer

我建议做类似的事情：

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 的项目。

c++ - 在 C++ 中创建稀疏数组的最佳方法是什么？

11 回答 11

Related

Reference