Square root¶

A matrix \(B\) satisfying

\[B^2 = A\]

is referred to as the square-root of the matrix \(A\). Such a matrix can easily be seen to exist when \(A\) is diagonalizable: if \(A = X \Lambda X^{-1}\), then we may put

\[B = X \sqrt{\Lambda} X^{-1},\]

where each eigenvalue \(\lambda = r e^{i\theta}\) maps to \(\sqrt{\lambda} = \sqrt{r} e^{i\theta/2}\).

General¶

The current implementation for \(n \times n\) matrices uses the Newton iteration

\[X_{k+1} := \frac{1}{2} ( X_k + X_k^{-1} A )\]

and convergence is determined to occur when

\[\| X_{k+1} - X_k \|_1 \le \| X_{k+1} \|_1^{q+1} \text{tol},\]

where the exponent \(q\) is typically one, and the relative tolerance, tol, defaults to \(n \epsilon\), where \(\epsilon\) is the machine precision. Please see Nicholas J. Higham and Awad H. Al-Mohy’s Computing Matrix Functions for more details.

C++ API¶

void SquareRoot(Matrix<F> &A, const SquareRootCtrl<Base<F>> ctrl = SquareRootCtrl<Base<F>>())¶

void SquareRoot(AbstractDistMatrix<F> &A, const SquareRootCtrl<Base<F>> ctrl = SquareRootCtrl<Base<F>>())¶

template<> type SquareRootCtrl<Real>¶

Int maxIts¶: The maximum number of Newton iterations

Real tol¶: The relative tolerance for convergence of the Newton iteration

Real power¶: The power of the one norm of the new iterate that should be used to scale the error when checking for convergence

bool progress¶: Whether or not to print convergence progress

SquareRootCtrl()¶: Sets maxIts=100, tol=0 (which signals that matrix-dependent value will be computed), power=1, and progress=false.

C API¶

ElError ElSquareRoot_s(ElMatrix_s A)¶

ElError ElSquareRoot_d(ElMatrix_d A)¶

ElError ElSquareRoot_c(ElMatrix_c A)¶

ElError ElSquareRoot_z(ElMatrix_z A)¶

ElError ElSquareRootDist_s(ElDistMatrix_s A)¶

ElError ElSquareRootDist_d(ElDistMatrix_d A)¶

ElError ElSquareRootDist_c(ElDistMatrix_c A)¶

ElError ElSquareRootDist_z(ElDistMatrix_z A)¶

TODO: Expert versions

Hermitian Positive Semi-Definite¶

Computes the Hermitian EVD, square-roots the eigenvalues, and then reforms the matrix. If any of the eigenvalues were sufficiently negative, a NonHPSDMatrixException is thrown.

C++ API¶

void HPSDSquareRoot(UpperOrLower uplo, Matrix<F> &A)¶

void HPSDSquareRoot(UpperOrLower uplo, AbstractDistMatrix<F> &A)¶

C API¶

ElError ElHPSDSquareRoot_s(ElUpperOrLower uplo, ElMatrix_s A)¶

ElError ElHPSDSquareRoot_d(ElUpperOrLower uplo, ElMatrix_d A)¶

ElError ElHPSDSquareRoot_c(ElUpperOrLower uplo, ElMatrix_c A)¶

ElError ElHPSDSquareRoot_z(ElUpperOrLower uplo, ElMatrix_z A)¶

ElError ElHPSDSquareRootDist_s(ElUpperOrLower uplo, ElDistMatrix_s A)¶

ElError ElHPSDSquareRootDist_d(ElUpperOrLower uplo, ElDistMatrix_d A)¶

ElError ElHPSDSquareRootDist_c(ElUpperOrLower uplo, ElDistMatrix_c A)¶

ElError ElHPSDSquareRootDist_z(ElUpperOrLower uplo, ElDistMatrix_z A)¶

TODO: HermitianSquareRoot