Matrixklasse - kleiner Fehler

[RobertB]

cooky451 schrieb:

<vector> nicht #included. https://ideone.com/p6RKY0
Und meinen Beitrag oben nicht gelesen. Bei der Matrixgröße kannst du diesen Ansatz sicher vergessen.

Deinen Beitrag habe ich gelesen. Der Ansatz war nur aus Interesse. Aber dieser ist dann allen Anschein nicht geeignet.

[Werner Salomon]

Hallo RobertB,

bei der angestrebten Größe der Matrix kommst Du nur mit einer Sparse Matrix weiter - vorausgesetzt natürlich, Dein Problem, welches Du lösen willst, erzeugt auch nur dünn besetzte Matrizen.
boost.ublas bietet sparse matrix an. Ich selber habe aber keine Erfahrungen damit.

Wenn Du an einer einfach gestrickten Lösung interessiert bist, so kannst Du Dir auch selber was bauen. Im Prinzip lässt sich so eine Matrix mit einem vector von maps lösen. So wäre zumindest die Implementierung relativ straight forward.

#include <iostream>
#include <map>
#include <vector>
#include <cassert>
#include <iomanip> // setw
#include <numeric> // accumulate

template< typename T >
class SparseMatrix
{
    typedef std::map< std::size_t, T > Row; // key -> column
    typedef std::vector< Row > Field; // index -> row

    struct Ref
    {
        Ref( Row& row, std::size_t columnIdx )
            : row_( &row )
            , column_( columnIdx )
        {}
        Ref& operator=( const T& element )
        {
            if( element == T() )
            {
                auto iElement = row_->find( column_ );
                if( iElement != end(*row_) )
                    row_->erase( iElement );
            }
            else
                (*row_)[column_] = element;
            return *this;
        }
        operator T() const
        {
            auto iElement = row_->find( column_ );
            return iElement == end(*row_)? T(): iElement->second;
        }
    private:
        Row* row_;
        std::size_t column_;
    };

public:
    SparseMatrix( std::size_t rows, std::size_t columns )
        : columns_( columns )
        , field_( rows )
    {}
    std::size_t Columns() const { return columns_; }
    std::size_t Rows() const { return field_.size(); }

    // --  Element Zugriff
    T operator()( std::size_t rowIdx, std::size_t columnIdx ) const
    {
        assert( rowIdx < Rows() );
        assert( columnIdx < Columns() );
        auto iElement = field_[rowIdx].find( columnIdx );
        return iElement == end(field_[rowIdx])? T(): iElement->second;
    }
    Ref operator()( std::size_t rowIdx, std::size_t columnIdx )
    {
        assert( rowIdx < Rows() );
        assert( columnIdx < Columns() );
        return Ref( field_[rowIdx], columnIdx );
    }

    // --  Arithmetik
    SparseMatrix& operator+=( const SparseMatrix& b );
    SparseMatrix& operator-=( const SparseMatrix& b );

    // --  Ausgabe
    template< typename E, typename Traits > friend
        std::basic_ostream< E, Traits >& operator<<( std::basic_ostream< E, Traits >& out, const SparseMatrix& mx )
    {
        std::streamsize width = out.width();
        if( width == 0 )
            width = 4;
        for( auto iRow = begin(mx.field_); iRow != end(mx.field_); ++iRow )
        {
            auto iElement = begin(*iRow);
            for( std::size_t columnIdx = 0; columnIdx < mx.columns_; ++columnIdx )
            {
                out << " " << std::setw( width );
                if( iElement != end(*iRow) && iElement->first == columnIdx )
                {
                    out << iElement->second;
                    ++iElement;
                }
                else
                    out << T();
            }
            out << "\n";
        }
        return out << std::flush;
    }

    std::size_t NoElements() const // Anzahl belegter Elemente; für Debugzwecke
    {
        return std::accumulate( begin(field_), end(field_), std::size_t(0),
            []( std::size_t sum, const Row& row )-> std::size_t { return sum + row.size(); } );
    }

private:
    std::size_t columns_;
    Field field_;
};

Was bei dieser Matrixgröße ein No-GO ist, wenn Du jetzt zwei Schleifen baust, die über alle Zeilen und Spalte laufen und dann irgendwelche Operationen ausführen. Bei einer Kantenlänge von 240.000 Elementen einer quadratischen Matrix schlägt jede NanoSekunde der Operation mit einer knappen Minute Laufzeit zu Buche.
D.h. man muss sich auch bei der Operation die dünne Besetzung der Matrix zu Nutze machen.

Als Demonstration habe ich die Addition zweier Matrizen implementiert:

