Matrix Addition effizient mit Sparse

[bucada]

Hi !

Ich habe noch nicht soviel Programmiererfahrung und habe nach einer Anleitung zur Sparse-Darstellung von Matrizen mit Hilfe von map, versucht diese umzusetzen.
Das ist dabei herausgekommen:

#ifndef SPARSE_H
#define SPARSE_H

#include <assert.h>
#include <map>

typedef map<int,double>::iterator iterator;
typedef map<int,double>::const_iterator const_iterator;

class Sparse : public Matrix
{
public:
  Sparse(int M = 0);
  Sparse(int M, int N, int Z = 0);
  Sparse(const Sparse& A);
  ~Sparse();

  double& operator() (int i, int j)       { return s[i][j]; }
  double  operator() (int i, int j) const { return s[i][j]; }
  map<int,double>& operator[] (int i) const { return s[i];  }

  const Sparse& operator= (const Sparse& A);

  Sparse* operator+ (const Sparse& A) const;
private:
  // Sizes of Matrix
  int n;
  int m;
  int z;
  // Matrix elements
  map<int,double> *s;
};
#endif

#include "sparse.h"
Sparse::Sparse(int M) : m(M) 
{ 
  s=new map<int,double>[m]; 
}

Sparse::Sparse(int M, int N, int Z): m(M), n(N), z(Z)
{ 
  if ((n != 0) && m!= 0 )
    if( n > m )
      s=new map<int,double>[n]; 
    else
      s=new map<int,double>[m]; 
}

Sparse::Sparse(const Sparse& A)
{
  s=new map<int,double>[z=A.size()];

  for (int i=0; i<z; i++) 
    s[i]=A[i];
}

Sparse::~Sparse()
{
  delete[] s;
}

const Sparse& Sparse::operator= (const Sparse& A)
{
    map<int,double>* tmp = new map<int,double>[z=A.size()];
    for (int i=0; i<z; i++)
      tmp[i]=A[i];

    delete [] s;
    m = A.m;
    n = A.n;
    z = A.z;
    s = tmp;
    return *this;
}

Sparse* Sparse::operator+ (const Sparse& A) const
{ 
  assert( (m == A.m) && (n == A.n) );

  Sparse* result;
  result = new Sparse(m, n, 0);
  for ( int i = 0; i < m; i++ ) 
    for ( int j = 0; j < n; j++ )
    {
      if((*this)(i,j) != 0. || A(i,j) != 0. )
        (*result)(i,j) = (*this)(i,j) + A(i,j); 
    }
  return result;
}

Doch hier ist die Addtion nicht effizient. Es wird nicht ausgenutzt, dass die Matrix dünnbesetzt ist, weil ich in der for-Schleife die kompletten Zeilen und Spalten durchgehe.

Was ich mir überlegt habe nachdem ich mir die STL "map" genauer angeschaut habe, so sollte irgendwie da der Zugriff mittels first und second möglich sein. So könnte ich dann die Schleife für die Anzahl der Elemente (also z) durchgehen und dann irgendwie nachschauen, welche Zeilen und Spalten zueinander passen und diese dann gegebenfalls addieren, oder so in result reinschreiben.
Doch so kriege ich das momentan nicht hin, weil ich bei der Kompilierung Fehler bekomme wenn ich folgende Zeile reinschriebe:
int a = (*this)->first;

"error: base operand of ->' has non-pointer typeconst Sparse'"

Kann mir da jemand behilflich sein, der sich mit map etwas auskennt?

[Neku]

Also Sparse kenne ich nur bei Dateien: http://de.wikipedia.org/wiki/Sparse

'this' ist ein Zeiger auf ein Objekt. Wenn du diesen dereferenzierst, ist es kein Zeiger mehr, sondern eine Referenz.

Das was du meinst ist vermutlich this->s->first (bzw. this->s->first(), wenn 'first' eine Methode ist - weiß ich gerade nicht genau).

Viellciht kann dir jemand besser helfen, wenn du Sparse erklärst oder die Anleitung verlinkst.

[bucada]

Mit einer Sparse Darstellung der Matrix ist folgendes gemeint:
Es werden nur die Elemente ungleich Null gespeichert.
Zudem muss dann zu jedem Element die Information mitgegeben werden, in welcher Zeile und in welcher Spalte sich das Element befindet.
Beispiel für eine Matrix:
2 0 3
0 0 1
0 1 0

gespeichert wird:

Werte: 2, 3, 1, 1
Zeile: 1, 1, 2, 3
Spalte:1, 3, 3, 2

Dabei stehen die zueinander gehörenden Zahlen hier untereinander.

[Neku]

