In numpy i have the following

```
X1 = np.arange(1, 6).reshape(-1, 1)
X1 = np.hstack([X1, X1])
# Output: X1
# array([ [1, 1],
# [2, 2],
# [3, 3],
# [4, 4],
# [5, 5] ])
X2 = X1[:, np.newaxis, :]
X3 = X1[np.newaxis, :, :]
# Output: Shapes
# X2.shape = (5, 1, 2)
# X3.shape = (1, 5, 2)
X4 = (X2 - X3) ** 2
# Output: X4
# X4.shape = (5, 5, 2)
# X4[:, :, 0] =
# array([ [ 0, 1, 4, 9, 16],
# [ 1, 0, 1, 4, 9],
# [ 4, 1, 0, 1, 4],
# [ 9, 4, 1, 0, 1],
# [16, 9, 4, 1, 0] ], dtype=int32)
# X4[:, :, 1] =
# array([ [ 0, 1, 4, 9, 16],
# [ 1, 0, 1, 4, 9],
# [ 4, 1, 0, 1, 4],
# [ 9, 4, 1, 0, 1],
# [16, 9, 4, 1, 0] ], dtype=int32)
```

Is this type of operations do-able in Eigen C++, if so how?

I would need this because I will later have to divide `X4`

by another (Row) vector with the same number of columns as `X1`

. I understand Eigen has its broadcasting operations (https://eigen.tuxfamily.org/dox/group__TutorialReductionsVisitorsBroadcasting.html) however my numpy implementation works with 3D matrices.

Source: Windows Questions C++