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)