[R] geodesic distance (solution)

[stefano iacus]

 > Hi,
 > has anyone ever seen implemented in R the following "geodesic"
 > distance between positive definite pxp matrices A and B?
 >
 > d(A,B) = \sum_{i=1}^p (\log \lambda_i)^2
 >
 > were \lambda is the solution of det(A -\lambda B)  = 0
 >
 > thanks
 > stefano

as I received few private email on the claimed solution, I'm posting  
it to r-help.

when matrix B is invertible (which is always my case), one approach  
is to notice that
solving

det(A -\lambda * B)  = 0

is equivalent to solve

det(B^-1*A -\lambda *I)  = 0

which is a standard eigen value problem for the matrix B^-1 * A, hence

eigen(solve(B) %*% A)$values

is the answer.

I'm pretty sure that the problem can also be solved using some svd  
decomposition when B is not invertible.

hope it helps
stefano