Wikipedia gives the expansion of the square root matrix as follows:
$$\frac{A^{1/2}}{\|A\|^{1/2}} = I – \sum_{n = 1}^{\infty}\left\lvert {1/2 \choose n}\right\rvert \left(I-\frac{A}{\|A\|}\right)^n = I – \frac 1 2\left(I-\frac{A}{\|A\|}\right) – \frac 1 8 \left(I-\frac{A}{\|A\|}\right)^2 \,…$$
In this post, I'm taking the matrix norm to be the max of the absolute values of the column-sums of matrix $A$.
I thought I might be able to get an estimate of square root matrix by taking the first two terms and using the third term as an error bound. But a simple numerical example appears to not work. I'm not sure if I'm using the expansion incorrectly, or there's some hypothesis I'm missing. In particular for $A = \begin{pmatrix}
4 & 0 \\
0 & 16
\end{pmatrix}, \|A\| = 16$
I get an error:$$\|A^{1/2}/\|A\|^{1/2} – I + (I-A/\|A\|)/2 \| = 1/8$$
Error bound using the third term: $$\|(I-A/\|A\|)^2/8 \| = 9/128$$
Best Answer
This must work like an ordinary Taylor series because it is an ordinary Taylor series (or at least, an ordinary Taylor series per component). In particular, consider the (matrix-valued) function $$ F(\lambda)=\left(I-\lambda\left(I-\frac{A}{\|A\|}\right)\right)^{1/2}, $$ which is equal to $I$ at $\lambda=0$ and to $\left(\frac{A}{\|A\|}\right)^{1/2}$ at $\lambda=1$. This equation can be expanded in a Taylor series around $\lambda=0$ by taking its derivatives w.r.t. $\lambda$. The first few are $$ \begin{eqnarray} F(\lambda) &=& &&\left(I-\lambda\left(I-\frac{A}{\|A\|}\right)\right)^{1/2} &\implies& \qquad F(0) &=& I \\ F'(\lambda) &=& -\frac{1}{2}\left(I-\frac{A}{\|A\|}\right)&&\left(I-\lambda\left(I-\frac{A}{\|A\|}\right)\right)^{-1/2} &\implies& \qquad F'(0) &=& -\frac{1}{2}\left(I-\frac{A}{\|A\|}\right) \\ F''(\lambda) &=& -\frac{1}{4}\left(I-\frac{A}{\|A\|}\right)^2&&\left(I-\lambda\left(I-\frac{A}{\|A\|}\right)\right)^{-3/2} &\implies& \qquad F''(0) &=& -\frac{1}{4}\left(I-\frac{A}{\|A\|}\right)^2, \end{eqnarray} $$ and so on. And as usual, the expansion $$ F(\lambda)=\sum_{n=0}^{\infty}\frac{F^{(n)}(0)}{n!}\lambda^n $$ converges to the correct value (within its radius of convergence, which is just the smallest of the component-wise radii of convergence), with the error in the truncated series bounded by the maximum absolute value of the first omitted term. So we expect that the error in the linear approximation should be bounded by the maximum absolute value between $\lambda=0$ and $\lambda=1$ of the quadratic term, which is $$D(\lambda)=\frac{1}{8}\left(I-\frac{A}{\|A\|}\right)^2\left(I-\lambda\left(I-\frac{A}{\|A\|}\right)\right)^{-3/2}.$$ (Indeed, the bound should hold for each component of the matrix.) In OP's case, $$I-\frac{A}{\|A\|}=I-\begin{pmatrix} 1/4 & 0 \\ 0 & 1 \end{pmatrix}=\begin{pmatrix} 3/4 & 0 \\ 0 & 0 \end{pmatrix},$$ so $$ D(\lambda)=\frac{1}{8}\begin{pmatrix} 3/4 & 0 \\ 0 & 0 \end{pmatrix}^2 \begin{pmatrix} 1-(3/4)\lambda & 0 \\ 0 & 1 \end{pmatrix}^{-3/2}=\begin{pmatrix} \frac{9}{128}\left(1-(3/4)\lambda\right)^{-3/2} & 0 \\ 0 & 0 \end{pmatrix}. $$ This assumes its maximum value at $\lambda=1$, where the upper-left component is $\frac{9}{16}$ (and the other components are $0$). So the error in the linear term, which as OP calculated is $$ F(0)+F'(0)-\sqrt{\frac{A}{\|A\|}}=I-\frac{1}{2}\left(I-\frac{A}{\|A\|}\right) -\begin{pmatrix} 1/2 & 0 \\ 0 & 1 \end{pmatrix} =\begin{pmatrix} 1/8 & 0 \\ 0 & 0 \end{pmatrix}, $$ fits into the theoretical bounds.