What, precisely, is the derivative function? I would like to illustrate my confusion about this topic by building up the issue from more or less first principles. Let $U \subseteq \mathbb{R}^n$ be an open subset and let $f:U \to \mathbb{R}^m$. We say that $f$ is totally differentiable at $a \in U$ if there exists a linear function $D_{a}f:\mathbb{R}^n \to \mathbb{R}^m$ such that the following holds.
$\lim_{x \to a} \frac{||f(x)-f(a)-D_{a}f(x-a)||}{||x-a||} = 0$
We call $D_{a}f$ the total derivative of $f$ at $a$. My question concerns how we then define the derivative function. For the familiar case in which $n=m=1$, this is simple since the total derivative reduces to just
$(D_{a}f)(b) = \big(\frac{df}{dx}\big{\rvert}_a\big)b$
for all $b \in \mathbb{R}$. Here, $\frac{df}{dx}\big{\rvert}_a$ denotes the usual definition of the derivative of $f$ at $a$. We then define the derivative function $f':U \to \mathbb{R}$ as
$f'(y) := \frac{df}{dx}\big{\rvert}_y$
for all $y \in U$. Let us now consider the case when $m=1$. Representing the action of $D_{a}f$ via the Jacobian matrix, we have
$(D_{a}f)(b) = \begin{bmatrix}
\frac{\partial{f}}{\partial{x_1}}\big{\rvert}_a & \cdots & \frac{\partial{f}}{\partial{x_n}}\big{\rvert}_a
\end{bmatrix}
\begin{bmatrix}
b_1 \\\
\vdots \\\
b_n
\end{bmatrix}$
for all $b=(b_1,\cdots,b_n)\in\mathbb{R}^n$. This suggests that we define the derivative function $Df:U \to \mathbb{R}^n$ as
$(Df)(y) := \big(\frac{\partial{f}}{\partial{x_1}}\big{\rvert}_y, \cdots, \frac{\partial{f}}{\partial{x_n}}\big{\rvert}_y\big)$
for all $y \in U$. We could perhaps make a similar argument for the case in which $n=1$. How does this extend to cases in which the Jacobian matrix is not a simple column or row vector? What definition of the derivative function is most natural in those cases? By natural, I mean that the definition should allow important identities like the product rule and chain rule to retain their obvious forms.
 A: The answer is that the derivative function is the function that inputs the point $a \in U \subset \mathbb R^n$ and that outputs the linear transformation $D_a \in L(\mathbb R^n,\mathbb R^m)$, where $L(V,W)$ is the set of linear transformations from a vector space $V$ to a vector space $W$. In other words, it is a function 
$$D : \mathbb U \to L(\mathbb R^n,\mathbb R^n)
$$
There is, of course, a natural bijection $L(\mathbb R^n,\mathbb R^m) \approx M_{m,n}(\mathbb R)$ where $\mathbb M_{n,m}(\mathbb R)$ is the set of $n \times m$ matrices of real numbers, so that one can also think of the derivative function 
$$D : \mathbb U \to M_{n,m}
$$
where the output is the Jacobian matrix of $f$ at $a$.
A: *

*The derivative $Df$ of a function $f$ is, most generally, the best linear approximation to that function at each point.


*This means that any definition of derivative ought to satisfy
$$f(x+\Delta) \approx f(x) + Df(x)\cdot \Delta $$
That is, when you move a small distance away from a point $x$, to first approximation the change in $f$ is proportional to the value of the derivative (i.e. the change in $f$ is linear, and $Df(x)$ is the constant of proportionality.)


*In detail, if $f:\mathbb{R}^m \rightarrow \mathbb{R}^n$, then $f$ has $m$ inputs with $n$ outputs.  If we know how to take derivatives of functions $\mathbb{R}\rightarrow \mathbb{R}$, we can ask how each of the individual $n$ outputs varies as you vary each of the individual $m$ inputs (while keeping the rest fixed.)
The resulting $m\times n$ matrix of numbers is the derivative of $f$. Note that it's a matrix, i.e. a linear transformation, in keeping with being the best linear approximation to $f$ at each point.


*This definition obeys rules such as the sum rule and multiplication by a constant (because $Df$ is linear). It also, importantly, obeys the chain rule: $D(f\circ g) = (Df\circ g) \bullet Dg$. You can use the chain rule to prove the other rules such as the product and quotient rules.
And this definition reduces to the usual definition of derivatives when $m=1$ and $n=1$, because the components of the matrix are just ordinary limits of difference quotients in each coordinate-axis direction.
