Analytical formula for numerical derivative of the matrix pseudo-inverse?

Is there a simple numerical procedure for obtaining the derivative (with respect to $x$) of the pseudo-inverse of a matrix $A(x)$, without approximations (except for the usual floating-point limitations)? The matrix $\frac{\mathrm{d}}{\mathrm{d}x}A(x)$ is supposed to be known.

In other words, are there analytical formulas that could be numerically evaluated so as to obtain the derivative of the pseudo-inverse? or, what formula would generalize $$\frac{\mathrm{d}}{\mathrm{d}x}A^{-1}(x) = -A^{-1}(x) \left(\frac{\mathrm{d}}{\mathrm{d}x}A(x)\right) A^{-1}(x)$$ for the pseudo-inverse?

I would be happy if this were possible, as this would allow my uncertainty calculation programming package to precisely calculate uncertainties on the pseudo-inverse of matrices whose elements have uncertainties (currently, a numerical differentiation is performed, which may yield imprecise results in some cases).

Any idea would be much appreciated!

The answer is known since at least 1973: a formula for the derivative of the pseudo-inverse of a matrix $A(x)$ of constant rank can be found in

The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate. Author(s): G. H. Golub and V. Pereyra. Source: SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432

References 29 and 30 in the above paper contain an earlier formula that can also be used to obtain the same result (papers by P.A. Wedin).

The case of non-constant rank is simple: the pseudo-inverse is not continuous, in this case (see Corollary 3.5 in On the Perturbation of Pseudo-Inverses, Projections and Linear Least Squares Problems. G. W. Stewart. SIAM Review, Vol. 19, No. 4. (Oct., 1977), pp. 634-662).

Here is the formula for a matrix of constant rank (equation (4.12), in the Golub paper):

$$\frac{\mathrm d}{\mathrm d x} A^+(x) = -A^+ \left( \frac{\mathrm d}{\mathrm d x} A \right) A^+ +A^+ A{^+}^T \left( \frac{\mathrm d}{\mathrm d x} A^T \right) (1-A A^+) + (1-A^+ A) \left( \frac{\mathrm d}{\mathrm d x} A^T \right) A{^+}^T A^+$$

(for a real matrix).

For complex matrices, the above formula works if Hermitian conjugates are used instead of transposes. I don't have any reference on this (anyone?), but this is verified by all the numerical tests I did (with matrices of various shapes and ranks).

This is not a complete answer.

According to the Wikipedia page you linked, the pseudoinverse $A^+$ is not a continuous function of $A$, as it jumps around when $A$ is ill-conditioned. Therefore, you can't expect $A^+(x)$ to always have a derivative in terms of the matrix derivative of $A(x)$.

I suppose it may be reasonable to ask for a formula that works when the trajectory of $A(x)$ is restricted to constant rank strata, but I don't know such a formula. If it exists, it should blow up as you approach a place where the rank jumps down.

This is a shot in the dark, but you could try working out a limit of Tikhonov regularizations, taking the derivative of $(A^\ast(x)A(x) - \epsilon I)^{-1}A^\ast(x)$ as $\epsilon \to 0$.

• +1: Thank you for the discussion about the matrix rank. I was hoping to avoid the Tikhonov regularization approach because of the numerical complications it implies, but it is the only solution I could think of. – Eric O Lebigot May 24 '10 at 21:02

Since your goal is to take into account the effect of perturbations in the matrix elements on the least squares solution you may find the following useful:

http://www.jstor.org/pss/2156807

http://en.wikipedia.org/wiki/Total_least_squares

Depending on what assumptions you make about the nature of the perturbations, these results might save you some work on the analysis.