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.

  • 1
    $\begingroup$ 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. $\endgroup$
    – Stefano M
    Aug 18 '12 at 16:50
  • $\begingroup$ >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. $\endgroup$
    – MLE
    Aug 18 '12 at 17:21
  • $\begingroup$ 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)? $\endgroup$ 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$.


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.

  • $\begingroup$ 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. $\endgroup$
    – Stefano M
    Aug 19 '12 at 0:06
  • 1
    $\begingroup$ 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. $\endgroup$
    – Jed Brown
    Aug 19 '12 at 2:40
  • $\begingroup$ 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). $\endgroup$
    – Stefano M
    Aug 19 '12 at 10:21
  • $\begingroup$ @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. $\endgroup$ Aug 19 '12 at 11:39
  • $\begingroup$ @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. $\endgroup$
    – Stefano M
    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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.