1

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?

justHelloWorld
  • 6,478
  • 8
  • 58
  • 138

3 Answers3

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.

justHelloWorld
  • 6,478
  • 8
  • 58
  • 138
  • Looks good. Shuffle is O(n), better than sort. But you solution need to generate p random numbers rather than t. Could you replace the fixed numbers like 10 and 3 with variable names like p and t so it is more readable? You don't need std::vector as VectorXi can do the same thing. – kangshiyin Jul 13 '16 at 04:09
  • Sorry I replaced the numbers with `p` and `t`. Anyway I'm pretty sure that I can't use `std::shuflle` with `Eigen::VectorXi` (I need standard iterators) or use `std::fill_n`. Am I wrong? – justHelloWorld Jul 13 '16 at 04:15
  • Pointer is a iterator. See how I call std::sort – kangshiyin Jul 13 '16 at 04:24
0

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;
}
Ryan
  • 314
  • 1
  • 10
0

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;
}
kangshiyin
  • 9,681
  • 1
  • 17
  • 29