# Computing $\sqrt{A+B}$ where A,B are matrices using Taylor expansion at $X_0 = A$

I read an answer to a question in Physics.StackExchange about how to perturbatively compute $$\sqrt{A+B}$$ where $$A,B$$ are matrices, $$A$$ is positive definite and $$B$$ is "sufficiently small and smaller than $$A$$". In the original question, the poster assumed $$A$$ is diagonal, the answer does not assume this.

Using Taylor expansion the answer's poster wrote the following expression as an answer: $$\sqrt{A+B} = \sqrt{A} + \sum_{n=1}^{\infty}{\binom {1/2} n} \left( (L_B + \mathrm{ad}A )^{n-1} B \right) A^{\frac{1}{2}-n}$$ where $$L_B(X):= BX$$ and $$\mathrm{ad}A(B) := [A,B] = AB-BA$$. Also, $$\binom \alpha k = \frac{\alpha^{\underline k}}{k!} = \frac{\alpha(\alpha-1)(\alpha-2)\cdots(\alpha-k+1)}{k(k-1)(k-2)\cdots 1} \quad\text{for } k\in\mathbb{N} \text{ and arbitrary } \alpha.$$

Other than stating it is based on a Taylor expansion at $$X_0 = A$$ the poster doesn't provide details how to derive this formula (e.g. why the ad operator appears exactly in that form). Without proof, I suspect the above formula may be incorrect.

Moreover, after understanding this formula I wish to implement it in a computer code, particularly in Python. This is not trivial since $$L_B$$ and $$\mathrm{ad}A$$ do not commute. So I will be happy to get tips on that as well.

Edit: I found a way to implement this formula for general $$n$$ using recursion (or equivalently, recursive loop computing one term after another, each iteration applying $$(L_B + \mathrm{ad}A)$$ on the previous iteration's result).

However, for some matrices (in particular, when $$A$$ and $$B$$ do not commute) the expression does not converge, and even worse the error term diverges (and even explodes). It occurs even for matrices $$A$$ with condition number of 36 and determinant of around 3. The matrices are generated randomly as follow (in Python):

N = 2
np.random.seed(seed=1)       # fix a seed
X = np.random.rand(N,N)      # generate random matrix ~ U(0,1)
A = X @ X.T                  # multiplying X by X transposed
A = 10*A/np.linalg.norm(A)   # normalize A and make it "big enough"
Y = np.random.rand(N,N)      # generate random matrix ~ U(0,1)
B = Y @ Y.T                  # multiplying Y by Y transposed
B = 0.01*B/np.linalg.norm(B) # normalize B and rescale to be sufficiently small


The seed is set in order to have the same matrices A,B each time for debugging. These matrices seem to be well-conditioned and positive-definite, so there should be no reason for the algorithm (Taylor expansion around $$A$$) to diverge. Why does it happen?

Edit 2: I think the formula $$\sqrt{A+B} = \sqrt{A} + \sum_{n=1}^{\infty}{\binom {1/2} n} \left( (L_B + \mathrm{ad}A )^{n-1} B \right) A^{\frac{1}{2}-n}$$ is incorrect. First, the total power of $$A$$ is $$1/2 - n + n-1 = -1/2$$ while in the scalar Taylor expansion is $$1/2 - n$$. Second, I checked it numerically in Python.

So what is the correct formula?

• One thing to note is that inverting (non-diagonal) matrices is generally pretty slow, and can also be numerically badly behaved if the matrices are close to singular. A second comment is that since A is positive semi-definite it is Hermitian, so without loss of (much) generality you can assume it is diagonal by doing the calculation in a basis in which it is. Apr 30, 2020 at 10:42
• I think the $ad$ operator appears to account the fact that $A$ and $B$ do not commute in general. However, it makes the expansion of $(B + \mathrm{ad}A)^n$ untractable analytically for general $n$ (i.e. have an explicit formula/expansion). Apr 30, 2020 at 11:33
• It turns out that for my ad-hoc application, taking $n=4$ is enough, and when I implemented the formula in code it worked out just fine. However, the general answer still may be interesting. Apr 30, 2020 at 13:13
• After trying to explicitly develop the first two non-trivial terms, I suspect that the poster's formula is incorrect. What strengthen this suspicion is that it adds powers of $A$ true $ad A$ where they shouldn't be more powers of $A$ other than $A^{1/2-n}$. When $A$ and $B$ commute $(adA)B = [A,B]=0$ and these extra powers of $A$ cancel out. May 2, 2020 at 18:44
• Note that when $XY \ne YX$ in general $(XY)^{1/2} \ne X^{1/2} Y^{1/2}$, so the trick of pulling $A$ outside and expand $(I + A^{-1}B)^{1/2}$ would probably not work. May 3, 2020 at 16:41