Derivative of a linear vector function by matrix - 3d-Jacobian or a vector? If I differentiate
$$
Y = 
\pmatrix{a & b \\ c & d}
\pmatrix {x \\ y}
$$
with respect to the matrix on the left side, what do I get?
On one hand, I think I should get $\pmatrix{x \\ y}$ because it's a linear function.
On the other hand, I think I should get a 3-dimensional Jacobian matrix (if it's a real thing) where the each 4 entries are 2-d vectors
$$
{\partial Y \over \partial a} = \pmatrix{x \\ 0},
{\partial Y \over \partial b} = \pmatrix{y \\ 0} \\
{\partial Y \over \partial c} = \pmatrix{0 \\ x},
{\partial Y \over \partial d} = \pmatrix{0 \\ y} \\
$$
Which one is right?
 A: You can look at your function as $Y : \Bbb{R}^4 \to \Bbb{R}^2$ given by
$$Y(a,b,c,d) = \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} x \\ y\end{pmatrix} = \begin{pmatrix} ax+by \\ cx+dy\end{pmatrix}$$
so clearly the second interpretation is correct. Te resulting $2\times 4$ matrix is
$$\begin{pmatrix} \frac{\partial Y}{\partial a} & \frac{\partial Y}{\partial b} & \frac{\partial Y}{\partial c} & \frac{\partial Y}{\partial d} \end{pmatrix} = \begin{pmatrix} x & y & 0 & 0 \\ 0 & 0 & x & y\end{pmatrix}$$
as you calculated. Now, recall that the differential $DY(a,b,c, d)$ at a point $(a,b,c,d) \in \Bbb{R}^4$ is actually a linear map $\Bbb{R}^4\to\Bbb{R}^2$ whose matrix is precisely the above Jacobian.
Therefore it acts as
$$DY(a,b,c,d)\begin{pmatrix} a' \\ b' \\ c' \\ d'\end{pmatrix} = \begin{pmatrix} x & y & 0 & 0 \\ 0 & 0 & x & y\end{pmatrix}\begin{pmatrix} a' \\ b' \\ c' \\ d'\end{pmatrix} = \begin{pmatrix} a'x+b'y \\ c'x+d'y\end{pmatrix} = Y(a',b',c',d')$$
so the linear map $DY(a,b,c,d)$ is precisely $Y$ itself. This is because $Y$ is a linear map, as you noticed in your first interpretation.
A: You have a map $Y:\mathbb R^{2\times2}\to\mathbb R^2$ given by
$$\begin{pmatrix}a&b\\c&d\end{pmatrix}\mapsto\begin{pmatrix}a&b\\c&d\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}.$$
This map is linear, and the derivative of a linear map at any arbitrary point is always the linear map itself. So you have $\mathrm D Y(A)=Y$ for all $A\in\mathbb R^{2\times2}$.
You can also find a matrix representation of said map if you choose a basis. Let's say we choose the standard basis for $\mathbb R^2$, and the basis consisting of
$$\begin{pmatrix}1&0\\0&0\end{pmatrix},\begin{pmatrix}0&1\\0&0\end{pmatrix},\begin{pmatrix}0&0\\1&0\end{pmatrix},\begin{pmatrix}0&0\\0&1\end{pmatrix}$$
for $\mathbb R^{2\times2}$. I will call these matrices $B_1,\dots,B_4$.
Now we can simply calculate
$$Y(B_1)=\begin{pmatrix}x\\0\end{pmatrix},~Y(B_2)=\begin{pmatrix}y\\0\end{pmatrix},~Y(B_3)=\begin{pmatrix}0\\x\end{pmatrix},~Y(B_4)=\begin{pmatrix}0\\y\end{pmatrix}.$$
So the matrix representation with respect to this basis is
$$\begin{pmatrix}x&y&0&0\\0&0&x&y\end{pmatrix}.$$
The dimensions of the matrix shouldn't be surprising, since $Y$ maps from a 4d space $\mathbb R^{2\times2}$ to a 2d space $\mathbb R^2$, so it needs 4 columns and 2 rows.
