In this paper, we investigate the problem of strong approximation of the solutions of stochastic differential equations (SDEs) when the drift coefficient is given in integral form. We investigate its upper error bounds, in terms of the discretization parameter $n$ and the size $M$ of the random sample drawn at each step of the algorithm, in different subclasses of coefficients of the underlying SDE presenting various rates of convergence. Integral-form drift often appears when analyzing stochastic dynamics of optimization procedures in machine learning (ML) problems. Hence, we additionally discuss connections of the defined randomized Euler approximation scheme with the perturbed version of the stochastic gradient descent (SGD) algorithm. Finally, the results of numerical experiments performed using GPU architecture are also reported, including a comparison with other popular optimizers used in ML.