Derivado de $\sqrt{AA^T}$ con respecto a $A$

Question

Cómo hallar la derivada de la función $f: M_n(\mathbb{R})\to M_n(\mathbb{R}), A\mapsto \sqrt{AA^T},$ donde $A^T$ es el transpuesto de la matriz $A$ ?

\begin{align} Df_V(A) & = \lim_{h\to 0}\dfrac{f(A+hV)-f(A)}{h}\\ & = \lim_{h\to 0}\dfrac{\sqrt{(A+hV)(A+hV)^T}-\sqrt{AA^T}}{h}\\ & = \lim_{h\to 0} \dfrac{\sqrt{AA^T+hAV^T+hVA^T+h^2VV^T}-\sqrt{AA^T}}{h} \end{align}

Ahora, ¿qué debo hacer?

Answer 1

Definir la matriz $F$ tal que $$F = \sqrt{AA^T} \;\implies F^2 = AA^T $$ En operación vec para convertir estas matrices en vectores. $$a={\rm vec}(A),\quad f={\rm vec}(F)$$ El gradiente solicitado puede calcularse del siguiente modo. $$\eqalign{ &F\,F &= AA^T \cr &F\,dF\,(I)+(I)\,dF\,F &= A\,dA^T\,(I)+(I)\,dA\,A^T \cr &(I^T\otimes F+F^T\otimes I)\,{\rm vec}(dF) &= (I^T\otimes A)\,{\rm vec}(dA^T)+(A\otimes I)\,{\rm vec}(dA) \cr &\Big(I\otimes F+F\otimes I\Big)\,df &= \Big((I\otimes A)K+(A\otimes I)\Big)\,da \cr &\frac{\partial f}{\partial a} &= \Big(I\otimes F+F\otimes I\Big)^+ \Big((I\otimes A)K+(A\otimes I)\Big) \cr\cr }$$ donde $M^+$ denota el pseudoinverso de $M$ , $I$ es la matriz identidad, y $K$ es el matriz de conmutación asociado al producto de Kronecker. La solución también aprovecha el hecho de que $I$ y $F$ son simétricos.

El verdadero gradiente es un tensor de cuarto orden, mientras que el resultado anterior es un aplanado versión de la misma. Si quieres el resultado completo del tensor, los componentes no son demasiado difíciles de calcular, ya que hay una correspondencia uno a uno entre los elementos del tensor y los elementos de la matriz aplanada.

Al releer la pregunta, parece que está interesado en la derivada direccional, donde los cambios en $A$ se restringen a la dirección de la matriz $V$ . $$\eqalign{ v &= {\rm vec}(V) \\ df_v &= \bigg(\frac{\partial f}{\partial a}\bigg)\cdot v \\ }$$ Una vez más, desinflamar este vector en una matriz es sólo un mapeo uno a uno, es decir. $$\eqalign{ F &\in {\mathbb R}^{m\times n} \implies f \in {\mathbb R}^{mn\times 1} \\ F_{ij} &= f_{\alpha} \\ \alpha &= i+(j-1)\,m \\ i &= 1+(\alpha-1)\,{\rm mod}\,m \\ j &= 1+(\alpha-1)\,{\rm div}\,m \\ }$$

Answer 2

El mayor problema de esta derivada es que las matrices no son necesariamente conmutables. Este puesto debería darle una idea de las dificultades que uno puede encontrar.

No obstante, vamos a intentarlo. Pero primero, generalicemos un poco el problema considerando el mapa

$$ f\colon\mathbb S^n_+\longrightarrow\mathbb S^n_+,\, X\longmapsto X^{\alpha} \overset{\text{def}}{=} \exp(\alpha\log(X))$$

que asigna una matriz simétrica definida positiva a su potencia matricial, definida mediante el exponencial matricial y el logaritmo matricial. Se puede comprobar fácilmente que $\sqrt{X} = X^{1/2}$ diagonalizando.

1. El caso conmutativo

Si $\Delta X$ conmuta con $X$ entonces todo es fácil, porque entonces $\log(X\cdot\Delta X)=\log(X) + \log(\Delta X)$ y $\exp(X+\Delta X) = \exp(X)\cdot\exp(\Delta X)$ que no se cumple en el caso general

En el caso conmutativo tenemos

