Outer addition with Eigen

Question

I am attempting to translate the following python function into C++:

import numpy as np
from scipy.linalg import blas
def scaled_dist(a, b, ls):
    al = a/ls
    bl = b/ls
    tmp1 = np.sum(al**2, axis=1)
    tmp2 = np.sum(bl**2, axis=1)
    tmp3 = np.add.outer(tmp1, tmp2, order='F')
    tau = blas.dgemm(a=al, b=bl, alpha=-2.0, c=tmp3, beta=1, trans_b=1)
    np.clip(tau, 0, np.inf, out=tau)
    return tau

However I have hit a stumbling block with the line:

tmp3 = np.add.outer(tmp1, tmp2)

My c++ code compiles but encounters a runtime error when executed. The code (up to that line) is:

Eigen::MatrixXd test2(const Eigen::MatrixXd &x1, const Eigen::MatrixXd &x2,const Eigen::VectorXd &vec)
{
  Eigen::MatrixXd r = Eigen::MatrixXd::Zero(x1.rows(), x2.rows());
  Eigen::MatrixXd al = x1.array().rowwise() / vec.transpose().array();
  Eigen::VectorXd tmp1 = al.array().square().rowwise().sum();

  Eigen::MatrixXd bl = x2.array().rowwise() / vec.transpose().array(); 
  Eigen::VectorXd tmp2 = bl.array().square().rowwise().sum();

  r = tmp1.transpose().array() + tmp2.array();
  return r;
}

I am able to make sense of the runtime error, which is (I believe) an assertion error, complaining that the left and right hand sides of the addition expression do not match in size. My approach was motivated by the fact that tmp1.transpose() * tmp2 does appear to produce the expected result.

My question is as follows:

Given two vectors, vec1 and vec2, what is the idiomatic way using Eigen of achieving the same functionality as numpy.add.outer(vec1, vec2), namely an "outer" addition whereby a matrix is obtained by adding the (broadcast) rows of one vector to the (broadcast) columns of the other? i.e., if

vec1 = [1,2,3]
vec2 = [3,4,5]

then

outer_add(vec1, vec2) = 
[4, 5, 6]
[5, 6, 7]
[6, 7, 8]

ggael ggael · Accepted Answer · 2019-05-26T20:23:18

You can use replicate for that, e.g.:

Vector3f v1(1,2,3), v2(3,4,5);
MatrixXf r = v1.rowwise().replicate(v2.size())
           + v2.transpose().colwise().replicate(v1.size());

Outer addition with Eigen

1 Answers