template< typename T >
SparseMatrix< T >& SparseMatrix< T >::operator+=( const SparseMatrix< T >& b ) {
    assert( Rows() == b.Rows() );
    assert( Columns() == b.Columns() );
    auto iRowDst = begin(field_);
    for( auto iRow = begin(b.field_); iRow != end(b.field_); ++iRow, ++iRowDst )
    {
        auto iElementDst = begin(*iRowDst);
        for( auto iElement = begin(*iRow); iElement != end(*iRow); ++iElement )
        {
            for( ; iElementDst != end(*iRowDst) && iElementDst->first < iElement->first; ++iElementDst )
                ;
            if( iElementDst != end(*iRowDst) && iElementDst->first == iElement->first )
            {
                if( (iElementDst->second += iElement->second) == T() ) // Auslöschung beachten
                    iElementDst = iRowDst->erase( iElementDst );
                else
                    ++iElementDst;
            }
            else
                iRowDst->insert( iElementDst, *iElement );
        }
    }
    return *this;
}

Die Subtraktion ist identisch - nur in Zeile 15 das += durch ein -= ersetzen.

Bei mir (VS10) benötigt eine leere map<size_t,double> 20 Bytes Speicher. Damit hat eine leere SparseMatrix mit 240.000 Zeilen einem Speicherbedarf von knapp 5MByte. Für jedes Element kommen dann in der Größenordnung 20Byte hinzu (ein Knoten in der Map).
Wahrscheinlich ist es auch schlauer die Map durch eine sortierte Liste von pairstd::size_t,T (evt. forward_list) zu ersetzen. Aber mit der Map war es einfacher zu codieren.

Abschließend eine Demo, die die Funktion zeigt.

int main()
{
    using namespace std;
    SparseMatrix< double > mx( 8, 5 );
    mx( 0, 1 ) = 3.4;
    mx( 7, 4 ) = 99.;
    cout << fixed << setprecision( 1 );
    cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl;
    for( std::size_t c = 0; c < mx.Columns(); ++c )
        mx( 6, c ) = (c+1)*1.1;
    cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl;
    mx( 6, 1 ) = 0.;
    cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl;
    cout << "Der Wert von mx(6,2) ist: " << mx(6,2) << endl;

    SparseMatrix< double > m2( 8, 5 );
    m2(2,0) = 2.4;
    m2(6,2) = 3.3;
    m2(6,4) = -5.5;
    mx += m2;
    cout << "\n" << mx << "Anzahl Elemente !=0: " << mx.NoElements() << endl;

    return 0;
}

Wie viele Elemente!=0 enthält Deine Matrix so in etwa? Und welche Operationen benötigst Du für die Matrix?

Gruß
Werner

[RobertB]

Wow, ich bin beeindruckt vom Code! Vielen Dank!

Dünn besetze Matrizen habe ich. Muss dann noch ein Gleichungssystem lösen. Aber dafür werde ich wahrscheinlich Matlab nutzen und das rüberziehen. Oder ich benutze mit dem Mex-Interface eine Build-in Funktion von Matlab in C++ für das Lösen.

Weitere Operationen brauche ich eig. nicht. Es geht dann doch hauptsächlich um das Füllen und Bearbeiten einzelner Einträge dieser großen Matrix. Das geht ja alles mit dem Code, vielen Dank. Ich muss diesen nur noch ganz verstehen :).

Und leider brauche ich eine ganz große Schleife. Die erste geht nämlich über die Dimension der Länge bzw. Breite (habe eine quadratische Matrix) der Matrix, die zweite dann von 1 bis 3. Also im schlimmsten Fall habe ich 240.000*3 Durchläufe.

[ScottZhang]

Wenn du nen matlab interface hast bietet sich Werners Implementierung aus Performance Gründen nicht an. Soweit ich weiss brauchst du für Matlab alles in einem zusammenhängenden Speicher. D.h. ein haufen Kopieroperation werden fällig.
(nanosekunden und so ... :D)

Google einfach mal nach "Compressed row format" bzw. "Compressed column format".

[cooky451]

Du hast uns immer noch nicht verraten warum du C++ benutzen willst wenn du doch matlab hast?

[RobertB]

cooky451 schrieb:

Du hast uns immer noch nicht verraten warum du C++ benutzen willst wenn du doch matlab hast?

Aus Performancegründen soll ich C++ benutzen. Man sagte mir, dass wegen den Schleifen C++ schneller ist. Aber wegen den Kopiervorgängen habe ich mir auch schon "Sorgen" gemacht, wie das beim MEX-Interface abläuft.

Ich soll das Programm so schnell wie möglich machen, das ist das Ziel. Ich würde lieber bei Matlab bleiben :).

[cooky451]

Mein Tipp: Bleib einfach bei Matlab. Nicht, dass ich groß Ahnung davon hätte, aber ich weiß ein paar Dinge:
- Die Leute die Matlab implementiert haben, wissen bestimmt was sie tun.
- Du hast keine Ahnung von C++.
-> Es dauert wahrscheinlich mindestens einige Wochen bis du gut genug bist, um Matlab performancemäßig zu schlagen, auch mit Bibliotheken die für dich das Rechnen übernehmen.

