Here is a Mixed Integer Linear Problem (MILP) formulation that may be solved in practice for some instances with MILP-solvers like CPLEX.
For every integer vector $(i_1,\dots,i_k)\in [1,a_1]\times\cdots\times[1,a_k]$, let us introduce two variables: a nonnegative real (rational) $x_{i_1,\dots,i_k}$ and a binary $y_{i_1,\dots,i_k}$.
The $x$'s represent the elements of $M$ (with some elements being zero -- the more such elements, the better). They satisfy the following equalities: $$\forall j\in[1,k]\quad\forall t\in[1,a_j]\ :\qquad\sum_{i_1,\dots,i_k\atop i_j=t} x_{i_1,\dots,i_k} = \frac{1}{a_j}.$$
In an optimal solution, we want as many as possible $x$'s be zero and this is why we need $y$'s. In an optimal solution, they will represent indicator values for the positivity of $x$'s, i.e., $y_{i_1,\dots,i_k}=1$ iff $x_{i_1,\dots,i_k}>0$. This can be achieved with the inequalities: $$x_{i_1,\dots,i_k} \leq y_{i_1,\dots,i_k}$$ and the objective function: $$\min \sum_{i_1,\dots,i_k} y_{i_1,\dots,i_k}.$$
I have implemented this in Sage (with CPLEX) and here is a couple of examples of computed optimal $M$'s:
(3,4,5): [1/60, 1/30, 1/20, 1/12, 7/60, 2/15, 1/6, 1/5, 1/5]
(4,5,6): [1/60, 1/60, 1/30, 1/30, 1/15, 1/12, 7/60, 2/15, 1/6, 1/6, 1/6]