1

I have this matlab code that I want to write in Eigen:

[V_K,D_K] = eig(K);
d_k = diag(D_K);
ind_k = find(d_k > 1e-8);
d_k(ind_k) = d_k(ind_k).^(-1/2);
K_half = V_K*diag(d_k)*V_K';

So far the code that I've written is (the last part is based following this example):

Eigen::EigenSolver<Eigen::MatrixXf> es (K,true);
Eigen::MatrixXcf v = es.eigenvalues();
v = (v.array()).select(1/std::sqrt(v),v);

This is obviously wrong (first of all, you cannot apply std::sqrt(v) on Eigen::MatrixXcf), but I didn't figure out how to write it. Can you help me?

Community
  • 1
  • 1
justHelloWorld
  • 6,478
  • 8
  • 58
  • 138

1 Answers1

3

As in the example you linked to, you need to compare the .array() (.real()?) to > 1e-8 and use Eigens functions (or unaryExpr) to manipulate the values of the matrix:

v = (v.array().real() 
//            ^^^
//        can't compare complex numbers

    > 1e-8)
//   ^^^
//  Actually compare the number 
    .select(v.cwiseSqrt().cwiseInverse(), v);
//            ^^^
//         Use Eigens functions or alternatively unaryExpr

// In one line:
v = (v.array().real() > 1e-8).select(v.cwiseSqrt().cwiseInverse(), v);
Avi Ginsburg
  • 10,323
  • 3
  • 29
  • 56
  • Thanks so much! Can you please tell me the difference with `unaryExpr`? – justHelloWorld Jul 12 '16 at 09:58
  • 1
    You could write your own function that would be applied to each value. Usually, it's better to use Eigens functions when available, as this allows vectorization and expression templates. – Avi Ginsburg Jul 12 '16 at 10:01
  • Thanks for the explanation! Could you please give a look at [this](http://stackoverflow.com/questions/38327550/different-eigenvector-and-eigenvalues-in-eigen-and-matlab-could-generate-errors) question? – justHelloWorld Jul 12 '16 at 11:29