derivative of linsolve

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

The question asks for the derivative of the solution of $X$ of $A^*(AX-S)=0$ with respect to a parameter $t=g_i$ on which $A$ (but not $S$) depends.

Differentiation of the defining equation with respect to $t$ (denoted by a dot) gives $\dot A^*(AX-S) +A^*(\dot A X+A\dot X)=0$, hence $\dot X$ is the solution of the system of equations $A^*A\dot X=B:=\dot A^*(S-AX)-A^*\dot A X$.

The right-hand side $B$ can be computed after $X$ has been found, typically from an orthogonal factorization $A=QR$ (but a Cholesky factor $R$ of the normal equations matrix $A^*R$ would also do). Then $\dot X$ can be found without refactorization by solving $R^*R\dot X=B$.

In case a $QR$ factorization is used, a slightly more numerically stable version first solves $R^*Y=\dot A^*(S-AX)$ and then $R\dot X=Y-Q^*\dot A X$.

Derivatives with respect to all $g_i$ can be obtained in both cases more efficiently by stacking the right hand sides, using BLAS 3 routines.