Q:

derivative of linsolve

Consider a vector $\mathbf{g} \in \mathbb{R}^{m}$ and a matrix $\mathbf{A} \equiv \mathbf{A(g)} \in \mathcal{M}_{p\times q} [\mathbb{R}]$, a function of $\mathbf{g}$. 
Furthermore, let $\mathbf{S} \in \mathcal{M}_{p \times r} [\mathbb{R}]$, independent of $\mathbf{g}$.
Is it a way to calculate $\displaystyle \frac{\partial L(\mathbf{A}, \mathbf{S})}{\partial g_i}$, where $L(\cdot)$ is the linsolve Matlab function? 
One can suppose that $\displaystyle \frac{\partial \mathbf{A}}{\partial g_i}$ is known.
Note: $X =$ linsolve$(A,B)$ solves the linear system $AX = B$ using LU factorization with partial pivoting when $A$ is square and QR factorization with column pivoting otherwise.

A:

Let's simplify slightly. Suppose the overdetermined $m\times n$ system $A$ admits the factorization $Q R = A$ with $R$ full rank. Let
$$A^\dagger = R^{-1} Q^T$$
be the linsolve solution operator and observe that
$$I = R^{-1} Q^T Q R  = A^\dagger A.$$
Differentiating the above, we have
$$ 0 = \frac{\partial (A^\dagger A)}{\partial g} = \frac{\partial A^\dagger}{\partial g} A + A^\dagger \frac{\partial A}{\partial g} . $$
We cannot solve for $\frac{\partial A^\dagger}{\partial g}$ algebraically because it is not uniquely determined in this form, but we can compute
$$ \frac{\partial A^\dagger}{\partial g} Q Q^T = - A^\dagger \frac{\partial A}{\partial g} R^{-1} Q^T = - R^{-1} Q^T \frac{\partial A}{\partial g} R^{-1} Q^T = - A^\dagger \frac{\partial A}{\partial g} A^\dagger $$
where $QQ^T$ is the identity on the range of $A$.
Edit
If we need the rest, $\frac{\partial A^\dagger}{\partial g} (I - Q Q^T)$, we can formally differentiate the QR factorization,
$$ \frac{\partial A^\dagger}{\partial g} = \frac{\partial R^{-1}}{\partial g} Q^T + R^{-1} \frac{\partial Q^T}{\partial g},$$
and decompose into the parts in the range of $A$ (which yields expressions equivalent to that derived above) and in the left null space of $A$,
$$\begin{align} \frac{\partial A^\dagger}{\partial g} (I - Q Q^T) &= \left[ \frac{\partial R^{-1}}{\partial g} Q^T + R^{-1} \frac{\partial Q^T}{\partial g} \right] (I - Q Q^T) \\
&= R^{-1} \frac{\partial Q^T}{\partial g} (I - Q Q^T) .
\end{align} $$
Note that only perturbations $dg$ that change $Q$ can contribute to this second term. We can rewrite in terms of $\frac{\partial A}{\partial g}$ by
$$\begin{align}
  \frac{\partial A^\dagger}{\partial g} (I - Q Q^T) &= R^{-1} R^{-T} \left[ R^T \frac{\partial Q^T}{\partial g} + \frac{\partial R^T}{\partial g} Q^T \right] (I - Q Q^T) \\
&= R^{-1} R^{-T} \frac{\partial A^T}{\partial g} (I - Q Q^T) \\
&= (A^T A)^{-1} \frac{\partial A^T}{\partial g} (I - A A^\dagger)
\end{align}$$
where the last line is the second term in Stefano's derivation, though the QR form
$$\frac{\partial A^\dagger}{\partial g} = R^{-1} \left[ -Q^T \frac{\partial A}{\partial g} R^{-1} Q^T + R^{-T} \frac{\partial A^T}{\partial g} (I - Q Q^T) \right]$$
is preferable to compute with for numerical stability and because it reuses more computation and operates in a smaller space.
You can apply similar arguments to rank-deficient and under-determined cases.

A:

X = linsolve(A,B) is a convenience function of MATLAB that tries to find a sensible numerical solution to the linear system $$AX = B$$ in the general case, $A\in\mathbb{C}^{m\times n}$ with $m \gtreqqless n$. 
For (numerical) full rank A, linsolve operates in such a way that the result is linear in the second argument, so we can assume that the solution to the above problem is defined as $$X = A^\dagger B,$$ where $A^\dagger$ is a function of $A$ alone. The MATLAB approximation to $A^\dagger$ can be computed as 
Ad = linsolve(A, eye(size(A,1)));

If the question at hand is to compute $\frac{\partial X}{\partial g_i}$, the the proble reduces to computing $$\frac{\partial A^\dagger}{\partial g_i}.$$
I will discuss only the case of real $A$, with $m>n$ and $\mathrm{rank}(A) = n$, i.e. an overdetermined full rank system. In this case MATLAB computes the least square solution, and we have $$A^\dagger = (A^TA)^{-1}A^T.$$ Now 
\begin{multline}
\frac{\partial A^\dagger}{\partial g_i} = 
\frac{\partial (A^TA)^{-1}}{\partial g_i} A^T +
(A^TA)^{-1} \frac{\partial A^T}{\partial g_i} 
= \\ =
- (A^TA)^{-1} \left[  \frac{\partial A^T}{\partial g_i}A +  A^T \frac{\partial A}{\partial g_i} \right] (A^TA)^{-1} A^T
+ (A^TA)^{-1} \frac{\partial A^T}{\partial g_i}
= \\ =
-A^\dagger \frac{\partial A}{\partial g_i} A^\dagger + (A^TA)^{-1} \frac{\partial A^T}{\partial g_i} (I-AA^\dagger).
\end{multline} 
Analogous calculations can be performed in the remaing cases, with reference to the actual problem MATLAB is solving.

