C++ find eigenvalues and eigenvectors of matrix

ceigenvalueeigenvectormatrix

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?

Best Answer

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 ]

Setting a bit

Use the bitwise OR operator (|) to set a bit.

number |= 1UL << n;

That will set the nth bit of number. n should be zero, if you want to set the 1st bit and so on upto n-1, if you want to set the nth bit.

Use 1ULL if number is wider than unsigned long; promotion of 1UL << n doesn't happen until after evaluating 1UL << n where it's undefined behaviour to shift by more than the width of a long. The same applies to all the rest of the examples.

Clearing a bit

Use the bitwise AND operator (&) to clear a bit.

number &= ~(1UL << n);

That will clear the nth bit of number. You must invert the bit string with the bitwise NOT operator (~), then AND it.

Toggling a bit

The XOR operator (^) can be used to toggle a bit.

number ^= 1UL << n;

That will toggle the nth bit of number.

Checking a bit

You didn't ask for this, but I might as well add it.

To check a bit, shift the number n to the right, then bitwise AND it:

bit = (number >> n) & 1U;

That will put the value of the nth bit of number into the variable bit.

Changing the nth bit to x

Setting the nth bit to either 1 or 0 can be achieved with the following on a 2's complement C++ implementation:

number ^= (-x ^ number) & (1UL << n);

Bit n will be set if x is 1, and cleared if x is 0. If x has some other value, you get garbage. x = !!x will booleanize it to 0 or 1.

To make this independent of 2's complement negation behaviour (where -1 has all bits set, unlike on a 1's complement or sign/magnitude C++ implementation), use unsigned negation.

number ^= (-(unsigned long)x ^ number) & (1UL << n);

unsigned long newbit = !!x;    // Also booleanize to force 0 or 1
number ^= (-newbit ^ number) & (1UL << n);

It's generally a good idea to use unsigned types for portable bit manipulation.

number = (number & ~(1UL << n)) | (x << n);

(number & ~(1UL << n)) will clear the nth bit and (x << n) will set the nth bit to x.

It's also generally a good idea to not to copy/paste code in general and so many people use preprocessor macros (like the community wiki answer further down) or some sort of encapsulation.

C++ – What are the differences between a pointer variable and a reference variable in C++

A pointer can be re-assigned:

int x = 5;
int y = 6;
int *p;
p = &x;
p = &y;
*p = 10;
assert(x == 5);
assert(y == 10);

A reference cannot be re-bound, and must be bound at initialization:

int x = 5;
int y = 6;
int &q; // error
int &r = x;

A pointer variable has its own identity: a distinct, visible memory address that can be taken with the unary & operator and a certain amount of space that can be measured with the sizeof operator. Using those operators on a reference returns a value corresponding to whatever the reference is bound to; the reference’s own address and size are invisible. Since the reference assumes the identity of the original variable in this way, it is convenient to think of a reference as another name for the same variable.
```
int x = 0;
int &r = x;
int *p = &x;
int *p2 = &r;

assert(p == p2); // &x == &r
assert(&p != &p2);
```

You can have arbitrarily nested pointers to pointers offering extra levels of indirection. References only offer one level of indirection.

int x = 0;
int y = 0;
int *p = &x;
int *q = &y;
int **pp = &p;

**pp = 2;
pp = &q; // *pp is now q
**pp = 4;

assert(y == 4);
assert(x == 2);

A pointer can be assigned nullptr, whereas a reference must be bound to an existing object. If you try hard enough, you can bind a reference to nullptr, but this is undefined and will not behave consistently.

/* the code below is undefined; your compiler may optimise it
 * differently, emit warnings, or outright refuse to compile it */

int &r = *static_cast<int *>(nullptr);

// prints "null" under GCC 10
std::cout
    << (&r != nullptr
        ? "not null" : "null")
    << std::endl;

bool f(int &r) { return &r != nullptr; }

// prints "not null" under GCC 10
std::cout
    << (f(*static_cast<int *>(nullptr))
        ? "not null" : "null")
    << std::endl;

You can, however, have a reference to a pointer whose value is nullptr.

Pointers can iterate over an array; you can use ++ to go to the next item that a pointer is pointing to, and + 4 to go to the 5th element. This is no matter what size the object is that the pointer points to.
A pointer needs to be dereferenced with * to access the memory location it points to, whereas a reference can be used directly. A pointer to a class/struct uses -> to access its members whereas a reference uses a ..
References cannot be put into an array, whereas pointers can be (Mentioned by user @litb)
Const references can be bound to temporaries. Pointers cannot (not without some indirection):
```
const int &x = int(12); // legal C++
int *y = &int(12); // illegal to take the address of a temporary.
```
This makes const & more convenient to use in argument lists and so forth.