10

I'm writing an algorithm with a lot of steps (PCA), and two of them are finding eigenvalues and eigenvectors of a given matrix.

I do not wish to write the whole code for it because I know it is a long job, so I searched for some adhoc code for that but just found 1 or 2 libraries and at first I prefer not to include libraries and I don't want to move to matlab.

Is there any algorithm/tutorial/code that doesn't seem very hard to follow?

Daniel
  • 7,357
  • 7
  • 32
  • 84
  • 4
    Eigen (https://bitbucket.org/eigen/eigen/)? Armadillo (http://arma.sourceforge.net/) ? – Severin Pappadeux May 22 '18 at 01:57
  • 2
    `Eigen` is a headers only package ( so technically, you won't need to link to any library ) and [here](https://eigen.tuxfamily.org/dox/classEigen_1_1EigenSolver.html#title13) is an example from the documentation for calculating eigenvalues. –  May 23 '18 at 00:37

1 Answers1

7

If someone needs this, here how I did it

    Eigen::EigenSolver<Eigen::MatrixXf> eigensolver;
    eigensolver.compute(covmat);
    Eigen::VectorXf eigen_values = eigensolver.eigenvalues().real();
    Eigen::MatrixXf eigen_vectors = eigensolver.eigenvectors().real();
    std::vector<std::tuple<float, Eigen::VectorXf>> eigen_vectors_and_values; 

    for(int i=0; i<eigen_values.size(); i++){
        std::tuple<float, Eigen::VectorXf> vec_and_val(eigen_values[i], eigen_vectors.row(i));
        eigen_vectors_and_values.push_back(vec_and_val);
    }
    std::sort(eigen_vectors_and_values.begin(), eigen_vectors_and_values.end(), 
        [&](const std::tuple<float, Eigen::VectorXf>& a, const std::tuple<float, Eigen::VectorXf>& b) -> bool{ 
            return std::get<0>(a) <= std::get<0>(b); 
    });
    int index = 0;
    for(auto const vect : eigen_vectors_and_values){
        eigen_values(index) = std::get<0>(vect);
        eigen_vectors.row(index) = std::get<1>(vect);
        index++;
    }

Here covmat in which eigenvectors and eigenvalues to be found out. Also, I sort them according to descending order which we do the most of the times. One important thing is that when you selecting which eigendecomposition technique is to be used to be careful because those are not doing the same way. You may find out more information here [https://eigen.tuxfamily.org/dox/group__Eigenvalues__Module.html ]

GPrathap
  • 7,336
  • 7
  • 65
  • 83
  • Whether the solution is real or complex depends entirely on the matrix that you feed. The solver, `Eigen::EigenSolver` admits general matrices, so using ".real()" to get rid of the imaginary part will give the wrong result (also, eigenvectors may have an arbitrary complex phase!). Judging from the name `covmat`, I'm assuming you are feeding a covariance matrix, which is symmetric (or hermitian/self adjoint). In that case you are better off using `Eigen::SelfAdjointEigenSolver` instead, for both performance and accuracy. I don't know what you are trying to do after line 4. – DavidAce Jun 06 '19 at 21:01
  • @DavidAce "admits general matrices, so using ".real()" to get rid of the imaginary part will give the wrong result " since I know the result will be real and needed to get Eigen::VectorXf instead of complex vector – GPrathap Jun 07 '19 at 08:44
  • @DavidAce "I don't know what you are trying to do after line 4" read the below comment: "I sort them according to descending order which we do the most of the times. " – GPrathap Jun 07 '19 at 08:46
  • What is it's alternative for C? – mohammadsdtmnd Jun 25 '23 at 11:50