[R] Re: Help : generating correlation matrix with a particular

Martin Maechler maechler at stat.math.ethz.ch
Mon Dec 13 17:03:36 CET 2004

>>>>> "Herbert" == Herbert Desson <Herbert_Desson at jltgroup.com>
>>>>>     on Mon, 13 Dec 2004 12:24:10 -0000 writes:

    Herbert> Here is some code we have used.

    Herbert> a<-array(c(1,.9,.7,.9,1,.3,.7,.3,1),dim=c(3,3))
    Herbert> a
    Herbert> s<-eigen(a)$vectors
    Herbert> l<-diag(eigen(a)$values)
    Herbert> l[l<0]<-0
    Herbert> b<-s%*%sqrt(l)
    Herbert> for(i in 1:nrow(b)){b[i,]<-b[i,]/sqrt(sum(b[i,]^2))}
    Herbert> ap<-b%*%t(b)
    Herbert> ap

This code does the same thing as my (simplistic, but slightly more
general) function posdefify()  in package "sfsmisc" :

  a <- matrix(c(1,.9,.7,.9,1,.3,.7,.3,1), 3)



          [,1]      [,2]      [,3]
[1,] 1.0000000 0.8940242 0.6963190
[2,] 0.8940242 1.0000000 0.3009691
[3,] 0.6963190 0.3009691 1.0000000

    Herbert> It is based on a paper by Rebonato etal that formerly was at
    Herbert> www.rebonato.com/correlationmatrix.pdf.
    Herbert> Unfortunately the website has disappeared.

The idea is very simple and has been re-invented many times as
far as I know.

More sophisticated methods for "posdefiying" a matrix exist in
other places. Given symmetrix matrix A, they try to find the
matrix Ap, positive definite, such  ||A - Ap||  is minimal.
The eigen-value based simple solution that you've used above
and I've also coded in posdefify(), is not the same one would
get for `usual' matrix norms  || . || 

[[NB:  posdefify() also has a 2nd method the implementation of
       which has an embarassing bug. The next version of
       sfsmisc, due in a day or two, will have it fixed.

Does anyone know of rigorous mathematical results in this

Martin Maechler, ETH Zurich

More information about the R-help mailing list