# Frechet derivatives and Taylor expansion of a rational power of a matrix

Let $$A$$ and $$B$$ be two real square non-commuting matrices ($$A,B \in \mathbb{R}^{N \times N}$$ with $$[A,B] = AB-BA \ne 0$$). I also assume that $$A$$ is positive-definite. Consider the power function $$f : \mathbb{R}^{N \times N} \to \mathbb{R}^{N \times N}, \quad X \mapsto X^q.$$ When $$q \in \mathbb{N}$$ it is simply multiplying $$X$$ by itself $$q$$ times. I am interested in the case when $$q \in \mathbb{Q}$$ is fractional, e.g. $$q = \frac{1}{2}$$.

Although probably the easiest way to compute $$(A+B)^q$$ is by eigen-decomposition and applying $$\lambda_i \mapsto \lambda_i^q$$ on the eigenvalues, I want to use a different approach to calculate it by using Taylor expansion around $$A$$ where $$B$$ is sufficiently small (assume I know $$A^{1/2}$$). In this paper the following definition of Taylor expansion is given: $$f(A+B) = \sum_{n=0}^{\infty} \frac{1}{n!} D_f^{[n]}(A,B)$$ where $$D_f^{[n]}(A,B) = \left. \frac{d^n}{dt^n}\right|_{t=0} f(A+tB)$$ is the Frechet derivative.

This is the first time I encounter Frechet derivative. I try to read some about it (e.g. Wikipedia, or this paper, but they lack explicit examples). The last source offers an algorithm to compute them order by order, but I think it is beyond the the level of desired solution (which I want to implement in Python).

However I don't how to explicitly compute these derivatives for $$f(X)=X^q$$ when $$q \in \mathbb{Q}-\mathbb{N}$$, in particular for $$q=1/2$$. For the 1st derivative I tried to use the definition: $$D_f^{[1]}(A,B) = \left( (A+tB)^q - A^q \right) + o(t)$$ but I don't know how to expand $$(A+tB)^q$$ around $$A$$ when $$A$$ and $$B$$ do not commute. When $$A$$ and $$B$$ commute then $$D_f^{[n]}(A,B) = f^{(n)}(A) B^n$$ where $$f^{(n)}(x) = \frac{d^n f(x)}{d x^n}$$ is the scalar n-th derivative of $$f$$ (I assume $$f$$ is infinitely differentiable in $$X = A$$). I wonder if there is a closed form or formula to these derivatives when $$f(X) = X^q$$.

For the sake of completeness, when $$q=1/2$$ and $$a,b \in \mathbb{R}^+$$ are scalars then the Taylor expansion is $$(a+b)^{1/2} = a^{1/2} + \sum_{n=1}^{\infty} \binom{1/2}{n} a^{\frac{1}{2}-n} b^n$$

Edit: the first Frechet derivative of $$Y = X^{1/2}$$ can be computed as follows: the Frechet derivative of $$X = Y^2$$ in $$E$$ is obtained via the definition $$L_{y^2}(Y,E) = (Y+E)^2 - Y^2 = Y^2 + YE + EY + E^2 - Y^2 = YE + EY + o(\| E \|)$$ and since the Frechet derivative of the inverse of $$X = Y^2$$ is the inverse of $$L_{y^2}(Y,E)$$, that is: $$L_{y^2}(Y,L_{x^{1/2}}(X,E)) = E$$ one need to solve the Sylvester equation $$X^{1/2} L + L X^{1/2} = E$$ where $$L = L_{x^{1/2}}(X,E)$$ is the desired 1st Frechet derivative. There is a known algorithm to solve it and it is even implemented in Python's SciPy package/extension.

But how do I compute higher Frechet derivatives in this case ($$X \mapsto X^q$$ in particular for $$q=1/2$$)?

So to summarize my question(s):

1. How to compute the Taylor expansion $$(A+B)^q$$ where $$A > 0$$ and $$B$$ is small but $$[A,B] \ne 0$$?
2. If there is no closed form in the general case, is there one at least for $$q=1/2$$?
3. How to compute the Frechet derivative of this Taylor expansion?
• Here ( physics.stackexchange.com/questions/196448/… ) an answer suggests a formula for the Taylor series $(A+B)^{1/2}$ - Eq.(2) - but this formula is incorrect (I also check it numerically via Python code). Commented May 3, 2020 at 16:26
• See also this question math.stackexchange.com/questions/3651527/… for the specific matrix form of the Taylor expansion of $(A+B)^{1/2}$ when $[A,B] = AB-BA \ne 0$. Commented May 5, 2020 at 9:32
• This question math.stackexchange.com/questions/3371209/… also deals with Taylor expansion for non-commutative matrices but for $f(X) = e^X$ there is a known answer (though not a closed form) in the literature. Commented May 6, 2020 at 11:33

Assume that $$A\in M_n$$ is $$>0$$ real symmetric and $$H\in M_n$$ is a small real symmetric matrix. Let $$q$$ be a positive integer and $$A^{1/q}=B$$. We search an approximation of $$(A+H)^{1/q}$$ and we know $$A,B,H$$. Let $$f:X\mapsto X^{1/q}$$.

Then $$f(A+H)=B+Df_A(H)+O(||H||^2)$$, where $$df_A(H)=K$$ satisfies

$$H=KB^{q-1}+BKB^{q-2}+\cdots+B^{q-1}K$$.

If $$q>2$$, we are dealing with a generalized Sylvester equation in the unknown $$K$$. If $$q=2$$, this is the standard Sylvester equation -it admits a sole solution $$K$$ that is symmetric-. Now, we write about this last case

Several algo. can solve the equation $$H=KB+BK$$ with complexity $$\approx 20n^3$$.

$$\textbf{Remark.}$$ We can also diagonalize $$A+H=PDP^T$$; then

$$f(A+H)=PD^{1/2}P^T$$; the complexity is approximately the same as above.