Processing method for high order tensor data

Abstract

A processing method for high-order tensor data, which can avoid that the vectorization process of the image observation sample set damage the internal structure of the data, simplify the large amount of redundant information in the high-order tensor data in the image observation sample set, and improve the image processing speed. In this processing method for high-order tensor data, the high-order tensor data are divided into three parts: the shared subspace component, the personality subspace component and the noise part; the shared subspace component and the personality subspace component respectively represent the high-order tensor data as a group of linear combination of the tensor base and the vector coefficient; the variational EM method is used to solve the base tensor and the vector coefficient; design a classifier to classify the test samples by comparing the edge distribution of samples.

Claims

1. A processing method for high-order tensor data, comprising: dividing the high-order tensor data into three parts comprising a shared subspace component, a personality subspace component and a noise part; wherein, the shared subspace component and the personality subspace component respectively represent the high-order tensor data as a group of linear combination of a tensor base and a vector coefficient; using a variational EM method to solve the tensor base and the vector coefficient; and designing a classifier to classify a plurality of test samples by comparing an edge distribution of the plurality of samples.

2. The processing method for high-order tensor data according to claim 1, further comprising forming a two dimensional model assuming an observed data set y={Y.sub.ij}.sub.i=1,j=1.sup.I,N.sup.i; wherein the observed data set consist of N images of I individuals and each individual has Ni images; Yij represents a jth image of an ith individual and Yij is modeled as formula (1):
Y.sub.ij=×.sub.3h.sub.i.sup.T+×.sub.3g.sub.ij.sup.T+M+E.sub.ij  (1) in the formula (1), ×.sub.3 h.sub.i.sup.T or ×.sub.3 g.sub.ij.sup.T is a generalization of vector representation for a given sample, and are factor tensors consisting of a plurality of base matrices in a shared subspace and an individual subspace, respectively; ×.sub.3 represents a product of tensor and matrix (vector), h.sub.i∈.sup.K.sup.1 is a shared latent variable crossing all of the plurality of images Y.sub.i1, . . . , Y.sub.iN.sub.i of object/class i and g.sub.ij∈.sup.K.sup.2 represents the vector coefficient of the jth sample of the ith object in the individual space; M is a mean matrix of the high-order data; E.sub.ij is a stochastic noise term, each element in the stochastic noise term follows the normal distribution (0, σ.sup.2).

3. The processing method for high-order tensor data according to the claim 2, wherein, the vector of the formula (1) is formula (2)
vec(Y.sub.ij)=Fh.sub.i+Wg.sub.ij+vec(M)+vec(E.sub.ij)  (2) wherein, columns f.sub.k.sub.1 and w.sub.k.sub.2 in the matrix F and W is respectively k1, k2 face of tensor and , are used to restrict the base tensor by CP decomposition of tensor; wherein, $ℱ=〚\lambda ;F^{\left(1\right)},F^{\left(2\right)},F^{\left(h\right)}〛=\underset{r=1}{\overset{R_1}{.Math.}}{\lambda }_r{f}_{1,r}^{\left(1\right)}\circ f_{1,r}^{\left(2\right)}\circ f_{1,r}^{\left(h\right)}\mathrm{and}$ $\mathrm{��}=〚\eta ;W^{\left(1\right)},W^{\left(2\right)},W^{\left(g\right)}〛=\underset{r=1}{\overset{R_2}{.Math.}}{\eta }_r{w}_{1,r}^{\left(1\right)}\circ w_{1,r}^{\left(2\right)}\circ w_{1,r}^{\left(g\right)}.$

4. The processing method for high-order tensor data according to the claim 3, wherein, Q(h, g, ρ) is any joint distribution over a plurality of hidden variables, a lower bound on a likelihood of the observed data is formula (2), $\begin{array}{cc}\begin{array}{c}ℒ\left(\Theta \right)=\log p\left(y.Math.\Theta \right)\\ \ge {\int }_h{\int }_g{\int }_{\rho }Q\left(h,g,\rho \right)\log \frac{p\left(y,h,g,\rho \right)}{Q\left(h,g,\rho \right)}\mathrm{dhdgd}\rho \\ =E_{h,g,\rho }\left[\log p\left(y.Math.h,g,\rho \right)\right]+E_h\left[\log \frac{p\left(h\right)}{Q\left(h\right)}\right]+\\ E_g\left[\log \frac{p\left(g\right)}{Q\left(g\right)}\right]+E_{\rho }\left[\log \frac{p\left(\rho \right)}{Q\left(\rho \right)}\right]\\ \stackrel{\Delta }{=}\widehat{ℒ}\left(Q\left(h\right),Q\left(g\right),Q\left(\rho \right),\Theta \right)\end{array}& \left(3\right)\end{array}$ wherein, Q(h, g, ρ)=Q(h)Q(g)Q(ρ).

5. The processing method for high-order tensor data according to claim 1, wherein, a high-order observed data set is {y.sub.ij|y.sub.ij∈.sup.D.sup.1.sup.× . . . ×D.sup.M; i=1, . . . , I, j=1 . . . , N.sub.i, Σ.sub.iN.sub.i=N}; each sample is a tensor of order M, and a linear discriminant analysis model of the higher-order tensor is formula (4)
y.sub.ij=×.sub.(M+1)h.sub.i.sup.T+×.sub.(M+1)g.sub.ij.sup.T+.sub.ij  (4); wherein h.sub.i∈.sup.K.sup.1, g.sub.ij∈.sup.K.sup.2, and are a tensor of order (M+1), wherein a size of and is respectively D.sub.1× . . . ×D.sub.M×K.sub.1 and D.sub.1× . . . ×D.sub.M×K.sub.2; wherein, the base tensor has a structure of CP decomposition, $\begin{array}{c}ℱ=〚F^{\left(1\right)},F^{\left(2\right)},\mathrm{.Math.},F^{\left(M\right)},F^{\left(h\right)}〛\\ =\underset{r=1}{\overset{R_1}{.Math.}}{f}_{1,r}^{\left(1\right)}\circ f_{1,r}^{\left(2\right)}\circ \mathrm{.Math.}\circ f_{1,r}^{\left(M\right)}\circ f_{1,r}^{\left(h\right)}\end{array}$ $\mathrm{and}$ $\begin{array}{c}\mathrm{��}=〚W^{\left(1\right)},W^{\left(2\right)},\mathrm{.Math.},W^{\left(M\right)},W^{\left(g\right)}〛\\ =\underset{r=1}{\overset{R_2}{.Math.}}{w}_{1,r}^{\left(1\right)}\circ w_{1,r}^{\left(2\right)}\circ \mathrm{.Math.}\circ w_{1,r}^{\left(M\right)}\circ w_{1,r}^{\left(g\right)}.\end{array}$

