Taking a diagonal matrix to a power

Question

I have a Vector X X <-rnorm(n). I've turned it into a diagonal matrix V <- matrix(diag(abs(X),ncol = n). I want to raise this to the power of -1/2. I've tried to use expm::expm with K <- V %^% (-1/2), But when I do this it just gets turned into a diagonal matrix of 1's. When it should be (v_i)^(-1/2). How would I fix this? Out put:

[1,] 0.08378436 0.0000000 0.000000 0.0000000 0.0000000
[2,] 0.00000000 0.9829437 0.000000 0.0000000 0.0000000
[3,] 0.00000000 0.0000000 1.875067 0.0000000 0.0000000
[4,] 0.00000000 0.0000000 0.000000 0.1861447 0.0000000
[5,] 0.00000000 0.0000000 0.000000 0.0000000 0.6334857
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    0    1    0    0
[4,]    0    0    0    1    0
[5,]    0    0    0    0    1

`help("%^%")` says that the exponent needs to be a non-negative integer. `X %^% -1` is an error, so the "problem" is likely just that the exponent is truncated towards zero. — Mikael Jagan, Apr 08 '23 at 22:05
more specifically, the source code for the function is `.Call(R_matpow, x, as.integer(k))`; `as.integer()` will do the truncation. It does seem too bad that no warning is given ... — Ben Bolker, Apr 08 '23 at 23:48

score 4 · Answer 1 · edited Apr 09 '23 at 06:17

The other answers here are good, and general, but for your particular use case (taking the inverse square root of a diagonal matrix), you can take advantage of the fact that if D is a diagonal matrix, then computing D^m is equivalent to raising the diagonal elements of D to the mth power.

diag(V) will extract the diagonal (although in your specific example you already know this is just abs(X))
diag(V)^(-1/2) will take the inverse square root elementwise
diag(diag(V)^(-1/2)) will re-embed the results as the diagonal of a matrix

D <- diag(diag(V)^(-1/2))
all.equal(D, expm::sqrtm(solve(V)))  ## TRUE
all.equal(D, chol(solve(V)))   ## TRUE

Taking advantage of the special case is also much faster than the general solutions:

microbenchmark(diag(diag(V)^(-1/2)), expm::sqrtm(solve(V)), chol(solve(V)))
Unit: microseconds
                 expr      min       lq       mean    median       uq      max
 diag(diag(V)^(-1/2))    5.701    7.093    9.61355    9.7385   11.251   39.554
      sqrtm(solve(V)) 1412.884 1503.633 1559.03027 1518.5315 1544.601 3366.459
       chol(solve(V))   21.150   24.070   26.30623   26.3740   27.702   45.565
 neval cld
   100  a 
   100   b
   100  a

Thanks - I didn't notice the word "diagonal". I would maybe still suggest (or mention) `diag(as.complex(diag(V))^-0.5)` to allow negative diagonal entries. — Mikael Jagan, Apr 09 '23 at 00:03

score 3 · Answer 2 · edited Apr 09 '23 at 06:14

I think you misused %^% for your purpose. When you type ?`%^%` in the console, you will see

Compute the k-th power of a matrix. Whereas x^k computes element wise powers, x %^% k corresponds to k - 1 matrix multiplications, x %% x %% ... %*% x.

and

If you think you need x^k for k<0, then consider instead solve(x %^% (-k)).

A possible workaround: You can try below for the equivalence to V^(-1/2)

m <- chol(solve(V))

and verify it by

crossprod(V, crossprod(m))

Mikael Jagan · Answer 3 · 2023-04-08T22:54:51.393

You can use the sqrtm function, also from expm. It gives the "principal" square root, where the square root (which need not be unique) exists and can be computed. It does not require its argument to be positive semidefinite or even symmetric and so may return a matrix of type complex.

library(expm)
set.seed(0)
n <- 6L
x <- matrix(rnorm(n * n), n, n)
inv.sqrt.x <- sqrtm(solve(x))
stopifnot(all.equal(x %*% inv.sqrt.x %*% inv.sqrt.x, diag(1+0i, n, n)))

The Wikipedia entry on the matrix square root is quite thorough and probably worth a read. This answer directly addresses some of the "bad" cases.

Taking a diagonal matrix to a power

3 Answers3