gallon li
2007-Jan-24 07:04 UTC
[R] Matrix question: obtaining the square root of a positive definite matrix?
I want to compute B=A^{1/2} such that B*B=A. For example a=matrix(c(1,.2,.2,.2,1,.2,.2,.2,1),ncol=3) so> a[,1] [,2] [,3] [1,] 1.0 0.2 0.2 [2,] 0.2 1.0 0.2 [3,] 0.2 0.2 1.0> a%*%a[,1] [,2] [,3] [1,] 1.08 0.44 0.44 [2,] 0.44 1.08 0.44 [3,] 0.44 0.44 1.08> b=a%*%ai have tried to use singular value decomposion> c=svd(b)> c$u%*%diag(sqrt(c$d))[,1] [,2] [,3] [1,] -0.8082904 2.043868e-18 0.6531973 [2,] -0.8082904 -5.656854e-01 -0.3265986 [3,] -0.8082904 5.656854e-01 -0.3265986 this does not come close to the original a. Can anybody on this forum enlight me on how to get a which is the square root of b? [[alternative HTML version deleted]]
Prof Brian Ripley
2007-Jan-24 08:02 UTC
[R] Matrix question: obtaining the square root of a positive definite matrix?
On Wed, 24 Jan 2007, gallon li wrote:> I want to compute B=A^{1/2} such that B*B=A.According to your subject line A is positive definite and hence symmetric? The usual definition of a matrix square root involves a transpose, e.g. B'B = A. There are many square roots: were you looking for a symmetric one? For such an A,> e <- eigen(A) > V <- e$vectors > V %*% diag(e$values) %*% t(V)recovers A (up to rounding errors), and> B <- V %*% diag(sqrt(e$values)) %*% t(V)is such that B %*% B = A. Even that is not unique, e.g. -B is an equally good answer.> For example(with A = b and B = a, it seems)> a=matrix(c(1,.2,.2,.2,1,.2,.2,.2,1),ncol=3) > > so >> a > [,1] [,2] [,3] > [1,] 1.0 0.2 0.2 > [2,] 0.2 1.0 0.2 > [3,] 0.2 0.2 1.0 >> a%*%a > [,1] [,2] [,3] > [1,] 1.08 0.44 0.44 > [2,] 0.44 1.08 0.44 > [3,] 0.44 0.44 1.08 >> b=a%*%a > > i have tried to use singular value decomposion > >> c=svd(b) > >> c$u%*%diag(sqrt(c$d)) > [,1] [,2] [,3] > [1,] -0.8082904 2.043868e-18 0.6531973 > [2,] -0.8082904 -5.656854e-01 -0.3265986 > [3,] -0.8082904 5.656854e-01 -0.3265986 > > this does not come close to the original a. Can anybody on this forum > enlight me on how to get a which is the square root of b? > > [[alternative HTML version deleted]] > > ______________________________________________ > R-help at stat.math.ethz.ch 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. >-- Brian D. Ripley, ripley at stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UK Fax: +44 1865 272595