내가 열 확률 행렬 A를 가지고 ++ C에서 다음 식 해결 할

Question

C에서 열 확률 매트릭스 ++의 고유 벡터 계산 방법 : 도끼 = X 나는이를 찾을 필요가 있다고 가정하고

을 eigenvector x eigenvalue는 1 (오른쪽?)으로 설정되어 있지만 C++에서는 이해할 수 없습니다. 지금까지 셀던 (Seldon), CPPScaLapack, 아이겐 (Eigen)과 같은 몇 가지 수학 라이브러리를 살펴 보았습니다. 그 중 Eigen은 좋은 옵션 인 것처럼 보이지만 위의 방정식을 풀 때 그 중 하나를 활용하는 방법을 이해할 수 없었습니다.

방정식을 해결하기위한 제안/코드 스 니펫이나 아이디어를 제공해 줄 수 있습니까? 도움을 주시면 감사하겠습니다.

감사합니다.

편집 : A는 nxn 실제, 음수가 아닌 행렬입니다.

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

두 솔버 클래스는 고유 값과 고유 벡터 컬렉션에 대한 다른 유형을 가지고 있지만, 그들이있는 한 두 클래스 매트릭스를 기반으로 행렬이 위의 작동합니다, 컨버터블이다.

또는,은 상 삼각 행렬 A-I _N 변환에 의해 해결 될 수 homogeneous linear equation (A-I _N) X = 0 같은 문제를 접근 할 수있다. Gaussian elimination은 그렇게 할 것입니다 (그러나 선행 계수가 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-correctness 및 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(); 
};