pracma (version 1.8.8)

sqrtm,rootm: Matrix Square and p-th Roots

Description

Computes the matrix square root and matrix p-th root of a nonsingular real matrix.

Usage

sqrtm(A, kmax = 20, tol = .Machine$double.eps^(1/2))
signm(A, kmax = 20, tol = .Machine$double.eps^(1/2))
rootm(A, p, kmax = 20, tol = .Machine$double.eps^(1/2))

Arguments

A

numeric, i.e. real, matrix.

p

p-th root to be taken.

kmax

maximum number of iterations.

tol

absolut tolerance, norm distance of A and B^p.

Value

• sqrtm(A) returns a list with components
• B
  square root matrix of A.
• Binv
  inverse of the square root matrix.
• k
  number of iterations.
• acc
  accuracy or absolute error.
• rootm(A) returns a list with components
• B
  square root matrix of A.
• k
  number of iterations.
• acc
  accuracy or absolute error.
• If k is negative the iteration has not converged.

  signm just returns one matrix, even when there was no convergence.

Details

A real matrix may or may not have a real square root; if it has no real negative eigenvalues. The number of square roots can vary from two to infinity. A positive definite matric has one distinguished square root, called the principal one, with the property that the eigenvalues lie in the segment {z | -pi/p < arg(z) < pi/p} (for the p-th root).

The matrix square root sqrtm(A) is computed here through the Denman-Beavers iteration (see the references) with quadratic rate of convergence, a refinement of the common Newton iteration determining roots of a quadratic equation.

The matrix p-th root rootm(A) is computed as a complex integral $$A^{1/p} = \frac{p \sin(\pi/p)}{\pi} A \int_0^{\infty} (x^p I + A)^{-1} dx$$ applying the trapezoidal rule along the unit circle.

One application is the computation of the matrix logarithm as $$\log A = 2^k log A^{1/2^k}$$ such that the argument to the logarithm is close to the identity matrix and the Pade approximation can be applied to $\log(I + X)$.

The matrix sector function is defined as sectm(A,m)=(A^m)^(-1/p)%*%A; for p=2 this is the matrix sign function.

S=signm(A) is real if A is and has the following properties: S^2=Id; S A = A S singm([0 A; B 0])=[0 C; C^-1 0] where C=A(BA)^-1

These functions arise in control theory.

References

N. J. Higham (1997). Stable Iterations for the Matrix Square Root. Numerical Algorithms, Vol. 15, pp. 227--242.

D. A. Bini, N. J. Higham, and B. Meini (2005). Algorithms for the matrix pth root. Numerical Algorithms, Vol. 39, pp. 349--378.

See Also

expm, expm::sqrtm

Examples

A1 <- matrix(c(10,  7,  8,  7,
                7,  5,  6,  5,
                8,  6, 10,  9,
                7,  5,  9, 10), nrow = 4, ncol = 4, byrow = TRUE)

X <- sqrtm(A1)$B    # accuracy: 2.352583e-13
X
A2 <- matrix(c(90.81, 8.33, 0.68, 0.06, 0.08, 0.02, 0.01, 0.01,
                0.70, 90.65, 7.79, 0.64, 0.06, 0.13, 0.02, 0.01,
                0.09, 2.27, 91.05, 5.52, 0.74, 0.26, 0.01, 0.06,
                0.02, 0.33, 5.95, 85.93, 5.30, 1.17, 1.12, 0.18,
                0.03, 0.14, 0.67, 7.73, 80.53, 8.84, 1.00, 1.06,
                0.01, 0.11, 0.24, 0.43, 6.48, 83.46, 4.07, 5.20,
                0.21, 0, 0.22, 1.30, 2.38, 11.24, 64.86, 19.79,
                0, 0, 0, 0, 0, 0, 0, 100
              ) / 100, nrow = 8, ncol = 8, byrow = TRUE)

X <- rootm(A2, 12)  # k = 6, accuracy: 2.208596e-14

##  Matrix sign function
signm(A1)                               # 4x4 identity matrix
B <- rbind(cbind(zeros(4,4), A1),
           cbind(eye(4), zeros(4,4)))
signm(B)                                # [0, signm(A1)$B; signm(A1)$Binv 0]