$$\begin{aligned} \log(X+\Delta X)-\log(X) &= \log(X+\Delta X) + \log(X^{-1}) \\&= \log((X+\Delta X)X^{-1}) \\&= \log(I+X^{-1}\Delta X) \\&= \sum_{k=1}^{\infty}(-1)^{k+1} \frac{(X^{-1}\Delta X)^{k}}{k} \sim X^{-1}\Delta X \\ \exp(X+\Delta X) - \exp(X) &= \exp(X)\exp(\Delta X) - \exp(X) \\&= \exp(X)\big(\exp(\Delta X) -I\big)\\ &= \exp(X)\sum_{k=1}^{\infty} \frac{1}{k!}\Delta X^k \sim \exp(X)\Delta X \end{aligned}$$ Por lo tanto $\partial_{\Delta X}\log(X) = X^{-1}$ y $\partial_{\Delta X} \exp(X) = \exp(X)$ en el caso de conmutación $[X, \Delta X]=0$ . En particular, mediante la regla de la cadena se deduce que tenemos $\partial_{\Delta X} X^\alpha = \alpha X^{\alpha-1}$ que satisface nuestras expectativas.

2. El caso no conmutativo

Este es mucho más difícil, se puede tratar de utilizar el Fórmula Baker-Campbell-Hausdorff y/o el Fórmula Zassenhaus pero el resultado final no será bonito. En general, la bonita fórmula de antes ya no vale. Por ejemplo, si $\alpha=2$ entonces una vez puede comprobar fácilmente que

$$(X+\Delta X)^2-X^2 \sim X\cdot\Delta X + \Delta X \cdot X \neq 2X\cdot\Delta X$$

También podemos comprobarlo numéricamente, por ejemplo con el programa python (3.7) que se muestra a continuación. Los residuos son estables al aumentar el tamaño de la muestra $N$ en el caso conmutativo, pero serán cada vez mayores en el caso general. (Tomar muestras aleatorias de matrices con un conmutador extremadamente grande es bastante raro...)

import numpy as np
from scipy.linalg import norm, fractional_matrix_power as powm
from scipy.stats import ortho_group, wishart  # to sample orthogonal/ spd matrices

alphas = np.linspace(0.5, 5, 10)
N = 100 # sample size
eps = 10**-8
n=6  # matrix size

print("Commutative case")
# using simultaneous diagonalizable => commuting
for a in alphas:
    r = 0
    for _ in range(N):
        D = np.diag(np.random.rand(n))
        S = np.diag(np.random.rand(n))
        U = ortho_group.rvs(n)
        X = U.T @ D @ U
        DX = eps* U.T@ S @ U
        num = powm(X+DX, a) - powm(X, a)  # numerical derivative
        ana = a*powm(X, a-1) @ DX         # formula
        r =max(r, norm( num-ana))
    print(F"alpha: {a:.2f}, max residual {r}")
# residuals should be numerically close to zero   

print("General case")
for a in alphas:
    r = 0
    for _ in range(N):
        X = wishart(scale=np.eye(n)).rvs()
        DX= eps*wishart(scale=np.eye(n)).rvs()
        num = powm(X+DX, a) - powm(X, a) # numerical derivative
        ana = a*powm(X, a-1) @ DX        # formula
        r =max(r, norm( num-ana))
    print(F"alpha: {a:.2f}, max residual {r}")
# residuals should be much larger

Answer 3

Para empezar, $$\frac{\partial}{\partial A_{ij}}(AA^T)_{kl}=\frac{\partial}{\partial A_{ij}}(A_{km}A_{lm})=\delta_{ik}\delta_{jm}A_{lm}+A_{km}\delta_{il}\delta_{jm}=\delta_{ik}A_{lj}+\delta_{il}A_{kj}.$$ A continuación $X:=\sqrt{AA^T}$ así que $(AA^T)_{kl}=X_{kp}X_{pl}$ Así que $$\frac{\partial X_{kp}}{\partial A_{ij}}X_{pl}+\frac{\partial X_{pl}}{\partial A_{ij}}X_{kp}=\delta_{ik}A_{lj}+\delta_{il}A_{kj}.$$ Para evaluar la derivada $D_{ijrs}:=\frac{\partial X_{rs}}{\partial A_{ij}}$ Obsérvese que $$D_{ijrs}(\delta_{kr}X_{sl}+\delta_{sl}X_{kr})=\delta_{ik}A_{lj}+\delta_{il}A_{kj}.$$ Por tanto, el resultado es de la forma $B^{-1}C$ pero con $B,\,C$ como rank- $4$ tensores en lugar de matrices, y su multiplicación e inversión convenientemente definidas.