Pergunta

Estou trabalhando em um projeto que requer a manipulação de matrizes enormes, especificamente somatório piramidal para cálculo de cópula.

Resumindo, preciso acompanhar um número relativamente pequeno de valores (geralmente um valor de 1 e, em casos raros, mais de 1) em um mar de zeros na matriz (matriz multidimensional).

Uma matriz esparsa permite ao usuário armazenar um pequeno número de valores e assumir que todos os registros indefinidos são um valor predefinido.Como não é fisicamente possível armazenar todos os valores na memória, preciso armazenar apenas alguns elementos diferentes de zero.Isso pode representar vários milhões de entradas.

A velocidade é uma grande prioridade e eu também gostaria de escolher dinamicamente o número de variáveis ​​na classe em tempo de execução.

Atualmente trabalho em um sistema que usa uma árvore de pesquisa binária (árvore b) para armazenar entradas.Alguém conhece um sistema melhor?

Foi útil?

Solução

Para C++, um mapa funciona bem.Vários milhões de objetos não serão um problema.10 milhões de itens demoraram cerca de 4,4 segundos e cerca de 57 megas no meu computador.

Minha aplicação de teste é a seguinte:

#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;
}

Agora, para escolher dinamicamente o número de variáveis, a solução mais fácil é representar índice como uma string, e, em seguida, use string como chave para o mapa.Por exemplo, um item localizado em [23][55] pode ser representado pela string "23,55".Também podemos estender esta solução para dimensões superiores;tal como para três dimensões, um índice arbitrário será semelhante a "34,45,56".Uma implementação simples desta técnica é a seguinte:

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;

Outras dicas

Como conselho geral, um método que usa strings como índices é na verdade muito lento.Uma solução muito mais eficiente, mas equivalente, seria usar vetores/matrizes.Não há absolutamente nenhuma necessidade de escrever os índices em uma string.

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;

No entanto, usando um map na prática, muitas vezes não é muito eficiente devido à implementação em termos de uma árvore de pesquisa binária balanceada.Uma estrutura de dados com melhor desempenho neste caso seria uma tabela hash, conforme fornecido por std::unordered_map.

Boost tem uma implementação modelada de BLAS chamada uBLAS que contém uma matriz esparsa.

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

Pequeno detalhe na comparação do índice.Você precisa fazer uma comparação lexicográfica, caso contrário:

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

Editar:Portanto, a comparação provavelmente deveria ser:

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

As tabelas hash têm uma inserção e consulta rápida.Você poderia escrever uma função hash simples, pois sabe que estaria lidando apenas com pares inteiros como chaves.

Eigen é uma biblioteca de álgebra linear C++ que possui um implementação de uma matriz esparsa.Ele ainda suporta operações de matrizes e solucionadores (fatoração LU, etc.) que são otimizados para matrizes esparsas.

A lista completa de soluções pode ser encontrada na Wikipedia.Por conveniência, citei as seções relevantes a seguir.

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

Dicionário de chaves (DOK)

A DOK consiste em um dicionário que mapeia (linha, coluna)-suba para o valor dos elementos.Os elementos que estão faltando no dicionário são considerados zero.O formato é bom para a construção de uma matriz esparsa em ordem aleatória, mas ruim para iterar sobre valores diferentes de zero na ordem lexicográfica.Normalmente, constrói uma matriz neste formato e depois se converte para outro formato mais eficiente para o processamento. [1

Lista de listas (LIL)

O LIL armazena uma lista por linha, com cada entrada contendo o índice da coluna e o valor.Normalmente, essas entradas são mantidas classificadas pelo índice de coluna para uma pesquisa mais rápida.Este é outro formato bom para a construção da matriz incremental. [2

Lista de coordenadas (COO)

COO armazena uma lista de tuplas (linha, coluna, valor).Idealmente, as entradas são classificadas (por índice de linha, depois Índice de coluna) para melhorar os tempos de acesso aleatório.Este é outro formato que é bom para a construção da matriz incremental. [3

Linha esparsa compactada (formato CSR, CRS ou Yale)

A linha esparsa compactada (CSR) ou o formato de armazenamento de linha comprimido (CRS) representa uma matriz M por três matrizes (unidimensionais), que contêm respectivamente valores de zero, as extensões das linhas e os índices de coluna.É semelhante ao COO, mas comprime os índices da linha, daí o nome.Este formato permite acesso rápido à linha e multiplicações de vetor de matriz (MX).

A melhor maneira de implementar matrizes esparsas é não implementá-las - pelo menos não por conta própria.Eu sugeriria ao BLAS (que acho que faz parte do LAPACK), que pode lidar com matrizes realmente grandes.

Como apenas valores com [a][b][c]...[w][x][y][z] são importantes, armazenamos apenas os próprios índices, não o valor 1 que está em quase todos os lugares - sempre o mesmo + não há como fazer hash.Observando que a maldição da dimensionalidade está presente, sugiro usar alguma ferramenta estabelecida NIST ou Boost, pelo menos ler as fontes para evitar erros desnecessários.

Se o trabalho precisar capturar as distribuições de dependência temporal e tendências paramétricas de conjuntos de dados desconhecidos, então um Mapa ou Árvore B com raiz de valor único provavelmente não será prático.Podemos armazenar apenas os próprios índices, com hash se a ordenação (sensibilidade para apresentação) puder subordinar à redução do domínio do tempo em tempo de execução, para todos os valores 1.Como os valores diferentes de zero além de um são poucos, um candidato óbvio para eles é qualquer estrutura de dados que você possa encontrar prontamente e compreender.Se o conjunto de dados for realmente do tamanho de um universo vasto, sugiro algum tipo de janela deslizante que gerencie você mesmo o arquivo/disco/io persistente, movendo partes dos dados para o escopo conforme necessário.(escrever código que você possa entender) Se você está comprometido em fornecer uma solução real para um grupo de trabalho, deixar de fazê-lo o deixa à mercê de sistemas operacionais de consumo que têm o único objetivo de tirar seu almoço de você.

Aqui está uma implementação relativamente simples que deve fornecer uma pesquisa razoavelmente rápida (usando uma tabela hash), bem como uma iteração rápida sobre elementos diferentes de zero em uma linha/coluna.

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

Para simplificar, é immutable, mas você pode torná-lo mutável;certifique-se de mudar std::vector para std::set se você deseja "inserções" razoavelmente eficientes (alterando zero para diferente de zero).

Eu sugeriria fazer algo como:

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

Para ajudar a manter seus dados esparsos, você pode querer escrever uma subclasse de unorderd_map, cujos iteradores ignoram (e apagam) automaticamente quaisquer itens com valor 0.

Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top