# Re: [R] Coding matrix equation

From: Petr Savicky <savicky_at_praha1.ff.cuni.cz>
Date: Mon, 11 Apr 2011 14:17:09 +0200

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

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

but i am not sure. The package mvtnorm

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

Hope this helps.

Petr Savicky.

R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Mon 11 Apr 2011 - 12:21:29 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 11 Apr 2011 - 12:50:28 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.