Doubly-enhanced EM algorithm

Description

Doubly-enhanced EM algorithm for tensor clustering

Usage

DEEM(X, nclass, niter = 100, lambda = NULL, dfmax = n, pmax = nvars, pf = rep(1, nvars),
eps = 1e-04, maxit = 1e+05, sml = 1e-06, verbose = FALSE, ceps = 0.1,
initial = TRUE, vec_x = NULL)

Arguments

Details

The DEEM function implements the Doubly-Enhanced EM algorithm (DEEM) for tensor clustering. The observations \mathbf{X}_i are assumed to be following the tensor normal mixture model (TNMM) with common covariances across different clusters:

\mathbf{X}_i\sim\sum_{k=1}^K\pi_k \mathrm{TN}(\bm{\mu}_k;\bm{\Sigma}_1,\ldots,\bm{\Sigma}_M),\quad i=1,\dots,n,

where 0<\pi_k<1 is the prior probability for \mathbf{X} to be in the k-th cluster such that \sum_{k=1}^{K}\pi_k=1, \bm{\mu}_k is the cluster mean of the k-th cluster and \bm{\Sigma}_1,\ldots,\bm{\Sigma}_M) are the common covariances across different clusters. Under the TNMM framework, the optimal clustering rule can be showed as

\widehat{Y}^{opt}=\arg\max_k\{\log\pi_k+\langle\mathbf{X}-(\bm{\mu}_1+\bm{\mu}_k)/2,\mathbf{B}_k\rangle\},

where \mathbf{B}_k=[\![\bm{\mu}_k-\bm{\mu}_1;\bm{\Sigma}_1^{-1},\ldots,\bm{\Sigma}_M^{-1}]\!]. In the enhanced E-step, DEEM imposes sparsity directly on the optimal clustering rule as a flexible alternative to popular low-rank assumptions on tensor coefficients \mathbf{B}_k as

\min_{\mathbf{B}_2,\dots,\mathbf{B}_K}\bigg[\sum_{k=2}^K(\langle\mathbf{B}_k,[\![\mathbf{B}_k,\widehat{\bm{\Sigma}}_1^{(t)},\ldots,\widehat{\bm{\Sigma}}_M^{(t)}]\!]\rangle-2\langle\mathbf{B}_k,\widehat{\bm{\mu}}_k^{(t)}-\widehat{\bm{\mu}}_1^{(t)}\rangle) +\lambda^{(t+1)}\sum_{\mathcal{J}}\sqrt{\sum_{k=2}^Kb_{k,\mathcal{J}}^2}\bigg],

where \lambda^{(t+1)} is a tuning parameter. In the enhanced M-step, DEEM employs a new estimator for the tensor correlation structure, which facilitates both the computation and the theoretical studies.

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Mai, Q., Zhang, X., Pan, Y. and Deng, K. (2021). A Doubly-Enhanced EM Algorithm for Model-Based Tensor Clustering. Journal of the American Statistical Association.

See Also

tune_lamb, tune_K

Examples

dimen = c(5,5,5)
nvars = prod(dimen)
K = 2
n = 100
sigma = array(list(),3)

sigma[[1]] = sigma[[2]] = sigma[[3]] = diag(5)

B2=array(0,dim=dimen)
B2[1:3,1,1]=2

y = c(rep(1,50),rep(2,50))
M = array(list(),K)
M[[1]] = array(0,dim=dimen)
M[[2]] = B2

vec_x=matrix(rnorm(n*prod(dimen)),ncol=n)
X=array(list(),n)
for (i in 1:n){
  X[[i]] = array(vec_x[,i],dim=dimen)
  X[[i]] = M[[y[i]]] + X[[i]]
}

myfit = DEEM(X, nclass=2, lambda=0.05)

Fit the Tensor Envelope Mixture Model (TEMM)

Description

Fit the Tensor Envelope Mixture Model (TEMM)

Usage

TEMM(Xn, u, K, initial = "kmeans", iter.max = 500, 
stop = 1e-3, trueY = NULL, print = FALSE)

Arguments

Details

