What is the best way to create a sparse array in C++?

The accepted answer recommends using strings to represent multi-dimensional indices.

However, constructing strings is needlessly wasteful for this. If the size isn’t known at compile time (and thus std::tuple doesn’t work), std::vector works well as an index, both with hash maps and ordered trees. For std::map, this is almost trivial:

#include <vector> #include <map> using index_type = std::vector<int>; template <typename T> using sparse_array = std::map<index_type, T>; 

For std::unordered_map (or similar hash table-based dictionaries) it’s slightly more work, since std::vector does not specialise 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>; 

Either way, the usage is the same:

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 } 

Boost has a templated implementation of BLAS called uBLAS that contains a sparse matrix.

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

Complete list of solutions can be found in the wikipedia. For convenience, I have quoted relevant sections as follows.

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

Dictionary of keys (DOK)

DOK consists of a dictionary that maps (row, column)-pairs to the value of the elements. Elements that are missing from the dictionary are taken to be zero. The format is good for incrementally constructing a sparse matrix in random order, but poor for iterating over non-zero values in lexicographical order. One typically constructs a matrix in this format and then converts to another more efficient format for processing.[1]

List of lists (LIL)

LIL stores one list per row, with each entry containing the column index and the value. Typically, these entries are kept sorted by column index for faster lookup. This is another format good for incremental matrix construction.[2]

Coordinate list (COO)

COO stores a list of (row, column, value) tuples. Ideally, the entries are sorted (by row index, then column index) to improve random access times. This is another format which is good for incremental matrix construction.[3]

Compressed sparse row (CSR, CRS or Yale format)

The compressed sparse row (CSR) or compressed row storage (CRS) format represents a matrix M by three (one-dimensional) arrays, that respectively contain nonzero values, the extents of rows, and column indices. It is similar to COO, but compresses the row indices, hence the name. This format allows fast row access and matrix-vector multiplications (Mx).

Eigen is a C++ linear algebra library that has an implementation of a sparse matrix. It even supports matrix operations and solvers (LU factorization etc) that are optimized for sparse matrices.

Small detail in the index comparison. You need to do a lexicographical compare, otherwise:

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

Edit: So the comparison should probably be:

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 

Hash tables have a fast insertion and look up. You could write a simple hash function since you know you'd be dealing with only integer pairs as the keys.

The best way to implement sparse matrices is to not to implement them - atleast not on your own. I would suggest to BLAS (which I think is a part of LAPACK) which can handle really huge matrices.

Since only values with [a][b][c]...[w][x][y][z] are of consequence, we only store the indice themselves, not the value 1 which is just about everywhere - always the same + no way to hash it. Noting that the curse of dimensionality is present, suggest go with some established tool NIST or Boost, at least read the sources for that to circumvent needless blunder.

If the work needs to capture the temporal dependence distributions and parametric tendencies of unknown data sets, then a Map or B-Tree with uni-valued root is probably not practical. We can store only the indice themselves, hashed if ordering ( sensibility for presentation ) can subordinate to reduction of time domain at run-time, for all 1 values. Since non-zero values other than one are few, an obvious candidate for those is whatever data-structure you can find readily and understand. If the data set is truly vast-universe sized I suggest some sort of sliding window that manages file / disk / persistent-io yourself, moving portions of the data into scope as need be. ( writing code that you can understand ) If you are under commitment to provide actual solution to a working group, failure to do so leaves you at the mercy of consumer grade operating systems that have the sole goal of taking your lunch away from you.

Here is a relatively simple implementation that should provide a reasonable fast lookup (using a hash table) as well as fast iteration over non-zero elements in a row/column.

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

For simplicity, it's immutable, but you can can make it mutable; be sure to change std::vector to std::set if you want a reasonable efficient "insertions" (changing a zero to a non-zero).

I would suggest doing something like:

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; /* ... */ } 

To help keep your data sparse, you might want to write a subclass of unorderd_map, whose iterators automatically skip over (and erase) any items with a value of 0.