Assume $m=1$ and $f$ is of bounded variation on $[0,1]$. The problem is to estimate $$ \|f\|^2_{2}-\frac{1}{n}\sum\limits_{i=1}^n\left(f({p_i})\right)^2=\int_{0}^{1}f^{2}(t)dt -\frac{1}{n}\sum\limits_{i=1}^n f^{2}({p_i}), $$ as the number of points grows. Setting $g=f^{2}$, which is also of bounded variation, the question is about the rate of convergence of the Riemann sums of $g$ to its integral.

For a regular mesh $\{1/n,2/n,\ldots,1\}$, one has $$ \left|\int_{0}^{1}g(t)dt-\frac1n\sum_{k=1}^{n}g(k/n)\right|\leq\int_{0}^{1/n}\sum_{k=1}^{n} |g(t+(k-1)/n)-g(k/n)|dt\leq\frac{V(g)}{n}, $$ where $V(g)$ denotes the variation of $g$.

More generally, consider a tagged mesh $T=\{\sigma_{k},[s_{k-1},s_{k}],~k=1,\ldots,n\}$ of $[0,1]$ such that $T\ll\delta$ meaning that $\max_{k}(s_{k}-s_{k-1})<\delta$, and set $$ g(T)=\sum_{k=1}^{n}g(\sigma_{k})(s_{k}-s_{k-1}),\qquad \psi_{\delta}(g)=\sup _{T \ll \delta}\left|g(T)-\int_{0}^{1} g(t) dt\right| $$ Then, the following holds, for any function $g$, $$ \sup _{\delta>0} \frac{\psi_{\delta}(g)}{\delta} \leq V(g) \leq 2\liminf _{\delta \to 0}\frac{\psi_{\delta}(g)}{\delta}, $$ see J.A. Alewine, Rates of uniform convergence for Riemann integrals. Missouri J. Math. Sci. 26 (2014), 48-56.

Hence, for a function of bounded variation, its Riemann sums converge to its integral at a rate of $O(\delta)$, and that rate cannot be improved.

In several variables ($m>1$), $f$ being of bounded variation means that its distributional derivative is a finite Radon measure. It is not clear what follows from that assumption.

Assume $m=1$ and $f$ is of bounded variation on $[0,1]$. The problem is to estimate $$ \|f\|^2_{2}-\frac{1}{n}\sum\limits_{i=1}^n\left(f({p_i})\right)^2=\int_{0}^{1}f^{2}(t)dt -\frac{1}{n}\sum\limits_{i=1}^n f^{2}({p_i}), $$ as the number of points grows. Setting $g=f^{2}$, which is also of bounded variation, the question is about the rate of convergence of the Riemann sums of $g$ to its integral.

For a regular mesh $\{1/n,2/n,\ldots,1\}$, one has $$ \left|\int_{0}^{1}g(t)dt-\frac1n\sum_{k=1}^{n}g(k/n)\right|\leq\int_{0}^{1/n}\sum_{k=1}^{n} |g(t+(k-1)/n)-g(k/n)|dt\leq\frac{V(g)}{n}, $$ where $V(g)$ denotes the variation of $g$.

More generally, consider a tagged mesh $T=\{\sigma_{k},[s_{k-1},s_{k}],~k=1,\ldots,n\}$ of $[0,1]$ such that $T\ll\delta$ meaning that $\max_{k}(s_{k}-s_{k-1})<\delta$, and set $$ g(T)=\sum_{k=1}^{n}g(\sigma_{k})(s_{k}-s_{k-1}),\qquad \psi_{\delta}(g)=\sup _{T \ll \delta}\left|g(T)-\int_{0}^{1} g(t) dt\right| $$ Then, the following holds, for any function $g$, $$ \sup _{\delta>0} \frac{\psi_{\delta}(g)}{\delta} \leq V(g) \leq 2\liminf _{\delta \to 0}\frac{\psi_{\delta}(g)}{\delta}, $$ see J.A. Alewine, Rates of uniform convergence for Riemann integrals. Missouri J. Math. Sci. 26 (2014), 48-56.

Hence, for a function of bounded variation, its Riemann sums converge to its integral at a rate of $O(\delta)$, and that rate cannot be improved.

For several variables, you may have a look here.