Wie der Eigenvektor einer Spalte stochastische Matrix in C++

Question

Ich habe eine Spalte stochastische Matrix A berechnen und will die folgende Gleichung in C++ lösen: Ax = x

ich davon aus bin, dass ich ein finden muß aus Eigenvektor x wobei Eigenwert auf 1 gesetzt ist (richtig?), aber ich konnte es in C++ nicht herausfinden. Bisher habe ich einige Math-Bibliotheken wie Seldon, CPPScaLapack, Eigen ... ausgecheckt. Unter ihnen scheint Eigen eine gute Option zu sein, aber ich konnte nicht verstehen, wie ich sie nutzen könnte, um die obige Gleichung zu lösen.

Können Sie mir einige Vorschläge/Code-Schnipsel oder Ideen geben, um die Gleichung zu lösen? Jede Hilfe wird sehr geschätzt.

Danke.

edit: A ist n-mal-n echte, nicht negative Matrix.

Answer 1

Denken Sie daran, ich habe nicht Eigen verwendet, aber seine EigenSolver und SelfAdjointEigenSolver aussehen, um dies zu lösen. Diese werden als "experimentell" aufgeführt, sodass möglicherweise Fehler auftreten und sich die API in Zukunft ändern kann.

// simple, not very efficient 
template <typename _M> 
bool isSelfAdjoint(const _M& a) { 
    return a == a.adjoint(); 
} 

template <typename _M> 
std::pair<Eigen::EigenSolver<_M>::EigenvalueType Eigen::EigenSolver<_M>::EigenvectorType> 
eigenvectors(const _M& a) { 
    if (isSelfAdjoint(a)) { 
     Eigen::EigenSolver<M> saes(a); 
     return pair(saes.eigenvalues(), saes.eigenvectors()); 
    } else { 
     Eigen::EigenSolver<M> es(a); 
     return pair(es.eigenvalues, es.eigenvectors()); 
    } 
} 

Die beiden Solver Klassen haben verschiedene Arten für den Eigenwert und Eigenvektor Sammlungen, aber da sie beide basierend auf der Klasse Matrix sind und Matrices umwandelbar sind, sollte die oben funktioniert.

Alternativ könnte man das Problem als homogeneous linear equation (A-I _n) x = , Ansatz, der durch Umwandeln eines I- _n zu einer oberen Dreiecksmatrix gelöst werden kann. Gaussian elimination wird das tun (obwohl Sie den Normalisierungsschritt für jede Zeile überspringen müssen, wo Sie sicherstellen, dass der führende Koeffizient 1 ist, da ganze Zahlen kein Feld sind). Eine schnelle Durchsicht der oben genannten Projekte ergab keine Unterstützung für die Konvertierung der Zeilenschaltung, was wahrscheinlich bedeutet, dass ich sie verpasst habe. In jedem Fall ist es nicht zu schwer mit wenigen Hilfsklassen zu implementieren (RowMajor, RowMajor::iterator, RowWithPivot im Folgenden). Ich habe nicht einmal getestet, ob dies kompiliert wird oder nicht, also nimm es eher als Beispiel für den Algorithmus als eine vollständige Drop-In-Lösung. Obwohl das Beispiel Funktionen verwendet, könnte es sinnvoller sein, eine Klasse (à EigenSolver) zu verwenden.

/* Finds a row with the lowest pivot index in a range of matrix rows. 
* Arguments: 
* - start: the first row to check 
* - end: row that ends search range (not included in search) 
* - pivot_i (optional): if a row with pivot index == pivot_i is found, search 
*  no more. Can speed things up if the pivot index of all rows in the range 
*  have a known lower bound. 
* 
* Returns an iterator p where p->pivot_i = min([start .. end-1]->pivot_i) 
* 
*/ 
template <typename _M> 
RowMajor<_M>::iterator 
find_lead_pivot (RowMajor<_M>::iterator start, 
       const RowMajor<_M>::iterator& end, 
       int pivot_i=0) 
{ 
    RowMajor<_M>::iterator lead=start; 
    for (; start != end; ++start) { 
     if (start->pivot() <= pivot_i) { 
      return start; 
     } 
     if (start->pivot() < lead->pivot()) { 
      lead = start; 
     } 
    } 
    return end; 
} 

