View source for Square root of a matrix

{{Use American English|date = March 2019}}
{{Short description|Mathematical operation}}
In [[mathematics]], the '''square root of a matrix''' extends the notion of [[square root]] from numbers to [[Matrix (mathematics)|matrices]]. A matrix {{mvar|B}} is said to be a square root of {{mvar|A}} if the [[matrix product]] {{math|''BB''}} is equal to  {{mvar|A}}.<ref name="Higham">{{Citation|last=Higham|first=Nicholas J.|author-link=Nicholas Higham|date=April 1986|title=Newton's Method for the Matrix Square Root|journal=[[Mathematics of Computation]]|volume=46|issue=174|pages=537–549|doi=10.2307/2007992|url=http://www.ams.org/journals/mcom/1986-46-174/S0025-5718-1986-0829624-5/S0025-5718-1986-0829624-5.pdf|jstor=2007992}}</ref>

Some authors use the name ''square root'' or the notation {{math|''A''<sup>1/2</sup>}} only for the specific case when {{mvar|A}} is [[positive semidefinite matrix|positive semidefinite]], to denote the unique matrix {{mvar|B}} that is positive semidefinite and such that {{math|1=''BB'' = ''B''<sup>T</sup>''B'' = ''A''}} (for real-valued matrices, where {{math|''B''<sup>T</sup>}} is the [[matrix transpose|transpose]] of {{mvar|B}}).

Less frequently, the name ''square root'' may be used for any factorization of a positive semidefinite matrix {{mvar|A}} as {{mvar|1=''B''<sup>T</sup>''B'' = ''A''}}, as in the [[Cholesky factorization]], even if {{math|''BB'' ≠ A}}. This distinct meaning is discussed in ''{{section link|Positive definite matrix#Decomposition}}''.

==Examples==

In general, a matrix can have several square roots. In particular, if <math>A = B^2</math> then <math>A=(-B)^2</math> as well.

For example, the 2×2 [[identity matrix]] <math>\textstyle\begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix} </math> has infinitely many square roots. They are given by
:<math>\begin{bmatrix} \pm 1 &  ~~0\\  ~~0 & \pm 1\end{bmatrix}</math> and <math>\begin{bmatrix} a &  ~~b\\  c & -a\end{bmatrix}</math>
where <math>(a, b, c)</math> are any numbers (real or complex) such that <math> a^2+ bc = 1 ~.</math>
In particular if <math>\ (r,s,t)\ </math> is any [[Pythagorean triple]]{{efn|group=Note|
A [[Pythagorean triple]] is any triplet of positive integers such that <math>\ r^2 + s^2 = t^2 ~.</math>
}} then
<math>\ \frac{1}{t}\begin{bmatrix}r & ~~s\\ s & -r\end{bmatrix}\ </math> is one of the matrix square roots of <math>\ I_{2}\ </math> which happens to be symmetric and has rational entries.<ref>{{cite journal |last=Mitchell |first=Douglas W. |date=November 2003 |title=Using Pythagorean triples to generate square roots of {{mvar|I}}{{sub|2}} |journal=[[The Mathematical Gazette]] |volume=87 |issue=510 |pages=499–500 |doi=10.1017/s0025557200173723 |doi-access=free }}</ref>
Thus
: <math>\ \begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix} = \begin{bmatrix}0 & 1\\ 1 & 0\end{bmatrix}^2 = \begin{bmatrix}\frac{4}{5} & ~~\frac{3}{5}\\ \frac{3}{5} & -\frac{4}{5}\end{bmatrix}^2 ~.</math>

Minus {{mvar|I}}{{sub|2}} also has a square root, for example:
: <math>\ -\begin{bmatrix} 1 & 0\\ 0 & 1 \end{bmatrix} = \begin{bmatrix} 0 & -1 \\ 1 & ~~0 \end{bmatrix}^2\ ,</math>
which can be used to represent the [[imaginary unit]] <math>\ i\ </math> and hence all [[complex numbers]] using 2×2 real matrices, see [[Complex numbers#Matrix representation of complex numbers|matrix representation of complex numbers]].

Just as with the [[real numbers]], a real matrix may fail to have a real square root, but have a square root with [[Complex number|complex]]-valued entries.
Some matrices have no square root. An example is the matrix <math>\ \begin{bmatrix}0 & 1\\ 0 & 0\end{bmatrix} ~.</math>

Notice that some ideas from [[number theory]] do not carry over to matrices: The square root of a nonnegative [[integer]] must either be another integer or an [[irrational number]], excluding non-integer [[rational number|rationals]]. Contrast that to a matrix of integers, which can have a square root whose entries are all non-integer rational numbers, as demonstrated in some of the above examples.

== Positive semidefinite matrices==
{{See also|Positive definite matrix#Decomposition}}
A symmetric real ''n'' × ''n'' matrix <math>A</math> is called ''[[positive semidefinite matrix|positive semidefinite]]'' if <math>x^\textsf{T} A x \geq 0</math> for all <math>x \in \mathbb{R}^n</math> (here <math>x^\textsf{T}</math> denotes the [[matrix transpose|transpose]], changing a column vector {{mvar|x}} into a row vector).
A square real matrix <math>A</math> is positive semidefinite if and only if <math>A = B^\textsf{T} B</math> for some matrix {{mvar|B}}.
There can be many different such matrices {{mvar|B}}.
A positive semidefinite matrix {{mvar|A}} can also have many matrices {{mvar|B}} such that <math>A = B B</math>.
However, {{mvar|A}} always has precisely one square root {{mvar|B}} that is ''both'' positive semidefinite and symmetric.
In particular, since {{mvar|B}} is required to be symmetric, <math>B=B^\textsf{T}</math>, so the two conditions <math>A = B B</math> or <math>A = B^\textsf{T} B</math> are equivalent.

For complex-valued matrices, the [[conjugate transpose]] <math>B^*</math> is used instead and positive semidefinite matrices are [[Hermitian matrix|Hermitian]], meaning <math>B^*=B</math>. 
{{Math theorem
 |name=Theorem<ref name="HJ-7.2.6">{{harvtxt|Horn|Johnson|2013}}, p. 439, Theorem 7.2.6 with <math>k = 2</math></ref>
 | Let {{mvar|A}} be a positive semidefinite matrix that is also symmetric.{{efn|group=Note|
Note that positive semidefinite matrices can be asymmetric, so the symmetric condition is a necessary, distinct requirement.
}} Then there is exactly one positive semidefinite and symmetric matrix {{mvar|B}} such that {{nowrap|<math>A = BB</math>}}.{{efn|group=Note|There can be more than one non-symmetric and positive semidefinite matrix <math>B</math> such that <math>A = B^TB ~.</math>}}
}}

This unique matrix is called the '''principal''', '''non-negative''', or '''positive square root''' (the latter in the case of [[positive definite matrix|positive definite matrices]]).

The principal square root of a real positive semidefinite matrix is real.<ref name="HJ-7.2.6"/>
The principal square root of a positive definite matrix is positive definite; more generally, the rank of the principal square root of {{mvar|A}} is the same as the rank of {{mvar|A}}.<ref name="HJ-7.2.6"/>

The operation of taking the principal square root is continuous on this set of matrices.<ref>{{cite book|last1=Horn|first1=Roger A.|last2=Johnson|first2=Charles R.|title=Matrix analysis|date=1990|publisher=Cambridge Univ. Press|location=Cambridge|isbn=9780521386326|page=411}}</ref> These properties are consequences of the [[holomorphic functional calculus]] applied to matrices.<ref>For analytic functions of matrices, see
*{{harvnb|Higham|2008}} 
*{{harvnb|Horn|Johnson|1994}}</ref><ref>
For the holomorphic functional calculus, see:
*{{harvnb|Rudin|1991}}
*{{harvnb|Bourbaki|2007}}
*{{harvnb|Conway|1990}}
</ref>
The existence and uniqueness of the principal square root can be deduced directly from the [[Jordan normal form]] (see below).

== Matrices with distinct eigenvalues==
An {{math|''n''×''n''}} matrix with {{mvar|n}} ''distinct nonzero eigenvalues'' has 2<sup>''n''</sup> square roots. Such a matrix, {{mvar|A}},  has an [[Eigendecomposition of a matrix|eigendecomposition]] {{math|''VDV''<sup>−1</sup>}} where {{mvar|V}} is the matrix whose columns are eigenvectors of {{mvar|A}} and {{mvar|D}} is the diagonal matrix whose diagonal elements are the corresponding {{mvar|n}} eigenvalues {{math|''λ<sub>i</sub>''}}. Thus the square roots of {{mvar|A}} are given by {{math|''VD''<sup>1/2</sup> ''V''<sup>−1</sup>}}, where {{math|''D''<sup>1/2</sup>}} is any square root matrix of {{mvar|D}}, which,  for distinct eigenvalues, must be diagonal with diagonal elements equal to square roots of the diagonal elements of  {{mvar|D}};  since there are two possible choices for a square root of each diagonal element of {{mvar|D}}, there are 2<sup>''n''</sup> choices for the matrix {{math|''D''<sup>1/2</sup>}}.

This also leads to a proof of the above observation, that a positive-definite matrix has precisely one positive-definite square root: a positive definite matrix has only positive eigenvalues, and each of these eigenvalues has only one positive square root; and since the eigenvalues of the square root matrix are the diagonal elements of {{math|''D''<sup>1/2</sup>}}, for the square root matrix to be itself positive definite necessitates the use of only the unique positive square roots of the original eigenvalues.

== Solutions in closed form ==
{{See also|Square root of a 2 by 2 matrix}}

If a matrix is [[idempotent matrix|idempotent]], meaning <math>A^2 = A</math>, then by definition one of its square roots is the matrix itself.

=== Diagonal and triangular matrices===
If {{mvar|D}} is a [[diagonal matrix|diagonal]] ''n'' × ''n'' matrix <math>D = \operatorname{diag}(\lambda_1,\dots,\lambda_n)</math>,
then some of its square roots are diagonal matrices <math>\operatorname{diag}(\mu_1,\dots,\mu_n)</math>, where <math>\mu_i = \pm \sqrt{\lambda_i}</math>.
If the diagonal elements of {{mvar|D}} are real and non-negative then it is positive semidefinite, and if the square roots are taken with the (+) sign (i.e. all non-negative), the resulting matrix is the principal root of {{mvar|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 {{mvar|U}} is an [[upper triangular matrix]] (meaning its entries are <math>u_{i,j} = 0</math> for <math>i > j</math>) and at most one of its diagonal entries is zero, then one upper triangular solution of the equation <math>B^2 = U</math> can be found as follows.
Since the equation <math>u_{i,i} = b_{i,i}^2</math> should be satisfied, let <math>b_{i,i}</math> be the [[Square root#Principal square root of a complex number|principal square root]] of the complex number <math>u_{i,i}</math>.
By the assumption <math>u_{i,i} \neq 0</math>, this guarantees that <math>b_{i,i} + b_{j,j} \neq 0</math> for all {{mvar|i,j}} (because the principal square roots of complex numbers all lie on one half of the complex plane).
From the equation
: <math>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}</math>
we deduce that <math>b_{i,j}</math> can be computed recursively for <math>j-i</math> increasing from 1 to ''n''-1 as:
: <math>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).</math>

If {{mvar|U}} is upper triangular but has multiple zeroes on the diagonal, then a square root might not exist, as exemplified by <math>\left(\begin{smallmatrix}0 & 1\\ 0 & 0\end{smallmatrix}\right)</math>.
Note the diagonal entries of a triangular matrix are precisely its [[eigenvalue]]s (see [[Triangular matrix#Properties]]).

=== By diagonalization ===
An ''n'' × ''n'' matrix {{mvar|A}} is [[diagonalizable matrix|diagonalizable]] if there is a matrix {{mvar|V}} and a diagonal matrix {{mvar|D}} such that {{math|1=''A'' = ''VDV''<sup>−1</sup>}}. This happens if and only if {{mvar|A}}  has ''n'' [[eigenvector]]s which constitute a basis for {{math|'''C'''<sup>''n''</sup>}}.  In this case, {{mvar|V}} can be chosen to be the matrix with the ''n'' eigenvectors as columns, and thus a square root of {{mvar|A}}  is
: <math>R = V S V^{-1}~,</math>

where {{mvar|S}} is any square root of {{mvar|D}}. Indeed,
: <math>\left(V D^\frac{1}{2} V^{-1}\right)^2 = V D^\frac{1}{2} \left(V^{-1} V\right) D^\frac{1}{2} V^{-1} = V D V^{-1} = A ~.</math>

For example, the matrix <math>A = \left(\begin{smallmatrix} 33 & 24\\ 48 & 57\end{smallmatrix}\right)</math> can be diagonalized as  {{math|''VDV''<sup>−1</sup>}}, where 
: <math>V = \begin{pmatrix} 1 & 1\\ 2 & -1\end{pmatrix}</math> and <math>D = \begin{pmatrix} 81 & 0\\ 0 & 9\end{pmatrix}</math>.

{{mvar|D}} has principal square root 
: <math>D^\frac{1}{2} = \begin{pmatrix} 9 & 0\\ 0 & 3\end{pmatrix}</math>,

giving the square root
: <math>A^\frac{1}{2} = V D^\frac{1}{2} V^{-1} = \begin{pmatrix} 5 & 2\\ 4 & 7\end{pmatrix}</math>.

When {{mvar|A}} is symmetric, the diagonalizing matrix {{mvar|V}} can be made an [[orthogonal matrix]] by suitably choosing the eigenvectors (see [[Eigendecomposition#Decomposition for special matrices|spectral theorem]]).  Then the inverse of  {{mvar|V}} is simply the transpose, so that
:<math>R = V S V^\textsf{T}~.</math>

=== By Schur decomposition ===
Every complex-valued square matrix <math>A</math>, regardless of diagonalizability, has a [[Schur decomposition]] given by <math>A=QUQ^*</math> where <math>U</math> is upper triangular and <math>Q</math> is [[unitary matrix|unitary]] (meaning {{nowrap|<math>Q^* = Q^{-1}</math>).}}
The [[eigenvalue]]s of <math>A</math> [[Triangular matrix#Properties|are exactly]] the diagonal entries of <math>U</math>;
if at most one of them is zero, then the following is a square root<ref>{{Citation|last=Deadman|first=Edvin|title=Blocked Schur Algorithms for Computing the Matrix Square Root|date=2013|work=Applied Parallel and Scientific Computing|pages=171–182|publisher=Springer Berlin Heidelberg|isbn=978-3-642-36802-8|last2=Higham|first2=Nicholas J.|last3=Ralha|first3=Rui|doi=10.1007/978-3-642-36803-5_12|url=http://eprints.maths.manchester.ac.uk/1951/1/paper16.pdf}}</ref>
: <math>A^\frac{1}{2} = Q U^\frac{1}{2} Q^*.</math>
where a square root <math>U^\frac{1}{2}</math> of the upper triangular matrix <math>U</math> can be found as described above.

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

=== By Jordan decomposition ===
Similarly as for the Schur decomposition, every square matrix <math>A</math> can be decomposed as <math>A = P^{-1} J P</math> where {{mvar|P}} is [[invertible matrix|invertible]] and {{mvar|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 matrix|nilpotent]]. If {{math|1=(1 + ''z'')<sup>1/2</sup> = 1 + ''a''<sub>1</sub> ''z'' + ''a''<sub>2</sub>  ''z''<sup>2</sup> + ⋯}} 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 
{{math|1=''S'' = √λ (''I'' + ''a''<sub>1</sub> ''N'' + ''a''<sub>2</sub>  ''N''<sup>2</sup> + ⋯)}} gives a square root of the Jordan block with eigenvalue {{math|1=√λ}}.

It suffices to check uniqueness for a Jordan block with λ = 1. The square constructed above has the form {{nowrap|1=''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 {{mvar|M}} nilpotent, commuting with ''N'' and hence ''L''. But then {{math|1=0 = ''S''<sup>2</sup> − ''T''<sup>2</sup> = 2(''L'' − ''M'')(''I'' + (''L'' + ''M'')/2)}}. Since {{mvar|L}} and {{mvar|M}} commute, the matrix {{math|1=''L'' + ''M''}} is nilpotent and {{math|1=''I'' + (''L'' + ''M'')/2}} is invertible with inverse given by a [[Neumann series]]. Hence {{math|1=''L'' = ''M''}}.

If  {{mvar|A}} is a matrix with positive eigenvalues and [[Minimal polynomial (linear algebra)|minimal polynomial]] {{math|''p''(''t'')}}, then the Jordan decomposition into generalized eigenspaces of {{mvar|A}} can be deduced from the partial fraction expansion of {{math|''p''(''t'')<sup>−1</sup>}}. The corresponding projections onto the generalized eigenspaces are given by real polynomials in  {{mvar|A}}. On each eigenspace, {{mvar|A}} has the form {{math|''λ''(''I'' + ''N'')}} as above. The [[power series]] expression for the square root on the eigenspace show that the principal square root of {{mvar|A}} has the form ''q''(''A'') where ''q''(''t'') is a polynomial with real coefficients.

=== Power series ===
Recall the formal power series <math display="inline">(1 - z)^\frac{1}{2} = \sum_{n=0}^\infty (-1)^n \binom{1/2}{n} z^n</math>, which converges provided <math>\|z\| \leq 1</math> (since the coefficients of the power series are summable). Plugging in <math>z = I - A</math> into this expression yields
: <math>A^\frac{1}{2} := \sum_{n=0}^\infty (-1)^n {\frac{1}{2} \choose n} (I - A)^n </math>

provided that <math display="inline">\limsup_n\|(I - A)^n\|^\frac{1}{n} < 1</math>. By virtue of [[Spectral radius#Gelfand's formula|Gelfand formula]], that condition is equivalent to the requirement that the spectrum of <math> A </math> is contained within the disk <math> D(1,1) \subseteq \mathbb{C} </math>. This method of defining or computing <math>A^\frac{1}{2}</math> is especially useful in the case where <math> A </math> is positive semi-definite. In that case, we have <math display="inline">\left\|I - \frac{A}{\|A\|}\right\| \leq 1</math> and therefore <math display="inline">\left\|\left(I - \frac{A}{\|A\|}\right)^n\right\| \leq \left\|I - \frac{A}{\|A\|}\right\|^n \leq 1</math>, so that the expression <math display="inline">\|A\|^\frac{1}{2} = \left(\sum_{n=0}^\infty (-1)^n \binom{1/2}{n} \left(I - \frac{A}{\|A\|}\right)^n\right)</math> defines a square root of <math> A </math> 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.

== Iterative solutions ==

=== By Denman–Beavers iteration ===

Another way to find the square root of an {{nowrap|''n'' × ''n''}} matrix ''A'' is the Denman–Beavers square root iteration.<ref>{{harvnb|Denman|Beavers|1976}}; {{harvnb|Cheng|Higham|Kenney|Laub|2001}}</ref>

Let ''Y''<sub>0</sub> = ''A'' and ''Z''<sub>0</sub> = ''I'', where ''I''  is the {{nowrap|''n'' × ''n''}} [[identity matrix]]. The iteration is defined by
:<math> \begin{align}
  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{align} </math>

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 [[Invertible matrix#Newton's method|computing inverses]],
:<math>X_{n+1} = 2X_n - X_n B X_n.</math>

With this, for later values of {{mvar|k}} one would set <math>X_0 = Z_{k-1}^{-1}</math> and <math>B = Z_k,</math> and then use <math>Z_k^{-1} = X_n</math> for some small <math>n</math> (perhaps just 1), and similarly for <math>Y_k^{-1}.</math>

Convergence is not guaranteed, even for matrices that do have square roots, but if the process converges, the matrix <math>Y_k</math> converges quadratically to a square root {{mvar|A}}<sup>1/2</sup>, while <math>Z_k</math> converges to its inverse, {{mvar|A}}<sup>−1/2</sup>.

=== By the Babylonian method ===

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''<sub>0</sub> = ''I'', where ''I''  is the [[identity matrix]]. The iteration is defined by
:<math>X_{k+1} = \frac{1}{2} \left(X_k + A X_k^{-1}\right)\,.</math>

Again, convergence is not guaranteed, but if the process converges, the matrix <math>X_k</math> converges quadratically to a square root ''A''<sup>1/2</sup>. 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 [[Invertible matrix#Newton's method|computing inverses]] (see [[#Denman–Beavers iteration|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.<ref name="Higham"/>

The Babylonian method follows from [[Newton's method]] for the equation <math>X^2-A=0</math> and using <math>AX_k=X_k A</math> for {{nowrap|all <math>k.</math><ref>{{Cite journal|last=Higham|first=Nicholas J.|date=1997|title=Stable iterations for the matrix square root|journal=Numerical Algorithms|volume=15|issue=2|pages=227–242|doi=10.1023/A:1019150005407|bibcode=1997NuAlg..15..227H}}</ref>}}

== Square roots of positive operators ==

{{More citations needed|date=July 2010}}

In [[linear algebra]] and [[operator theory]], given a [[Bounded operator|bounded]] [[positive semidefinite matrix|positive semidefinite operator]] (a non-negative operator) ''T'' on a complex Hilbert space, ''B'' is a square root of ''T'' if ''T'' = ''B* B'', where ''B*'' denotes the [[Hermitian adjoint]] of ''B''. {{Citation needed|date=July 2010}} According to the [[spectral theorem]], the [[continuous functional calculus]] can be applied to obtain an operator ''T''<sup>1/2</sup> such that
''T''<sup>1/2</sup> is itself positive and (''T''<sup>1/2</sup>)<sup>2</sup> = ''T''. The operator ''T''<sup>1/2</sup> is the '''unique non-negative square root''' of ''T''.  {{Citation needed|date=July 2010}}

A bounded non-negative operator on a complex Hilbert space is self adjoint by definition. So {{nowrap|1=''T'' = (''T''<sup>1/2</sup>)* ''T''<sup>1/2</sup>.}} Conversely, it is trivially true that every operator of the form ''B* B'' is non-negative. Therefore, an operator ''T'' is non-negative [[if and only if]] ''T'' = ''B* B'' for some ''B'' (equivalently, ''T'' = ''CC*'' for some ''C'').

The [[Cholesky factorization]] provides another particular example of square root, which should not be confused with the unique non-negative square root.

=== Unitary freedom of square roots ===

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 matrix|unitary]] ''U'' such that ''A'' = ''UB''.

Indeed, take ''B'' = ''T''<sup>{{sfrac|1|2}}</sup> to be the unique non-negative square root of ''T''. If ''T'' is strictly positive, then ''B'' is invertible, and so {{nowrap|1=''U'' = ''AB''<sup>−1</sup>}} is unitary:
:<math>\begin{align}
  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{align}</math>

If ''T'' is non-negative without being strictly positive, then the inverse of ''B'' cannot be defined, but the [[Moore–Penrose pseudoinverse]] ''B''<sup>+</sup> can be.  In that case, the operator {{nowrap|1=''B''<sup>+</sup>''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 theorem|closed range]]. In general, if ''A'', ''B'' are [[unbounded operator|closed and densely defined operators]] on a Hilbert space ''H'', and ''A* A'' = ''B* B'', then ''A = UB'' where ''U'' is a partial isometry.

=== Some applications ===

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

==== Polar decomposition ====
{{Main|Polar decomposition}}
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
:<math>A = UP;</math>
this is the polar decomposition of ''A''.  The positive operator ''P'' is the unique positive square root of the positive operator ''A''<sup>∗</sup>''A'', and ''U'' is defined by {{nowrap|1=''U'' = ''AP''<sup>−1</sup>}}.

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''<sup>+</sup> is a unitary operator from the range of ''A'' to itself, which can be extended by the identity on the kernel of ''A''<sup>∗</sup>.  The resulting unitary operator ''U'' then yields the polar decomposition of ''A''.

==== Kraus operators ====
{{Main|Choi's theorem on completely positive maps}}

By Choi's result, a linear map

:<math>\Phi : C^{n \times n} \to C^{m \times m}</math>

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

:<math>\Phi(A) = \sum_i ^k V_i A V_i^*</math>

where ''k'' ≤ ''nm''. Let {''E<sub>pq</sub>''} ⊂ '''C'''<sup>''n'' × ''n''</sup> be the ''n''<sup>2</sup> elementary matrix units. The positive matrix

:<math>M_\Phi = \left(\Phi \left(E_{pq}\right)\right)_{pq} \in C^{nm \times nm}</math>

is called the ''Choi matrix'' of Φ. The Kraus operators correspond to the, not necessarily square, square roots of ''M''<sub>Φ</sub>: For any square root ''B'' of ''M''<sub>Φ</sub>, one can obtain a family of Kraus operators ''V<sub>i</sub>'' by undoing the Vec operation to each column ''b<sub>i</sub>'' of ''B''. Thus all sets of Kraus operators are related by partial isometries.

==== Mixed ensembles ====
{{Main|Density matrix}}

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

:<math>\rho = \sum_i p_i v_i v_i^*</math>

where <math> p_i > 0 </math> and Σ ''p<sub>i</sub>'' = 1, the set

:<math>\left\{p_i,  v_i\right\}</math>

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

:<math>\rho = \sum_j a_j a_j^*.</math>

The trace 1 condition means

:<math>\sum_j a_j ^* a_j = 1.</math>

Let

:<math>p_i = a_i ^* a_i,</math>

and ''v<sub>i</sub>'' be the normalized ''a<sub>i</sub>''. We see that

:<math>\left\{p_i,  v_i\right\}</math>

gives the mixed state ''ρ''.

==Footnotes==
{{notelist|group=Note}}

==See also==
* [[Matrix function]]
* [[Holomorphic functional calculus]]
* [[Logarithm of a matrix]]
* [[Sylvester's formula]]
* [[Square root of a 2 by 2 matrix]]

==Citations==
{{reflist|25em}}

== References ==

* {{citation|title=Théories spectrales, chapitres 1 et 2|first= Nicolas|last= Bourbaki|publisher= Springer|year= 2007|isbn= 978-3540353317}}
* {{citation|title=A Course in Functional Analysis|series=Graduate Texts in Mathematics|volume=96|first=	John B.|last= Conway|publisher=Springer|year= 1990|pages=199–205|isbn=978-0387972459}}, Chapter IV, Reisz functional calculus
* {{citation | first1 = Sheung Hun | last1 = Cheng | first2 = Nicholas J. | last2 = Higham | author2-link = Nicholas Higham | first3 = Charles S. | last3 = Kenney | first4 = Alan J. | last4 = Laub | title = Approximating the Logarithm of a Matrix to Specified Accuracy | url = http://www.eeweb.ee.ucla.edu/publications/journalAlanLaubajlaub_simax22(4)_2001.pdf | journal = SIAM Journal on Matrix Analysis and Applications | volume = 22 | year = 2001 | pages = 1112&ndash;1125 | doi = 10.1137/S0895479899364015 | issue = 4 | url-status = dead | archive-url = https://web.archive.org/web/20110809202647/https://eeweb.ee.ucla.edu/publications/journalAlanLaubajlaub_simax22(4)_2001.pdf | archive-date = 2011-08-09 | citeseerx = 10.1.1.230.912 }}
* {{citation | first1 = Donald R. | last1 = Burleson | title = Computing the square root of a Markov matrix: eigenvalues and the Taylor series | url = http://www.blackmesapress.com/TaylorSeries.htm}}
* {{citation | first1 = Eugene D. | last1 = Denman | first2 = Alex N. | last2 = Beavers | title = The matrix sign function and computations in systems | journal = Applied Mathematics and Computation | volume = 2 | year = 1976 | pages = 63&ndash;94 | doi = 10.1016/0096-3003(76)90020-5 | issue = 1 }}
* {{Citation | last1=Higham | first1=Nicholas | author1-link=Nicholas Higham | title=Functions of Matrices. Theory and Computation | publisher=[[Society for Industrial and Applied Mathematics|SIAM]] | isbn=978-0-89871-646-7 | year=2008}}
* {{Cite book | last1=Horn | first1=Roger A. | last2=Johnson | first2=Charles R. | title=Matrix Analysis | publisher=[[Cambridge University Press]] | isbn=978-0-521-54823-6 | year=2013 |edition=2nd}}
* {{citation|title=Topics in Matrix Analysis|first1=Roger A.|last1= Horn|first2=Charles R. |last2=Johnson|
publisher=Cambridge University Press|year= 1994|isbn=978-0521467131}}
* {{Rudin Walter Functional Analysis|edition=2}} <!-- {{sfn | Rudin | 1991 | p=}} -->

[[Category:Matrix theory]]