6. The processing method for high-order tensor data according to claim 5, wherein, for a higher-order model, an EM algorithm is used to solve a problem, a posterior of h.sub.i is a multivariate normal distribution; a mean matrix and a covariance matrix are, respectively, $\begin{array}{c}{\stackrel{\_}{h}}_i=\left(\underset{ℱ}{.Math.}+\frac{1}{\stackrel{\_}{\rho }{N}_i}{I}_{K_1}\right)a_i\\ \underset{h}{.Math.}={\left(I_{K1}+\frac{\stackrel{\_}{\rho }N}{I}{.Math.}_{ℱ}\right)}^{-1}\end{array}$ where,
=F.sup.(h)[(F.sup.(m)TF.sup.(m))⊙(F.sup.(m)TF.sup.(m))]F.sup.(h)T
a.sub.i=F.sup.(h)diag(F.sup.(m)TY.sub.(m).sup.iF.sup.(m)), where ⊙ represents the Hadamard elementwise product; ${\stackrel{\_}{y}}_i=\frac{1}{N_i}\underset{j=1}{\overset{N_i}{.Math.}}\left(y_{\mathrm{ij}}-\mathrm{��}{\times }_{\left(M+1\right)}{g}_{\mathrm{ij}}\right),$ wherein, Y.sub.(m).sup.i represents a mode-m metricized version of y.sub.i, and
F.sup.(m)=F.sup.(M) . . . F.sup.(m+1)F.sup.(m−1) . . . F.sup.(1) where represents the Khatri-Rao product; by computing, the posterior of g.sub.ij is a multivariate Gauss distribution, and a mean and a covariance of a latent variable g.sub.ij are
g.sub.ij=+1/ρI.sub.K).sup.−1b.sub.ij
Σ.sub.g=(I.sub.K+ρ).sup.−1
where
b.sub.ij=W.sup.(g)diag(W.sup.(m)TŶ.sub.(m).sup.ijW.sup.(m))
=W.sup.(g)[(W.sup.(m)TW.sup.(m))⊙(W.sup.(m)TW.sup.(m))]W.sup.(g)T wherein, ŷ.sub.ij=y.sub.ij−×.sub.(M+1)h.sub.i.sup.T, and Ŷ.sub.(m).sup.ij represents a mode-m metricized version of ŷ.sub.ij; and
W.sup.(m)=W.sup.(M) . . . W.sup.(m+1)W.sup.(m−1) . . . W.sup.(1) an optimal Q*(ρ) is a Gamma distribution; in M-step, taking derivatives of with respect to F.sup.(n) (n=1, . . . , N) and F.sup.(h), to be zero, obtain $F^{\left(m\right)}=\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\tilde{Y}}_{\left(m\right)}^{\mathrm{ij}}{\stackrel{\_}{F}}^{\left(m\right)}\mathrm{Diag}\left(h_i^T{F}^{\left(h\right)}\right)\right]\hspace{1em}{\left[\left({\stackrel{\_}{F}}^{\left(m\right)T}{\stackrel{\_}{F}}^{\left(m\right)}\right)\odot \left({\mathrm{NF}}^{\left(h\right)T}\underset{h}{.Math.}{F}^{\left(h\right)}+F^{\left(h\right)T}\tilde{H}{\tilde{H}}^T{F}^{\left(h\right)}\right)\right]}^{-1}{F}^{\left(h\right)}={\left[\tilde{H}{\tilde{H}}^T+N\underset{h}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{h_i\left[\mathrm{diag}\left({\stackrel{\_}{F}}^{\left(m\right)T}{\tilde{Y}}_{\left(m\right)}^{\left(\mathrm{ij}\right)T}{F}^{\left(m\right)}\right)\right]}^T\right]\hspace{1em}\left[\left(F^{\left(m\right)T}{F}^{\left(m\right)}\right)\odot \left({\stackrel{\_}{F}}^{\left(m\right)T}{\stackrel{\_}{F}}^{\left(m\right)}\right)\right]\text{-}1.$ wherein, {tilde over (Y)}.sub.(m).sup.ij represents a mode-m metricized version of {tilde over (y)}.sub.ij; and {tilde over (y)}.sub.ij=y.sub.ij−×.sub.(M+1) g.sub.ij, taking derivatives of with respect to W.sup.(m) (m=1, . . . , M) and W.sup.(h), to be zero, obtain $\begin{array}{cc}{W}^{\left(m\right)}=\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\hat{Y}}_{\left(m\right)}^{\mathrm{ij}}{\stackrel{\_}{W}}^{\left(m\right)}\mathrm{Diag}\left(g_{\mathrm{ij}}^T{W}^{\left(g\right)}\right)\right]\hspace{1em}{\left[\left({\stackrel{\_}{W}}^{\left(m\right)T}{\stackrel{\_}{W}}^{\left(m\right)}\right)\odot \left({\mathrm{NW}}^{\left(g\right)T}\underset{g}{.Math.}{W}^{\left(g\right)}+W^{\left(g\right)T}{\mathrm{GG}}^T{W}^{\left(g\right)}\right)\right]}^{-1}& \\ \mathrm{and}& \\ W^{\left(h\right)}={\left[{\mathrm{GG}}^T+N\underset{g}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{g_{\mathrm{ij}}\left[\mathrm{diag}\left({\stackrel{\_}{W}}^{\left(m\right)T}{\hat{Y}}_{\left(m\right)}^{\left(\mathrm{ij}\right)T}{W}^{\left(m\right)}\right)\right]}^T\right]\hspace{1em}{\left[\left(W^{\left(m\right)T}{W}^{\left(m\right)}\right)\odot \left({\stackrel{\_}{W}}^{\left(m\right)T}{W}^{\left(m\right)}\right)\right]}^{-1}.& \end{array}$

7. The processing method for high-order tensor data according to claim 6, wherein, for a sample Y.sub.p to be classified, in a category by calculating a likelihood under a model; a key point of recognition is the plurality of samples in a same class share a same identity variable; assuming Y.sub.1 . . . N is a plurality of observed samples, if Y.sub.p and one particular Yi belong to a same category, Y.sub.p and one particular Yi share a same h.

8. The processing method for high-order tensor data according to claim 7, wherein, Y1 and Yp belongs to the same category and both of identity variables are h1; a likelihood of the plurality of observed samples and Yp as p.sub.1 (Y.sub.1, Y.sub.2, Y.sub.p), Yp and Y2 came from the same category with the identity variable h2 and the likelihood is represented as p.sub.2 (Y.sub.1, Y.sub.2, Y.sub.p); then comparing p.sub.1 and p.sub.2, and finding out a largest likelihood value; where Yp belongs to a sample's label corresponding to the largest likelihood; the likelihood of the plurality of samples is calculated as follows; Y.sub.ij follows a matrix-variate normal conditioned on a latent variables h and g; proposing a model
y.sub.ij=Fh.sub.i+Wg.sub.ij+e.sub.ij wherein, y.sub.ij:=vec(Y.sub.ij), F=(F.sup.(2)⊙F.sup.(1))F.sup.(h)T and W=(W.sup.(2)⊙W.sup.(1))W.sup.(g)T; a plurality of factor matrices have been learned in a training phrase, the plurality of factor matrices are known to a recognition processing; the likelihood of N equals $\begin{array}{cc}{p}_n\left(y_{1\mathrm{.Math.}N,p}\right)=\underset{i=1,i\ne n}{\overset{N}{.Math.}}p\left(y_i\right)p\left(y_p,y_n\right);& \left(20\right)\end{array}$ illustrating a marginal distribution of n images y.sub.1 . . . n that share the same identity variable h; Forming generative equations $\begin{array}{cc}\left[\begin{array}{c}{y}_1\\ y_2\\ \mathrm{.Math.}\\ y_n\end{array}\right]=\left[\begin{array}{ccccc}F& W& 0& \mathrm{.Math.}& 0\\ F& 0& W& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ F& 0& 0& \mathrm{.Math.}& W\end{array}\right]\left[\begin{array}{c}h\\ g_1\\ g_2\\ \mathrm{.Math.}\\ g_n\end{array}\right]+\left[\begin{array}{c}{e}_1\\ e_2\\ \mathrm{.Math.}\\ e_n\end{array}\right],& \left(5\right)\end{array}$ wherein,
y′=Az′+m′+e′ taking the following probabilistic assumption
p(z′)=(0,I)
and
p(e′)=(0,σ.sup.2I) converting the model (5) into the form of probabilistic PCA to obtain the marginal distribution over the observed y′ by integrating out the latent variables as
p(y′)=(m′,AA.sup.T+σ.sup.2I)  (21); using (21), calculating both p(y.sub.i) and p(y.sub.p, y.sub.n) and in (20) and then obtaining p.sub.n(y.sub.1, . . . , N, p); comparing all p.sub.n(y.sub.1, . . . , N, p) and finding out a largest value, then dividing the sample Yp the sample's label corresponding to the largest p.sub.n(y.sub.1, . . . , N, p).

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows a schematic diagram of formula (1) according to the invention.

(2) FIG. 2 shows a graphic illustration of probabilistic recognition method according to the present invention.

(3) FIG. 3 shows a visual effect of the reconstruction of the invention on the extended Yale database.

(4) FIG. 4 shows reconstructed images of the invention on AR database.

DETAILED DESCRIPTION OF THE EMBODIMENTS

(5) As shown as FIG. 1, in this processing method for high-order tensor data, the high-order tensor data are divided into three parts: the shared subspace component, the personality subspace component and the noise part; the shared subspace component and the personality subspace component respectively represent the high-order tensor data as a group of linear combination of the tensor base and the vector coefficient; the variational EM method is used to solve the base tensor and the vector coefficient; design a classifier to classify the test samples by comparing the edge distribution of samples.

(6) Compared with the traditional linear discriminant method, the present invention extracts the common feature and the individual feature directly from the tensor data in two spaces respectively. In each space, by the structural feature of the tensor, it constructs a representation way from the tensor data to the vector data, so as to extract the discriminant feature of the tensor sample, thus avoiding that the vectorization process of the image observation sample set damage the internal structure of the data, which is flexible and simple, which can extract the discriminant features of the tensor sample, simplify the large amount of redundant information in the high-order tensor data in the image observation sample set, and improve the image processing speed.

(7) Preferably, two dimensional model is assuming an observed data set y={Y.sub.ij}.sub.i=1,j=1.sup.I,N.sup.i, which consist of N images of I individuals and each individual has Ni images. Yij represents the jth image of the ith individual and it will be modeled as formula (1):
Y.sub.ij=×.sub.3h.sub.i.sup.T+×.sub.3g.sub.ij.sup.T+M+E.sub.ij  (1)
In this model, ×.sub.3 h.sub.i.sup.T or ×.sub.3 g.sub.ij.sup.T can be regarded as the generalization of vector representation for a given sample, and are the factor tensors consisting of base matrices in shared subspace and individual subspace, respectively. ×.sub.3 represents the product of tensor and matrix (vector), h.sub.i∈.sup.K.sup.1 is a shared latent variable crossing all the images Y.sub.i1, . . . , Y.sub.iN.sub.i of object/class i and g.sub.ij∈.sup.K.sup.2 represents the coefficient of the jth sample of the ith object in individual space. M is the mean matrix of the data, suppose M=0. E.sub.ij is a stochastic noise term, each element in which follows the normal distribution (0, σ.sup.2).

(8) Preferably, the vector of model (1) is formula (2)
vec(Y.sub.ij)=Fh.sub.i+Wg.sub.ij+vec(M)+vec(E.sub.ij)  (2)
in which, the columns and f.sub.k.sub.1 and w.sub.k.sub.2 in matrix F and W is respectively k1, k2 face of tensor and , are used to restrict the base tensor by CP decomposition of tensor. Assuming

(9) $ℱ=〚\lambda ;F^{\left(1\right)},F^{\left(2\right)},F^{\left(h\right)}〛=\underset{r=1}{\overset{R_1}{.Math.}}{\lambda }_r{f}_{1,r}^{\left(1\right)}{\mathrm{•f}}_{1,r}^{\left(2\right)}{\mathrm{•f}}_{1,r}^{\left(h\right)}$$\mathrm{and}$$\mathrm{��}=〚\eta ;W^{\left(1\right)},W^{\left(2\right)},W^{\left(g\right)}〛=\underset{r=1}{\overset{R_2}{.Math.}}{\eta }_r{w}_{1,r}^{\left(1\right)}•w_{1,r}^{\left(2\right)}•w_{1,r}^{\left(g\right)}.$

(10) Preferably, assuming Q(h, g, ρ) is any joint distribution over the hidden variables, find that a lower bound on likelihood of the observed data is formula (2),

(11) $\begin{array}{cc}\begin{array}{c}ℒ(\ominus )=\log p(\mathrm{��}.Math.\ominus \\ \ge \hspace{0.17em}{\int }_h{\int }_g{\int }_{\rho }\mathrm{��}\left(h,g,\rho \right)\log \frac{p\left(\mathrm{��},h,g,\rho \right)}{\mathrm{��}\left(h,g,\rho \right)}\mathrm{dhdgd}\rho \\ =E_{h,g,\rho }\left[\log p\left(\mathrm{��}.Math.h,g,\rho \right)\right]+E_h\left[\log \frac{p\left(h\right)}{\mathrm{��}\left(h\right)}\right]\\ +E_g\left[\log \frac{p\left(g\right)}{\mathrm{��}\left(g\right)}\right]+E_{\rho }\left[\log \frac{p\left(\rho \right)}{\mathrm{��}\left(\rho \right)}\right]\\ \stackrel{\Delta }{=}\widehat{ℒ}\left(\mathrm{��}\left(h\right),\mathrm{��}\left(g\right),\mathrm{��}\left(\rho \right),\ominus \right)\end{array}\hspace{1em}& \left(3\right)\end{array}$
The second equation of the formula (3) is based on the feature of separation of Q(h, g, ρ), Q(h, g, ρ)=Q(h)Q(g)Q(ρ).

(12) Preferably, assuming a high-order observed data set {y.sub.ij|y.sub.ij∈.sup.D.sup.1.sup.× . . . ×D.sup.M; i=1, . . . , I, j=1, . . . , N.sub.i, Σ.sub.iN.sub.i=N}.

(13) each sample is a tensor of order M, and the linear discriminant analysis model of higher-order tensor is formula (4)
y.sub.ij=×.sub.(M+1)h.sub.i.sup.T+×.sub.(M+1)g.sub.ij.sup.T+.sub.ij  (4)
Wherein h.sub.i∈.sup.K.sup.1, g.sub.ij∈.sup.K.sup.2, and is a tensor of order (M+1), whose size is respectively D.sub.1× . . . ×D.sub.M×K.sub.1 and D.sub.1× . . . ×D.sub.M×K.sub.2. Assuming base tensor has the structure of CP decomposition,

(14) $\begin{array}{c}ℱ=〚F^{\left(1\right)},F^{\left(2\right)},\mathrm{.Math.},F^{\left(M\right)},F^{\left(h\right)}〛\\ =\underset{r=1}{\overset{R_1}{.Math.}}{f}_{1,r}^{\left(1\right)}•f_{1,r}^{\left(2\right)}•\mathrm{.Math.}•f_{1,r}^{\left(M\right)}•f_{1,r}^{\left(h\right)}\end{array}$$\mathrm{and}$$\begin{array}{c}\mathrm{��}=〚W^{\left(1\right)},W^{\left(2\right)},\mathrm{.Math.},W^{\left(M\right)},W^{\left(g\right)}〛\\ =\underset{r=1}{\overset{R_2}{.Math.}}{w}_{1,r}^{\left(1\right)}•w_{1,r}^{\left(2\right)}•\mathrm{.Math.}•w_{1,r}^{\left(M\right)}•w_{1,r}^{\left(g\right)}.\end{array}$

(15) Preferably, for higher-order model, EM algorithm is used to solve the problem, therefore the posterior of h.sub.i is still a multivariate normal distribution. The mean and covariance matrix are, respectively,

(16) $\begin{array}{c}{\stackrel{\_}{h}}_i=\left(\underset{ℱ}{.Math.}+\frac{1}{\stackrel{\_}{\rho }{N}_i}{I}_{K_1}\right)a_i\\ \underset{h}{.Math.}={\left(I_{K1}+\frac{\stackrel{\_}{\rho }N}{I}\underset{ℱ}{.Math.}\right)}^{-1}\end{array}$
where,
=F.sup.(h)[(F.sup.(m)TF.sup.(m))⊙(F.sup.(m)TF.sup.(m)]F.sup.(h)T
a.sub.i=F.sup.(h)diag(F.sup.(m)TY.sub.(m).sup.iF.sup.(m)).
Where ⊙ represents the Hadamard elementwise product.
Define

(17) ${\stackrel{\_}{\mathrm{��}}}_i=\frac{1}{N_i}\underset{j=1}{\overset{N_i}{.Math.}}\left({\mathrm{��}}_{\mathrm{ij}}-W{\times }_{\left(M+1\right)}{g}_{\mathrm{ij}}\right),$
then Y.sub.(m).sup.i represents the mode-m metricized version of y.sub.i, and
F.sup.(m)=F.sup.(M) . . . F.sup.(m+1)F.sup.(m−1) . . . F.sup.(1)
Where represents the Khatri-Rao product.
By computing, the posterior of g.sub.ij is still a multivariate Gauss distribution, and the mean and covariance of latent variable g.sub.ij have become
g.sub.ij=(+1/ρI.sub.K).sup.−1b.sub.ij
Σ.sub.g=(I.sub.K+ρ).sup.−1
where
b.sub.ij=W.sup.(g)diag(W.sup.(m)TŶ.sub.(m).sup.ijW.sup.(m))
Σ.sub.W=W.sup.(g)[(W.sup.(m)TW.sup.(m))⊙(W.sup.(m)TW.sup.(m))]W.sup.(g)T
Let ŷ.sub.ij=y.sub.ij−×.sub.(M+1)h.sub.i.sup.T, and Ŷ.sub.(m).sup.ij represents the mode-m metricized version of ŷ.sub.ij. And
W.sup.(m)=W.sup.(M) . . . W.sup.(m+1)W.sup.(m−1) . . . W.sup.(1)
The optimal Q*(ρ) is still a Gamma distribution.
In M-step, taking derivatives of with respect to F.sup.(n) (n=1, . . . , N) and F.sup.(h), and letting them be zero, obtain

(18) $F^{\left(m\right)}={\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\tilde{Y}}_{\left(m\right)}^{\mathrm{ij}}{\stackrel{\_}{F}}^{\left(m\right)}\mathrm{Diag}\left(h_i^T{F}^{\left(h\right)}\right)\right]\left[\left({\stackrel{\_}{F}}^{\left(m\right)T}{\stackrel{\_}{F}}^{\left(m\right)}\right)\odot \left({\mathrm{NF}}^{\left(h\right)T}\underset{h}{.Math.}{F}^{\left(h\right)}+F^{\left(h\right)T}\tilde{H}{\tilde{H}}^T{F}^{\left(h\right)}\right)\right]}^{-1}$$F^{\left(h\right)}={{\left[\tilde{H}{\tilde{H}}^T+N\underset{h}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{h_i\left[\mathrm{diag}\left({\stackrel{\_}{F}}^{\left(m\right)T}{\tilde{Y}}_{\left(m\right)}^{\left(\mathrm{ij}\right)T}{F}^{\left(m\right)}\right)\right]}^T\right]\left[\left(F^{\left(m\right)T}{F}^{\left(m\right)}\right)\odot \left({\stackrel{\_}{F}}^{\left(m\right)T}{\stackrel{\_}{F}}^{\left(m\right)}\right)\right]}^{-1}.$
Denote {tilde over (Y)}.sub.(m).sup.ij presents the mod-m metricized version of {tilde over (y)}.sub.ij. And {tilde over (y)}.sub.ij=y.sub.ij−×.sub.(M+1)g.sub.ij,
Taking derivatives of with respect to W.sup.(m) (m=1, . . . , M) and W.sup.(h), and letting them be zero, obtain

(19) $W^{\left(m\right)}={\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\tilde{Y}}_{\left(m\right)}^{\mathrm{ij}}{\stackrel{\_}{W}}^{\left(m\right)}\mathrm{Diag}\left(g_{\mathrm{ij}}^T{W}^{\left(g\right)}\right)\right]\left[\left({\stackrel{\_}{W}}^{\left(m\right)T}{\stackrel{\_}{W}}^{\left(m\right)}\right)\odot \left({\mathrm{NW}}^{\left(g\right)T}\underset{g}{.Math.}{W}^{\left(g\right)}+W^{\left(g\right)T}{\mathrm{GG}}^T{W}^{\left(g\right)}\right)\right]}^{-1}\mathrm{and}$$W^{\left(h\right)}={{\left[{\mathrm{GG}}^T+N\underset{g}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{g_{\mathrm{ij}}\left[\mathrm{diag}\left({\stackrel{\_}{W}}^{\left(m\right)T}{\tilde{Y}}_{\left(m\right)}^{\left(\mathrm{ij}\right)T}{W}^{\left(m\right)}\right)\right]}^T\right]\left[\left(W^{\left(m\right)T}{W}^{\left(m\right)}\right)\odot \left({\stackrel{\_}{W}}^{\left(m\right)T}{\stackrel{\_}{W}}^{\left(m\right)}\right)\right]}^{-1}.$

(20) Preferably, for a sample Y.sub.p to be classified, its category by calculating the likelihood under a model. The key point of recognition is that samples in the same class share the same identity variable. Assuming Y.sub.1 . . . N is observed samples, if Y.sub.p and one particular Yi belong to the same category, they shall share the same h.

(21) Preferably, Y1 and Yp belongs to the same category and both of identity variables are h1. In this case, the likelihood of observed samples and Yp as p.sub.1(Y.sub.1, Y.sub.2, Y.sub.p), Yp and Y2 came from the same category with the identity variable h2 and the likelihood can be represented as p.sub.2(Y.sub.1, Y.sub.2, Y.sub.p). Then we compare p.sub.1 and p.sub.2, and find out the largest likelihood value. Thus, Yp belongs to the sample's label corresponding to the largest likelihood. The likelihood of samples can be calculated as follows.

(22) Because Y.sub.ij follows a matrix-variate normal conditioned on the latent variables h and g. For convenience, we work with p(y.sub.ij) instead of p(Y.sub.ij). Here, we note y.sub.ij:=vec(Y.sub.ij). Thus, the proposed model (1) can become
y.sub.ij=Fh.sub.i+Wg.sub.ij+e.sub.ij
Where F=(F.sup.(2)⊙F.sup.(1))F.sup.(h)T and W=(W.sup.(2)⊙W.sup.(1))W.sup.(g)T
These factor matrices have been learned in training phrase, which are known to the recognition processing. Thus, the likelihood of N equals

(23) $\begin{array}{cc}{p}_n\left(y_{1\mathrm{.Math.}N,p}\right)=\underset{i=1,i\ne n}{\overset{N}{.Math.}}p\left(y_i\right)p\left(y_p,y_n\right)& \left(20\right)\end{array}$
Next, we illustrate the marginal distribution of n images y.sub.1 . . . n that share the same identity variable h. The generative equations become

(24) $\begin{array}{cc}\left[\begin{array}{c}{y}_1\\ y_2\\ \mathrm{.Math.}\\ y_n\end{array}\right]=\left[\begin{array}{ccccc}F& W& 0& \mathrm{.Math.}& 0\\ F& 0& W& \mathrm{.Math.}& 0\\ \mathrm{.Math.}& \mathrm{.Math.}& \mathrm{.Math.}& \ddots & \mathrm{.Math.}\\ F& 0& 0& \mathrm{.Math.}& W\end{array}\right]\left[\begin{array}{c}h\\ g_1\\ g_2\\ \mathrm{.Math.}\\ g_n\end{array}\right]+\left[\begin{array}{c}{e}_1\\ e_2\\ \mathrm{.Math.}\\ e_n\end{array}\right]& \left(5\right)\end{array}$
noting,
y′=Az′+m′+e′
Taking the following probabilistic assumption
p(z′)=(0,I)
and
p(e′)=(0,σ.sup.2I)
Converts the model (5) into the form of probabilistic PCA. Thus, the marginal distribution over the observed y′ can be easily obtained by integrating out the latent variables as
p(y′)=(m′,AA.sup.T+σ.sup.2I)  (21)
By (21), both p(y.sub.i) and p(y.sub.p, y.sub.n) in (20) can be calculated and then we can obtain P.sub.n(y.sub.1, . . . , N,p). Comparing all p.sub.n(y.sub.1, . . . , N,p) and finding out the largest value, then the sample Yp is divided into the sample's label corresponding to the largest p.sub.n(y.sub.1, . . . , N,p).

(25) The invention is described in more detail below.

(26) The technical solution adopted by the invention is a probabilistic linear discriminant analysis method based on tensor data vector expression, and the specific implementation process of the method is as follows.

(27) Two Dimensional Model Construction

(28) Assuming an observed data set y={Y.sub.ij}.sub.i=1,j=1.sup.I,N.sup.i, which consist of N images of I individuals and each individual has Ni images (Σ.sub.i=1.sup.IN.sub.i=N). Yij represents the jth image of the ith individual and it will be modeled as formula (1):
Y.sub.ij=×.sub.3h.sub.i.sup.T+×.sub.3g.sub.ij.sup.T+M+E.sub.ij  (1)

(29) In this model, ×.sub.3 h.sub.i.sup.T or ×.sub.3 g.sub.ij.sup.T can be regarded as the generalization of vector representation for a given sample, in which coefficient is vector, as shown as FIG. 1. and are the factor tensors consisting of base matrices in shared subspace and individual subspace, respectively. ×.sub.3 represents the product of tensor and matrix (vector). h.sub.i∈.sup.K.sup.1 is a shared latent variable crossing all the images Y.sub.i1, . . . , Y.sub.iN.sub.i of object/class i and g.sub.ij∈.sup.K.sup.2 represents the coefficient of the jth sample of the ith object in individual space. M is the mean matrix of the data. Without loss of generality, we suppose M=0 (a zero matrix). E.sub.ij is a stochastic noise term, each element in which follows the normal distribution (0, σ.sup.2).

(30) Furthermore, we impose a prior of multivariate normal distribution on latent variables to generate a Bayesian model

(31) 0$p\left(h_i\right)=\mathrm{��}\left(h_i.Math.0,I_{K_1}\right)={\left(\frac{1}{2\pi }\right)}^{K_1/2}\exp \left\{-\frac{1}{2}{h}_i^T{h}_i\right\}$$\mathrm{and}$$p\left(g_{\mathrm{ij}}\right)=\mathrm{��}\left(g_{\mathrm{ij}}.Math.0,I_{K_2}\right)={\left(\frac{1}{2\pi }\right)}^{K_2/2}\exp \left\{-\frac{1}{2}{g}_{\mathrm{ij}}^T{g}_{\mathrm{ij}}\right\}$

(32) Letting

(33) $\rho =\frac{1}{{\sigma }^2},$
we assume σ satisfies the following Gamma prior:

(34) $p_{\sigma }\left(\rho \right)=\Gamma \left(\rho .Math.a,b\right)=\frac{b^a}{\Gamma \left(a\right)}{\rho }^{a-1}\exp \left\{-b\rho \right\}$

(35) After vectorizing model (1), we can obtain the following model:
vec(Y.sub.ij)=Fh.sub.i+Wg.sub.ij+vec(M)+vec(E.sub.ij)  (2)

(36) in which, the columns f.sub.k.sub.1 and w.sub.k.sub.2 in matrix F and W is respectively k1, k2 face of tensor and , are used to restrict the base tensor by CP decomposition of tensor. Therefore, the model is transformed into the traditional probabilistic linear discriminant analysis. However, a major problem is that the parameters in the model will increase with the increase of the original data. For example, for N-order tensor data, the projections and are (N+1)-order tensors with K.sub.1Π.sub.nD.sub.n and K.sub.2Π.sub.nD.sub.n components, respectively. They will quickly reaches billions when the mode dimensionalities K.sub.1, K.sub.2, D.sub.1:N and N are moderate. These make it extremely difficult to learn parameters for base tensors. Hence, we propose to introduce tensor CP decompositions in constructing the multiplicative interactions in tensors and , respectively. Assuming

(37) $ℱ=〚\lambda ;F^{\left(1\right)},F^{\left(2\right)},F^{\left(h\right)}〛=\underset{r=1}{\overset{R_1}{.Math.}}{\lambda }_r{f}_{1,r}^{\left(1\right)}•f_{1,r}^{\left(2\right)}•f_{1,r}^{\left(h\right)}$$\mathrm{and}$$\mathrm{��}=〚\eta ;W^{\left(1\right)},W^{\left(2\right)},W^{\left(g\right)}〛=\underset{r=1}{\overset{R_2}{.Math.}}{\eta }_r{w}_{1,r}^{\left(1\right)}•w_{1,r}^{\left(2\right)}•w_{1,r}^{\left(g\right)}.$

(38) Solution of the Model

(39) For the above-mentioned model, h, g and ρ are the model latent variables and θ={F.sup.(1), F.sup.(2), F.sup.(h), W.sup.(1), W.sup.(2), W.sup.(g)} represents all the model parameters. Therefore, the joint distribution of a set of given observations y and all the latent variables can be expressed as

(40) $p\left(\mathrm{��},h,g,\rho .Math.\ominus \right)=\underset{i,j}{.Math.}\mathrm{��}\left(Y_i-ℱ\times {}_{3}h_i^T-\mathrm{��}\times {}_{3}w_{\mathrm{ij}}^T.Math.0,\sigma I,\sigma I\right)$$\mathrm{��}\left(h_i.Math.0,I_{K_1}\right)\mathrm{��}\left(g_{\mathrm{ij}}.Math.0,I_{K_2}\right)p_{\sigma }\left(\rho \right).$

(41) So the likelihood of the sample is
(θ)=log p(y|θ)=log ∫∫∫p(y,h,g,ρ|θ)dhdgdρ

(42) Furthermore, we investigate the variational electromagnetic (EM) algorithm for the proposed TPLDA model.

(43) As the proposed model depends on unobserved latent variables, we need to find a lower bound on likelihood of the observed data for any variational distribution Q(h, g, ρ) over the hidden variables.

(44) $\begin{array}{cc}ℒ(\ominus )=\log p\left(\mathrm{��}.Math.\ominus \right)\ge {\int }_h{\int }_g{\int }_{\rho }\mathrm{��}\left(h,g,\rho \right)\log \frac{p\left(\mathrm{��},h,g,\rho \right)}{\mathrm{��}\left(h,g,\rho \right)}\mathrm{dhdgdp}=E_{h,g,\rho }\left[\log p\left(\mathrm{��}.Math.h,g,\rho \right)\right]+E_h\left[\log \frac{p\left(h\right)}{\mathrm{��}\left(h\right)}\right]+E_g\left[\log \frac{p\left(g\right)}{\mathrm{��}\left(g\right)}\right]+E_p\left[\log \frac{p\left(\rho \right)}{\mathrm{��}\left(\rho \right)}\right]\stackrel{\Delta }{=}\widehat{ℒ}\left(\mathrm{��}\left(h\right),\mathrm{��}\left(g\right),\mathrm{��}\left(\rho \right),\ominus \right)& \left(3\right)\end{array}$

(45) By taking a separable Q(h, g, ρ), that is Q(h, g, ρ)=Q(h)Q(g)Q(ρ), the second equation can be deduced. The variational EM algorithm can be implemented as follows. 1) Variational E-Step: with fixing the current parameter value θ, Q-distributions of latent variables are calculated in E-Step. (i). Updating Q-Distribution of h.sub.i: From model (1), we can see that and g.sub.ij are constant with respect to h.sub.i; hence, we can rewrite the model (1) as

