The Lagrangian (or variational) formulation of the Euler-Lagrange equations and the Hamiltonian formulation are equivalent. This equivalence can be made quite explicit and goes a bit deeper than the standard treatments show. The equivalence can be established in several steps, which I'll try to outline below with references.
The Hamiltonian formalism is a special case of the Lagrangian formalism.
There is a generic operation that can be performed on variational problems: adjunction and elimination of auxiliary fields or variables. Given a Lagrangian $L(x,y)$, the variable $y$ (which could be vector valued) is called auxiliary if the Euler-Lagrange equations obtained from the variation of $y$ can be algebraically solved for solved for in terms of the remaining variables and their derivatives $y=y(x)$. The important point here is that we need not solve any differential equations to obtain $y(x)$. The Lagrangian $L'(x) = L(x,y(x))$ gives a new variational principle with the auxiliary field $y$ eliminated. The critical points of $L'(x)$ are in one-to-one correspondence with those of the original $L(x,y)$. The adjunction of an auxiliary field is the reverse operation. Given $L(x)$, we look for another Lagrangian $L'(x,y)$ where $y$ is auxiliary and its elimination gives $L'(x,y(x))=L(x)$.
It is straightforward to check that given a Lagrangian $L(x)$, with Hamiltonian $H(x,p)$, the new Lagrangian $L'(x,p) = p\dot{x} - H(x,p)$, associated to Hamilton's Least Action Principle, is a special case of adjoining some auxiliary fields, namely the momenta $p$. The elimination $p=p(x)$ is precisely the inverse Legendre transform.
The moral here is that the Legendre transform is not sacred. I learned this point of view from the following paper of Barnich, Henneaux and Schomblond (PRD, 1991). This point of view is particularly helpful in higher order field theories (multiple independent variables, second or higher order derivatives in the Lagrangian) where a unique notion of Legendre transform is lacking.
The phase space is the space of solutions of the Euler-Lagrange equations.
When the Euler-Lagrange equations have a well-posed initial value problem. Then solutions can be put into one-to-one correspondence with initial data. The initial data are uniquely specified by the canonical position and momentum variables that are commonly used to define the canonical phase space in the Hamiltonian picture. This is a rather common identification nowadays and can be found in many places, so I won't give a specific reference.
If the Euler-Lagrange equations do not have a well-posed initial value problem, the definitions change somewhat on either side, but an equivalence can still be established. See the references in the next step for more details.
Any Lagrangian defines a (pre)symplectic form (current).
This is, unfortunately, less well known than it should be. The (pre)symplectic structure of classical mechanics and field theory can be defined straightforwardly directly from the Lagrangian. In case the equations of motion are degenerate, the form is degenerate and hence only presymplectic, otherwise symplectic. There is more than one way to do this, but a particularly transparent one is referred to as the covariant phase space method. A very nice (though not original) reference is Lee & Wald (JMP, 1990). See also this nLab page, which also has a more extensive reference list.
Applying this method to the Lagrangian of Hamilton's Least Action Principle gives the standard symplectic form in terms of the canonical position and momentum coordinates.
Briefly, and without going into the details of differential forms on jet spaces where this is most easily formalized, the construction is as follows. Denote by $d$ the space-time (which is just 1-dimensional if you only have time) exterior differential and by $\delta$ the field variation exterior differential. Without dropping boundary or total divergence terms, the total variation of the Lagrangian can be expressed as $\delta L(x) = \mathrm{EL}\delta x + d\theta$. Here the Lagrangian is a space-time volume form, $\mathrm{EL}$ denotes the Euler-Lagrange equations and $d\theta$ consists of all the terms that are usually dropped during partial integration. Clearly $\theta$ is of 1-form in terms of field variations and as a spacetime form of one degree lower than a volume form (aka a current). Taking another exterior field variation, we get $\omega=\delta\theta$, which is the desired (pre)symplectic current. If $\Sigma$ is a Cauchy surface (in particular it is codimension-1), the (pre)symplectic form is defined as $\Omega=\int_\Sigma \omega$, now a 2-form in terms of field variations, as expected. It can be shown that the (pre)symplectic form is closed, $d\omega=0$, when evaluated on solutions of the Euler-Lagrange equations. Hence, by Stokes' theorem, $\Omega$ is independent of the choice of $\Sigma$. One space-time is 1-dimensional, $\Omega$ is just $\omega$ evaluated at a particular time.
A (pre)symplectic form (current) defines a Lagrangian.
As alluded to in the question, Hamilton's equations of motion are often expressed in a special form that highlights a certain geometrical structure that is not obvious in the original Lagrangian form. It is well known from the symplectic formulation of classical mechanics that this structure can be seen as a consequence of the fact that they correspond time evolution generated by a Hamiltonian via the Poisson bracket. The symplectic form is then preserved by the evolution. There are analogous statements for field theory. In fact, no variational formulation is necessary to discuss this geometrical structure of the equations.
What is quite remarkable is that a kind of converse to this statement holds as well. Namely, given a system of (partial) differential equations, if there is a conserved (pre)symplectic current $\omega$ on the space of solutions ($\omega$ is field dependent and $d\omega=0$ when evaluated on solutions, see previous step), then a subsystem of the equations is derived from a variational principle. There is a subtlety here. Even if there exists a conserved $\omega$, if it is degenerate (not symplectic) other independent equations need to be added to the Euler-Lagrange equations of the corresponding variational principle to obtain a system of equations equivalent to the original one. Onece the Lagrangian of this variational principle is known, the Hamiltonian and symplectic forms could be defined in the usual way and the original system of equations recast in the canonical Hamilton form.
To my knowledge, the above observation first appeared in Henneaux (AnnPhys, 1982) for ODEs and in Bridges, Hydon & Lawson (MathProcCPS, 2010) for PDEs. The calculation demonstrating this observation is given in a bit more detail on this nLab page.
Another way to look at this result is to consider a conserved (pre)symplectic form as a certificate for the solution of the inverse problem of the calculus of variations.
A final note about the usefulness of the Hamiltonian formulation. Despite the fact that as a consequence of the above discussion it's not strictly necessary. Any symplectic manifold has local coordinates in which the symplectic form is canonical (Darboux's theorem). The Legendre transform identifies this choice of coordinates explicitly.
