[R] Coding matrix equation

[Petr Savicky]

On Mon, Apr 11, 2011 at 08:43:03AM +0100, matthew.r.robinson at sheffield.ac.uk wrote:
> Hi all,
> 
> I have two matrices:
> 
> G<-matrix(c(2.0, 0.5, 0.5, 0.5, 2.0, 0.5, 0.5, 0.5,2.0),3,3)
> P<-matrix(c(1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.5, 0.5,1.0),3,3)
> 
> and I want to run this equation to get a new matrix F:
> 
> F = [P+2G]^-1/2 P [P+2G]^-1/2
> 
> Could someone please tell me how to code this in R?

Hi.

Try the following.

  sqrtSymMat <- function(A)
  {
      out <- eigen(A)
      D <- diag(out$values)
      U <- out$vectors
      U %*% sqrt(D) %*% t(U)
  }

  A <- sqrtSymMat(solve(P + 2*G))

  F <- A %*% P %*% A

See also the function svd() and its help ?svd.

The operator A^(1/2) works component-wise. There may be a function
computing the square root of a positive semidefinite matrix in some
of the CRAN extension packages

  http://cran.at.r-project.org/web/packages/index.html

but i am not sure. The package mvtnorm

  http://cran.at.r-project.org/web/packages/mvtnorm/index.html

computes the square root of a matrix internally. See the help
of the function ?rmvnorm.

Hope this helps.

Petr Savicky.