(46) $\underset{j=1}{\overset{N_i}{.Math.}}\left(Y_{\mathrm{ij}}-\mathrm{��}\times {}_{3}g_{\mathrm{ij}}^T\right)=N_iℱ\times {}_{3}h_i^T+\underset{j=1}{\overset{N_i}{.Math.}}{E}_{\mathrm{ij}}$

(47) Define

(48) ${\stackrel{\_}{Y}}_i=\frac{1}{N_i}\underset{j=1}{\overset{N_i}{.Math.}}\left(Y_{\mathrm{ij}}-W\times {}_{3}g_{\mathrm{ij}}\right)\mathrm{and}{\stackrel{\_}{E}}_i=\frac{1}{N_i}\underset{j=1}{\overset{N_i}{.Math.}}{E}_{\mathrm{ij}}.$
Then, the above-mentioned equation becomes
Y.sub.i=×.sub.3h.sub.i.sup.T+E.sub.i

(49) And each element in Ē.sub.i follows the normal distribution

(50) $\mathrm{��}\left(0,\frac{{\sigma }^2}{N_i}\right)$
with zero mean.

(51) In (3), the last two expectation terms are constant with respect to h.sub.i, hence we only need consider computing the first two expectations. The likelihood of observed data p(y|h, g, ρ) in (3) can be obtained using the following equation.

