<p>I am working on a project that requires the manipulation of enormous matrices, specifically pyramidal summation for a copula calculation.  </p>

<p>In short, I need to keep track of a relatively small number of values (usually a value of 1, and in rare cases more than 1) in a sea of zeros in the matrix (multidimensional array).  </p>

<p>A sparse array allows the user to store a small number of values, and assume all undefined records to be a preset value.  Since it is not physically possibly to store all values in memory, I need to store only the few non-zero elements. This could be several million entries.</p>

<p>Speed is a huge priority, and I would also like to dynamically choose the number of variables in the class at runtime.</p>

<p>I currently work on a system that uses a binary search tree (b-tree) to store entries. Does anyone know of a better system?</p>

[{'answer_id': 4311, 'author': 'Mark Harrison', 'author_id': 116, 'author_profile': 'https://Stackoverflow.com/users/116', 'pm_score': 6, 'selected': True, 'text': '<p>For C++, a map works well.  Several million objects won\'t be a problem.  10 million items took about 4.4 seconds and about 57 meg on my computer.</p>\n\n<p>My test application is as follows:</p>\n\n<pre><code>#include &lt;stdio.h&gt;\n#include &lt;stdlib.h&gt;\n#include &lt;map&gt;\n\nclass triple {\npublic:\n    int x;\n    int y;\n    int z;\n    bool operator&lt;(const triple &amp;other) const {\n        if (x &lt; other.x) return true;\n        if (other.x &lt; x) return false;\n        if (y &lt; other.y) return true;\n        if (other.y &lt; y) return false;\n        return z &lt; other.z;\n    }\n};\n\nint main(int, char**)\n{\n    std::map&lt;triple,int&gt; data;\n    triple point;\n    int i;\n\n    for (i = 0; i &lt; 10000000; ++i) {\n        point.x = rand();\n        point.y = rand();\n        point.z = rand();\n        //printf("%d %d %d %d\\n", i, point.x, point.y, point.z);\n        data[point] = i;\n    }\n    return 0;\n}\n</code></pre>\n\n<p>Now to dynamically choose the number of variables, the easiest solution is to represent <strong>index as a string</strong>, and then use string as a key for the map. For instance, an item located at [23][55] can be represented via "23,55" string. We can also extend this solution for higher dimensions; such as for three dimensions an arbitrary index  will look like "34,45,56". A simple implementation of this technique is as follows:</p>\n\n<pre><code>std::map data&lt;string,int&gt; data;\nchar ix[100];\n\nsprintf(ix, "%d,%d", x, y); // 2 vars\ndata[ix] = i;\n\nsprintf(ix, "%d,%d,%d", x, y, z); // 3 vars\ndata[ix] = i;\n</code></pre>\n'}, {'answer_id': 4354, 'author': 'nlucaroni', 'author_id': 157, 'author_profile': 'https://Stackoverflow.com/users/157', 'pm_score': 2, 'selected': False, 'text': "<p>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.</p>\n"}, {'answer_id': 5722, 'author': 'JSN', 'author_id': 108465, 'author_profile': 'https://Stackoverflow.com/users/108465', 'pm_score': 1, 'selected': False, 'text': '<p>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.</p>\n'}, {'answer_id': 17051, 'author': 'Mat Noguchi', 'author_id': 1799, 'author_profile': 'https://Stackoverflow.com/users/1799', 'pm_score': 2, 'selected': False, 'text': '<p>Small detail in the index comparison. You need to do a lexicographical compare, otherwise:</p>\n\n<pre><code>a= (1, 2, 1); b= (2, 1, 2);\n(a&lt;b) == (b&lt;a) is true, but b!=a\n</code></pre>\n\n<p>Edit: So the comparison should probably be:</p>\n\n<pre><code>return lhs.x&lt;rhs.x\n    ? true \n    : lhs.x==rhs.x \n        ? lhs.y&lt;rhs.y \n            ? true \n            : lhs.y==rhs.y\n                ? lhs.z&lt;rhs.z\n                : false\n        : false\n</code></pre>\n'}, {'answer_id': 21450, 'author': 'Nic Strong', 'author_id': 2281, 'author_profile': 'https://Stackoverflow.com/users/2281', 'pm_score': 4, 'selected': False, 'text': '<p>Boost has a templated implementation of BLAS called uBLAS that contains a sparse matrix.</p>\n\n<p><a href="https://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm" rel="nofollow noreferrer">https://www.boost.org/doc/libs/release/libs/numeric/ublas/doc/index.htm</a></p>\n'}, {'answer_id': 39052, 'author': 'Konrad Rudolph', 'author_id': 1968, 'author_profile': 'https://Stackoverflow.com/users/1968', 'pm_score': 5, 'selected': False, 'text': "<p>The accepted answer recommends using strings to represent multi-dimensional indices.</p>\n<p>However, constructing strings is needlessly wasteful for this. If the size isn’t known at compile time (and thus <code>std::tuple</code> doesn’t work), <code>std::vector</code> works well as an index, both with hash maps and ordered trees. For <code>std::map</code>, this is almost trivial:</p>\n<pre><code>#include &lt;vector&gt;\n#include &lt;map&gt;\n\nusing index_type = std::vector&lt;int&gt;;\n\ntemplate &lt;typename T&gt;\nusing sparse_array = std::map&lt;index_type, T&gt;;\n</code></pre>\n<p>For <code>std::unordered_map</code> (or similar hash table-based dictionaries) it’s slightly more work, since <code>std::vector</code> does not specialise <code>std::hash</code>:</p>\n<pre><code>#include &lt;vector&gt;\n#include &lt;unordered_map&gt;\n#include &lt;numeric&gt;\n\nusing index_type = std::vector&lt;int&gt;;\n\nstruct index_hash {\n    std::size_t operator()(index_type const&amp; i) const noexcept {\n        // Like boost::hash_combine; there might be some caveats, see\n        // &lt;https://stackoverflow.com/a/50978188/1968&gt;\n        auto const hash_combine = [](auto seed, auto x) {\n            return std::hash&lt;int&gt;()(x) + 0x9e3779b9 + (seed &lt;&lt; 6) + (seed &gt;&gt; 2);\n        };\n        return std::accumulate(i.begin() + 1, i.end(), i[0], hash_combine);\n    }\n};\n\ntemplate &lt;typename T&gt;\nusing sparse_array = std::unordered_map&lt;index_type, T, index_hash&gt;;\n</code></pre>\n<p>Either way, the usage is the same:</p>\n<pre><code>int main() {\n    using i = index_type;\n\n    auto x = sparse_array&lt;int&gt;();\n    x[i{1, 2, 3}] = 42;\n    x[i{4, 3, 2}] = 23;\n\n    std::cout &lt;&lt; x[i{1, 2, 3}] + x[i{4, 3, 2}] &lt;&lt; '\\n'; // 65\n}\n</code></pre>\n"}, {'answer_id': 1484067, 'author': 'Nicholas Jordan', 'author_id': 177505, 'author_profile': 'https://Stackoverflow.com/users/177505', 'pm_score': 0, 'selected': False, 'text': '<p>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. </p>\n\n<p>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.</p>\n'}, {'answer_id': 24370334, 'author': 'eold', 'author_id': 395744, 'author_profile': 'https://Stackoverflow.com/users/395744', 'pm_score': 0, 'selected': False, 'text': '<p>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.</p>\n\n<pre><code>// Copyright 2014 Leo Osvald\n//\n// Licensed under the Apache License, Version 2.0 (the "License");\n// you may not use this file except in compliance with the License.\n// You may obtain a copy of the License at\n//\n//     http://www.apache.org/licenses/LICENSE-2.0\n//\n// Unless required by applicable law or agreed to in writing, software\n// distributed under the License is distributed on an "AS IS" BASIS,\n// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n// See the License for the specific language governing permissions and\n// limitations under the License.\n\n#ifndef UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_\n#define UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_\n\n#include &lt;algorithm&gt;\n#include &lt;limits&gt;\n#include &lt;map&gt;\n#include &lt;type_traits&gt;\n#include &lt;unordered_map&gt;\n#include &lt;utility&gt;\n#include &lt;vector&gt;\n\n// A simple time-efficient implementation of an immutable sparse matrix\n// Provides efficient iteration of non-zero elements by rows/cols,\n// e.g. to iterate over a range [row_from, row_to) x [col_from, col_to):\n//   for (int row = row_from; row &lt; row_to; ++row) {\n//     for (auto col_range = sm.nonzero_col_range(row, col_from, col_to);\n//          col_range.first != col_range.second; ++col_range.first) {\n//       int col = *col_range.first;\n//       // use sm(row, col)\n//       ...\n//     }\ntemplate&lt;typename T = double, class Coord = int&gt;\nclass SparseMatrix {\n  struct PointHasher;\n  typedef std::map&lt; Coord, std::vector&lt;Coord&gt; &gt; NonZeroList;\n  typedef std::pair&lt;Coord, Coord&gt; Point;\n\n public:\n  typedef T ValueType;\n  typedef Coord CoordType;\n  typedef typename NonZeroList::mapped_type::const_iterator CoordIter;\n  typedef std::pair&lt;CoordIter, CoordIter&gt; CoordIterRange;\n\n  SparseMatrix() = default;\n\n  // Reads a matrix stored in MatrixMarket-like format, i.e.:\n  // &lt;num_rows&gt; &lt;num_cols&gt; &lt;num_entries&gt;\n  // &lt;row_1&gt; &lt;col_1&gt; &lt;val_1&gt;\n  // ...\n  // Note: the header (lines starting with \'%\' are ignored).\n  template&lt;class InputStream, size_t max_line_length = 1024&gt;\n  void Init(InputStream&amp; is) {\n    rows_.clear(), cols_.clear();\n    values_.clear();\n\n    // skip the header (lines beginning with \'%\', if any)\n    decltype(is.tellg()) offset = 0;\n    for (char buf[max_line_length + 1];\n         is.getline(buf, sizeof(buf)) &amp;&amp; buf[0] == \'%\'; )\n      offset = is.tellg();\n    is.seekg(offset);\n\n    size_t n;\n    is &gt;&gt; row_count_ &gt;&gt; col_count_ &gt;&gt; n;\n    values_.reserve(n);\n    while (n--) {\n      Coord row, col;\n      typename std::remove_cv&lt;T&gt;::type val;\n      is &gt;&gt; row &gt;&gt; col &gt;&gt; val;\n      values_[Point(--row, --col)] = val;\n      rows_[col].push_back(row);\n      cols_[row].push_back(col);\n    }\n    SortAndShrink(rows_);\n    SortAndShrink(cols_);\n  }\n\n  const T&amp; operator()(const Coord&amp; row, const Coord&amp; col) const {\n    static const T kZero = T();\n    auto it = values_.find(Point(row, col));\n    if (it != values_.end())\n      return it-&gt;second;\n    return kZero;\n  }\n\n  CoordIterRange\n  nonzero_col_range(Coord row, Coord col_from, Coord col_to) const {\n    CoordIterRange r;\n    GetRange(cols_, row, col_from, col_to, &amp;r);\n    return r;\n  }\n\n  CoordIterRange\n  nonzero_row_range(Coord col, Coord row_from, Coord row_to) const {\n    CoordIterRange r;\n    GetRange(rows_, col, row_from, row_to, &amp;r);\n    return r;\n  }\n\n  Coord row_count() const { return row_count_; }\n  Coord col_count() const { return col_count_; }\n  size_t nonzero_count() const { return values_.size(); }\n  size_t element_count() const { return size_t(row_count_) * col_count_; }\n\n private:\n  typedef std::unordered_map&lt;Point,\n                             typename std::remove_cv&lt;T&gt;::type,\n                             PointHasher&gt; ValueMap;\n\n  struct PointHasher {\n    size_t operator()(const Point&amp; p) const {\n      return p.first &lt;&lt; (std::numeric_limits&lt;Coord&gt;::digits &gt;&gt; 1) ^ p.second;\n    }\n  };\n\n  static void SortAndShrink(NonZeroList&amp; list) {\n    for (auto&amp; it : list) {\n      auto&amp; indices = it.second;\n      indices.shrink_to_fit();\n      std::sort(indices.begin(), indices.end());\n    }\n\n    // insert a sentinel vector to handle the case of all zeroes\n    if (list.empty())\n      list.emplace(Coord(), std::vector&lt;Coord&gt;(Coord()));\n  }\n\n  static void GetRange(const NonZeroList&amp; list, Coord i, Coord from, Coord to,\n                       CoordIterRange* r) {\n    auto lr = list.equal_range(i);\n    if (lr.first == lr.second) {\n      r-&gt;first = r-&gt;second = list.begin()-&gt;second.end();\n      return;\n    }\n\n    auto begin = lr.first-&gt;second.begin(), end = lr.first-&gt;second.end();\n    r-&gt;first = lower_bound(begin, end, from);\n    r-&gt;second = lower_bound(r-&gt;first, end, to);\n  }\n\n  ValueMap values_;\n  NonZeroList rows_, cols_;\n  Coord row_count_, col_count_;\n};\n\n#endif  /* UTIL_IMMUTABLE_SPARSE_MATRIX_HPP_ */\n</code></pre>\n\n<p>For simplicity, it\'s <code>immutable</code>, but you can can make it mutable; be sure to change <code>std::vector</code> to <code>std::set</code> if you want a reasonable efficient "insertions" (changing a zero to a non-zero).</p>\n'}, {'answer_id': 27696214, 'author': 'Emile Cormier', 'author_id': 245265, 'author_profile': 'https://Stackoverflow.com/users/245265', 'pm_score': 3, 'selected': False, 'text': '<p><a href="http://eigen.tuxfamily.org/" rel="noreferrer">Eigen</a> is a C++ linear algebra library that has an <a href="http://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html" rel="noreferrer">implementation</a> of a sparse matrix. It even supports matrix operations and solvers (LU factorization etc) that are optimized for sparse matrices.</p>\n'}, {'answer_id': 39989181, 'author': 'BenGoldberg', 'author_id': 3093194, 'author_profile': 'https://Stackoverflow.com/users/3093194', 'pm_score': 0, 'selected': False, 'text': '<p>I would suggest doing something like:</p>\n\n<pre><code>typedef std::tuple&lt;int, int, int&gt; coord_t;\ntypedef boost::hash&lt;coord_t&gt; coord_hash_t;\ntypedef std::unordered_map&lt;coord_hash_t, int, c_hash_t&gt; sparse_array_t;\n\nsparse_array_t the_data;\nthe_data[ { x, y, z } ] = 1; /* list-initialization is cool */\n\nfor( const auto&amp; element : the_data ) {\n    int xx, yy, zz, val;\n    std::tie( std::tie( xx, yy, zz ), val ) = element;\n    /* ... */\n}\n</code></pre>\n\n<p>To help keep your data sparse, you might want to write a subclass of <code>unorderd_map</code>, whose iterators automatically skip over (and erase) any items with a value of 0.</p>\n'}, {'answer_id': 41024354, 'author': 'Validus Oculus', 'author_id': 1516116, 'author_profile': 'https://Stackoverflow.com/users/1516116', 'pm_score': 3, 'selected': False, 'text': '<p>Complete list of solutions can be found in the wikipedia. For convenience, I have quoted relevant sections as follows.</p>\n\n<p><a href="https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_.28DOK.29" rel="noreferrer">https://en.wikipedia.org/wiki/Sparse_matrix#Dictionary_of_keys_.28DOK.29</a></p>\n\n<p><strong>Dictionary of keys (DOK)</strong></p>\n\n<blockquote>\n  <p>DOK consists of a dictionary that maps (row, column)-pairs to the\n  value of the elements. Elements that are missing from the dictionary\n  are taken to be zero. The format is good for incrementally\n  constructing a sparse matrix in random order, but poor for iterating\n  over non-zero values in lexicographical order. One typically\n  constructs a matrix in this format and then converts to another more\n  efficient format for processing.[1]</p>\n</blockquote>\n\n<p><strong>List of lists (LIL)</strong></p>\n\n<blockquote>\n  <p>LIL stores one list per row, with each entry containing the column\n  index and the value. Typically, these entries are kept sorted by\n  column index for faster lookup. This is another format good for\n  incremental matrix construction.[2]</p>\n</blockquote>\n\n<p><strong>Coordinate list (COO)</strong></p>\n\n<blockquote>\n  <p>COO stores a list of (row, column, value) tuples. Ideally, the entries\n  are sorted (by row index, then column index) to improve random access\n  times. This is another format which is good for incremental matrix\n  construction.[3]</p>\n</blockquote>\n\n<p><strong>Compressed sparse row (CSR, CRS or Yale format)</strong></p>\n\n<blockquote>\n  <p>The compressed sparse row (CSR) or compressed row storage (CRS) format\n  represents a matrix M by three (one-dimensional) arrays, that\n  respectively contain nonzero values, the extents of rows, and column\n  indices. It is similar to COO, but compresses the row indices, hence\n  the name. This format allows fast row access and matrix-vector\n  multiplications (Mx).</p>\n</blockquote>\n'}]
['https://Stackoverflow.com/questions/4306', 'https://Stackoverflow.com', 'https://Stackoverflow.com/users/522/']