This question is a slight generalization of this question.

Let M be a fixed rational m $\times$ n matrix, e.g., M = Rationalize[RandomReal[{-1, 1}, {m, n}], 0], 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.

Another way of stating the problem is to compute a matrix $A$ and a vector $b$ such that $\forall x: ~ x \in [0,1]^n \Leftrightarrow A M x \leq b$.