(52) $p\left(\mathrm{��}.Math.h,g,\rho \right)=\underset{i=1}{\overset{I}{.Math.}}\mathrm{��}\left({\stackrel{\_}{Y}}_i.Math.ℱ\times {}_{3}h_i^T,\frac{{\sigma }^2}{N_i}\right)$

(53) Therefore, the posterior of h.sub.i is still a multivariate normal distribution. The mean and covariance matrix are, respectively,

(54) 0${\stackrel{\_}{h}}_i=\left(\underset{ℱ}{.Math.}+\frac{1}{\stackrel{\_}{\rho }{N}_i}{I}_{K_1}\right)a_i$$\mathrm{and}$$\underset{h}{.Math.}={\left(I_{K1}+\frac{\stackrel{\_}{\rho }N}{I}{\Sigma }_{ℱ}\right)}^{-1}$

(55) where
a.sub.i=F.sup.(h)diag(F.sup.(1)TY.sub.iF.sup.(2))
and
=F.sup.(h)[(F.sup.(1)TF.sup.(1))⊙(F.sup.(2)TF.sup.(2))]F.sup.(h)T

(56) Where ⊙ represents the Hadamard elementwise product.

(57) ii) Updating Q-Distribution of g.sub.ij:

(58) we rewrite model (1) as,
Y.sub.ij−×.sub.3h.sub.i.sup.T=×.sub.3g.sub.ij.sup.T+E.sub.ij

