How to raise each element of each row of matrix M to the power of each element of vector x in C++ using Eigen?

Question

I am using the Eigen library to perform the following operation. Given a matrix M and a vector x, I would like to raise each row of M to the power of x element-wise. That is, if M = [[10, 20], [20, 30], [30, 40], [40, 50]] and x = [2, 3], I would like to get a result of [[100, 8000], [400, 27000], [900, 64000], [1600, 125000]]. Is there a way of achieving this using Eigen's SIMD features and without any explicit for loop?

I tried the following code

VectorXd func() { 
    // [[10, 20], 
    //  [20, 30], 
    //  [30, 40], 
    //  [40, 50]]
    Eigen::MatrixXd M(4, 2);
    M(0, 0) = 10; M(0, 1) = 20;
    M(1, 0) = 20; M(1, 1) = 30;
    M(2, 0) = 30; M(2, 1) = 40;
    M(3, 0) = 40; M(3, 1) = 50; 

    std::vector<double> x_numbers  = {2, 3}; 
    Eigen::Map<Eigen::VectorXd> x(x_numbers.data(), x_numbers.size());

    auto retval = M.array().pow(x.array()); 
    // I expect a 4x2 matrix returned but I am getting a vector of 
    // size two instead. Changing the return type to MatrixXd 
    // does not help either. 
    return retval;
}

However, the return value of this function for the above values of M and x is [100, 8000]. I also tried M.array().rowwise().pow(x.array()) but I get "error: no member named 'pow' ...". What am I missing? I am using Eigen-3.4.0.

Answer 1

For working with Eigen's array operations, either both operands need to be the same shape, or one of them is a scalar. Because you want an operation between matrix and a vector, you need to broadcast the vector to match the matrix first. This can be done with the replicate function.

Eigen::MatrixXd func1(const Eigen::MatrixXd& M, const Eigen::VectorXd& x) {
    return M.array().pow(x.transpose().array().replicate(M.rows(),1));
}

See Demo

If you want to avoid the extra copies of x then you can only work on a column by column basis with M and you will need to write a for loop, but Eigen's array operation can still be used:

Eigen::MatrixXd func2(const Eigen::MatrixXd& M, const Eigen::VectorXd& x) {
    Eigen::MatrixXd retval(M.rows(), M.cols());
    for (Eigen::Index col = 0; col < M.cols(); ++col) {
        retval.col(col) = M.col(col).array().pow(x(col));
    }
    return retval;
}

(In the demo is an extra version using NullaryExpr which is not using a for loop but does not exhibit Eigen's array operations)

Answer 2

When you call M.array().pow(x.array()), it’s raising each element of M to the corresponding element of x for each corresponding element of M and x. That’s why you’re only getting the first two elements of the result you expect. It is not doing the multiplication for each row of M with the vector x.

Here’s how you can modify your function:

MatrixXd func(const Ref<const Matrix<double,Dynamic,Dynamic,RowMajor>>& M, VectorXd x) { 
    auto retval = M.array().pow(x.transpose().array().replicate(M.rows(),1));
    return retval;
}

x.transpose().array().replicate(M.rows(),1) creates a matrix where each row is the vector x, and then M.array().pow(...) raises each row of M to the corresponding row of that matrix.

Edit: Sorry, there was an error in the code. It has been fixed and the answer is updated now.

Output:

Result:
   100   8000
   400  27000
   900  64000
  1600 125000