# 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.

• Why are you defining $L$ with reference to MATLAB? Strictly speaking MATLAB functions are defined for floating point (finite precision) matrices, so your question makes no sense to me. Aug 18 '12 at 16:50
• >X = linsolve(A,B) solves the linear system A*X = B using LU factorization with partial pivoting when A is square and QR factorization with column pivoting otherwise. This should be a comment.
– MLE
Aug 18 '12 at 17:21
• Are you just interested in differentiating something that involves the solution of a linear system? Or are you interested in the derivative of a floating-point algorithm (which would be very strange and subtle)? Aug 19 '12 at 7:23

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.

• Strange answer: linsolve uses the QR factorization only for numerical stability reasons. For the full rank real overdetermined least square problem wouldn't be more easy to define $A^\dagger = (A^TA)^{-1} A^T$ and differentiate this expression? Sorry if I repeat myself: linsolve does not define $A^\dagger$, but tries to compute a numerical approximation Adagger = A\eye(size(A,1)). I think that while MATLAB is a wonderful tool for learning numerical analysis, it is dangerous as a linear algebra testbed. Aug 19 '12 at 0:06
• 1. QR factorization is an extremely useful way to think about non-square matrices, not just a handy numerically stable method for computing things. 2. I don't understand what point you are trying to make. I don't think your suggestion would be simpler or provide more understanding. Aug 19 '12 at 2:40
• Jed, I checked your result and I suspect that dropping the projection is wrong. In fact the null space of $dA$ is different from the null space of $A$. Please see the result in my own answer. In fact $QQ^T = AA^\dagger$, but $I-AA^\dagger \neq 0$ while $(I-AA^\dagger)AA^\dagger = 0$. I agree with your remarks: QR gives an interesting interpretation of the least square formulation. The problem is that the $QR$ factorization is not unique, unless you somehow constrain $R$, and for the problem at hand the least squares formula is easier to handle (although tedious). Aug 19 '12 at 10:21
• @StefanoM: I agree that the result should be checked. $A^{\dagger}$ is a $\{1,2\}$-generalized inverse of $A$. ( See Ben-Israel and Greville, Generalized Inverses: Theory and Applications.) There are multiple such generalized inverses; your definition computes one of them. Using Jed's notation for $A = QR$, then your definition of $A^{\dagger}$ yields $A^{\dagger} = R^{-1}R^{-T}Q^{T}R^{T}$, which differs from Jed's definition by changing the basis used to express $Q$. Jed's definition is consistent with the linsolve algorithm. Aug 19 '12 at 11:39
• @GeoffOxberry I feel like beeing overly pedantic, but $(A^TA)^{-1}A^T = (R^TQ^TQR)^{-1} R^TQ^T = R^{-1} Q^T$, so the two $A^\dagger$ are the same. Aug 19 '12 at 15:35

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.

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.