This question is a slight generalization of this question.

Let M be a constant m $\times$ n matrix and X = Table[Symbol["x" <> ToString[i]], {i, n}] be a vector of n variables.

Assuming that 0 <= X <= 1, I would like to compute the linear map M.X. This can be expressed with Reduce:

Y = Table[Symbol["y" <> ToString[i]], {i, m}]
formula = Exists[Evaluate[X], And@@Thread[Y == M.X] && And@@Thread[0 <= X <= 1]]
Reduce[formula, Y, Reals]

The solution Y forms a polyhedron. How can this polyhedron be efficiently computed with Mathematica? The solution above does not scale.