Wie der Eigenvektor einer Spalte stochastische Matrix in C ++ berechnen

Question

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

Ich gehe davon aus, dass ich ein Eigenvektor, um herauszufinden, wo x Eigenwert 1 gesetzt ist (oder?), Aber ich konnte es nicht in C ++ herauszufinden. Bisher habe ich einige Mathe-Libs ausgecheckt wie Seldon, CPPScaLapack, Eigen ... Unter ihnen Eigen scheint eine gute Option, aber ich konnte nicht verstehen, wie einer von ihnen zu nutzen, über die Gleichung zu lösen.

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

Danke.

bearbeiten: a. Ist n-mal-n-real, nicht-negative Matrix

Answer 1

Beachten Sie habe ich nicht verwendet Eigen, aber seine EigenSolver und SelfAdjointEigenSolver Blick in der Lage sein, diese zu lösen. Diese sind aufgelistet als „experimentell“, so könnte es Fehler sein, und die API könnte sich in Zukunft ändern.

// 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 und Matrices sind umwandelbar sind, sollte die oben funktioniert.

Alternativ können Sie das Problem als homogene lineare Gleichung (AI _n ) x = 0 , die durch Umwandlung von AI _n , um eine obere Dreiecksmatrix gelöst werden kann. Gaußsche Eliminations wird das tun (obwohl Sie den Normalisierungsschritt für jede Zeile überspringen müssen, wo Sie sorgen für die Leitkoeffizient 1, als ganze Zahlen kein Feld sind). Eine schnelle Durchsicht der oben genannten Projekte aufdrehen nicht Unterstützung für Zeilenstufenumwandlung, was wahrscheinlich bedeutet, dass ich es verpasst. Auf jedem Fall ist es nicht zu hart mit einigen Helferklassen (RowMajor, RowMajor::iterator, RowWithPivot im folgenden) zu implementieren. Ich habe nicht einmal geprüft, ob diese kompiliert, nehmen so es eher als eine Darstellung des Algorithmus als eine vollständige, Drop-in-Lösung. Obwohl die Probe verwenden Funktionen, könnte es mehr Sinn machen, um eine Klasse zu verwenden (à la EigenSolver).

/* 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 nicht überprüft haben, die C ++ Arbeit zu machen.

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();
};