Eigen: random binary vector with t 1s

Question

I want to compute K*es where K is an Eigen matrix (dimension pxp) and es is a px1 random binary vector with 1s.

For example if p=5 and t=2 a possible es is [1,0,1,0,0]' or [0,0,1,1,0]' and so on...

How do I easily generate es with Eigen?

Answer 1

So you're using Eigen. I'm not sure what matrix type you're using, but I'll go off the class Eigen::MatrixXd.

What you need to do is:

• Create a 1xp matrix that's all 0
• Choose random spots to flip from 0 to 1 that are between 0-p, and make sure that spot is unique.

The following code should do the trick, although you could implement it other ways.

//Your p and t
int p = 5;
int t = 2;

//px1 matrix
MatrixXd es(1, p);

//Initialize the whole 1xp matrix
for (int i = 0; i < p; ++i)
        es(1, i) = 0;

//Get a random position in the 1xp matrix from 0-p
for (int i = 0; i < t; ++i)
{
    int randPos = rand() % p;

    //If the position was already a 1 and not a 0, get a different random position
    while (es(1, randPos) == 1)
        randPos = rand() % p;

    //Change the random position from a 0 to a 1
    es(1, randPos) = 1;
}

Answer 2

I came up with even a better solution, which is a combination of std::vector, Egien::Map and std::shuffle.

std::vector<int> esv(p,0);
std::fill_n(esv.begin(),t,1);
Eigen::Map<Eigen::VectorXi> es (esv.data(), esv.size());
std::random_device rd;
std::mt19937 g(rd());
std::shuffle(std::begin(esv), std::end(esv), g);

This solution is memory efficient (since Eigen::Map doesn't copy esv) and has the big advantage that if we want to permute es several times (like in this case), then we just need to repeat std::shuffle(std::begin(esv), std::end(esv), g);

Maybe I'm wrong, but this solution seems more elegant and efficient than the previous ones.

Answer 3

When t is close to p, Ryan's method need to generate much more than t random numbers. To avoid this performance degrade, you could solve your original problem

find t different numbers from [0, p) that are uniformly distributed

by the following steps

1. generate t uniformly distributed random numbers idx[t] from [0, p-t+1)

2. sort these numbers idx[t]

3. idx[i]+i, i=0,...,t-1 are the result

The code:

VectorXi idx(t);
VectorXd es(p);
es.setConstant(0);

for(int i = 0; i < t; ++i) {
  idx(i) = int(double(rand()) / RAND_MAX * (p-t+1));
}

std::sort(idx.data(), idx.data() + idx.size());

for(int i = 0; i < t; ++i) {
  es(idx(i)+i) = 1.0;
}