Square root of a matrix

If D is a diagonal n × n matrix ${\displaystyle D=\operatorname {diag} (\lambda _{1},\dots ,\lambda _{n})}$, then some of its square roots are diagonal matrices ${\displaystyle \operatorname {diag} (\mu _{1},\dots ,\mu _{n})}$, where ${\displaystyle \mu _{i}=\pm {\sqrt {\lambda _{i}}}}$. If the diagonal elements of D are real and non-negative then it is positive semidefinite, and if the square roots are taken with non-negative sign, the resulting matrix is the principal root of D. A diagonal matrix may have additional non-diagonal roots if some entries on the diagonal are equal, as exemplified by the identity matrix above.

If U is an upper triangular matrix (meaning its entries are ${\displaystyle u_{i,j}=0}$ for ${\displaystyle i>j}$) and at most one of its diagonal entries is zero, then one upper triangular solution of the equation ${\displaystyle B^{2}=U}$ can be found as follows. Since the equation ${\displaystyle u_{i,i}=b_{i,i}^{2}}$ should be satisfied, let ${\displaystyle b_{i,i}}$ be the principal square root of the complex number ${\displaystyle u_{i,i}}$. By the assumption ${\displaystyle u_{i,i}\neq 0}$, this guarantees that ${\displaystyle b_{i,i}+b_{j,j}\neq 0}$ for all i,j (because the principal square roots of complex numbers all lie on one half of the complex plane). From the equation

${\displaystyle u_{i,j}=b_{i,i}b_{i,j}+b_{i,i+1}b_{i+1,j}+b_{i,i+2}b_{i+2,j}+\dots +b_{i,j}b_{j,j}}$

we deduce that ${\displaystyle b_{i,j}}$ can be computed recursively for ${\displaystyle j-i}$ increasing from 1 to n as:

${\displaystyle b_{i,j}={\frac {1}{b_{i,i}+b_{j,j}}}\left(u_{i,j}-b_{i,i+1}b_{i+1,j}-b_{i,i+2}b_{i+2,j}-\dots -b_{i,j-1}b_{j-1,j}\right).}$

