Eigen: (AA')^-(1/2) を計算

Eigenで正定値対象行列AA'のinverse sqrtを計算する方法。 (A'はAの転置)
Matlabなら(A*A.t())^-(1/2)で一発。 Eigenでも比較的簡単にできた。

MatrixXf A;
SelfAdjointEigenSolver<MatrixXf> es(A*A.transpose());
MatrixXf ans = es.operatorInverseSqrt();
assert((ans * ans).inverse() == A*A.transpose());

このコードだと、AA'がフルランクじゃない場合(A.rows() > A.cols())におかしなことになる。(たぶん中で固有値ゼロに対して1/sqrt(0)としてしまっている)
ということで、ゼロ固有値は無視するように。

 SelfAdjointEigenSolver<MatrixXf> es(A * A.transpose());
 MatriXfx dvec = es.eigenvalues();
 for(int i=0; i < dvec.rows(); ++i){
     if(abs(dvec(i)) > 1e-5) dvec(i) = 1.0 / sqrt(dvec(i)); // ignore eigen value = zero
 }
 MatrixXf ans = es.eigenvectors() * dvec.asDiagonal() * es.eigenvectors().inverse();