The TEMM function fits the Tensor Envelope Mixture Model (TEMM) through a subspace-regularized EM algorithm. For mode m, let (\bm{\Gamma}_m,\bm{\Gamma}_{0m})\in R^{p_m\times p_m} be an orthogonal matrix where \bm{\Gamma}_{m}\in R^{p_{m}\times u_{m}}, u_{m}\leq p_{m}, represents the material part. Specifically, the material part \mathbf{X}_{\star,m}=\mathbf{X}\times_{m}\bm{\Gamma}_{m}^{T} follows a tensor normal mixture distribution, while the immaterial part \mathbf{X}_{\circ,m}=\mathbf{X}\times_{m}\bm{\Gamma}_{0m}^{T} is unimodal, independent of the material part and hence can be eliminated without loss of clustering information. Dimension reduction is achieved by focusing on the material part \mathbf{X}_{\star,m}=\mathbf{X}\times_{m}\bm{\Gamma}_{m}^{T}. Collectively, the joint reduction from each mode is

\mathbf{X}_{\star}=[\![\mathbf{X};\bm{\Gamma}_{1}^{T},\dots,\bm{\Gamma}_{M}^{T}]\!]\sim\sum_{k=1}^{K}\pi_{k}\mathrm{TN}(\bm{\alpha}_{k};\bm{\Omega}_{1},\dots,\bm{\Omega}_{M}),\quad \mathbf{X}_{\star}\perp\!\!\!\perp\mathbf{X}_{\circ,m},

where \bm{\alpha}_{k}\in R^{u_{1}\times\cdots\times u_{M}} and \bm{\Omega}_m\in R^{u_m\times u_m} are the dimension-reduced clustering parameters and \mathbf{X}_{\circ,m} does not vary with cluster index Y. In the E-step, the membership weights are evaluated as

\widehat{\eta}_{ik}^{(s)}=\frac{\widehat{\pi}_{k}^{(s-1)}f_{k}(\mathbf{X}_i;\widehat{\bm{\theta}}^{(s-1)})}{\sum_{k=1}^{K}\widehat{\pi}_{k}^{(s-1)}f_{k}(\mathbf{X}_i;\widehat{\bm{\theta}}^{(s-1)})},

where f_k denotes the conditional probability density function of \mathbf{X}_i within the k-th cluster. In the subspace-regularized M-step, the envelope subspace is iteratively estimated through a Grassmann manifold optimization that minimize the following log-likelihood-based objective function:

G_m^{(s)}(\bm{\Gamma}_m) = \log|\bm{\Gamma}_m^T \mathbf{M}_m^{(s)} \bm{\Gamma}_m|+\log|\bm{\Gamma}_m^T (\mathbf{N}_m^{(s)})^{-1} \bm{\Gamma}_m|,

where \mathbf{M}_{m}^{(s)} and \mathbf{N}_{m}^{(s)} are given by

\mathbf{M}_m^{(s)} = \frac{1}{np_{-m}}\sum_{i=1}^{n} \sum_{k=1}^{K}\widehat{\eta}_{ik}^{(s)} (\bm{\epsilon}_{ik}^{(s)})_{(m)}(\widehat{\bm{\Sigma}}_{-m}^{(s-1)})^{-1} (\bm{\epsilon}_{ik}^{(s)})_{(m)}^T,

\mathbf{N}_m^{(s)} = \frac{1}{np_{-m}}\sum_{i=1}^{n} (\mathbf{X}_i)_{(m)}(\widehat{\bm{\Sigma}}_{-m}^{(s-1)})^{-1}(\mathbf{X}_i)_{(m)}^T.

The intermediate estimators \mathbf{M}_{m}^{(s)} can be viewed the mode-m conditional variation estimate of \mathbf{X}\mid Y and \mathbf{N}_{m}^{(s)} is the mode-m marginal variation estimate of \mathbf{X}.

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Deng, K. and Zhang, X. (2021). Tensor Envelope Mixture Model for Simultaneous Clustering and Multiway Dimension Reduction. Biometrics.

See Also

TGMM, tune_u_sep, tune_u_joint

Examples

  A = array(c(rep(1,20),rep(2,20))+rnorm(40),dim=c(2,2,10))
  myfit = TEMM(A,u=c(2,2),K=2)

Fit the Tensor Gaussian Mixture Model (TGMM)

Description

Fit the Tensor Gaussian Mixture Model (TGMM)

Usage

TGMM(Xn, K, shape = "shared", initial = "kmeans", 
iter.max = 500, stop = 1e-3, trueY = NULL, print = FALSE)

Arguments

Details

The TGMM function fits the Tensor Gaussian Mixture Model (TGMM) through the classical EM algorithm. TGMM assumes the following tensor normal mixture distribution of M-way tensor data \mathbf{X}:

\mathbf{X}\sim\sum_{k=1}^K\pi_k \mathrm{TN}(\bm{\mu}_k,\mathcal{M}_k),\quad i=1,\dots,n,

