diff --git a/src/faust_linear_operator/CPU/faust_MatDense.h b/src/faust_linear_operator/CPU/faust_MatDense.h index ee224c37d842cc0a1dd7b66ec278dedbfee51683..ccef5e5d3d6025bc2e03d11457d6535f65d7b1a7 100644 --- a/src/faust_linear_operator/CPU/faust_MatDense.h +++ b/src/faust_linear_operator/CPU/faust_MatDense.h @@ -522,6 +522,8 @@ namespace Faust void best_low_rank(const int &r, MatDense<FPP,Cpu> &bestX, MatDense<FPP, Cpu> &bestY) const; + void best_low_rank2(MatDense<FPP,Cpu> &bestX, MatDense<FPP, Cpu> &bestY) const; + void initJacobiSVD(Eigen::JacobiSVD<Eigen::Matrix<FPP, Eigen::Dynamic, Eigen::Dynamic>>& svd); /** * \brief Returns the best rank-1 approximate this = bestX * bestY using the svd/eigensolver/power iteration. diff --git a/src/faust_linear_operator/CPU/faust_MatDense.hpp b/src/faust_linear_operator/CPU/faust_MatDense.hpp index eb76baaa4289f81949d00096a9328d4b716d52fe..328b4e4bc2ad839dfd2afe84eec0fa72c6ca8cdb 100644 --- a/src/faust_linear_operator/CPU/faust_MatDense.hpp +++ b/src/faust_linear_operator/CPU/faust_MatDense.hpp @@ -58,6 +58,8 @@ #include <Eigen/SVD> #include "faust_init_from_matio.h" +#include <Eigen/Eigenvalues> + namespace Faust { @@ -1404,7 +1406,7 @@ template<typename FPP> void MatDense<FPP, Cpu>::best_low_rank(const int &r, MatDense<FPP,Cpu> &bestX, MatDense<FPP, Cpu> &bestY) const { #if(EIGEN_WORLD_VERSION > 3 || EIGEN_WORLD_VERSION >= 3 && EIGEN_MAJOR_VERSION >= 3) - Eigen::JacobiSVD<Eigen::Matrix<FPP, Eigen::Dynamic, Eigen::Dynamic>> svd(this->mat, Eigen::ComputeThinU | Eigen::ComputeThinV); + Eigen::BDCSVD<Eigen::Matrix<FPP, Eigen::Dynamic, Eigen::Dynamic>> svd(this->mat, Eigen::ComputeThinU | Eigen::ComputeThinV); #else //#ifdef _MSC_VER // as far as I tested eigen3.4rc1 doesn't compile with VS 14 @@ -1414,6 +1416,27 @@ void MatDense<FPP, Cpu>::best_low_rank(const int &r, MatDense<FPP,Cpu> &bestX, M best_low_rank(r, bestX, bestY, svd); } + +template<typename FPP> +void MatDense<FPP, Cpu>::best_low_rank2(MatDense<FPP,Cpu> &bestX, MatDense<FPP, Cpu> &bestY) const +{ + // NOTE: this function is another try to do faster than best_low_rank, but it doesn't work better for rank == 1 + MatDense<FPP, Cpu> AtA; + bestX.resize(this->getNbRow(), 1); + bestY.resize(this->getNbCol(), 1); + gemm(*this, *this, AtA, FPP(1.0), FPP(0.0), 'N', 'H'); + Eigen::SelfAdjointEigenSolver<Eigen::Matrix<FPP, Eigen::Dynamic, Eigen::Dynamic>> es; + es.compute(AtA.mat); +// std::cout << "eigenvalues:" << es.eigenvalues().transpose() << std::endl; + int loc; + es.eigenvalues().maxCoeff(&loc); +// std::cout << "max of eigenvalues:" << loc << std::endl; + bestX.mat = es.eigenvectors().col(loc); + bestX.mat *= sqrt(es.eigenvalues()[loc]); + bestY.mat = bestX.mat.householderQr().solve(mat); + bestY.adjoint(); +} + template<typename FPP> void MatDense<FPP, Cpu>::approx_rank1(MatDense<FPP,Cpu> &bestX, MatDense<FPP, Cpu> &bestY) const {