-
-
Save pshriwise/67c2ae78e5db3831da38390a8b2a209f to your computer and use it in GitHub Desktop.
// method for calculating the pseudo-Inverse as recommended by Eigen developers | |
template<typename _Matrix_Type_> | |
_Matrix_Type_ pseudoInverse(const _Matrix_Type_ &a, double epsilon = std::numeric_limits<double>::epsilon()) | |
{ | |
Eigen::JacobiSVD< _Matrix_Type_ > svd(a ,Eigen::ComputeFullU | Eigen::ComputeFullV); | |
// For a non-square matrix | |
// Eigen::JacobiSVD< _Matrix_Type_ > svd(a ,Eigen::ComputeThinU | Eigen::ComputeThinV); | |
double tolerance = epsilon * std::max(a.cols(), a.rows()) *svd.singularValues().array().abs()(0); | |
return svd.matrixV() * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint(); | |
} |
Glad you found it helpful, @gokhansolak! And thanks for the info about non-square matrices. I've added a comment line with this info.
Hi, thanks for the share.
I am using it with an Eigen::Matrix<double, 6, 9>, but in return get the same shape instead of an Eigen::Matrix<double, 9, 6>.
Is it normal behavior?
Thanks
Hi @muttistefano! As @gokhansolak noted above, this implementation doesn't seem to work for a non-square matrix. Perhaps you could give the version linked above a try?
Hi @pshriwise , i tried both versions and they both have problems.
I actually need to use Eigen::ComputeFullU | Eigen::ComputeFullV
cause I need the full matrices of the svd.
If it's useful, I can post here a working version, based on the Eigen::BDCSVD class.
If I will have time I will dig more into this.
For now, I will stick to the working version.
I post it here cause if you look up for eigen pseudoinverse you might end up here.
Thanks
template<typename MatType>
PseudoInverseType<MatType> pseudoInverse(const MatType &a, double epsilon = std::numeric_limits<double>::epsilon())
{
using WorkingMatType = Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, Eigen::Dynamic, 0,
MatType::MaxRowsAtCompileTime, MatType::MaxColsAtCompileTime>;
Eigen::BDCSVD<WorkingMatType> svd(a, Eigen::ComputeFullU | Eigen::ComputeFullV);
svd.setThreshold(epsilon*std::max(a.cols(), a.rows()));
Eigen::Index rank = svd.rank();
Eigen::Matrix<typename MatType::Scalar, Eigen::Dynamic, MatType::RowsAtCompileTime,
0, Eigen::BDCSVD<WorkingMatType>::MaxDiagSizeAtCompileTime, MatType::MaxRowsAtCompileTime>
tmp = svd.matrixU().leftCols(rank).adjoint();
tmp = svd.singularValues().head(rank).asDiagonal().inverse() * tmp;
return svd.matrixV().leftCols(rank) * tmp;
}
Thanks for posting this @muttistefano! If I have time I'll try this version in the application where I applied the inverse in this gist to ensure it also has the behavior I'd expect.
Thank you for sharing.
It gives the following error for a non-square matrix:
Replacing
Eigen::ComputeFullU | Eigen::ComputeFullV
above withEigen::ComputeThinU | Eigen::ComputeThinV
solves this issue as suggested here.Non-square matrix fork