Die zeit ist recht grob geschätzt, da ich weder dich noch Matlab gut genug kenne, aber das sollte man trotzdem mal erwähnen. Abgesehen davon: Wenn du wirklich, wirklich Performance brauchst, warum lädst du dann wieder Sachen nach Matlab? Das bringt doch nichts. Wenn, dann rechne gleiche alles in C++ aus. Nimm Eigen oder MTL und lass die möglichst alles rechnen, was die können. Denn deren Algorithmen und Implementierungen sind mit an Sicherheit grenzender Wahrscheinlichkeit deutlich schneller als alles was du oder wir hier mal eben zusammenbasteln könnten.

[RobertB]

Das mit Matlab entscheide leider nicht ich allein, sondern meine "Chefs". Aber Danke für Deinen Rat.

Denn, ich habe mir etwas überlegt: ich schmeiße ganz am Anfang in Matlab erstellte Daten in C++ (mit build-in Funktionen von Matlab), und am Ende hole ich mir das Ergebnis und lasse es in Matlab plotten. Dann habe ich nur die beiden Schnittstellen.

Und, was hälts Du davon, dass ich mein ganzes Matlabprogramm am Ende einfach mit einem Coder in C++-Code übersetze? Matlab bietet dafür eigens eine Lösung. Ob das dann schneller ist?

Auf jeden Fall ziehe ich den Hut vor euch C++-Programmieren! Das ist doch die wahre Kunst.

[ScottZhang]

Das wollt ich gerade vorschlagen. Das moderne Matlab kann C++ Code exportieren (vorausgesetzt ihr habt die entsprechenden Lizenzen).

Ich habe zwar nur seeehhr wenig ahnung vom mex-Interface. Aber könntest du nicht rauskriegen wie Matlab die Sparse-Matrizen implementiert (ich tippe auf eine von den von mir vorgeschlagenden Formaten) und einfach für diese Implementierung die entsprechend Funktionen schreiben die du brauchst (Addition , Multiplikation, Assemblierung und was sonst noch anfällt. Bei den ersten beiden würd ich ma vermuten das sind schon build-ins).

[Werner Salomon]

ScottZhang schrieb:

Wenn du nen matlab interface hast bietet sich Werners Implementierung aus Performance Gründen nicht an. Soweit ich weiss brauchst du für Matlab alles in einem zusammenhängenden Speicher. D.h. ein haufen Kopieroperation werden fällig.
(nanosekunden und so ... :D)

Google einfach mal nach "Compressed row format" bzw. "Compressed column format".

Dass es dafür auch Lösungen gibt, die nur mit Arrays und Indizes funktionieren war ja klar; die gute alte Fortran-Welt lässt grüßen. Aber dass es eine Datenstruktur als Quasi-Standard bis in die Wikipedia geschafft; damit hätte ich jetzt nicht gerechnet.
Ich habe nicht herausbekommen, ob MatLab die Sparse matrix im CRS abspeichert, aber ich habe mir mal die Mühe gemacht, meine class SparseMatrix auf dieses Format umzustellen.
Es ist nicht wirklich mehr Code geworden, aber es war schwieriger zu implementieren. Ich habe bestimmt mehr als die doppelte Zeit gebraucht wie für die std::map-Lösung.

Keine Ahnung, ob das anderen auch so geht, aber mit den STL-Container kann man auf einem höheren Abstraktionsniveau codieren, und das geht bei mir wesentlich flüssiger.

Die class SparseMatrixc im CRS Datenformat:

#include <algorithm> // for_each, find_if,..
#include <cassert>
#include <iterator> // distance
#include <ostream>
#include <vector>
#include <utility> // std::pair<>

template< typename T >
class SparseMatrix
{
    struct Ref
    {
        Ref( SparseMatrix& mx, std::size_t rowIdx, std::size_t columnIdx )
            : mx_( &mx )
            , rowIdx_( rowIdx )
            , columnIdx_( columnIdx )
        {}
        Ref& operator=( const T& val )
        {
            auto element = mx_->find_element( rowIdx_, columnIdx_ );
            if( element.second )
            {
                if( (mx_->val_[element.first] = val) == T() )
                    mx_->erase_( rowIdx_, element.first );
            }
            else if( val != T() )
                mx_->insert_( rowIdx_, columnIdx_, element.first, val );
            return *this;
        }
        operator T() const
        {
            auto element = mx_->find_element( rowIdx_, columnIdx_ );
            return element.second? mx_->val_[element.first]: T();
        }
    private:
        SparseMatrix* mx_;
        std::size_t rowIdx_;
        std::size_t columnIdx_;
    };
    friend struct Ref;

public:
    SparseMatrix( std::size_t rows, std::size_t columns )
        : columns_( columns )
        , val_()
        , colInd_()
        , rowPtr_( rows+1 )
    {}
    std::size_t Columns() const { return columns_; }
    std::size_t Rows() const { return rowPtr_.size()-1; }

