Fwd: matrix polynomial

["Richard M. Heiberger" <rmh@temple.edu>]

matrix.poly <- function(X, cc, sym=TRUE) {
  X.e <- eigen(X)
  powers <- 0:dim(X)[1]
  
  ## I prefer the counting scheme where the power and coefficient number
  ## are the same.
  
  X.e.p <- outer(X.e$values, powers, "^")

  result.diag <- (X.e.p %*% cc)[,1]
  if(sym)
    t(X.e$vectors) %*% diag(result.diag) %*% X.e$vectors
  else
    solve(X.e$vectors) %*% diag(result.diag) %*% X.e$vectors
}

## trace(matrix.poly, exit=browser)

X <- matrix(rnorm(25),5,5)
cc <- sample(1:6)

matrix.poly(X, cc, sym=FALSE) ## Not symmetric

Xs <- X + t(X)
matrix.poly(Xs, cc, sym=TRUE) ## symmetric


## with the reversed numbering scheme for the coefficients
matrix.poly(Xs, rev(cc), sym=TRUE) ## symmetric


---- Original message ----
>Date: Sun, 16 Oct 2005 18:29:23 -0700
>From: "Seth O'Beeth" <seth3116@gmail.com>  
>Subject: [S] matrix polynomial  
>To: s-news@lists.biostat.wustl.edu
>
>   Dear S-users:
>    
>   Can anybody suggest an elegant vectorized way of creating a square matrix, 
> B,  which is
>   the polynomial of order n of the square matix X?
>   Namely, B=c0*X%*%....n times....%*%X + c1*X%*%...n-1 times...%*%X +.....+cn
>   The coefficient c1,...., cn and matrix X are considered known
>    
>    
>   Thank you in advance
>   Seth O'Beeth
>