If U is upper triangular but has multiple zeroes on the diagonal, then a square root might not exist, as exemplified by ${\displaystyle \left({\begin{smallmatrix}0&1\\0&0\end{smallmatrix}}\right)}$. Note the diagonal entries of a triangular matrix are precisely its eigenvalues (see Triangular matrix#Properties).

An n × n matrix A is diagonalizable if there is a matrix V and a diagonal matrix D such that A = VDV⁻¹. This happens if and only if A has n eigenvectors which constitute a basis for Cⁿ. In this case, V can be chosen to be the matrix with the n eigenvectors as columns, and thus a square root of A is

${\displaystyle R=VSV^{-1}~,}$

where S is any square root of D. Indeed,

${\displaystyle \left(VD^{\frac {1}{2}}V^{-1}\right)^{2}=VD^{\frac {1}{2}}\left(V^{-1}V\right)D^{\frac {1}{2}}V^{-1}=VDV^{-1}=A~.}$

For example, the matrix ${\displaystyle A=\left({\begin{smallmatrix}33&24\\48&57\end{smallmatrix}}\right)}$ can be diagonalized as VDV⁻¹, where

${\displaystyle V={\begin{pmatrix}1&1\\2&-1\end{pmatrix}}}$ and ${\displaystyle D={\begin{pmatrix}81&0\\0&9\end{pmatrix}}}$.

D has principal square root

${\displaystyle D^{\frac {1}{2}}={\begin{pmatrix}9&0\\0&3\end{pmatrix}}}$,

giving the square root

${\displaystyle A^{\frac {1}{2}}=VD^{\frac {1}{2}}V^{-1}={\begin{pmatrix}5&2\\4&7\end{pmatrix}}}$.

When A is symmetric, the diagonalizing matrix V can be made an orthogonal matrix by suitably choosing the eigenvectors (see spectral theorem). Then the inverse of V is simply the transpose, so that

${\displaystyle R=VSV^{\textsf {T}}~.}$

Every complex-valued square matrix ${\displaystyle A}$, regardless of diagonalizability, has a Schur decomposition given by ${\displaystyle A=QUQ^{*}}$ where ${\displaystyle U}$ is upper triangular and ${\displaystyle Q}$ is unitary (meaning ${\displaystyle Q^{*}=Q^{-1}}$). The eigenvalues of ${\displaystyle A}$ are exactly the diagonal entries of ${\displaystyle U}$; if at most one of them is zero, then the following is a square root[7]

${\displaystyle A^{\frac {1}{2}}=QU^{\frac {1}{2}}Q^{*}}$.

where a square root ${\displaystyle U^{\frac {1}{2}}}$ of the upper triangular matrix ${\displaystyle U}$ can be found as described above.

If ${\displaystyle A}$ is positive definite, then the eigenvalues are all positive reals, so the chosen diagonal of ${\displaystyle U^{\frac {1}{2}}}$ also consists of positive reals. Hence the eigenvalues of ${\displaystyle QU^{\frac {1}{2}}Q^{*}}$ are positive reals, which means the resulting matrix is the principal root of ${\displaystyle A}$.

Similarly as for the Schur decomposition, every square matrix ${\displaystyle A}$ can be decomposed as ${\displaystyle A=P^{-1}JP}$ where P is invertible and J is in Jordan normal form.

To see that any complex matrix with positive eigenvalues has a square root of the same form, it suffices to check this for a Jordan block. Any such block has the form λ(I + N) with λ > 0 and N nilpotent. If (1 + z)^1/2 = 1 + a₁ z + a₂ z² + ⋅⋅⋅⋅⋅ is the binomial expansion for the square root (valid in |z| < 1), then as a formal power series its square equals 1 + z. Substituting N for z, only finitely many terms will be non-zero and S = √λ (I + a₁ N + a₂ N² + ⋅⋅⋅⋅⋅) gives a square root of the Jordan block with eigenvalue √λ.

It suffices to check uniqueness for a Jordan block with λ = 1. The square constructed above has the form S = I + L where L is polynomial in N without constant term. Any other square root T with positive eigenvalues has the form T = I + M with M nilpotent, commuting with N and hence L. But then 0 = S² − T² = 2(L − M)(I + (L + M)/2). Since L and M commute, the matrix L + M is nilpotent and I + (L + M)/2 is invertible with inverse given by a Neumann series. Hence L = M.

If A is a matrix with positive eigenvalues and minimal polynomial p(t), then the Jordan decomposition into generalized eigenspaces of A can be deduced from the partial fraction expansion of p(t)⁻¹. The corresponding projections onto the generalized eigenspaces are given by real polynomials in A. On each eigenspace, A has the form λ(I + N) as above. The power series expression for the square root on the eigenspace show that the principal square root of A has the form q(A) where q(t) is a polynomial with real coefficients.

Recall the formal power series ${\displaystyle (1-z)^{1/2}=1-\sum _{n=1}^{\infty }\left|{1/2 \choose n}\right|z^{n}}$, which converges provided ${\displaystyle \|z\|\leq 1}$ (since the coefficients of the power series are summable). Plugging in ${\displaystyle z=I-A}$ into this expression yields

${\displaystyle A^{1/2}:=I-\sum _{n=1}^{\infty }\left|{1/2 \choose n}\right|(I-A)^{n}~}$

provided that ${\displaystyle \limsup _{n}\|(I-A)^{n}\|^{1/n}<1}$. By virtue of Gelfand formula, that condition is equivalent to the requirement that the spectrum of ${\displaystyle A}$ is contained within the disk ${\displaystyle D(1,1)\subseteq \mathbb {C} }$. This method of defining or computing ${\displaystyle A^{1/2}}$ is especially useful in the case where ${\displaystyle A}$ is positive semi-definite. In that case, we have ${\displaystyle \|I-A/\|A\|\|\leq 1}$ and therefore ${\displaystyle \|(I-A/\|A\|)^{n}\|\leq \|I-A/\|A\|\|^{n}\leq 1}$, so that the expression ${\displaystyle \|A\|^{1/2}\left(I-\sum _{n=1}^{\infty }\left|{1/2 \choose n}\right|(I-A/\|A\|)^{n}\right)}$ defines a square root of ${\displaystyle A}$ which moreover turns out to be the unique positive semi-definite root. This method remains valid to define square roots of operators on infinite-dimensional Banach or Hilbert spaces or certain elements of (C*) Banach algebras.

Another way to find the square root of an n × n matrix A is the Denman–Beavers square root iteration.[8]

Let Y₀ = A and Z₀ = I, where I is the n × n identity matrix. The iteration is defined by

${\displaystyle {\begin{aligned}Y_{k+1}&={\frac {1}{2}}\left(Y_{k}+Z_{k}^{-1}\right),\\Z_{k+1}&={\frac {1}{2}}\left(Z_{k}+Y_{k}^{-1}\right).\end{aligned}}}$

As this uses a pair of sequences of matrix inverses whose later elements change comparatively little, only the first elements have a high computational cost since the remainder can be computed from earlier elements with only a few passes of a variant of Newton's method for computing inverses,

${\displaystyle X_{n+1}=2X_{n}-X_{n}BX_{n}.}$

With this, for later values of k one would set ${\displaystyle X_{0}=Z_{k-1}^{-1}}$ and ${\displaystyle B=Z_{k},}$ and then use ${\displaystyle Z_{k}^{-1}=X_{n}}$ for some small n (perhaps just 1), and similarly for ${\displaystyle Y_{k}^{-1}.}$

Convergence is not guaranteed, even for matrices that do have square roots, but if the process converges, the matrix ${\displaystyle Y_{k}}$ converges quadratically to a square root A^1/2, while ${\displaystyle Z_{k}}$ converges to its inverse, A^−1/2.

Yet another iterative method is obtained by taking the well-known formula of the Babylonian method for computing the square root of a real number, and applying it to matrices. Let X₀ = I, where I is the identity matrix. The iteration is defined by

${\displaystyle X_{k+1}={\frac {1}{2}}\left(X_{k}+AX_{k}^{-1}\right)\,.}$

Again, convergence is not guaranteed, but if the process converges, the matrix ${\displaystyle X_{k}}$ converges quadratically to a square root A^1/2. Compared to Denman–Beavers iteration, an advantage of the Babylonian method is that only one matrix inverse need be computed per iteration step. On the other hand, as Denman–Beavers iteration uses a pair of sequences of matrix inverses whose later elements change comparatively little, only the first elements have a high computational cost since the remainder can be computed from earlier elements with only a few passes of a variant of Newton's method for computing inverses (see Denman–Beavers iteration above); of course, the same approach can be used to get the single sequence of inverses needed for the Babylonian method. However, unlike Denman–Beavers iteration, the Babylonian method is numerically unstable and more likely to fail to converge.[1]

The Babylonian method follows from Newton's method for the equation ${\displaystyle X^{2}-A=0}$ and using ${\displaystyle AX_{k}=X_{k}A}$ for all ${\displaystyle k}$.[9]

If T is a non-negative operator on a finite-dimensional Hilbert space, then all square roots of T are related by unitary transformations. More precisely, if T = A*A = B*B, then there exists a unitary U such that A = UB.

Indeed, take B = T^½ to be the unique non-negative square root of T. If T is strictly positive, then B is invertible, and so U = AB⁻¹ is unitary:

${\displaystyle {\begin{aligned}U^{*}U&=\left(\left(B^{*}\right)^{-1}A^{*}\right)\left(AB^{-1}\right)=\left(B^{*}\right)^{-1}T\left(B^{-1}\right)\\&=\left(B^{*}\right)^{-1}B^{*}B\left(B^{-1}\right)=I.\end{aligned}}}$

If T is non-negative without being strictly positive, then the inverse of B cannot be defined, but the Moore–Penrose pseudoinverse B⁺ can be. In that case, the operator B⁺A is a partial isometry, that is, a unitary operator from the range of T to itself. This can then be extended to a unitary operator U on the whole space by setting it equal to the identity on the kernel of T. More generally, this is true on an infinite-dimensional Hilbert space if, in addition, T has closed range. In general, if A, B are closed and densely defined operators on a Hilbert space H, and A* A = B* B, then A = UB where U is a partial isometry.

Square roots, and the unitary freedom of square roots, have applications throughout functional analysis and linear algebra.

If A is an invertible operator on a finite-dimensional Hilbert space, then there is a unique unitary operator U and positive operator P such that

${\displaystyle A=UP;\,}$

this is the polar decomposition of A. The positive operator P is the unique positive square root of the positive operator A^∗A, and U is defined by U = AP⁻¹.

If A is not invertible, then it still has a polar composition in which P is defined in the same way (and is unique). The unitary operator U is not unique. Rather it is possible to determine a "natural" unitary operator as follows: AP⁺ is a unitary operator from the range of A to itself, which can be extended by the identity on the kernel of A^∗. The resulting unitary operator U then yields the polar decomposition of A.

By Choi's result, a linear map

${\displaystyle \Phi :C^{n\times n}\rightarrow C^{m\times m}}$

is completely positive if and only if it is of the form

${\displaystyle \Phi (A)=\sum _{i}^{k}V_{i}AV_{i}^{*}}$

where k ≤ nm. Let {E_pq} ⊂ C^{n × n} be the n² elementary matrix units. The positive matrix

${\displaystyle M_{\Phi }=\left(\Phi \left(E_{pq}\right)\right)_{pq}\in C^{nm\times nm}}$

is called the Choi matrix of Φ. The Kraus operators correspond to the, not necessarily square, square roots of M_Φ: For any square root B of M_Φ, one can obtain a family of Kraus operators V_i by undoing the Vec operation to each column b_i of B. Thus all sets of Kraus operators are related by partial isometries.

In quantum physics, a density matrix for an n-level quantum system is an n × n complex matrix ρ that is positive semidefinite with trace 1. If ρ can be expressed as

${\displaystyle \rho =\sum _{i}p_{i}v_{i}v_{i}^{*}}$

where ${\displaystyle p_{i}>0}$ and ∑ p_i = 1, the set

${\displaystyle \left\{p_{i},v_{i}\right\}}$

is said to be an ensemble that describes the mixed state ρ. Notice {v_i} is not required to be orthogonal. Different ensembles describing the state ρ are related by unitary operators, via the square roots of ρ. For instance, suppose

${\displaystyle \rho =\sum _{j}a_{j}a_{j}^{*}.}$

The trace 1 condition means

${\displaystyle \sum _{j}a_{j}^{*}a_{j}=1.}$

Let

${\displaystyle p_{i}=a_{i}^{*}a_{i},}$

and v_i be the normalized a_i. We see that

${\displaystyle \left\{p_{i},v_{i}\right\}}$

gives the mixed state ρ.

In the Unscented Kalman Filter (UKF),[10] the square root of the state error covariance matrix is required for the unscented transform which is the statistical linearization method used. A comparison between different matrix square root calculation methods within a UKF application of GPS/INS sensor fusion was presented, which indicated that the Cholesky decomposition method was best suited for UKF applications.[11]