Как вычислить собственный вектор стохастической матрицы столбца в C++
У меня есть стохастическая матрица столбцов A и я хочу решить следующее уравнение в C++: Ax=x
Я предполагаю, что мне нужно найти собственный вектор x, где собственное значение установлено равным 1(верно?), Но я не смог понять это в C++. До сих пор я проверял некоторые математические библиотеки, такие как Seldon, CPPScaLapack, Eigen... Среди них Eigen кажется хорошим вариантом, но я не мог понять, как использовать любую из них для решения уравнения выше.
Можете ли вы дать мне некоторые предложения / фрагменты кода или идеи для решения уравнения? Любая помощь высоко ценится.
Благодарю.
редактировать: A - это n-на-n вещественная неотрицательная матрица.
1 ответ
Имейте в виду, я не использовал 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;
}
Интерфейсы для вспомогательных классов, которые не были проверены на предмет правильности const и других необходимых деталей, которые заставляют 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();
};