I feel your pain, and I've been working on a multidimensional array library
that's designed to let Eigen do all the matrix math:
http://code.google.com/p/ndarray/
While that isn't the complete solution solution you're looking for, at least
it's a clear division of labor.
Essentially, while you can make an ndarray::Array of any dimension, you can
view 1- or 2-dimensional subarrays of those as Eigen-conforming objects, or
construct ndarray views into Eigen::Matrix objects.
This is a good approach, I'd like to encourage you in this direction.
I almost wrote back to the list to suggest that this might be the best
approach to take balancing reasonable performance with reasonable ease
of implementation.
Like with matrices, the hard part is products. The name of the game is
memory locality / cache friendliness. By leveraging Eigen's optimized
matrix product implementation and applying it to 2-dimensional slices
in your N-dim arrays, you should be able to get somewhere. Depending
on the storage order/ cases, you may have to copy slices into
temporaries to have them contiguously stored in memory.