(59) Define Ŷ.sub.ij=Y.sub.ij−×.sub.3 h.sub.i.sup.T. Then the above-mentioned model has become
Ŷ.sub.ij=×.sub.3g.sub.ij.sup.T+E.sub.ij

(60) Thus, the posterior of g.sub.ij is still a multivariate normal distribution. The mean and covariance matrix are, respectively,
g.sub.ij=(+1/ρI.sub.K).sup.−1b.sub.ij
and
Σ.sub.g=(I.sub.K+ρ).sup.−1
Where
b.sub.ij=W.sup.(g)diag(W.sup.(1)TŶ.sub.ijW.sup.(2))
and
Σ.sub.W=W.sup.(g)[(W.sup.(1)TW.sup.(1))⊙(W.sup.(2)TW.sup.(2))]W.sup.(g)T

(61) iii) Updating the posterior of ρ:

(62) utilizing Bayesian Theory, the log of optimal Q*(ρ) distribution can be calculated as expectation of the log joint distribution over all the hidden and visible variables with respect to all the other latent variables. Thus, we can get the following equation:

(63) $\begin{array}{c}\log Q^{*}\left(\rho \right)=E_{h,g}\left[\log p\left(\mathrm{��},h,g,\rho .Math.\ominus \right)\right]\\ =\frac{D_1{D}_{2}N}{2}\log \rho -\frac{\rho }{2}\psi +\left(a-1\right)\log \rho -b\rho \end{array}\hspace{1em}$

