Matrix Square and p-th Roots
Computes the matrix square root and matrix p-th root of a nonsingular real matrix.
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))
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 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(π/p)}{π} A \int_0^{∞} (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 Ssingm([0 A; B 0])=[0 C; C^-1 0] where C=A(BA)^-1
These functions arise in control theory.
sqrtm(A) returns a list with components
B |
square root matrix of |
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 |
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.
The p-th root of a positive definite matrix can also be computed from its eigenvalues as
E <- eigen(A)V <- E\$vectors; U <- solve(V)D <- diag(E\$values)B <- V %*% D^(1/p) %*% U
or by applying the functions expm, logm in package ‘expm’:
B <- expm(1/p * logm(A))
As these approaches all calculate the principal branch, the results are identical (but will numerically slightly differ).
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.
expm, expm::sqrtm
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]Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.