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.
There are two different views about the semiclassical limit in quantum mechanics, the first is based on a somewhat shaky ground due to the fact that the existence of the Feynman integral is not proved yet. On the other side, Wiener integral, its imaginary time counterpart does exist and one could pretend to work things out from this and then move to the Feynman integral. The other approach relies on substantial mathematical theorems due to Elliott Lieb and Barry Simon in the '70 and is essentially valid for many-body physics. These latter results make the limit $\hbar\rightarrow 0$ and $N\rightarrow\infty$ equivalent while the former is not really a physical limit due to the fact that Planck constant is never zero.
Starting from Feynman path integral, the standard formulation applies to a mechanial problem described from a Lagrangian $L$, normally $L=\frac{\dot x^2}{2}-V(x)$ but one can extend this to more general cases, and then the postulate is that, given a path $x(t)$, this must contribute to the full quantum mechanical amplitude of a particle going from the point $x_a$ to $x_b$ with a term $e^{\frac{i}{\hbar}S}$ being $S=\int_{t_a}^{t_b}dtL(\dot x,x,t)$ the action. All the possible paths contribute and so, the full amplitude will be given by the formal writing
$$
A(x_a,x_b)\sim\int[dx(t)]e^{\frac{i}{\hbar}\int_{t_a}^{t_b}dtL(\dot x,x,t)}.
$$
Be warned that this integral is not proved to exist yet, but the Wiener counterpart, that can be obtained changing $t\rightarrow it$, exists and describes Brownian motion. Now, if you take the formal limit $\hbar\rightarrow 0$ to this integral you will immediately recognize the conditions to apply the stationary phase method to it. This implies that the functional must have an extremum and this can be obtained by pretending that
$$
\delta S=\delta \int_{t_a}^{t_b}dtL(\dot x,x,t)=0
$$
that is, the paths that give the greatest contribution are the classical ones and one recover the classical limit as a variational principle as learned from standard textbooks.
While this is a quite common approach, to extend what really happens to a macroscopic system that we can see to respect all the laws of classical mechanics, we have to turn our attention to the limit of a large number of particles $N\rightarrow\infty$. In this case one has more rigorous results. These are due to Lieb and Simon as already said above. They published two papers about
Lieb E. H. and Simon B. 1973 Phys. Rev. Lett. 31, 681.
Lieb E. H. and Simon B. 1977 Adv. in Math. 23, 22.
In the first paper, their theorem 4 states
Theorem: For $\lambda < Z$, let $E_N^0$ and $\rho_N^0(x)$ denote the ground-state energy and one-electron distribution function for N spin-$\frac{1}{2}$ electrons obeying the Pauli principle and interacting with $k$ nuclei as described above. Then (a) $N^{-\frac{7}{3}}E_N^0\rightarrow E_1$, as $N\rightarrow\infty$; (b) $N^{-2}\rho_N^0(N^{-\frac{1}{3}}x)\rightarrow\rho_1(x)$ as $N\rightarrow\infty$, where convergence in (b) means that for any domain $D\subset R^3$, the expected fraction of electrons in $N^{-\frac{1}{3}}D$ approaches $\int_D\rho_1 (x)d^3x$.
Where $\rho_1(x)$ and $E_1$ refer to the Thomas-Fermi distribution and the corresponding energy. This theorem states that the limit $N\rightarrow\infty$ for a quantum system, under some mild conditions, is the Thomas-Fermi distribution. A system with this distribution is a classical system. The fact that a system with a Thomas-Fermi distribution is a classical one can be seen through the following two references:
W. Thirring(Ed.), The Stability of Matter: From Atoms to Stars - Selecta of E. Lieb, Springer-Verlag (1997).
L. Hörmander, Comm. Pure. Appl. Math. 32, 359 (1979).
The second paper just gives the mathematical support to derive Thomas-Fermi approximation as the leading order of a classical expansion for $\hbar\rightarrow 0$ that I will not present here.
Best Answer
Well, there is always the trivially enforced solution $$\tag{1} S[x,\lambda]~=~\int\! dt \sum_{i=1}^3\lambda_i(t) \left(\ddot{x}^i(t)+\alpha \dot{x}^i(t) \right),$$ where $\lambda_i(t)$ are three Lagrange multiplier variables. From now on we assume that we are not allowed to use other variables than $x^i(t)$.
Whether the 3D ODE $\ddot{\bf x}+\alpha\dot{\bf x}=0$ has a Lagrangian formulation or not is in general a hard problem, see e.g. Douglas' theorem and the Helmholtz conditions on this Wikipedia page. (We are unaware if these conditions have been completely analyzed in 3D. On the other hand, it is easy to see that the corresponding 1D problem has a Hamiltonian action formulation, at least locally, cf. e.g. this Phys.SE answer.)
The 3D problem simplifies enormously if we assume that the form of the kinetic part $T$ of the Lagrangian $L=T-U$ is just $T=\frac{1}{2}\dot{\bf x}^2$. (This assumption is of course heavily motivated from physics, and is undoubtedly necessary if we would like to make quantum mechanical sense of the system.) Then OP's question essentially boils down to if there exists a velocity-dependent potential $U({\bf x}, \dot{\bf x})$ for the friction force ${\bf F}=-\alpha \dot{\bf x}$. This is equivalent to asking if the Helmholtz conditions [1] $$\tag{2} \frac{\partial F_i}{\partial x^j} -\frac{1}{2}\frac{d}{dt}\frac{\partial F_i}{\partial \dot{x}^j} ~=~[i \longleftrightarrow j], \qquad \frac{\partial F_i}{\partial \dot{x}^j}~=~-[i \longleftrightarrow j].$$ are satisfied. It is a straightforward to see that this is not the case. (Already the second condition (2b) is not met.) See also this related Phys.SE post. For a general discussion of the notion of conservative force for velocity-dependent forces, see e.g. this Phys.SE answer.
References: