| Title: | Covariance Estimation for Matrix Data with the Kronecker-Core Decomposition |
|---|---|
| Description: | Matrix-variate covariance estimation via the Kronecker-core decomposition. Computes the Kronecker and core covariance matrices corresponding to an arbitrary covariance matrix, and provides an empirical Bayes covariance estimator that adaptively shrinks towards the space of separable covariance matrices. For details, see Hoff, McCormack and Zhang (2022) <arXiv:2207.12484> "Core Shrinkage Covariance Estimation for Matrix-variate data". |
| Authors: | Peter Hoff [aut, cre] |
| Maintainer: | Peter Hoff <[email protected]> |
| License: | GPL-3 |
| Version: | 0.1 |
| Built: | 2025-01-19 02:43:34 UTC |
| Source: | https://github.com/pdhoff/covkcd |
Reshape a covariance array to a covariance matrix.
ca2cm(A)ca2cm(A)
A |
a covariance array of dimension p1*p2*p1*p2. |
a p1*p2 by p1*p2 covariance matrix.
Peter Hoff
p1<-4 ; p2<-7 ; p<-p1*p2 S<-rWishart(1,p,diag(p))[,,1] A<-cm2ca(S,p1,p2) range(S-ca2cm(A))p1<-4 ; p2<-7 ; p<-p1*p2 S<-rWishart(1,p,diag(p))[,,1] A<-cm2ca(S,p1,p2) range(S-ca2cm(A))
Reshape a covariance matrix to a covariance array.
cm2ca(S, p1, p2)cm2ca(S, p1, p2)
S |
a covariance matrix of dimension (p1p2)*(p1p2). |
p1 |
the row dimension. |
p2 |
the column dimension. |
a four-way array where entry i1,j1,i2,j2 gives the covariance between element i1,j1 and element i2,j2 of a random matrix.
Peter Hoff
p1<-4 ; p2<-7 ; p<-p1*p2 S<-rWishart(1,p,diag(p))[,,1] A<-cm2ca(S,p1,p2) range(S-ca2cm(A))p1<-4 ; p2<-7 ; p<-p1*p2 S<-rWishart(1,p,diag(p))[,,1] A<-cm2ca(S,p1,p2) range(S-ca2cm(A))
Estimate a covariance matrix by adaptively shrinking the core.
covCSE(data, n = NULL, p1 = NULL, p2 = NULL, tol = 1e-08)covCSE(data, n = NULL, p1 = NULL, p2 = NULL, tol = 1e-08)
data |
either a numeric n*p1*p2 array consisting of n data matrices each of dimension p1*p2, or a p1*p2 covariance matrix of data of this type. If the latter, the values of n, p1 and p2 must be specified. |
n |
the sample size. |
p1 |
the row dimension of the data matrices. |
p2 |
the column dimension of the data matrices. |
tol |
the convergence tolerance of the iterative algorithm. |
a covariance matrix of the same dimension as S.
The attribute w of S gives the shrinkage weight on the
Kronecker covariance of S.
Peter Hoff
p1<-4 ; p2<-3 ; n<-20 # create a matrix Y with separable covariance Sig1<-rWishart(1,p1,diag(p1))[,,1] Sig2<-rWishart(1,p2,diag(p2))[,,1] Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2)) Y<-aperm( apply(Y,c(1,3),function(y){ msqrt(Sig1)%*%y } ),c(2,1,3)) Y<-aperm( apply(Y,c(1,2),function(y){ msqrt(Sig2)%*%y } ),c(2,3,1)) # covariance S<-mcov(Y) covCSE(S,n,p1,p2) # now an unstructured covariance S<-rWishart(1,p1*p2,diag(p1*p2))[,,1] covCSE(S,n,p1,p2)p1<-4 ; p2<-3 ; n<-20 # create a matrix Y with separable covariance Sig1<-rWishart(1,p1,diag(p1))[,,1] Sig2<-rWishart(1,p2,diag(p2))[,,1] Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2)) Y<-aperm( apply(Y,c(1,3),function(y){ msqrt(Sig1)%*%y } ),c(2,1,3)) Y<-aperm( apply(Y,c(1,2),function(y){ msqrt(Sig2)%*%y } ),c(2,3,1)) # covariance S<-mcov(Y) covCSE(S,n,p1,p2) # now an unstructured covariance S<-rWishart(1,p1*p2,diag(p1*p2))[,,1] covCSE(S,n,p1,p2)
Computes the Kronecker-core decomposition of a covariance matrix.
covKCD(S, p1, p2, tol = 1e-08)covKCD(S, p1, p2, tol = 1e-08)
S |
a covariance matrix of dimension (p1p2)*(p1p2). |
p1 |
the row dimension. |
p2 |
the column dimension. |
tol |
the convergence tolerance of the iterative algorithm. |
The Kronecker-core decomposition is a representation of an arbitrary covariance matrix S in terms of a separable Kronecker covariance matrix K and a complementary non-separable core covariance matrix C. The Kronecker covariance is the separable covariance matrix that is closest to S in terms of the divergence function
The core covariance matrix C is computed from S and K via
covKCD returns a list with the following elements:
the Kronecker covariance matrix;
the core covariance matrix;
the row covariance matrix;
the column covariance matrix;
the divergence between S and K across iterations of the algorithm.
Peter Hoff
p1<-4 ; p2<-3 ; n<-200 # create a matrix Y with separable covariance A<-matrix(rnorm(p1*p1),p1,p1) B<-matrix(rnorm(p2*p2),p2,p2)/3 Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2)) Y<-aperm( apply(Y,c(1,3),function(y){ A%*%y } ),c(2,1,3)) Y<-aperm( apply(Y,c(1,2),function(y){ B%*%y } ),c(2,3,1)) # covariance S<-mcov(Y) KCD<-covKCD(S,p1,p2) plot(A%*%t(A), KCD$K1) plot(B%*%t(B), KCD$K2)p1<-4 ; p2<-3 ; n<-200 # create a matrix Y with separable covariance A<-matrix(rnorm(p1*p1),p1,p1) B<-matrix(rnorm(p2*p2),p2,p2)/3 Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2)) Y<-aperm( apply(Y,c(1,3),function(y){ A%*%y } ),c(2,1,3)) Y<-aperm( apply(Y,c(1,2),function(y){ B%*%y } ),c(2,3,1)) # covariance S<-mcov(Y) KCD<-covKCD(S,p1,p2) plot(A%*%t(A), KCD$K1) plot(B%*%t(B), KCD$K2)
Compute the logarithm of the multivariate gamma function
.
lmvgamma(a, p)lmvgamma(a, p)
a |
a numeric scalar. |
p |
a positive integer. |
a scalar
Peter Hoff
Compute the covariance matrix of a sample of data matrices.
mcov(Y, use = "everything")mcov(Y, use = "everything")
Y |
a numeric n*p1*p2 data array corresponding to n data matrices of dimension p1*p2. |
use |
a character string giving method for dealing with missing
values, fed to the |
a p1*p2 by p1*p2 sample covariance matrix of the n vectorized data matrices.
Peter Hoff
p1<-4 ; p2<-3 ; n<-200 # create a matrix Y with separable covariance Sig1<-rWishart(1,p1,diag(p1))[,,1] Sig2<-rWishart(1,p2,diag(p2))[,,1] Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2)) Y<-aperm( apply(Y,c(1,3),function(y){ msqrt(Sig1)%*%y } ),c(2,1,3)) Y<-aperm( apply(Y,c(1,2),function(y){ msqrt(Sig2)%*%y } ),c(2,3,1)) # covariance S<-mcov(Y) image(S) plot(S,kronecker(Sig2,Sig1)) ; abline(0,1)p1<-4 ; p2<-3 ; n<-200 # create a matrix Y with separable covariance Sig1<-rWishart(1,p1,diag(p1))[,,1] Sig2<-rWishart(1,p2,diag(p2))[,,1] Y<-array(rnorm(n*p1*p2),dim=c(n,p1,p2)) Y<-aperm( apply(Y,c(1,3),function(y){ msqrt(Sig1)%*%y } ),c(2,1,3)) Y<-aperm( apply(Y,c(1,2),function(y){ msqrt(Sig2)%*%y } ),c(2,3,1)) # covariance S<-mcov(Y) image(S) plot(S,kronecker(Sig2,Sig1)) ; abline(0,1)
Compute the symmetric square root of a matrix.
msqrt(M)msqrt(M)
M |
a positive semidefinite matrix. |
a positive semidefinite matrix.
Peter Hoff
S<-rWishart(1,5,diag(5))[,,1] S Sh<-msqrt(S) Sh%*%ShS<-rWishart(1,5,diag(5))[,,1] S Sh<-msqrt(S) Sh%*%Sh
Compute the inverse of the symmetric square root of a matrix.
msqrtInv(M)msqrtInv(M)
M |
a positive definite matrix. |
a positive definite matrix.
Peter Hoff
S<-rWishart(1,5,diag(5))[,,1] solve(S) iSh<-msqrtInv(S) iSh%*%iShS<-rWishart(1,5,diag(5))[,,1] solve(S) iSh<-msqrtInv(S) iSh%*%iSh