############################################################################ ### Functions for the paper "Efficient recursive algorithms for functionals ### based on higher order derivatives of the multivariate Gaussian density" ### by J.E. Chacón and T. Duong ### Last modified: March 12, 2014 ############################################################################ ############################################################################ ### First the necessary libraries are loaded ### Some further auxiliary functions are put at the bottom of this script ############################################################################ library(combinat) library(multicool) library(ks) ############################################################################ ### Symmetrizer matrix computation ############################################################################ ############################################################################ ### Sdr.direct computes the symmetrizer matrix S_{d,r} based on Equation (4) ### as described in Section 3 of Chacón and Duong (2014) ############################################################################ Sdr.direct<-function(d,r){ S<-matrix(0,nc=d^r,nr=d^r) per<-permn(r) per.rep<-pinv.all(d,r) nper<-factorial(r) nper.rep<-d^r per<-matrix(unlist(per), byrow=TRUE, ncol=r, nrow=nper) pow<-0:(r-1) dpow<-d^pow if(nper.rep<=nper){ dpow.mat<-matrix(rep(dpow,nper),byrow=TRUE,ncol=r,nrow=nper) for(i in 1:nper.rep){ ## Loop over no. perms with reps (d^r) pinvi<-per.rep[i,] sigpinvi<-matrix(pinvi[per],byrow=FALSE,nrow=nrow(per),ncol=ncol(per)) psigpinvi<-drop(1+rowSums((sigpinvi-1)*dpow.mat)) S[i,]<-tabulate(psigpinvi,nbins=d^r) }} if(nper=2){ Id<-diag(d) T<-Id A<-Kmat(d,d) for(j in 2:r){ T<-((j-1)/j)*(A%*%(T%x%Id)%*%A)+A/j S<-(S%x%Id)%*%T if(j=2){for(i in 2:r){ hmnew<-mat.Kprod(hmold1,arg)-(i-1)*mat.Kprod(nvGinv,hmold0) hmnew<-matrix(Sdrv.recursive(d=d,r=i,v=hmnew), nrow=n) hmold0<-hmold1 hmold1<-hmnew } } dens<-dmvnorm(x,mean=rep(0,d),sigma=Sigma) result<-matrix(rep(dens,d^r),byrow=FALSE,nrow=n,ncol=d^r)*hmnew*(-1)^r return(drop(result)) } ############################################################################### ### dmvnorm.deriv.unique computes the whole vector derivative of the Gaussian ### density phi_Sigma(x) from its unique coordinates, based on Algorithm 3 as ### described in Section 5 of Chacón and Duong (2014) ############################################################################### dmvnorm.deriv.unique<-function(x,Sigma,deriv.order=0){ if(is.vector(x)){x<-matrix(x,nrow=1)} d<-ncol(x) n<-nrow(x) r<-deriv.order G<-Sigma Ginv<-chol2inv(chol(G)) arg<-x%*%Ginv hmold0 <- matrix(1,nrow=n,ncol=1) hmold1 <- arg hmnew <- hmold0 udind0<-matrix(rep(0,d),nrow=1,ncol=d) udind1<-diag(d) if(r==1){hmnew<-hmold1} if(r>=2){for(i in 2:r){ Ndi1<-ncol(hmold1) Ndi0<-ncol(hmold0) hmnew<-numeric() for(j in 1:d){ nrecj<-choose(d-j+i-1,i-1) hmnew.aux<-arg[,j]*hmold1[,Ndi1-(nrecj:1)+1] for(k in j:d){ udind0.aux<-matrix(udind1[Ndi1-(nrecj:1)+1,],ncol=d,byrow=FALSE) udind0.aux[,k]<-udind0.aux[,k]-1 valid.udind0<-as.logical(apply(udind0.aux>=0,1,min)) enlarged.hmold0<-matrix(0,ncol=nrow(udind0.aux),nrow=n) for(l in 1:nrow(udind0.aux)){ if(valid.udind0[l]){ pos<-which(rowSums((udind0-matrix(rep(udind0.aux[l,],nrow(udind0)), nrow=nrow(udind0),byrow=TRUE))^2)==0) enlarged.hmold0[,l]<-hmold0[,pos]} } hmnew.aux<-hmnew.aux-Ginv[j,k]*matrix(rep(udind1[Ndi1-(nrecj:1)+1,k],n), nrow=n,byrow=TRUE)*enlarged.hmold0 } hmnew<-cbind(hmnew,hmnew.aux) } hmold0 <- hmold1 hmold1 <- hmnew ##Compute the unique i-th derivative multi-indexes nudind1<-nrow(udind1) udindnew<-numeric() for(j in 1:d){ Ndj1i<-choose(d+i-1-j,i-1) udind.aux<-matrix(udind1[nudind1-(Ndj1i:1)+1,],ncol=d,byrow=FALSE) udind.aux[,j]<-udind.aux[,j]+1 udindnew<-rbind(udindnew,udind.aux) } udind0<-udind1 udind1<-udindnew }} if(r==0) result<-dmvnorm(x,mean=rep(0,d),sigma=Sigma) if(r==1) result<-(-1)*matrix(rep(dmvnorm(x,mean=rep(0,d),sigma=Sigma),d), nrow=n,byrow=FALSE)*hmnew if(r>=2){ per<-pinv.all(d=d,r=r) dind<-numeric() udind<-udind1 dind.base<-rep(0,d^r) udind.base<-rep(0,choose(d+r-1,r)) for(i in 1:d){ dind<-cbind(dind,rowSums(per==i)) ## Matrix of derivative indices dind.base<-dind.base+dind[,i]*(r+1)^(d-i) ## Transform each row to base r+1 udind.base<-udind.base+udind[,i]*(r+1)^(d-i) ## Transform each row to base r+1 } dlabs<-match(dind.base,udind.base) deriv.vector<-hmnew[,dlabs]*matrix(rep((-1)^rowSums(dind),n),nrow=n,byrow=TRUE) result<-matrix(rep(dmvnorm(x,mean=rep(0,d),sigma=Sigma),ncol(deriv.vector)), nrow=n,byrow=FALSE)*deriv.vector } return(drop(result)) } ########################################################################## ### Vector moments of the normal distribution ########################################################################## ############################################################################# ### mur.direct computes the vector moment E[X^{\otimes r}] for a random ### vector with N(mu,Sigma) distribution, on the basis of Equation (8) in ### Section 6 of Chacón and Duong (2014) ############################################################################# mur.direct<-function(r,mu,Sigma){ d<-ncol(Sigma) result<-as.vector(Kpow(mu,r)) vS<-vec(Sigma) vSj<-1 if(r>=2){ for(j in 1:floor(r/2)){ vSj<-as.vector(vSj%x%vS) cj<-prod(r:(r-2*j+1))/(prod(1:j)*2^j) mur2j<-as.vector(Kpow(mu,r-2*j)) result<-result+cj*as.vector(mur2j%x%vSj) } } return(drop(Sdrv.recursive(d=d,r=r,v=result))) } ############################################################################# ### mur.recursive computes the vector moment E[X^{\otimes r}] for a random ### vector with N(mu,Sigma) distribution, on the basis of Equation (9) in ### Section 6 of Chacón and Duong (2014), using Equation (7) in Section 5 ### to obtain the Hermite polynomial ############################################################################# mur.recursive<-function(r,mu,Sigma){ d<-ncol(Sigma) G<- -Sigma vG<-vec(G) arg<-mu hmold0 <- 1 hmold1 <- arg hmnew <- hmold0 if(r==1){hmnew<-hmold1} if(r>=2){for(i in 2:r){ hmnew<-as.vector(arg%x%hmold1-(i-1)*(vG%x%hmold0)) hmold0<-hmold1 hmold1<-hmnew } } return(drop(Sdrv.recursive(d=d,r=r,v=hmnew))) } ############################################################################### ### mur.unique computes the vector moment E[X^{\otimes r}] for a random vector ### with N(mu,Sigma) distribution, on the basis of Equation (9) in Section 6 ### of Chacón and Duong (2014), using Algorithm 3 in Section 5, based on the ### unique partial derivatives, to obtain the Hermite polynomial ############################################################################### mur.unique<-function(r,mu,Sigma){ d<-ncol(Sigma) G<- -Sigma arg<-mu hmold0 <- 1 hmold1 <- arg hmnew <- hmold0 udind0<-matrix(rep(0,d),nrow=1,ncol=d) udind1<-diag(d) if(r==1){hmnew<-hmold1} if(r>=2){for(i in 2:r){ Ndi1<-length(hmold1) Ndi0<-length(hmold0) hmnew<-numeric() for(j in 1:d){ nrecj<-choose(d-j+i-1,i-1) hmnew.aux<-arg[j]*hmold1[Ndi1-(nrecj:1)+1] for(k in j:d){ udind0.aux<-matrix(udind1[Ndi1-(nrecj:1)+1,],ncol=d,byrow=FALSE) udind0.aux[,k]<-udind0.aux[,k]-1 valid.udind0<-as.logical(apply(udind0.aux>=0,1,min)) enlarged.hmold0<-rep(0,nrow(udind0.aux)) for(l in 1:nrow(udind0.aux)){ if(valid.udind0[l]){ pos<-which(rowSums((udind0-matrix(rep(udind0.aux[l,],nrow(udind0)), nrow=nrow(udind0),byrow=TRUE))^2)==0) enlarged.hmold0[l]<-hmold0[pos]} } hmnew.aux<-hmnew.aux-G[j,k]*udind1[Ndi1-(nrecj:1)+1,k]*enlarged.hmold0 } hmnew<-c(hmnew,hmnew.aux) } hmold0 <- hmold1 hmold1 <- hmnew ## Compute the unique i-th derivative multi-indexes nudind1<-nrow(udind1) udindnew<-numeric() for(j in 1:d){ Ndj1i<-choose(d+i-1-j,i-1) udind.aux<-matrix(udind1[nudind1-(Ndj1i:1)+1,],ncol=d,byrow=FALSE) udind.aux[,j]<-udind.aux[,j]+1 udindnew<-rbind(udindnew,udind.aux) } udind0<-udind1 udind1<-udindnew }} if(r==0){ result<-1 } if(r==1){ result<-(-1)*hmnew } if(r>=2){ per<-pinv.all(d=d,r=r) dind<-numeric() udind<-udind1 dind.base<-rep(0,d^r) udind.base<-rep(0,choose(d+r-1,r)) for(i in 1:d){ dind<-cbind(dind,rowSums(per==i)) ## Matrix of derivative indices dind.base<-dind.base+dind[,i]*(r+1)^(d-i) ## Transform each row to base r+1 udind.base<-udind.base+udind[,i]*(r+1)^(d-i) ## Transform each row to base r+1 } dlabs<-match(dind.base,udind.base) result<-hmnew[dlabs]*(-1)^rowSums(dind) } return(drop(result*(-1)^r)) } ########################################################################## ### Moments of quadratic forms in normal variables ########################################################################## ############################################################################# ### nur.direct computes the moment E[(X^T AX)^r] of the quadratic form ### X^T AX where X is a random vector with N(mu,Sigma) distribution, using ### Equation (10) in Section 6 of Chacón and Duong (2014), and the direct ### implementation mur.direct of the normal moments ############################################################################# nur.direct<-function(r,A,mu,Sigma){ vA<-vec(A) result<-drop(Kpow(t(vA),r)%*%mur.direct(2*r,mu,Sigma)) return(result) } ############################################################################# ### nur.recursive computes the moment E[(X^T AX)^r] of the quadratic form ### X^T AX where X is a random vector with N(mu,Sigma) distribution, using ### Equation (10) in Section 6 of Chacón and Duong (2014), and the recursive ### implementation mur.recursive of the normal moments ############################################################################# nur.recursive<-function(r,A,mu,Sigma){ vA<-vec(A) result<-drop(Kpow(t(vA),r)%*%mur.recursive(2*r,mu,Sigma)) return(result) } ############################################################################# ### nur.unique computes the moment E[(X^T AX)^r] of the quadratic form ### X^T AX where X is a random vector with N(mu,Sigma) distribution, using ### Equation (10) in Section 6 of Chacón and Duong (2014), and the function ### mur.unique to compute the normal moments from its unique coordinates ############################################################################# nur.unique<-function(r,A,mu,Sigma){ vA<-vec(A) result<-sum(Kpow(vA,r)*mur.unique(2*r,mu,Sigma)) return(result) } ############################################################################# ### nur.cumulant computes the moment E[(X^T AX)^r] of the quadratic form ### X^T AX where X is a random vector with N(mu,Sigma) distribution, using ### the recursive formula relating moments and cumulants ############################################################################# nur.cumulant<-function(r,A,mu,Sigma){ if(r==0){result<-1} if(r==1){result<-sum(diag(A%*%Sigma))+t(mu)%*%A%*%mu} if(r>=2){ ASigma<-A%*%Sigma AS<-ASigma Amu<-A%*%mu kappas<-sum(diag(ASigma)+mu*Amu) nus<-kappas for(k in 2:r){ knew<-k*t(mu)%*%ASigma%*%Amu ASigma<-ASigma%*%AS knew<-(knew+sum(diag(ASigma)))*factorial(k-1)*2^(k-1) nnew<-knew+sum(choose(k-1,1:(k-1))*nus*rev(kappas)) kappas<-c(kappas,knew) nus<-c(nus,nnew) } result<-nnew } return(drop(result)) } ############################################################################# ### nurs.direct computes the joint moment E[(X^T AX)^r (X^T BX)^s] of the ### quadratic forms X^T AX and X^T BX, where X is a random vector with ### N(mu,Sigma) distribution, using Equation (10) in Section 6 of Chacón and ### Duong (2014), and the direct implementation mur.direct of the ### normal moments ############################################################################# nurs.direct<-function(r,s,A,B,mu,Sigma){ vA<-vec(A) vB<-vec(B) result<-(Kpow(t(vA),r)%x%Kpow(t(vB),s))%*%mur.direct(2*r+2*s,mu,Sigma) return(drop(result)) } ############################################################################# ### nurs.recursive computes the joint moment E[(X^T AX)^r (X^T BX)^s] of the ### quadratic forms X^T AX and X^T BX, where X is a random vector with ### N(mu,Sigma) distribution, using Equation (10) in Section 6 of Chacón and ### Duong (2014), and the recursive implementation mur.recursive of the ### normal moments ############################################################################# nurs.recursive<-function(r,s,A,B,mu,Sigma){ vA<-vec(A) vB<-vec(B) result<-(Kpow(t(vA),r)%x%Kpow(t(vB),s))%*%mur.recursive(2*r+2*s,mu,Sigma) return(drop(result)) } ############################################################################# ### nurs.unique computes the joint moment E[(X^T AX)^r (X^T BX)^s] of the ### quadratic forms X^T AX and X^T BX, where X is a random vector with ### N(mu,Sigma) distribution, using Equation (10) in Section 6 of Chacón and ### Duong (2014), and the function mur.unique to compute the normal moments ### from its unique coordinates ############################################################################# nurs.unique<-function(r,s,A,B,mu,Sigma){ vA<-vec(A) vB<-vec(B) result<-drop((Kpow(t(vA),r)%x%Kpow(t(vB),s))%*%mur.unique(2*r+2*s,mu,Sigma)) return(drop(result)) } ############################################################################# ### nurs.cumulant computes the joint moment E[(X^T AX)^r (X^T BX)^s] of the ### quadratic forms X^T AX and X^T BX, where X is a random vector with ### N(mu,Sigma) distribution, using the recursive formula (11) in Section 6 ### of Chacón and Duong (2014), relating moments and cumulants. The cumulants ### are computed using the function kappars, which is based on Theorem 3 ############################################################################# kappars<-function(r,s,A,B,mu,Sigma){ d<-ncol(A) if(r+s>1 & r>0 & s>0){ind<-allPerm(initMC(c(rep(1,r),rep(2,s))))} if(r+s==1 | r==0 | s==0){ind<-matrix(c(rep(1,r),rep(2,s)),nrow=1)} if(r+s==0){return(0)} nper<-nrow(ind) result<-0 Dmat<-solve(Sigma)%*%mu%*%t(mu) ASigma<-A%*%Sigma BSigma<-B%*%Sigma Id<-diag(d) for(i in 1:nper){ product<-Id for(j in 1:(r+s)){ if(ind[i,j]==1){product<-product%*%ASigma} else if(ind[i,j]==2){product<-product%*%BSigma} } result<-result+sum(diag(product%*%(Id/(r+s)+Dmat))) } result<-result*factorial(r)*factorial(s)*2^(r+s-1) return(drop(result)) } nurs.cumulant<-function(r,s,A,B,mu,Sigma){ if(r==0 & s>0){nurs<-nur.cumulant(s,B,mu,Sigma)} if(r>0 & s==0){nurs<-nur.cumulant(r,A,mu,Sigma)} if(r==0 & s==0){nurs<-1} if((r>0)&(s>0)){ K<-matrix(0,nrow=r+1,ncol=s) for(i in 0:r){for(j in 1:s){ K[i+1,j]<-kappars(i,j,A,B,mu,Sigma) }} N<-matrix(0,nrow=r+1,ncol=s) for(i in 0:r){ N[i+1,1]<-nur.cumulant(r=i,A=A,mu=mu,Sigma=Sigma) } if(s>1){ for(j in 1:(s-1)){for(i in 0:r){ Choose<-outer(choose(i,0:i),choose(j-1,0:(j-1))) N[i+1,j+1]<-sum(Choose*N[1:(i+1),1:j]*K[(i:0)+1,j:1]) }} } Choose<-outer(choose(r,0:r),choose(s-1,0:(s-1))) nurs<-sum(Choose*N[1:(r+1),1:s]*K[(r:0)+1,s:1]) } return(nurs) } ########################################################################## ### V-statistics with multivariate Gaussian derivatives kernel ########################################################################## ############################################################################ ### Qr.direct computes the V-statistic using the direct approach, as ### described in Section 6 of Chacón and Duong (2014) ############################################################################ Qr.direct <- function(x, Sigma, r=0, inc=1) { if (is.vector(x)) x <- matrix(x, nrow=1) eta <- nrow(x)^(-2)*drop(Kpow(vec(diag(d)), r/2) %*% ks:::dmvnorm.deriv.sum(x=x, Sigma=Sigma, deriv.order=r)) return(eta) } ############################################################################ ### Qr.cumulant computes the V-statistic using the relationship with nur ### shown in Theorem 4 of Chacón and Duong (2014) ############################################################################ Qr.cumulant <- function(x, Sigma, r=0, inc=1) { if (is.vector(x)) x <- matrix(x, nrow=1) d <- ncol(x) r <- r/2 y <- x nx <- as.numeric(nrow(x)) ny <- as.numeric(nrow(y)) G <- Sigma Ginv <- chol2inv(chol(G)) G2inv <- Ginv%*%Ginv G3inv <- G2inv%*%Ginv trGinv <- sum(diag(Ginv)) trG2inv <- sum(diag(G2inv)) detG <- det(G) ## indices for separating into blocks for double sum calculation n.seq <- ks:::block.indices(nx, ny, d=d, r=r, diff=FALSE) if (r==0) { xG <- x%*%Ginv a <- rowSums(xG*x) eta <- 0 for (i in 1:(length(n.seq)-1)) { nytemp <- n.seq[i+1] - n.seq[i] ytemp <- matrix(y[n.seq[i]:(n.seq[i+1]-1),], ncol=d) aytemp <- rowSums((ytemp %*% Ginv) *ytemp) M <- a%*%t(rep(1,nytemp)) + rep(1, nx)%*%t(aytemp) - 2*(xG%*%t(ytemp)) em2 <- exp(-M/2) eta <- eta + (2*pi)^(-d/2)*detG^(-1/2)*sum(em2) } } else if (r==1) { xG <- x%*%Ginv xG2 <- x%*%G2inv a <- rowSums(xG*x) a2 <- rowSums(xG2*x) eta <- 0 for (i in 1:(length(n.seq)-1)) { nytemp <- n.seq[i+1] - n.seq[i] ytemp <- matrix(y[n.seq[i]:(n.seq[i+1]-1),], nrow=nytemp) aytemp <- rowSums((ytemp %*% Ginv) *ytemp) aytemp2 <- rowSums((ytemp %*% G2inv) *ytemp) M <- a%*%t(rep(1,nytemp))+rep(1,nx)%*%t(aytemp)-2*(xG%*%t(ytemp)) M2 <- a2%*%t(rep(1,nytemp))+rep(1,nx)%*%t(aytemp2)-2*(xG2%*%t(ytemp)) eta <- eta + (2*pi)^(-d/2)*detG^(-1/2)*sum(exp(-M/2)*(M2-trGinv)) } } else if (r==2) { xG <- x%*%Ginv xG2 <- x%*%G2inv xG3 <- x%*%G3inv a <- rowSums(xG*x) a2 <- rowSums(xG2*x) a3 <- rowSums(xG3*x) eta <- 0 for (i in 1:(length(n.seq)-1)) { nytemp <- n.seq[i+1] - n.seq[i] ytemp <- matrix(y[n.seq[i]:(n.seq[i+1]-1),], ncol=d) aytemp <- rowSums((ytemp %*% Ginv) *ytemp) aytemp2 <- rowSums((ytemp %*% G2inv) *ytemp) aytemp3 <- rowSums((ytemp %*% G3inv) *ytemp) M <- a%*%t(rep(1,nytemp))+rep(1,nx)%*%t(aytemp)-2*(xG%*%t(ytemp)) M2 <- a2%*%t(rep(1,nytemp))+rep(1,nx)%*%t(aytemp2)-2*(xG2%*%t(ytemp)) M3 <- a3%*%t(rep(1,nytemp))+rep(1,nx)%*%t(aytemp3)-2*(xG3%*%t(ytemp)) eta <- eta + (2*pi)^(-d/2)*detG^(-1/2)*sum(exp(-M/2)*(2*trG2inv-4*M3 +(-trGinv+M2)^2)) } } else if (r>2) { xG <- x%*%Ginv a <- rowSums(xG*x) eta <- 0 for (i in 1:(length(n.seq)-1)) { nytemp <- n.seq[i+1] - n.seq[i] ytemp <- matrix(y[n.seq[i]:(n.seq[i+1]-1),], ncol=d) aytemp <- rowSums((ytemp %*% Ginv) *ytemp) M <- a %*% t(rep(1,nytemp)) + rep(1,nx)%*%t(aytemp) - 2*(xG%*%t(ytemp)) edv2 <- exp(-M/2) P0<-Ginv kappas <- matrix(nrow=as.numeric(nx*nytemp), ncol=r) for (j in 1:r) { Gi1inv <- P0%*%Ginv trGi0inv <- sum(diag(P0)) xGi1inv <- x%*%Gi1inv xGi1invx <- rowSums(xGi1inv*x) aytemp <- rowSums((ytemp %*% Gi1inv) *ytemp) dvi1 <- xGi1invx%*%t(rep(1,nytemp))+rep(1,nx)%*%t(aytemp)-2*(xGi1inv%*%t(ytemp)) kappas[,j] <- (-2)^(j-1)*factorial(j-1)*(-trGi0inv+j*dvi1) P0 <- Gi1inv } nus <- matrix(nrow=as.numeric(nx*nytemp), ncol=r+1) nus[,1] <- 1 for (j in 1:r) { js<-0:(j-1) if (j==1) nus[,2] <- kappas[,1] else nus[,j+1] <- rowSums(kappas[,j:1]*nus[,1:j]/matrix(rep(factorial(js)* factorial(rev(js)),nx*nytemp),nrow=nx*nytemp,byrow=TRUE))*factorial(j-1) } eta <- eta + (2*pi)^(-d/2)*detG^(-1/2)*sum(edv2*nus[,r+1]) } } if (inc==0) eta <- (eta - (-1)^r*nx*nu.cumulant(r=r, A=Ginv)*(2*pi)^(-d/2)* detG^(-1/2))/(nx*(ny-1)) if (inc==1) eta <- eta/(nx*ny) return(eta) } ########################################################################## ### Auxiliary functions ########################################################################## ########################################################################## ### pinv.all generates all the permutations PR_{d,r} as described in ### Appendix B of Chacón and Duong (2014) ########################################################################## pinv.all<-function(d,r){ i<-1:d^r n<-i-1 dpow<-d^(0:r) n.mat<-matrix(rep(n,r+1),byrow=FALSE,nrow=d^r,ncol=r+1) dpow.mat<-matrix(rep(dpow,d^r),byrow=TRUE,nrow=d^r,ncol=r+1) ndf.mat<-floor(n.mat/dpow.mat) ans<-ndf.mat[,r:1]-d*ndf.mat[,(r+1):2] return(ans+1) } ########################################################################## ### Kmat computes the commutation matrix of orders m,n ########################################################################## Kmat<-function(m,n){ K<-matrix(0,ncol=m*n,nrow=m*n) i<-1:m;j<-1:n rows<-rowSums(expand.grid((i-1)*n,j)) cols<-rowSums(expand.grid(i,(j-1)*m)) positions<-cbind(rows,cols) K[positions]<-1 return(K) } ########################################################################## ### mat.Kprod computes row-wise Kronecker products of matrices ########################################################################## mat.Kprod<-function(U,V){ #### Returns a matrix with rows U[i,]%x%V[i,] n1<-nrow(U) n2<-nrow(V) if(n1!=n2)stop("U and V must have the same number of vectors") p<-ncol(U) q<-ncol(V) onep<-rep(1,p) oneq<-rep(1,q) P<-(U%x%t(oneq))*(t(onep)%x%V) return(P) } ########################################################################## ##### Kpow computes the Kronecker power of a matrix A ########################################################################## Kpow<-function(A,pow){ if(floor(pow)!=pow)stop("pow must be an integer") Apow<-A if(pow==0){Apow<-1} if(pow>1){ for(i in 2:pow) Apow<-Apow%x%A } return(Apow) } ########################################################################## ##### matrix.pow computes the nth power of a matrix A ########################################################################## matrix.pow <- function(A, n) { if (nrow(A)!=ncol(A)) stop("A must a a square matrix") if (floor(n)!=n) stop("n must be an integer") if (n==0) return(diag(ncol(A))) if (n < 0) return(matrix.pow(A=solve(A), n=-n)) # trap non-integer n and return an error if (n == 1) return(A) result <- diag(1, ncol(A)) while (n > 0) { if (n %% 2 != 0) { result <- result %*% A n <- n - 1 } A <- A %*% A n <- n / 2 } return(result) }