Lohnt sich das wirklich?
Das dürfte doch nur effektiv sein, wenn man wirklich oft viele Nullen hat? Das Suchen nach Nullen bzw. das Zuordnen der Werte zu Zeilen/Spalten dürfte doch mehr Rechenzeit verbrauchen, als das Rechnen?

[bucada]

Die kleinste Matrix, die ich anwende ist: 2048x2048 mit 2047 Einträgen ungleich Null. Da lohnt sich das schon.
Die Addtion zwei solcher Matrizen braucht bei mir auf dem Rechner über 8 Sekunden.
Ist ja auch kein Wunder, ich gehe ja dort 4.194.304 mal in die for-Schleife.
Wenn ich die anders durchgehen könnte, so wäre das nur 2047 mal. Es ist die kleinste Matrix (andere Matrizen haben sicher eine Grösse von über 5000x5000)

Ich hab mir auch schon einen Algo überlegt, wie ich die Einträge suchen kann, was mir fehlt ist der Zugriff jeweils auf den Zeilenwert, Spaltenwert und Wert-Wert .

Ah ja, hab letztens noch den Link vergessen: http://physik.uni-graz.at/~uxh/lineare-algebra/
Ich meine dort das Kapitel 7 über Sparse-Matrizen
Es ist da alles schon ganz gut erklärt, doch meine Kenntnisse reichen leider nicht aus, um das Sinnvoll anwenden zu können

Vielleicht komme ich selber drauf, wenn ich alles besser verstehe. Wahrscheinlich hätte ich die folgende Frage reinschreiben sollen.
Dieser Ausdruck:

map<int,double> *s;

ist mir nicht ganz klar.
Ich habe dort eine map aus int und double Werten. Jedem int Wert ist ein double Wert zugeordnet.
So wie ich aus anderen Quellen weiss, dürfen in einer map im Key (hier int) nicht zweimal die selben Elemente vorkommen. Da habe ich schon das erste Problem, denn sowohl in Zeilen als auch in Spalten können Elemente vorkommen, die gleich sind. Das müsste dann hier etwas mit dem *s zu tun haben. Doch das ist mir unklar.

[CStoll]

Du könntest per set_union() eine Liste aller Positionen erzeugen, die in beiden Matrizen besetzt sind und dann nur die Elemente addieren, die in dieser Liste enthalten sind. Damit kämst du auf etwa 4000 Schleifendurchläufe für die eigentliche Addition (und nochmal genauso viel Zeit für den Aufbau der Indexliste).

[Braunstein]

Hallo,

Vielleicht kannst du auch mal hier reinschauen
Komponenten entwerfen mit der C++ STL, Ulrich Breymann
Komponenten entwerfen mit der C++ STL | ISBN: 3827314747
Kapitel 9.4 beschäftigt sich mit dünn besetzten Matrizen.

[bucada]

Könnte mir jemand diesen Ausdruck erklären?

map<int,double> *s;

Wie genau werde die Elemente dort gespeichert?

[CStoll]

Das ist ein Zeiger auf eine map, die irgendwo im Speicher liegt - in der map werden die Objekte als Binärbaum aus pair<const int,double> (sortiert nach den "const int"-Anteil) verwaltet.

[bucada]

Ja, soweit ist das klar, doch wie kommen hier die Zeilen und Spalten der Matrizen rein?

[CStoll]

Du initialisierst das mit "s=new map<...>[n];" im Konstruktor - d.h. der Zeiger verweist auf ein Array von n einzelnen map's (vermutlich eine pro Matrixzeile).

[Jester]

Neku schrieb:

Lohnt sich das wirklich?
Das dürfte doch nur effektiv sein, wenn man wirklich oft viele Nullen hat? Das Suchen nach Nullen bzw. das Zuordnen der Werte zu Zeilen/Spalten dürfte doch mehr Rechenzeit verbrauchen, als das Rechnen?

Stell Dir vor, Du hast ne 1000x1000-Matrix und da sind nur 5000 Werte ungleich 0 drin. Da ist es schon allein vom Speicher her effizienter, lieber die Zahlen ungleich 0 und wo sie stehen zu speichern, als trotzdem alle Zahlen aufzulisten. Und bei Rechenoperationen sparste auch nochmal ne ganze Menge.

[Neku]

Jester schrieb:

Neku schrieb:

Lohnt sich das wirklich?
Das dürfte doch nur effektiv sein, wenn man wirklich oft viele Nullen hat? Das Suchen nach Nullen bzw. das Zuordnen der Werte zu Zeilen/Spalten dürfte doch mehr Rechenzeit verbrauchen, als das Rechnen?

