How to see the ${\rm SO}(4)$ symmetry of the classical Kepler problem? It is well-known that the hidden symmetry of the $1/r$-problem $$H=\frac{{\bf p}^2}{2m}-\frac{k}{r}\tag{1}$$ is ${\rm SO(4)}$ in the sense that the components of angular momentum $L_i$ and those of the Runge-Lenz vector $A_i$ satisfy the $\mathfrak{so}(4)$ algebra, and also $[H,L_i]=[H,A_i]=0$ in quantum mechanics (QM). Therefore, in QM the invariance of the Hamiltonian under the ${\rm SO}(4)$ is trivial to check in the sense $$H'=UHU^{-1}=H\tag{2}$$ where $U$ is the unitary operator implementing the ${\rm SO}(4)$ symmetry and involves the exponential of the linear combinations of $L_i$'s and $A_i$'s.  That part is fine!
Question How about the classical scenario? When we talk about the ${\rm SO}(3)$ symmetry of the Hamiltonian of a central force problem, the symmetry is manifest in both the potential and kinetic term.
${\rm SO}(4)$ transformations can be thought of transformations which leave quantities of the form $$\Delta\equiv \zeta_1^2+\zeta_2^2+\zeta_3^2+\zeta^2_4,~~ \zeta_i\in\mathbb{R}\tag{3}$$ invariant. Is there a way to reduce the Hamiltonian of $(1)$ in the form of $\Delta$ so that the ${\rm SO}(4)$ symmetry becomes manifest?
 A: The classical treatment is actually simpler than the quantum one, since there are no operator ordering problems. This is what Pauli relied on in his solution of the Hydrogen atom, before Schroedinger, as every decent physicist of his generation had worked out the Kepler problem like this.
The idea is to find as many independent constants of the motion Q as possible, and break them up into independent,  ideally PB-commuting, quantities.
Every such invariant must PB-commute with the Hamiltonian, since that is its time derivative,
$$ \{ Q,H   \}= 0.$$
As a result, the analog of your quantum (2) is
$$
H'= e^{ \theta \{Q ,} H= H+\theta \{ Q,H\}+ {\theta^2 \over 2!} \{Q,\{Q,H \}\}+... =H.
$$
The linear first order phase-space gradients $\{Q, \bullet \}={ \partial Q\over \partial x^i} {\partial \over  \partial p^i}- { \partial Q\over\partial p^i } {\partial \over \partial x^i} $ are the generators of Lie transformations in phase space. They close into  a Lie algebra    purely classically, via PB commutation,
$$
\{Q, \{S, - \{S, \{ Q ,= \{ \{Q,S\}, 
$$
by Jacobi's identity. They are thus Lie algebra elements, and can act adjointly.
For the specific Kepler problem, let me follow an older talk of mine.
$$
H = \frac{{\bf p}^2}{2}-\frac{1}{r} ~,  
$$
in simplified (rescaled) notation.
The invariants of the hamiltonian are the angular momentum vector,
$$
{\bf L} = {\bf r} \times {\bf p}~,
$$
and the Hermann-Bernoulli-Laplace vector,
now usually called the Pauli-Runge-Lenz vector,
$$
{\bf A}=  {\bf p} \times {\bf L} -\hat{{\bf r}}~.  
$$
(Dotting it by $\hat{{\bf r}}$ instantly yields Kepler's elliptical orbits by inspection,
$\hat{{\bf r}}\cdot {\bf A}+1 = {\bf L}^2 /r$.)
Since $~{\bf A}\cdot {\bf L}=0$, it follows that
$$
H=\frac{{\bf A}^2-1}{2{\bf L}^2 } ~~.
$$
However, to simplify the PB Lie-algebraic structure,
$$
\{ L_i, L_j \} = \epsilon^{ijk} L_k, \qquad
\{ L_i, A_j \} = \epsilon^{ijk} A_k, \qquad
\{ A_i, A_j \} =-2H  \epsilon^{ijk} L_k, 
$$
it is useful to redefine
${\bf D} \equiv \frac{{\bf A}}{\sqrt{- 2 H}}$, and further
$$  
{\cal R}\equiv {\bf L} + {\bf D}, \qquad {\cal L}\equiv {\bf L} - {\bf D}. 
$$
These six simplified invariants obey the standard
$SU(2)\times SU(2) \sim SO(4)$ symmetry algebra,
$$
\{ {\cal R}_i, {\cal R}_j \} = \epsilon^{ijk} {\cal R}_k, \qquad
\{ {\cal R}_i, {\cal L}_j \} = 0, \qquad
\{ {\cal L}_i, {\cal L}_j \} =\epsilon^{ijk} {\cal L}_k,  
$$
and depend on each other and the hamiltonian through
$$
\bbox[yellow]{H=\frac{-1}{2 {\cal R}^2}=\frac{-1}{2 {\cal L}^2} }~,  
$$
so only five of the invariants are actually algebraically independent.
It is then manifest that all six respective infinitesimal transformations leave their Casimirs invariant, the Casimirs of the other SU(2) invariant, and hence the Hamiltonian invariant; consequently the corresponding finite transformations do likewise,
$$
H'= e^{ \theta^i \{ {\cal R}_i,} ~ e^{\phi^j\{ {\cal L}_j ,}  ~ H= H ~.
$$
You could work out these explicit infinitesimal phase-space transformations involved, if you were really inclined, and Fock and Bargmann do that, as linked in the comments, but this is the whole point of systematic formalism, that it cuts through the confusion and forestalls missing the forest for the trees. It's Lie's major contribution to geometry.