/* Returns a matrix that's the row echelon form of the passed in matrix. 
*/ 
template <typename _M> 
_M form_of_echelon(const _M& a) { 
    _M a_1 = a-_M::Identity(); 
    RowMajor<_M> rma_1 = RowMajor<_M>(a_1); 
    typedef RowMajor<_M>::iterator RMIter; 
    RMIter lead; 
    int i=0; 

    /* 
     Loop invariant: row(i).pivot_i <= row(j).pivot_i, for j st. j>i 
    */ 
    for (RMIter row_i = rma_1.begin(); 
     row_i != rma_1.end() && row_i->pivot() != 0; 
     ++row_i, ++i) 
    { 
     lead = find_lead_pivot(row_i, rma_1.end(), i); 
     // ensure row(i) has minimal pivot index 
     swap(*lead, *row_i); 

     // ensure row(j).pivot_i > row(i).pivot_i 
     for (RMIter row_j = row_i+1; 
      row_j != rma_1.end(); 
      ++row_j) 
     { 
      *row_j = *row_j * row_i->pivot() - *row_i * row_j->pivot(); 
     } 
     /* the best we can do towards true row echelon form is reduce 
     * the leading coefficient by the row's GCD 
     */ 
     // *row_i /= gcd(*row_i); 
    } 
    return static_cast<_M>(rma_1); 
} 

/* Converts a matrix to echelon form in-place 
*/ 
template <typename _M> 
_M& form_of_echelon(_M& a) { 
    a -= _M::Identity(); 
    RowMajor<_M> rma_1 = RowMajor<_M>(a); 
    typedef RowMajor<_M>::iterator RMIter; 
    RMIter lead; 
    int i=0; 

    /* 
     Loop invariant: row(i).pivot_i <= row(j).pivot_i, for j st. j>i 
    */ 
    for (RMIter row_i = rma_1.begin(); 
     row_i != rma_1.end() && row_i->pivot() != 0; 
     ++row_i, ++i) 
    { 
     lead = find_lead_pivot(row_i, rma_1.end(), i); 
     // ensure row(i) has minimal pivot index 
     swap(*lead, *row_i); 

     for (RMIter row_j = row_i+1; 
      row_j != rma_1.end(); 
      ++row_j) 
     { 
      *row_j = *row_j * row_i->pivot() - *row_i * row_j->pivot(); 
     } 
     /* the best we can do towards true row echelon form is reduce 
     * the leading coefficient by the row's GCD 
     */ 
     // *row_i /= gcd(*row_i); 
    } 
    return a; 
} 

Die Schnittstellen für die Hilfsklassen, die für const-Korrektheit und andere notwendige Details überprüft haben, die C++ funktioniert nicht.

template <typename _M> 
class RowWithPivot { 
public: 
    typedef _M::RowXpr Wrapped; 
    typedef _M::Scalar Scalar; 

    RowWithPivot(_M& matrix, size_t row); 

    Wrapped base(); 
    operator Wrapped(); 

    void swap(RowWithPivot& other); 

    int cmp(RowWithPivot& other) const; 
    bool operator <(RowWithPivot& other) const; 

    // returns the index of the first non-zero scalar 
    // (best to cache this) 
    int pivot_index() const; 
    // returns first non-zero scalar, or 0 if none 
    Scalar pivot() const; 
}; 

template <typename _M, typename _R = RowWithPivot<_M> > 
class RowMajor { 
public: 
    typedef _R value_type; 

    RowMajor(_M& matrix); 

    operator _M&(); 
    _M& base(); 

    value_type operator[](size_t i); 

    class iterator { 
    public: 
     // standard random access iterator 
     ... 
    }; 

    iterator begin(); 
    iterator end(); 
};