Description Usage Arguments Details Value Note References See Also Examples

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

1 2 3 4 |

`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 S`

`singm([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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | ```
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]
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

Please suggest features or report bugs in the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.