    // --  Element Zugriff
    T operator()( std::size_t rowIdx, std::size_t columnIdx ) const
    {
        assert( rowIdx < Rows() );
        assert( columnIdx < Columns() );
        auto element = find_element( rowIdx, columnIdx );
        return element->second? val_[element->first]: T();
    }
    Ref operator()( std::size_t rowIdx, std::size_t columnIdx )
    {
        assert( rowIdx < Rows() );
        assert( columnIdx < Columns() );
        return Ref( *this, rowIdx, columnIdx );
    }

    // --  Ausgabe
    template< typename E, typename Traits > friend
        std::basic_ostream< E, Traits >& operator<<( std::basic_ostream< E, Traits >& out, const SparseMatrix& mx )
    {
        std::streamsize width = out.width();
        if( width == 0 )
            width = 4;
        std::size_t i = 0; // Index in mx.colInd_[] & mx.val_[]
        for( auto iRow = begin(mx.rowPtr_); iRow != end(mx.rowPtr_)-1; ++iRow )
        {
            for( std::size_t columnIdx = 0; columnIdx < mx.columns_; ++columnIdx )
            {
                (out << " ").width( width );
                if( i < *(iRow+1) && mx.colInd_[i] == columnIdx )
                    out << mx.val_[i++];
                else
                    out << T();
            }
            out << "\n";
        }
        return out << std::flush;
    }

    std::size_t NoElements() const // Anzahl belegter Elemente; für Debugzwecke
    {
        return val_.size();
    }

private:
    void insert_( std::size_t rowIdx, std::size_t columnIdx, std::size_t i, const T& val )
    {
        assert( rowPtr_[rowIdx] <= i );
        assert( i <= rowPtr_[rowIdx+1] );
        assert( std::none_of( begin(colInd_) + rowPtr_[rowIdx], begin(colInd_) + rowPtr_[rowIdx+1],
            [columnIdx]( std::size_t c )->bool { return c == columnIdx; } ) );
        if( val != T() )
        {
            val_.insert( begin(val_) + i, val );
            colInd_.insert( begin(colInd_) + i, columnIdx );
            std::for_each( begin(rowPtr_) + rowIdx+1, end(rowPtr_), []( std::size_t& idx ) { ++idx; } );
        }
    }
    void erase_( std::size_t rowIdx, std::size_t i )
    {
        assert( i < val_.size() );
        assert( rowPtr_[rowIdx] <= i );
        assert( i < rowPtr_[rowIdx+1] );
        val_.erase( begin(val_) + i );
        colInd_.erase( begin(colInd_) + i );
        std::for_each( begin(rowPtr_) + rowIdx+1, end(rowPtr_), []( std::size_t& idx ) { --idx; } );
    }
    std::pair< std::size_t, bool > find_element( std::size_t rowIdx, std::size_t columnIdx ) const
    {
        auto iColumn = std::find_if( begin(colInd_) + rowPtr_[rowIdx], begin(colInd_) + rowPtr_[rowIdx+1],
            [columnIdx]( std::size_t c )->bool { return columnIdx <= c; } );
        return std::make_pair( std::size_t( std::distance( begin(colInd_), iColumn ) )
            , iColumn != begin(colInd_) + rowPtr_[rowIdx+1] && *iColumn == columnIdx );
    }

    // --   Members
    std::size_t columns_;
    // --   Datenstruktur der Sparse Matrix im CRS -> <http://de.wikipedia.org/wiki/Compressed_Row_Storage>
    std::vector< T > val_;
    std::vector< std::size_t > colInd_;
    std::vector< std::size_t > rowPtr_;
};

Der Speicherbedarf dürfte niedriger sein (ca. bei 20%-30% der vorherigen Lösung) dafür ist das Einfügen neuer Elemente nicht mehr performant. Zumindest nicht wenn man sie einzeln zuweist.

.. es fehlt noch die Addition. Wenn ihr das Demo-main aus meinem vorherigen Posting benutzt, so lasst die Addition einfach weg, oder schreibt Euch eine - viel Spaß damit.

Gruß
Werner

[ScottZhang]

Sieht wohl so aus als würde Matlab die columns compressen: Matlab: mxCreateSparse

[otze]

CSR ist auch in ublas implmentiert. das Problem an CSR ist wie gesagt das Einfügen bei matrix-vektor und matrix-matrix produkten.

Ein netter Trick dafür ist etwas mehr Speicher als benötigt zu reservieren und die Elemente bei der Berechnung des Produkts am Ende unsortiert einzufügen. Am Schluss wird dann gemerged.