Properties of the argmin function (continuity, differentiability..) Given vectors $x_1,\ldots,x_N \in \mathbb{R}^d$ and a function, say $\psi \colon \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}$, one is interested in the properties of the function $$\Phi\colon (x_1,\ldots, x_N) \mapsto  \underset{y \in \mathbb{R}^d}{\operatorname{argmin}} \sum_{i=1}^N \psi(y,x_i) $$
such as continuity or results regarding differentiability in any sense. 
Does someone have an idea if this was already the focus of some research and can probably refer to literatur or knows some results? 
Thanks in advance :-) 
 A: Let's consider a slightly more general case where we have a continuous $\Psi:\mathbb{R}^d\times\mathbb{R}^d \rightarrow \mathbb{R}$ and $\Phi: \mathbb{R}^d \rightarrow \mathbb{R}$ is given by:
\begin{equation}
\Phi(\underline{x}) := \operatorname{argmin}\limits_{\underline{y}} \Psi(\underline{x},\underline{y})
\end{equation}
To make things well-defined, let's assume that for all $\underline{x}$ there is a unique global minimum $\underline{y}$ for the function $\underline{y} \mapsto \Psi(\underline{x},\underline{y})$, which we will write $\underline{y}_{\underline{x}}$ for short. We will make our life easier by further requiring that $\Psi$ be twice-differentiable in the neighbourhood of any $(\underline{x},\underline{y}_{\underline{x}})$, with $\frac{\partial}{\partial \underline{y}}\frac{\partial}{\partial \underline{y}}\Psi(\underline{x},\underline{y}_{\underline{x}})$ always invertible. We want to show that, in this case, $\Phi$ is differentiable, and work out its Jacobian.
The idea is as follows: if we are standing at the point $(\underline{x},\underline{y}_{\underline{x}})$ and we change $\underline{x}$ to $\underline{x}+d\underline{x}$, then we should be able to "track" the corresponding displacement of $\underline{y}_{\underline{x}}$ to $\underline{y}_{\underline{x}+d\underline{x}} = \underline{y}_{\underline{x}}+d\underline{y}_{\underline{x},d\underline{x}}$, where $d\underline{y}_{\underline{x},d\underline{x}} = J({\underline{x}})\cdot  d\underline{x}$ for some matrix $J({\underline{x}})$. The matrix $J({\underline{x}})$, which we need to determine, is nothing else than the Jacobian of $\Phi$ at $\underline{x}$.
The reason why we expect to be able to follow the movement of $\underline{y}_{\underline{x}}$ is that uniqueness of global minimum for fixed $\underline{x}$ implies uniqueness of local minimum around $\underline{y}_{\underline{x}}$ for fixed $\underline{x}$, and hence we can track the change in $\underline{y}_{\underline{x}}$ by tracking the movement of the unique point of vanishing $\frac{\partial}{\partial \underline{y}}$ in the neighbourhood. Global continuity of $\Psi$ ensures that the local minimum we are tracking remains the global minimum throughout.
In order to find $J({\underline{x}})$, we solve the following vector equation, stating that our displaced $\underline{y}_{\underline{x}+d\underline{x}}$ is indeed again the global minimum (aka the unique local minimum in that neighbourhood):
\begin{equation}
\frac{\partial}{\partial \underline{y}}\Psi(\underline{x}+d\underline{x},\underline{y}_{\underline{x}}+J({\underline{x}})\cdot d\underline{x}) = \underline{0}
\end{equation}
We can move the differentials outside and re-write as follows:
\begin{equation}
\frac{\partial}{\partial \underline{y}}\Psi(\underline{x},\underline{y}_{\underline{x}}) +\frac{\partial}{\partial \underline{y}}\Big( \big(\frac{\partial}{\partial \underline{x}}\Psi(\underline{x},\underline{y}_{\underline{x}})\big)\cdot d\underline{x}+ \big(\frac{\partial}{\partial \underline{y}}\Psi(\underline{x},\underline{y}_{\underline{x}})\big)\cdot  J({\underline{x}})\cdot d\underline{x}\Big) = \underline{0}
\end{equation}
Because $\underline{y}_{\underline{x}}$ is a minimum, we have that $\frac{\partial}{\partial \underline{y}}\Psi(\underline{x},\underline{y}_{\underline{x}})=\underline{0}$, so the equation simplifies to:
\begin{equation}
\frac{\partial}{\partial \underline{y}}\Big( \big(\frac{\partial}{\partial \underline{x}}\Psi(\underline{x},\underline{y}_{\underline{x}})\big)\cdot d\underline{x}+ \big(\frac{\partial}{\partial \underline{y}}\Psi(\underline{x},\underline{y}_{\underline{x}})\big)\cdot J({\underline{x}})\cdot d\underline{x}\Big) = \underline{0}
\end{equation}
Now we write $H^{\underline{y}}_{\underline{x}}(\underline{x},\underline{y}_{\underline{x}})$ for the following matrix:
\begin{equation}
H^{\underline{y}}_{\underline{x}}(\underline{x},\underline{y}_{\underline{x}}) \cdot\underline{v}:= 
\frac{\partial}{\partial \underline{y}}\Big( \big(\frac{\partial}{\partial \underline{x}}\Psi(\underline{x},\underline{y}_{\underline{x}})\big)\cdot\underline{v}\Big) 
\end{equation}
Similarly, we write $H^{\underline{y}}_{\underline{y}}(\underline{x},\underline{y}_{\underline{x}})$ for the following matrix:
\begin{equation}
H^{\underline{y}}_{\underline{y}}(\underline{x},\underline{y}_{\underline{x}}) \cdot\underline{v}:= 
\frac{\partial}{\partial \underline{y}}\Big( \big(\frac{\partial}{\partial \underline{y}}\Psi(\underline{x},\underline{y}_{\underline{x}})\big)\cdot J({\underline{x}})\cdot\underline{v}\Big) 
\end{equation}
Our vector equation is then re-written more concisely as follows:
\begin{equation}
H^{\underline{y}}_{\underline{x}}(\underline{x},\underline{y}_{\underline{x}})+H^{\underline{y}}_{\underline{y}}(\underline{x},\underline{y}_{\underline{x}})\cdot J({\underline{x}}) = 0
\end{equation}
Because we originally assumed that $H^{\underline{y}}_{\underline{y}}$ be invertible around the global minima, we can finally solve for the Jacobian:
\begin{equation}
J({\underline{x}}) = -\big(H^{\underline{y}}_{\underline{y}}(\underline{x},\underline{y}_{\underline{x}})\big)^{-1} \cdot H^{\underline{y}}_{\underline{x}}(\underline{x},\underline{y}_{\underline{x}})
\end{equation}