(64) where

(65) $\psi =\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{.Math.Y_{\mathrm{ij}}-ℱ\times {}_3\stackrel{\_}{h}_i^T-\mathrm{��}\times {}_3\stackrel{\_}{g}_{\mathrm{ij}}^T.Math.}_F^2+N\mathrm{tr}\left(\underset{ℱ}{.Math.}\underset{h}{.Math.}\right)+N\mathrm{tr}\left({\Sigma }_{\mathrm{��}}{\Sigma }_g\right).$

(66) Thus, Q-Distribution of ρ is proportional to

(67) ${\mathrm{��}}^{*}\left(\rho \right)\propto {\rho }^{a-1}+\frac{D_1{D}_{2}N}{2}\exp \left\{-b\rho -\frac{\psi }{2}\rho \right\}$

(68) Therefore, the best Q*(ρ) is actually a Gamma distribution, wherein the parameters become

(69) $\stackrel{\_}{a}=a+\frac{D_1{D}_{2}N}{2}\mathrm{and}\stackrel{\_}{b}=b+\frac{1}{2}\psi .$

(70) 2) Variational M-Step: With all the fixed Q-Distribution of latent variables, maximize all the parameters of (Q(h), Q(g), Q(ρ), θ). We select all terms containing parameters θ in (3) and get

