The idea behind the differential calculus is to study potentially complicated functions $f:\mathbb{R}^n \to \mathbb{R}^m$ by means of linear approximations. Everything flows from this single idea.
For $x\in \mathbb{R}^n$ "the" linear approximation to $f$ near $x$ (if a unique one exists) is called the "derivative" or "gradient" $Df$. (It typically changes from one point to another and so is a function of $x$.)
By definition, then, $$Df: T_x\mathbb{R}^n \to T_{f(x)}\mathbb{R}^m$$ is a linear map from the space of all vectors in $\mathbb{R}^n$ originating at $x$ to the space of all vectors in $\mathbb{R}^m$ originating at $y=f(x)$.
It is a theorem (Spivak, Theorem 2-7) that when we use the directions determined by the coordinates $(x_1, x_2, \ldots, x_n)$ for $\mathbb{R}^n$ and $(y_1, y_2, \ldots, y_m)$ for $\mathbb{R}^m$ as bases for $T_x\mathbb{R}^n$ and $T_{f(x)}\mathbb{R}^m$, respectively, then the entries in the $m\times n$ matrix for $Df$ are the partial derivatives
$$(Df)_{ij} = \frac{\partial y_i}{\partial x_j}.\tag{1}$$
Suppose the "errors" $u_{x_i}$ in the $x_i$ are constant multiples of the standard deviations of random variables $U_i$ describing uncertainties in the $x_i$. One way to express this is in terms of the squared errors: assume there is some positive number $\lambda$ (often $1$, sometimes $2$, occasionally something else) for which
$$u_{x_i}^2 = \lambda^2 \operatorname{Var}(U_i)$$
for each $i$. Applying the linear approximation $Df(x)$, and taking $m=1$ for simplicity (although the general case is scarcely any more difficult), we know from $(1)$ that (at least approximately) the random variable governing uncertainties in $y$ is equal to
$$V = (Df(x)) (U_1, U_2, \ldots, U_n)^\prime = \frac{\partial y}{\partial x_i}U_1 + \frac{\partial y}{\partial x_2}U_2 + \cdots + \frac{\partial y}{\partial x_n}U_n.$$
The formula quoted in the question arises when it is assumed there is no correlation among the $U_i$, whence all the covariances in the calculation of $\operatorname{Var}(V)$ vanish, yielding
$$\eqalign{u_y^2 = \lambda^2\operatorname{Var}(V) &= \lambda^2\left(\operatorname{Var}\left(\frac{\partial y}{\partial x_1}U_1\right) + \cdots + \operatorname{Var}\left(\frac{\partial y}{\partial x_n}U_n\right)\right)\\
&= \left(\frac{\partial y}{\partial x_1}\right)^2 \lambda^2\operatorname{Var}(U_1)+ \cdots + \left(\frac{\partial y}{\partial x_n}\right)^2 \lambda^2\operatorname{Var}(U_n)\\
&= \left(\frac{\partial y}{\partial x_1}\right)^2 u_{x_1}^2 + \cdots + \left(\frac{\partial y}{\partial x_n}\right)^2 u_{x_n}^2
}$$
QED.
The assumptions needed to derive this conclusion provide insight into its meaning, interpretation, and scope:
$f$ must be differentiable at $x$: that is, $Df$ must exist at $x$. This means that within a neighborhood of $x$, $f$ can be approximated to second order in $|x|$ by a linear transformation.
The "errors" $u_{x_i}$ must be multiples of a standard deviation.
The random deviations (that model the uncertainty in $x$) must be uncorrelated.
The typical size of the deviations $U_i = x^{*}_i - x_i$, as measured by the errors $u_{x_i}$, must be small enough that $(Df(x))(x^{*}-x)$ remains a good approximation to $f(x^{*})$; and large deviations (where the approximation no longer holds) must be improbable.
Note that Taylor's Theorem is not needed: the result combines the most basic properties of differentiation with a fundamental property of variances.
Reference
Michael Spivak, Calculus on Manifolds, W. A. Benjamin (1965). Chapters 1 & 2.