where 0<\pi_k<1 is the prior probability for \mathbf{X} to be in the k-th cluster such that \sum_{k=1}^{K}\pi_k=1, \bm{\mu}_k is the mean of the k-th cluster, \mathcal{M}_k \equiv \{\bm{\Sigma}_{km}, m=1,\dots,M\} is the set of covariances of the k-th cluster. If \mathcal{M}_k's are the same for k=1,\dots,K, call TGMM with argument shape="shared".

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Deng, K. and Zhang, X. (2021). Tensor Envelope Mixture Model for Simultaneous Clustering and Multiway Dimension Reduction. Biometrics.

Tait, P. A. and McNicholas, P. D. (2019). Clustering higher order data: Finite mixtures of multidimensional arrays. arXiv:1907.08566.

See Also

TEMM

Examples

  A = array(c(rep(1,20),rep(2,20))+rnorm(40),dim=c(2,2,10))
  myfit = TGMM(A,K=2,shape="shared")

Select the number of clusters K in DEEM

Description

Select the number of clusters K along with tuning parameter lambda through BIC in DEEM.

Usage

tune_K(X, seqK, seqlamb, initial = TRUE, vec_x = NULL)

Arguments

Details

The tune_K function runs tune_lamb function length(seqK) times to choose the tuning parameter \lambda and number of clusters K simultaneously. Let \widehat{\bm{\theta}}^{\{\lambda,K\}} be the output of DEEM with the tuning parameter and number of clusters fixed at \lambda and K respectively, tune_K looks for the values of \lambda and K that minimizes

\mathrm{BIC}(\lambda,K)=-2\sum_{i=1}^n\log(\sum_{k=1}^K\widehat{\pi}^{\{\lambda,K\}}_kf_k(\mathbf{X}_i;\widehat{\bm{\theta}}_k^{\{\lambda,K\}}))+\log(n)\cdot |\widehat{\mathcal{D}}^{\{\lambda,K\}}|,

where \widehat{\mathcal{D}}^{\{\lambda,K\}}=\{(k, {\mathcal{J}}): \widehat b_{k,{\mathcal{J}}}^{\lambda} \neq 0 \} is the set of nonzero elements in \widehat{\bm{B}}_2^{\{\lambda,K\}},\ldots,\widehat{\bm{B}}_K^{\{\lambda,K\}}. The tune_K function intrinsically selects the initial point and return the optimal estimated labels.

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Mai, Q., Zhang, X., Pan, Y. and Deng, K. (2021). A Doubly-Enhanced EM Algorithm for Model-Based Tensor Clustering. Journal of the American Statistical Association.

See Also

DEEM, tune_lamb

Examples

dimen = c(5,5,5)
nvars = prod(dimen)
K = 2
n = 100
sigma = array(list(),3)

sigma[[1]] = sigma[[2]] = sigma[[3]] = diag(5)

B2=array(0,dim=dimen)
B2[1:3,1,1]=2

y = c(rep(1,50),rep(2,50))
M = array(list(),K)
M[[1]] = array(0,dim=dimen)
M[[2]] = B2

vec_x=matrix(rnorm(n*prod(dimen)),ncol=n)
X=array(list(),n)
for (i in 1:n){
  X[[i]] = array(vec_x[,i],dim=dimen)
  X[[i]] = M[[y[i]]] + X[[i]]
}

mytune = tune_K(X, seqK=2:4, seqlamb=seq(0.01,0.1,by=0.01))

Parameter tuning in enhanced E-step in DEEM

Description

Perform parameter tuning through BIC in DEEM.

Usage

tune_lamb(X, K, seqlamb, initial = TRUE, vec_x = NULL)

Arguments

Details

The tune_lamb function adopts a BIC-type criterion to select the tuning parameter \lambda in the enhanced E-step. Let \widehat{\bm{\theta}}^{\lambda} be the output of DEEM with the tuning parameter fixed at \lambda, tune_lamb looks for the value of \lambda that minimizes

\mathrm{BIC}(\lambda)=-2\sum_{i=1}^n\log(\sum_{k=1}^K\widehat{\pi}^{\lambda}_kf_k(\mathbf{X}_i;\widehat{\bm{\theta}}_k^{\lambda}))+\log(n)\cdot |\widehat{\mathcal{D}}^{\lambda}|,

