Come calcolare l'autovettore di una matrice colonna stocastico in C ++

Question

Ho una colonna stocastico matrice A e vuole risolvere la seguente equazione in C ++: Ax = x

Io parto dal presupposto che ho bisogno di scoprire un autovettore x dove autovalore è impostato per essere 1 (giusto?), Ma non riuscivo a capirlo in C ++. Finora ho controllato alcuni di matematica-libs quali Seldon, CPPScaLapack, Eigen ... Tra di loro, Eigen sembra una buona opzione, ma non riuscivo a capire come utilizzare qualcuno di loro per risolvere l'equazione di cui sopra.

Puoi darmi alcuni frammenti di codice suggerimenti / o idee per risolvere l'equazione? Ogni aiuto è molto apprezzato.

Grazie.

modifica:. A è n-da-n real, matrice non negativo

Answer 1

Tenete a mente che non ho usato Eigen, ma la sua EigenSolver e SelfAdjointEigenSolver aspetto di essere in grado di risolvere questo. Questi sono elencati come "sperimentale", quindi ci potrebbe essere bug e l'API potrebbe cambiare in futuro.

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

Le due classi risolutore hanno tipi differenti per le collezioni eigen valore e eigen vettoriali, ma come sono entrambi basati sulla classe Matrix e matrici sono convertibili, quanto sopra dovrebbe funzionare.

In alternativa, si potrebbe affrontare il problema come omogeneo lineare equazione (AI _n ) x = 0 , che può essere risolto convertendo aI _n ad una matrice triangolare superiore. gaussiana eliminazione lo farà (anche se avrete bisogno di saltare la fase di normalizzazione per ogni riga in cui si assicurare il primo coefficiente è 1, come numeri interi non sono un campo). Un rapido esame dei progetti di cui sopra non si presentò supporto per la conversione fila Echelon, che probabilmente significa che ho perso. In ogni caso, non è troppo difficile da attuare con un paio di classi di supporto (RowMajor, RowMajor::iterator, RowWithPivot nel seguito). Non ho ancora testato se questo verrà compilato, in modo da prendere come più di un'illustrazione dell'algoritmo di una completa, drop-in soluzione. Anche se le funzioni usi di esempio, potrebbe avere più senso utilizzare una classe (à l'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;
}

Le interfacce per le classi di supporto, che non sono stati controllati per const correttezza e altri dettagli necessari che rendono il lavoro C ++.

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