Comment calculer le vecteur propre d'une matrice stochastique de colonne en C ++

Question

Je une matrice stochastique colonne A et que vous voulez résoudre l'équation suivante en C ++: Ax = x

Je suppose que je dois savoir où x un vecteur propre valeur propre est fixé à être 1 (droit?), Mais je ne pouvais pas le comprendre en C ++. Jusqu'à présent, je l'ai vérifié quelques math-libs tels que Seldon, CPPScaLapack, Eigen ... Parmi eux, Eigen semble une bonne option, mais je ne pouvais pas comprendre comment utiliser l'un d'eux pour résoudre l'équation ci-dessus.

Pouvez-vous me donner quelques suggestions / extraits de code ou des idées pour résoudre l'équation? Toute aide est très appréciée.

Merci.

modifier:. A est n-par-n réelle, une matrice non négative

Answer 1

Gardez à l'esprit que je ne l'ai pas utilisé Eigen, mais son EigenSolver et SelfAdjointEigenSolver semblent être en mesure de résoudre ce problème. Ceux-ci sont répertoriés comme « expérimental », donc il pourrait y avoir des bugs et l'API pourrait changer à l'avenir.

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

Les deux classes de solveurs ont différents types pour la valeur eigen et eigen collections vecteur, mais comme ils sont tous les deux basés sur la classe Matrix et sont convertibles Matrices, qui précède, il devrait fonctionner.

Sinon, vous pouvez aborder le problème comme homogène équation linéaire (AI _n ) x = 0 , qui peut être résolu en convertissant AI _n à une matrice triangulaire supérieure. d'élimination gaussienne de se faire (bien que vous aurez besoin de sauter l'étape de normalisation pour chaque ligne où vous assurer le coefficient principal est 1, sous forme d'entiers ne sont pas un domaine). Une lecture rapide des projets ci-dessus ne se est pas un soutien pour la conversion de échelonnée, ce qui signifie sans doute je l'ai raté. Dans tous les cas, il est pas trop difficile à mettre en œuvre avec quelques classes d'aide (de RowMajor, RowMajor::iterator, RowWithPivot ci-après). Je n'ai même pas testé si oui ou non ce compilera, prenez donc comme plus d'une illustration de l'algorithme qu'un complet, goutte en solution. Bien que l'échantillon utilise les fonctions, il serait plus logique d'utiliser une classe (à 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;
}

Les interfaces pour les classes d'aide, qui ne sont pas vérifiés pour const-exactitude et d'autres détails nécessaires qui rendent le travail de 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();
};