where \widehat{\mathcal{D}}^{\lambda}=\{(k, {\mathcal{J}}): \widehat b_{k,{\mathcal{J}}}^{\lambda} \neq 0 \} is the set of nonzero elements in \widehat{\bm{B}}_2^{\lambda},\ldots,\widehat{\bm{B}}_K^{\lambda}. The tune_lamb function intrinsically selects the initial point and return the optimal estimated labels.

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Mai, Q., Zhang, X., Pan, Y. and Deng, K. (2021). A Doubly-Enhanced EM Algorithm for Model-Based Tensor Clustering. Journal of the American Statistical Association.

See Also

DEEM, tune_K

Examples

dimen = c(5,5,5)
nvars = prod(dimen)
K = 2
n = 100
sigma = array(list(),3)

sigma[[1]] = sigma[[2]] = sigma[[3]] = diag(5)

B2=array(0,dim=dimen)
B2[1:3,1,1]=2

y = c(rep(1,50),rep(2,50))
M = array(list(),K)
M[[1]] = array(0,dim=dimen)
M[[2]] = B2

vec_x=matrix(rnorm(n*prod(dimen)),ncol=n)
X=array(list(),n)
for (i in 1:n){
  X[[i]] = array(vec_x[,i],dim=dimen)
  X[[i]] = M[[y[i]]] + X[[i]]
}

mytune = tune_lamb(X, K=2, seqlamb=seq(0.01,0.1,by=0.01))

Tuning envelope dimension jointly by BIC in TEMM.

Description

Tuning envelope dimension jointly by BIC in TEMM.

Usage

tune_u_joint(u_candi, K, X, iter.max = 500, stop = 0.001, trueY = NULL)

Arguments

Details

The tune_u_joint function searches over all the combinations of u\equiv(u_1,\dots,u_M) in the neighborhood of \widetilde{u}, \mathcal{N}(\widetilde u)=\{u:\ \max(1,\widetilde u_m-2) \leq u_m \leq \min(\widetilde u_m+2,p_m),\ m=1,\dots,M\}, that minimizes

\mathrm{BIC}(u) = -2\sum_{i=1}^{n}\log(\sum_{k=1}^{K}\widehat{\pi}_k^u f_k(\mathbf{X}_i;\widehat{\bm{\theta}}^u)) + \log(n)\cdot K_u.

In the above BIC, K_u=(K-1)\prod_{m=1}^M u_m + \sum_{m=1}^{M}p_m(p_m+1)/2 is the total number of parameters in TEMM, \widehat{\pi}_k^u and \widehat{\bm{\theta}}^{u} are the estimated parameters with envelope dimension fixed at u. The tune_u_joint function intrinsically selects the initial point and return the optimal estimated labels.

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Deng, K. and Zhang, X. (2021). Tensor Envelope Mixture Model for Simultaneous Clustering and Multiway Dimension Reduction. Biometrics.

See Also

TEMM, tune_u_sep

Examples

  A = array(c(rep(1,20),rep(2,20))+rnorm(40),dim=c(2,2,10))
  mytune = tune_u_joint(u_candi=list(1:2,1:2),K=2,A)

Tuning envelope dimension separately by BIC in TEMM.

Description

Tuning envelope dimension separately by BIC in TEMM.

Usage

tune_u_sep(m, u_candi, K, X, C = 1, oneD = TRUE, 
iter.max = 500, stop = 0.001, trueY = NULL)

Arguments

Details

For tensor mode m=1,\dots,M, the tune_u_sep function selects the envelope dimension \widetilde{u}_m by minimizing the following BIC-type criterion over the set \{0,1,\dots,p_m\},

\mathrm{BIC}_m(u_m) = \log|\bm{\Gamma}_m^T \widehat{\mathbf{M}}_m \bm{\Gamma}_m|+\log|\bm{\Gamma}_{m}^T \widehat{\mathbf{N}}_m^{-1} \bm{\Gamma}_{m}| + C \cdot u_m \log(n)/n.

This separate selection over each mode m is less sensitive to the complex interrelationships of each mode of the tensor. The default constant C is set as 1 as suggested by Zhang and Mai (2018).

Value

Author(s)

Kai Deng, Yuqing Pan, Xin Zhang and Qing Mai

References

Deng, K. and Zhang, X. (2021). Tensor Envelope Mixture Model for Simultaneous Clustering and Multiway Dimension Reduction. Biometrics.

Zhang, X. and Mai, Q. (2018). Model-free envelope dimension selection. Electronic Journal of Statistics 12, 2193-2216.

See Also

TEMM, tune_u_joint

Examples

  A = array(c(rep(1,20),rep(2,20))+rnorm(40),dim=c(2,2,10))
  mytune = tune_u_sep(1,1:2,K=2,A)