The conceptually simplest mathematical proof I know uses a Coulomb's law integral over the charge distribution. This requires some knowledge of multi-dimensional integrals and matrices. It isn't the most condensed proof, but I think it is thorough and requires the least advanced math.
If you have some charge distribution $\rho(\vec{x})$ then the electric field is
$$\vec{E}(\vec{x}) = \int \, d^3 \vec{x}' \, \rho(x') \frac{\vec{x}-\vec{x'}}{|\vec{x}-\vec{x'}|^3}$$
Now suppose that the charge distribution has some symmetry that is a rotation or a reflection or a translation. That means there is a matrix $R$ and vector $\vec{b}$ that acts on $\vec{x}$ such that $\rho(R\vec{x}+\vec{b}) = \rho(\vec{x})$. Then we can use the integral to calculate the vector field at a transformed point:
$$\vec{E}(R\vec{x}+\vec{b}) = \int \, d^3 \vec{x}' \, \rho(x') \frac{R\vec{x}+\vec{b}-\vec{x'}}{|R\vec{x}+\vec{b}-\vec{x'}|^3}$$
We can change the variable of integration to $\vec{x}' = R\vec{y}+\vec{b}$. So
$$\vec{E}(R\vec{x}) = \int \, d^3 [Ry+b] \, \rho(Ry+\vec{b}) \frac{R(\vec{x}-\vec{y})}{|R(\vec{x}-y)|^3}$$
The charge distribution is symmetrical under $R$ so $\rho(Ry+b) = \rho(y)$. I already used the fact that rotations and reflections are linear in the line above, and that the translation cancels out of the vector difference. A reflection or rotation doesn't change the length of vectors so $|R(\vec{x}-\vec{y})| = |\vec{x}-\vec{y}|$. Since it's a linear transformation $d^3 [Ry+b] = |\det(R)| d^3 y$. Putting in all those relations and simplifying we get that
$$\vec{E}(R\vec{x}+\vec{b}) = |\det{R}|\cdot R \int \, d^3y \, \rho(y) \frac{\vec{x}-\vec{y}}{|\vec{x}-y|^3}$$
$$\vec{E}(R\vec{x}+\vec{b}) = |\det{R}|\cdot (R\vec{E}(\vec{x}))$$
So if the charge distribution is symmetrical under a rotation or reflection so is the electric field in the sense captured by the equation above. Rotations and reflections both have $|\det R| = 1$.
The equation above applies to any charge distribution with any set of rotation/reflection/translation symmetries. The rest of the proof uses the symmetries of the specific case to constrain the form of the electric field.
The infinite uniform cylinder has four symmetries: translation along z, reflection in z, reflection in x (or y), and rotation around z. We can go through each symmetry to further constrain the possible form of the electric field.
Translation: Using $\vec{b} = (0,0,z)$ we get $\vec{E}(x,y,z) = \vec{E}(x,y,0)$, so the entire z dependence of $\vec{E}$ is determined by just the $z=0$ plane.
Reflection in z: Using $R = diag(1,1,-1)$ we get $E_z(x,y,0) = -E_z(x,y,0) \Rightarrow E_z(x,y,0) = 0$, so the electric field has zero z-component everywhere.
Rotation around z: Here it helps to use cylindrical coordinates because they don't interchange under rotation. Using $R$ as the matrix for rotation by $\phi$ around the z axis we get $E_{r,\phi}(r,\phi,z) = E_{r,\phi}(r,0,0)$, so the electric field at any point with cylindrical radius $r$ is by the electric field at the point $(r,0,0)$.
Reflection in y: At the point $(r,0,0)$ a reflection the y axis sends $(r,0,0)$ to $(r,0,0)$. It flips the $\phi$ component of the electric field but not the $r$ component, so $E_\phi (r,0,0) = -E_\phi(r,0,0)\Rightarrow E_\phi (r,0,0) = 0$.
The only remaining component that hasn't been determined by the four symmetries is $E_r(r,0,0)$. We can then use this function and the symmetries to find out the electric field at any other point. Ultimately we have determined that the electric field points radially with no $\phi$ component (reflection in y), no $z$ component (reflection and translation in z), and is rotationally symmetric (rotation around z).