Stell Dir vor, Du hast ne 1000x1000-Matrix und da sind nur 5000 Werte ungleich 0 drin. Da ist es schon allein vom Speicher her effizienter, lieber die Zahlen ungleich 0 und wo sie stehen zu speichern, als trotzdem alle Zahlen aufzulisten. Und bei Rechenoperationen sparste auch nochmal ne ganze Menge.

An so große Matrizen habe ich gar nicht gedacht. Wozu braucht man die?

*mal rechnet*

Matrix: 1000x1000
Typ: double (~ 8 Byte)
Positionsspeicher: 2 ushort (x und y, ~ je 2 Byte)
Anzahl Elemente: uint (~ 4 Byte)

Normal: 1000 * 1000 * 8 = 8.000.000 Byte
Sparse (5000 Elemente != 0): 5000 * (8 + 4) + 4 = 60.004 Byte
Sparse (5000 Elemente == 0): 995.000 * (8 + 4) + 4 = 11.940.004 Byte

Stimmt, bis zu 666.666 Elementen != 0 (schöne Zahl :)) ist Sparse hier effizienter.

[camper]

den overhead von nodebasierten containern wie map sollte man dann aber nicht unterschlagen da map typischerweise auf einem red-black tree basiert, macht das zwei zeiger pro node, also nochmal jeweils 8 byte extra bei 32bit architektur. pi mal daumen würde ich schätzen, dass sich das ganze lohnt (bzgl. speicherplaz und rechenaufwand) wenn weniger als 10% der matrix besetzt sind. für richtige große matrizen ist eine unordered_map wahrscheinlich noch günstiger - eine hashfunktion dafür ist ja recht trivial.

CStoll schrieb:

Du initialisierst das mit "s=new map<...>[n];" im Konstruktor - d.h. der Zeiger verweist auf ein Array von n einzelnen map's (vermutlich eine pro Matrixzeile).

das macht wenig sinn. eine einzige map genügt, z.b. per:

struct Pos { unsigned short x,y; };
map<Pos,double> s;

mit geeignetem prädikat, oder man bildet gleich auf unsigneds ab per: index = zeile * spaltenzahl + spalte, was im grunde nichts anderes ist.

[Werner Salomon]

camper schrieb:

... oder man bildet gleich auf unsigneds ab per: index = zeile * spaltenzahl + spalte, was im grunde nichts anderes ist.

genau das würde ich auch vorschlagen. Weiter würde ich bucada empfehlen, keine Pointer als Member zu nutzen sondern die STL-Container unmittelbar zu aggregieren. Das spart viel Arbeit und verkleinert die Fehlerwahrscheinlichkeit.

Eine weitere Vereinfachung ist es statt eines

Sparse* operator+ (const Sparse& A) const;

lieber einen

Sparse& operator+=( const Sparse& A ) const;

zu schreiben und den fehlenden operator+ dann so zu erzeugen

Sparse operator+( Sparse a, const Sparse& b )
{
    return a += b;
}

Aus der Klasse Sparse wird dann etwa:

#include <cassert>
#include <iostream>
#include <iomanip>
#include <map>

class Sparse
{
public:
    typedef std::size_t size_type;

    Sparse( size_type zeilen, size_type spalten )
        : m_nZeilen( zeilen )
        , m_nSpalten( spalten )
    {}

    Sparse& operator+=( const Sparse& b )
    {
        assert( m_nZeilen == b.m_nZeilen );
        assert( m_nSpalten == b.m_nSpalten );
        std::map< size_type, double >::const_iterator i2 = b.m_data.begin();
        for( std::map< size_type, double >::iterator i1 = m_data.begin(), j = i1
            ; i1 != m_data.end() && i2 != b.m_data.end(); )
        {
            if( i1->first < i2->first )
                ++i1;   // b[i2] == 0, *i1 übernehmen
            else if( i1->first == i2->first )
            {
                // beide Werte sind != 0 -> also addieren
                i1->second += i2->second;
                ++i1;
                ++i2;
            }
            else // i2->first < i1->first
            {
                // *i1 == 0, *i2 übernehmen
                j = m_data.insert( j, *i2 );
                ++i2;
            }
        }
        m_data.insert( i2, b.m_data.end() );    // Rest hinzufügen
        return *this;
    }

private:
    size_type m_nZeilen;
    size_type m_nSpalten;
    std::map< size_type, double > m_data;
};

Gruß
Werner

[Werner Salomon]

Neku schrieb:

An so große Matrizen habe ich gar nicht gedacht. Wozu braucht man die?

Zum Beispiel hier für. Bei mittleren Problemen erhält man lineare Gleichungssysteme mit Matrizen von mehreren Tausend Elementen Kantenlänge. Dabei sind weniger als 1% der Werte != 0.

Gruß
Werner