(71) $\widehat{ℒ}\propto E_{h,g,\rho }\left[\log p\left(\mathrm{��},h,g,\rho \left[ℱ,\mathrm{��}\right)\right]\right]=\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{.Math.Y_{\mathrm{ij}}-ℱ\times {}_{3}h_i^T-\mathrm{��}\times {}_{3}g_{\mathrm{ij}}^T.Math.}_F^2+N\mathrm{tr}\left(\underset{ℱ}{.Math.}\underset{h}{.Math.}\right)+N\mathrm{tr}\left(\underset{\mathrm{��}}{.Math.}\underset{g}{.Math.}\right)+C\stackrel{\Delta }{=}\widetilde{ℒ}.$

(72) Where C represents a constant, which is independent with parameter θ. In addition, we can deduce
×.sub.3h.sub.i.sup.T=F.sup.(1)Diag(h.sub.i.sup.TF.sup.(h))F.sup.(2).sup.T
and
×.sub.3g.sub.ij.sup.T=W.sup.(1)Diag(g.sub.ij.sup.TW.sup.(g))W.sup.(2).sup.T.

(73) Thus, can be rewritten as

(74) $\widehat{ℒ}\propto \underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{.Math.{\tilde{Y}}_{\mathrm{ij}}-F^{\left(1\right)}\mathrm{Diag}\left(h_i^T{F}^{\left(h\right)}\right)F^{\left(2\right)T}.Math.}_F^2+N\mathrm{tr}\left(\left[\left(F^{\left(1\right)T}{F}^{\left(1\right)}\right)\odot \left(F^{\left(2\right)T}{F}^{\left(2\right)}\right)\right]F^{\left(h\right)T}\underset{h}{.Math.}{F}^{\left(h\right)}\right)$

(75) where
{tilde over (Y)}.sub.ij=Y.sub.ij−×.sub.3g.sub.ij.sup.T

(76) Taking derivatives of objective function with respect to F.sup.(1), F.sup.(2) and F.sup.(h), we have

(77) $F^{\left(1\right)}={\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\tilde{Y}}_{\mathrm{ij}}{F}^{\left(2\right)}\mathrm{Diag}\left(h_i^T{F}^{\left(h\right)}\right)\right]\left[\left(F^{\left(2\right)T}{F}^{\left(2\right)}\right)\odot \left({\mathrm{NF}}^{\left(h\right)T}\underset{h}{.Math.}{F}^{\left(h\right)}+F^{\left(h\right)T}\tilde{H}{\tilde{H}}^T{F}^{\left(h\right)}\right)\right]}^{-1},F^{\left(2\right)}={\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\tilde{Y}}_{\mathrm{ij}}^T{F}^{\left(1\right)}\mathrm{Diag}\left(h_i^T{F}^{\left(h\right)}\right)\right]\left[\left(F^{\left(1\right)T}{F}^{\left(1\right)}\right)\odot \left({\mathrm{NF}}^{\left(h\right)T}\underset{h}{.Math.}{F}^{\left(h\right)}+F^{\left(h\right)T}\tilde{H}{\tilde{H}}^T{F}^{\left(h\right)}\right)\right]}^{-1}\mathrm{and}$$F^{\left(h\right)}={{\left[\tilde{H}{\tilde{H}}^T+N\underset{h}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{h_i\left[d\mathrm{iag}\left(F^{\left(2\right)T}{\tilde{Y}}_{\mathrm{ij}}^T{F}^{\left(1\right)}\right)\right]}^T\right]\left[\left(F^{\left(1\right)T}{F}^{\left(1\right)}\right)\odot \left(F^{\left(2\right)T}{F}^{\left(2\right)}\right)\right]}^{-1}$

(78) Note that {tilde over (H)}∈.sup.K.sup.1.sup.×N is made up of all samples' h.sub.i and the columns corresponding to the same identity are equal. can be rewritten as

(79) $\widehat{ℒ}\propto \underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{.Math.{\tilde{Y}}_{\mathrm{ij}}-W^{\left(1\right)}\mathrm{Diag}\left(g_{\mathrm{ij}}^T{W}^{\left(g\right)}\right)W^{\left(2\right)T}.Math.}_F^2+N\mathrm{tr}\left(\left[\left(W^{\left(1\right)T}{W}^{\left(1\right)}\right)\odot \left(W^{\left(2\right)T}{W}^{\left(2\right)}\right)\right]W^{\left(g\right)T}{\Sigma }_g{W}^{\left(g\right)}\right)$

(80) Taking derivatives of objective function with respect to W.sup.(1), W.sup.(2) and W.sup.(h), and letting them be zero, we have

(81) $W^{\left(1\right)}={\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\hat{Y}}_{\mathrm{ij}}{W}^{\left(2\right)}\mathrm{Diag}\left(g_{\mathrm{ij}}^T{W}^{\left(g\right)}\right)\right]\left[\left(W^{\left(2\right)T}{W}^{\left(2\right)}\right)\odot \left({\mathrm{NW}}^{\left(g\right)T}\underset{g}{.Math.}{W}^{\left(g\right)}+W^{\left(g\right)T}{\mathrm{GG}}^T{W}^{\left(g\right)}\right)\right]}^{-1}$$W^{\left(2\right)}={\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\hat{Y}}_{\mathrm{ij}}^T{W}^{\left(1\right)}\mathrm{Diag}\left(g_{\mathrm{ij}}^T{W}^{\left(g\right)}\right)\right]\left[\left(W^{\left(1\right)T}{W}^{\left(1\right)}\right)\odot \left({\mathrm{NW}}^{\left(g\right)T}\underset{g}{.Math.}{W}^{\left(g\right)}+W^{\left(g\right)T}{\mathrm{GG}}^T{W}^{\left(g\right)}\right)\right]}^{-1}\mathrm{and}$$W^{\left(g\right)}={{\left[{\mathrm{GG}}^T+N\underset{g}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{g_{\mathrm{ij}}\left[\mathrm{diag}\left(W^{\left(2\right)T}{\hat{Y}}_{\mathrm{ij}}^T{W}^{\left(1\right)}\right)\right]}^T\right]\left[\left(W^{\left(1\right)T}{W}^{\left(1\right)}\right)\odot \left(W^{\left(2\right)T}{W}^{\left(2\right)}\right)\right]}^{-1}$

(82) 3) Overall Algorithm: The entire variational EM algorithm is composed of two steps. In E-Step, fixing all parameters in θ, hidden variables are updated as their Q-distributions. In M-Step, fixing all the hidden variable posteriors we update parameters. The two steps are conducted alternatively until a termination condition is satisfied.

(83) Define {tilde over (H)}.sub.t and G.sub.t as coefficient matrices in t-step. Thus, the reconstruction errors between the original and the reconstructed image sets can be represented as

(84) 0$e\left(t\right)=1-\frac{{.Math.\mathrm{��}-{ℱ}_t\times {}_3\tilde{H}_t^T-{\mathrm{��}}_t\times {}_{3}G_t^T.Math.}_F}{{.Math.\mathrm{��}.Math.}_F}$

(85) When the reconstruction errors are enough small or iterative number exceeds the defined value, the algorithm stops.

(86) High Dimensional Model Construction

(87) assuming a high-order observed data set {y.sub.ij|y.sub.ij∈.sup.D.sup.1.sup.× . . . ×D.sup.M; i=1, . . . , I, j=1, . . . , N.sub.i, Σ.sub.iN.sub.i=N}

(88) each sample is a tensor of order M, and the linear discriminant analysis model of higher-order tensor is
y.sub.ij=×.sub.(M+1)h.sub.i.sup.T+×.sub.(M+1)g.sub.ij.sup.T+ε.sub.ij  (4)

