An atomic-scale theory of the viscoelastic response of metallic glasses is derived from first principles, using a Zwanzig-Caldeira-Leggett system-bath Hamiltonian as a starting point within the framework of nonaffine linear response to mechanical deformation. This approach provides a Generalized-Langevin-Equation (GLE) as the average equation of motion for an atom or ion in the material, from which non-Markovian nonaffine viscoelastic moduli are extracted. These can be evaluated using the vibrational density of states (DOS) as input, where the boson peak plays a prominent role for the mechanics. To compare with experimental data of binary ZrCu alloys, numerical DOS was obtained from simulations of this system, which take also electronic degrees of freedom into account via the embedded atom method (EAM) for the interatomic potential. It is shown that the viscoelastic $\alpha$-relaxation, including the $\alpha$-wing asymmetry in the loss modulus, can be very well described by the theory if the memory kernel (the non-Markovian friction) in the GLE is taken to be a stretched-exponential decaying function of time. This finding directly implies strong memory effects in the atomic-scale dynamics, and suggests that the $\alpha$-relaxation time is related to the characteristic time-scale over which atoms retain memory of their previous collision history. This memory time grows dramatically below the glass transition.