Consider the heteroscedastic nonparametric regression model with random design $Y_i = f(X_i) + V^{1/2}(X_i)\varepsilon_i, \quad i=1,2,\ldots,n$, with $f(\cdot)$ and $V(\cdot)$ $\alpha$- and $\beta$-H\" older smooth, respectively. We show that the minimax rate of estimating $V(\cdot)$ under both local and global squared risks is of the order $n^{-\frac{8\alpha\beta}{4\alpha\beta + 2\alpha + \beta}} \vee n^{-\frac{2\beta}{2\beta+1}}$, where $a\vee b\define \max\{a,b\}$ for any two real numbers $a,b$. This result extends the fixed design rate $n^{-4\alpha} \vee n^{-2\beta/(2\beta+1)}$ derived in \cite{wang2008effect} in a non-trivial manner, as indicated by the appearance of both $\alpha$ and $\beta$ in the first term. In the special case of constant variance, we show that the minimax rate is $n^{-8\alpha/(4\alpha+1)}\vee n^{-1}$ for variance estimation, which further implies the same rate for quadratic functional estimation and thus unifies the minimax rate under the nonparametric regression model with those under the density model and the white noise model. Extensions to the estimation of higher order moments, multivariate models, and additive models will also be discussed.