  We establish a general theory of optimal strong error estimation for
numerical approximations of a second-order parabolic stochastic partial
differential equation with monotone drift driven by a multiplicative
infinite-dimensional Wiener process. The equation is spatially discretized by
Galerkin methods and temporally discretized by drift-implicit Euler and
Milstein schemes. By the monotone and Lyapunov assumptions, we use both the
variational and semigroup approaches to derive a spatial Sobolev regularity
under the $L_\omega^p L_t^\infty \dot H^{1+\gamma}$-norm and a temporal
H\"older regularity under the $L_\omega^p L_x^2$-norm for the solution of the
proposed equation with an $\dot H^{1+\gamma}$-valued initial datum for
$\gamma\in [0,1]$. Then we make full use of the monotonicity of the equation
and tools from stochastic calculus to derive the sharp strong convergence rates
$O(h^{1+\gamma}+\tau^{1/2})$ and $O(h^{1+\gamma}+\tau^{(1+\gamma)/2})$ for the
Galerkin-based Euler and Milstein schemes, respectively.