(89) Wherein h.sub.i∈.sup.K.sup.1, g.sub.ij∈.sup.K.sup.2, and is a tensor of order (M+1), whose size is respectively D.sub.1× . . . ×D.sub.M×K.sub.1 and D.sub.1× . . . ×D.sub.M×K.sub.2. Similar with the two order model, assuming base tensor has the structure of CP decomposition,

(90) $\begin{array}{c}ℱ=〚F^{\left(1\right)},F^{\left(2\right)},\mathrm{.Math.},F^{\left(M\right)},F^{\left(h\right)}〛\\ =\underset{r=1}{\overset{R_1}{.Math.}}{f}_{1,r}^{\left(1\right)}•f_{1,r}^{\left(2\right)}•\mathrm{.Math.}•f_{1,r}^{\left(M\right)}•f_{1,r}^{\left(h\right)}\end{array}\hspace{1em}\mathrm{and}\begin{array}{c}\mathrm{��}=〚W^{\left(1\right)},W^{\left(2\right)},\mathrm{.Math.},W^{\left(M\right)},W^{\left(g\right)}〛\\ =\underset{r=1}{\overset{R_2}{.Math.}}{w}_{1,r}^{\left(1\right)}•w_{1,r}^{\left(2\right)}•\mathrm{.Math.}•w_{1,r}^{\left(M\right)}•w_{1,r}^{\left(g\right)}.\end{array}\hspace{1em}$

(91) The Solution of High Dimensional Model Algorithm

(92) For higher-order model, EM algorithm is used to solve the problem, therefore the posterior of h.sub.i is still a multivariate normal distribution. The mean and covariance matrix are, respectively,

(93) $\begin{array}{c}{\stackrel{\_}{h}}_i=\left(\underset{ℱ}{.Math.}+\frac{1}{\stackrel{\_}{\rho }{N}_i}{I}_{K_1}\right)a_i\\ \underset{h}{.Math.}={\left(I_{K1}+\frac{\stackrel{\_}{\rho }N}{I}{.Math.}_{ℱ}\right)}^{-1}\end{array}$

(94) where,
=F.sup.(h)[(F.sup.(m)TF.sup.(m))⊙(F.sup.(m)TF.sup.(m))]F.sup.(h)T
a.sub.i=F.sup.(h)diag(F.sup.(m)TY.sub.(m).sup.iF.sup.(m)).

(95) Where ⊙ represents the Hadamard elementwise product.

(96) Define

(97) ${\stackrel{\_}{y}}_i=\frac{1}{N_i}\underset{j=1}{\overset{N_i}{.Math.}}\left(y_{\mathrm{ij}}-\mathrm{��}{\times }_{\left(M+1\right)}{g}_{\mathrm{ij}}\right),$
then Y.sub.(m).sup.i represents the mode-m metricized version of y.sub.i, and
F.sup.(m)=F.sup.(M) . . . F.sup.(m+1)F.sup.(m−1) . . . F.sup.(1)

(98) Where represents the Khatri-Rao product.

(99) By computing, the posterior of g.sub.ij is still a multivariate Gauss distribution, and the mean and covariance of latent variable g.sub.ij have become
g.sub.ij=(+1/ρI.sub.K).sup.−1b.sub.ij
Σ.sub.g=(I.sub.K+ρ).sup.−1
where
b.sub.ij=W.sup.(g)diag(W.sup.(m)TŶ.sub.(m).sup.ijW.sup.(m))
=W.sup.(g)[(W.sup.(m)TW.sup.(m))⊙(W.sup.(m)TW.sup.(m))]W.sup.(g)T

(100) Let ŷ.sub.ij=y.sub.ij−×.sub.(M+1)h.sub.i.sup.T, and Ŷ.sub.(m).sup.ij represents the mode-m metricized version of ŷ.sub.ij. And
W.sup.(m)=W.sup.(M) . . . W.sup.(m+1)W.sup.(m−1) . . . W.sup.(1)

(101) The optimal Q*(ρ) is still a Gamma distribution.

(102) In M-step, taking derivatives of with respect to F.sup.(n) (n=1, . . . , N) and F.sup.(h)and letting them be zero, obtain

(103) $F^{\left(m\right)}=\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\tilde{Y}}_{\left(m\right)}^{\mathrm{ij}}{\stackrel{\_}{F}}^{\left(m\right)}\mathrm{Diag}\left(h_i^T{F}^{\left(h\right)}\right)\right]\hspace{1em}{\left[\left({\stackrel{\_}{F}}^{\left(m\right)T}{\stackrel{\_}{F}}^{\left(m\right)}\right)\odot \left({\mathrm{NF}}^{\left(h\right)T}\underset{h}{.Math.}{F}^{\left(h\right)}+F^{\left(h\right)T}\tilde{H}{\tilde{H}}^T{F}^{\left(h\right)}\right)\right]}^{-1}{F}^{\left(h\right)}={\left[\tilde{H}{\tilde{H}}^T+N\underset{h}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{h_i\left[\mathrm{diag}\left({\stackrel{\_}{F}}^{\left(m\right)T}{\tilde{Y}}_{\left(m\right)}^{\left(\mathrm{ij}\right)T}{F}^{\left(m\right)}\right)\right]}^T\right]\hspace{1em}\left[\left(F^{\left(m\right)T}{F}^{\left(m\right)}\right)\odot \left({\stackrel{\_}{F}}^{\left(m\right)T}{\stackrel{\_}{F}}^{\left(m\right)}\right)\right]\text{-}1.$

(104) In which {tilde over (Y)}.sub.(m).sup.ij presents the mode-m metricized version of {tilde over (y)}.sub.ij. And {tilde over (y)}.sub.ij=y.sub.ij−×.sub.(M+1)g.sub.ij.

(105) Taking derivatives of with respect to W.sup.(m) (m=1, . . . , M) and W.sup.(h), and letting them be zero, obtain

(106) $\begin{array}{cc}{W}^{\left(m\right)}=\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{\hat{Y}}_{\left(m\right)}^{\mathrm{ij}}{\stackrel{\_}{W}}^{\left(m\right)}\mathrm{Diag}\left(g_{\mathrm{ij}}^T{W}^{\left(g\right)}\right)\right]\hspace{1em}{\left[\left({\stackrel{\_}{W}}^{\left(m\right)T}{\stackrel{\_}{W}}^{\left(m\right)}\right)\odot \left({\mathrm{NW}}^{\left(g\right)T}\underset{g}{.Math.}{W}^{\left(g\right)}+W^{\left(g\right)T}{\mathrm{GG}}^T{W}^{\left(g\right)}\right)\right]}^{-1}& \\ \mathrm{and}& \\ W^{\left(h\right)}={\left[{\mathrm{GG}}^T+N\underset{g}{.Math.}\right]}^{-1}\left[\underset{i=1}{\overset{I}{.Math.}}\underset{j=1}{\overset{N_i}{.Math.}}{g_{\mathrm{ij}}\left[\mathrm{diag}\left({\stackrel{\_}{W}}^{\left(m\right)T}{\hat{Y}}_{\left(m\right)}^{\left(\mathrm{ij}\right)T}{W}^{\left(m\right)}\right)\right]}^T\right]\hspace{1em}{\left[\left({\stackrel{\_}{W}}^{\left(m\right)T}{W}^{\left(m\right)}\right)\odot \left({\stackrel{\_}{W}}^{\left(m\right)T}{W}^{\left(m\right)}\right)\right]}^{-1}.& \end{array}$