Как вычислить собственный вектор стохастической матрицы столбца в C ++ - PullRequest
Новая
       60

Как вычислить собственный вектор стохастической матрицы столбца в C ++

1 голос
/ 05 декабря 2010

У меня есть стохастическая матрица столбца A, и я хочу решить следующее уравнение в C ++: Ax = x

Я предполагаю, что мне нужно найти собственный вектор x, где собственное значение установлено равным 1 (справа)?) но я не мог понять это в C ++.До сих пор я проверял некоторые математические библиотеки, такие как Seldon, CPPScaLapack, Eigen ... Среди них Eigen кажется хорошим вариантом, но я не мог понять, как использовать любую из них для решения приведенного выше уравнения.

Можете ли вы дать мне некоторые предложения / фрагменты кода или идеи для решения уравнения?Любая помощь высоко ценится.

Спасибо.

редактировать: A - это n-на-n действительная неотрицательная матрица.

1 Ответ

1 голос
/ 06 декабря 2010

Имейте в виду, что я не использовал Eigen, но его EigenSolver и SelfAdjointEigenSolver, похоже, могут решить эту проблему.Они перечислены как «экспериментальные», поэтому могут быть ошибки, и API может измениться в будущем. 

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

Два класса решателей имеют разные типы для собственных значений и коллекций собственных векторов, но поскольку ониВы оба на основе класса Matrix и Matrix являются конвертируемыми, тем выше это должно работать.

В качестве альтернативы, вы можете подойти к задаче как однородное линейное уравнение (AI n ) x = 0 , что можно решить путем преобразования AI n в верхнюю треугольную матрицу. Гауссово исключение сделает это (хотя вам нужно пропустить шаг нормализации для каждой строки, где вы убедитесь, что ведущий коэффициент равен 1, поскольку целые числа не являются полем).Быстрое ознакомление с вышеупомянутыми проектами не привело к поддержке преобразования эшелона строк, что, вероятно, означает, что я пропустил это.В любом случае, это не так сложно реализовать с помощью нескольких вспомогательных классов (RowMajor, RowMajor::iterator, RowWithPivot в следующем).Я даже не проверял, будет ли это компилироваться или нет, поэтому воспринимайте его как иллюстрацию алгоритма, а не как полное решение.Хотя в примере используются функции, возможно, имеет смысл использовать класс (например, 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;
}

Интерфейсы для вспомогательных классов, которые не были проверены на предмет правильности констант и других необходимыхдетали, которые заставляют 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();
};
...