Unsupervised feature selection using Bayesian Tucker decomposition Y-h. Taguchi1*† and Yoh-ichi Mototake2†
arXiv:2604.14949v1 [stat.ML] 16 Apr 2026
1* Department of Physics, Chuo University, Tokyo, 112-8551, Japan. 2 Graduate School of Social Data Science, Hitotsubashi University,
Tokyo, 186-8601, Japan.
*Corresponding author(s). E-mail(s): [email protected]; Contributing authors: [email protected]; † These authors contributed equally to this work. Abstract In this paper, we proposed Bayesian Tucker decomposition (BTuD) in which residual is supposed to obey Gaussian distribution analogous to linear regression. Although we have proposed an algorithm to perform the proposed BTuD, the conventional higher-order orthogonal iteration can generate Tucker decomposition consistent with the present implementation. Using the proposed BTuD, we can perform unsupervised feature selection successfully applied to various synthetic datasets, global coupled maps with randomized coupling strength, and gene expression profiles. Thus we can conclude that our newly proposed unsupervised feature selection method is promising. In addition to this, BTuD based unsupervised FE is expected to coincide with TD based unsupervised FE that were previously proposed and successfully applied to a wide range of problems. Keywords: tensor decomposition, feature selection, unsupervised learning
1 Introduction The purpose of this paper is to propose a Bayesian unsupervised feature selection. There is a way to perform unsupervised feature selection using tensor decomposition (TD) by TD based unsupervised feature extraction (FE) [1]. In this paper, we develop Bayesian Tucker decomposition (BTuD) to perform TD based unsupervised FE in the framework of Bayesian statistics.
1
TD is one of the old-fashioned embedding methods to process high dimensional data. It was invented long before the machine learning has become popular. In spite of that, it still remains one of the frequently used methods in data science. It is natural to modernize TD alongside the recent development. The proposal of Bayesian TD [2] is one of such effort. Nevertheless, most of those proposals are restricted to one of the popular variants of TD, CP decomposition. Besides Cheng et al’s book [2], although there were many proposals of Bayesian TD [3–5], most of them deal with only CP decomposition. Another frequently used TD, Tucker decomposition, was not frequently discussed along the framework of Bayes statistics [6–8]. The critical difference between these a few previous implementations and the one proposed in this study is that the previous methods require that the decomposed components themselves derived from tensor decomposition obey Gaussian distribution whereas the present implementation only requires that residuals follow a Gaussian distribution. The reason we extend Bayesian statistic to Tucker decomposition is that we would like to utilize Bayesian TD for PCA/TD based unsupervised FE [1] that we have proposed a long time ago. Especially, in TD based unsupervised FE, we definitely need Tucker decomposition since we have compared various TDs to establish TD based unsupervised FE and found that Tucker decomposition is the best one for TD based unsupervised FE. The common statistical basis of the above conventional Bayesian approaches is to assume a Gaussian distribution to components themselves following the tradition of probabilistic principal component analysis (PCA) [9], which is a bit problematic for our purposes. In PCA/TD based unsupervised feature extraction, feature selection is performed by selecting outliers under the null hypothesis that components are assumed to follow a Gaussian distribution. Thus, since there can be no outliers when assuming that components obey Gaussian as in the other implementations of other traditional Bayesian TD, the implementation where components are assumed to be Gaussian is not adequate for our purpose. To address this problem, instead of assuming the components themselves as Gaussian, employing the framework of linear regression, we assume that not the components themselves but the residuals obey Gaussian. Since the implementation of linear regression using Bayes statistics is well known [10], to get a representation of BTuD within this strategy is straightforward.
2 Previous feature selection approaches Most previous feature selection approaches are supervised ones. There are two branches of supervised feature selection. One is the filter approach and the other is the wrapper approach. In the filter approaches, individual features are evaluated based upon provided labels/classification. One of typical filter approaches is t test. In the t-test, P -values are attributed to individual features with assuming null hypothesis that all components of individual features are drawn from the same distribution and features that can reject the null hypothesis are selected. One of wrapper approaches is random forest [11]. In random forest, individual features are randomly selected and replaced until the desired performance is achieved. Anyway, in the supervised approaches, we need some (external) supervised criteria by which we can evaluate individual features. 2
Unsupervised approach is completely different from the supervised approaches. In unsupervised approaches, we have to have criteria by which we can evaluate individual features without any external information. One typical unsupervised approache is unimodality testing [12]. In unimodality testing, components in an individual feature are assumed to be unimodal (i.e., the distribution does not have multiple peaks but have only one peak). Although unimodality testing is powerful, it has limited ability to select features since features associated with single-peaked distribution cannot be selected even if it has some order. Thus, we need more powerful unsupervised approaches. PCA/TD based unsupervised FE is such a method, since the criterion for selecting features is independent of the evaluation of features (see Methods). Once singular value vectors are computed with attributing some of them to samples and others to features. Evaluation is based upon those attributed to samples and selection is based upon those attributed to features. Since the evaluation and selection is separated, we can select any kind of features coincident with external labels. Thus to perform TD based unsupervised FE in Bayesian statistics is promising approach and worthwhile trying.
3 Methods 3.1 Implementation of Bayesian Tucker decomposition with a framework of linear regression Assume that we have Tucker decomposition of third-order tensor, xijk ∈ RN ×M ×K as xijk ≃
L3 L2 X L1 X X
G(ℓ1 ℓ2 ℓ3 )uℓ1 i uℓ2 j uℓ3 k
(1)
ℓ1 =1 ℓ2 =1 ℓ3 =1
where 1 ≤ L1 ≤ N, 1 ≤ L2 ≤ M, 1 ≤ L3 ≤ K and G ∈ RL1 ×L2 ×L3 is a core tensor that represents the weight of the product uℓ1 i uℓ2 j uℓ3 k toward xijk , uℓ1 i ∈ RL1 ×N , uℓ2 j ∈ RL2 ×M , and uℓ3 k ∈ RL3 ×K are singular value matrices and orthogonal matrices. Top row of Fig. 1 shows the graph diagram to show the structure of Bayesian inference consistent with eq. (1). The likelihood that corresponds to the top row of Fig. 1 (or eq. (1)) is written as p(xijk |uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 ))
(2)
p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 ))
(3)
and prior is written as Thus the posterior is p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )|xijk ) =
p(xijk |uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 ))p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )) p(xijk ) (4)
3
u
l
x
1 i
u
l
2
u
ijk
j
l k 3
(
G l l l
(
1
2
∣u
p x
u
l
x
1 i
u
l
2
3
ijk
) u
l i 1
l
2
u
j
u
l k
(
G l l l
3
1
l
x
2 j
ijk
u
j
2
Φ
l k 3
(
G l l l 1
2
3
,l
(
G l l l
u
ik
,l
1
2
3
2
u
Φ
l
1
Φ l
2
(
x
j
ijk
l i
j
G l l l
l
l
u 2
)
1 i
u
x
1
Φ
l i
1
)
l
3k
u
1
jk
) u
l k
u
)
ijk
3
u
3
l
2
l
3
1
2
3
ij , l
3
)
ijk
, ijk
l k 3
(
G l l l 1
2
3
)
Fig. 1 Graphical diagram that shows the structure of Bayesian inference. Top row is consistent with Eq. (1). The middle row is consistent with eqs. (6), (8), and (10) (From left to right panels in the middle row). The bottom row is consistent with eq. (12). The cyclic computation among the middle and bottom rows is the approximation to the first row.
.
4
3.1.1 Bayesian Tucker decomposition interpreted as a set of linear regression Unfortunately, since it is not very practical to infer Eq. (1) directly, so we decompose it to a set of linear regressions as follows. To apply Bayesian linear regression to Tucker decomposition we re-write it as
X
G(ℓ1 ℓ2 ℓ3 )uℓ2 j uℓ3 k ∈ RM K×L1
(5)
Φjk,ℓ1 uℓ1 i
(6)
G(ℓ1 ℓ2 ℓ3 )uℓ1 i uℓ3 k ∈ RN K×L2
(7)
Φik,ℓ2 uℓ2 j
(8)
G(ℓ1 ℓ2 ℓ3 )uℓ1 i uℓ2 j ∈ RN M ×L3
(9)
Φij,ℓ3 uℓ3 k
(10)
Φℓ1 ℓ2 ℓ3 ,ijk ≡ uℓ1 i uℓ2 j uℓ3 k ∈ RL1 L2 L3 ×N M K X xijk ≃ Φℓ1 ℓ2 ℓ3 ,ijk G(ℓ1 ℓ2 ℓ3 )
(11)
Φjk,ℓ1 ≡
ℓ2 ,ℓ3
xijk ≃
X ℓ1
Φik,ℓ2 ≡
X ℓ1 ,ℓ3
xijk ≃
X ℓ2
Φij,ℓ3 ≡
X ℓ1 ,ℓ2
xijk ≃
X ℓ3
(12)
ℓ1 ℓ2 ℓ3
where Eqs. (6), (8), (10), and (12) are regarded as linear regression problems where uℓ1 i , uℓ2 j , uℓ3 k , and G(ℓ1 ℓ2 ℓ3 ) are the regression coefficients with the given Φjk,ℓ1 , Φik,ℓ2 , Φij,ℓ3 , and Φℓ1 ℓ2 ℓ3 ,ijk . Actual procedures are as follows. When applying linear regression to eq. (6), at first, i is fixed, uℓ1 i is assumed to be the set of regression coefficients, with L1 coefficients, and Φjk,l1 is supposed to be given (i.e., uℓ2 j and uℓ3 k are supposed to be known). Thus, there are L1 variables whereas the number of equations (conditions) is M K. This means that L1 can be larger than M K. To avoid the difficulty in this case, we employed Moore-Penrose pseudoinverse [13] matrix to perform linear regression (see below). This is also repeated toward uℓ2 j (eq. (8) or the middle panel of the middle row of Fig. 1) and uℓ3 k (eq. (10) or the right panel of the middle row of Fig. 1) with cycling indices among i, j, k and ℓ1 , ℓ2 , ℓ3 . The graphical diagram of this procedure is in the middle of Fig. 1. After that, the linear regression that corresponds to eq. (12) (or the bottom row of Fig. 1) is performed. The whole process that includes these four linear regression is an approximation to eq. (1) (or top row of Fig. 1) This means that MAP , uMAP , GMAP (ℓ ℓ ℓ )|x ) uMAP 1 2 3 ijk ℓa ia = arg max p(uℓa ia , uℓb ib ℓc ic uℓa ia
5
(13)
GMAP (ℓ1 ℓ2 ℓ3 ) = arg max p(uMAP , uMAP , uMAP , G(ℓ1 ℓ2 ℓ3 )|xijk ) ℓ1 i ℓ2 j ℓ3 k
(14)
G(ℓ1 ℓ2 ℓ3 )
where (ℓa , ℓb , ℓc ) = (ℓ1 , ℓ2 , ℓ3 ) and (ia , ib , ic ) = (i, j, k) in eq. (14) should be cyclically changed. Although likelihood, eq. (2), does not change, prior, Eq. (3), that corresponds to eqs. (6), (8), and (10) (or the left to right panels of the middle row of Fig. 1) are re-written as p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )) = p(uℓ1 i )p(uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )) MAP MAP = p(uℓa ia )δ(uℓb ib − uMAP (ℓ1 ℓ2 ℓ3 )) ℓb ib )δ(uℓc ic − uℓc ic )δ(G(ℓ1 ℓ2 ℓ3 ) − G
(15)
where (ℓa , ℓb , ℓc ) = (ℓ1 , ℓ2 , ℓ3 ) and (ia , ib , ic ) = (i, j, k) should be cyclically changed as well. and p(uℓa ia ), p(uℓb ib ), and p(uℓc ic ) are the uniform distributions. Thus, posterior, eqs. (4), also re-written as p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )|xijk ) p(xijk |uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 ))p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )) = p(xijk ) MAP MAP p(xijk |uℓa ia , uMAP , G , u (ℓ1 ℓ2 ℓ3 ))p(uℓa ia ) ℓb ib ℓc i c = p(xijk )
(16)
The optimization process for eqs. (6), (8), and (10) (or the left to right panels of the middle row of Fig. 1) is illustrated in Fig. 2. In summary, we decomposed Tucker decomposition eq. (1) into four linear regression subproblems, each of which estimates one of uℓ1 i , uℓ2 j , uℓ3 k , and G(ℓ1 ℓ2 ℓ3 ) as MAP solution. In each subproblem, the remaining three variables are treated as constants fixed at their current MAP estimates. Since this alternating optimization updates one subproblem at a time, it is not guaranteed to reach the global minimum and may instead converge to a local minimum (Fig. 2). Nevertheless, the obtained solution is expected to be close to the global minimum unless the optimization is trapped in a poor local minimum.
3.2 Likelihood Let xi ∈ RM K be the vectorization of {xijk }j,k , and let Φ ∈ RM K×L1 be the design matrix with Φ(j,k),ℓ1 = Φjk,ℓ1 in Eq.(2). Then Eq.(3) is written as εi ∼ N (0, β −1 I),
xi = Φ ui + εi , which yields the likelihood
p(xi | ui , Φ, β) = N xi | Φui , β −1 I .
6
u
l
b
i
,u b
l
c
i
,G
(l l l )
c
1
2
3
(
p u
l i
,u
1
u u
G
MAP
l i 1
,u
l i
(
,G l l l
1
1
2
3
)∣x ) ijk
MAP l
b
i
b
MAP l
c
i
c
(l l l ) 1
2
3
u
MAP l
a
i
u
l i a
a
a
Fig. 2 The process that optimizes posterior, p(uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 ℓ2 ℓ3 )|xijk ), which is repreMAP MAP (ℓ1 ℓ2 ℓ3 ) where (ℓa , ℓb , ℓc ) = sented by red contour, towards uℓa ia with fixed values uMAP ℓb ib , uℓc ic , G (ℓ1 , ℓ2 , ℓ3 ) and (ia , ib , ic ) = (i, j, k) should be cyclically changed. This does not clearly give us global minimum but local minimum.
Equivalently, for the whole tensor, p({xijk } | {uℓ1 i }, {uℓ2 j }, {uℓ3 k }, G, β) =
Y
X N xijk | Gℓ1 ℓ2 ℓ3 uℓ1 i uℓ2 j uℓ3 k , β −1 .
i,j,k
ℓ1 ,ℓ2 ,ℓ3
3.3 Prior probability We place independent zero-mean Gaussian priors on the regression coefficients: p(ui | α) = N (ui | 0, α−1 I),
(2)
p(uj
(2)
| α) = N (uj
| 0, α−1 I),
(3)
(3)
p(uk | α) = N (uk | 0, α−1 I),
and similarly for the vectorized core tensor g = vec(G): p(g | α) = N (g | 0, α−1 I).
3.4 Posterior Probability The posterior probability of uℓ1 i , uℓ2 j , uℓ3 k , G(ℓ1 , ℓ2 , ℓ3 ) are given as follows.
P {uℓa ia }ℓa |xijk , {G}ℓb ,ℓc , {muℓb ib }ℓb , {muℓc ic }ℓc = N {uℓa ia }ℓa |{muℓa ia }ℓa , {Sℓa ℓ′a }ℓa ,ℓ′a (17)
7
(P
Φ†ib ic ,ℓa xijk (α = 0) ∈ RL1 ×N ′ Φi i ,ℓ′ xijk (α ̸= 0) βS ′ ℓ ℓ a c b ib ic ℓa a a X −1 Sℓa ℓ′ = αI + β Φib ic ,ℓa Φib ic ,ℓ′a ∈ RL1 ×L1
Pib ic
muℓa ia =
(18) (19)
a
ib ic
where (ℓa , ℓb , ℓc ) = (ℓ1 , ℓ2 , ℓ3 ) and (ia , ib , ic ) = (i, j, k) should be cyclically changed.
P {G}ℓ |xijk , {muℓ1 i }ℓ1 , {muℓ2 j }ℓ2 , {muℓ3 k }ℓ3 = N ({G}ℓ |{mG }ℓ , {Sℓℓ′ }ℓℓ′ ) (P mG =
Φ†ℓ1 ℓ2 ℓ3 ,ijk xijk (α = 0) ∈ RL1 L2 L3 ×L1 L2 L3 ′ ijkℓ′ βSℓℓ Φℓ,ijk xijk (α ̸= 0)
Pijk
(20)
(21) (22)
−1 Sℓℓ ′ = αI + β
X
L1 L2 L3 ×L1 L2 L3
Φℓ,ijk Φℓ′ ,ijk ∈ R
(23)
ijk
where, for example, {· · · }ℓ1 means that it is composed of 1 ≤ ℓ1 ≤ L1 and so on whereas ℓ = (ℓ1 , ℓ2 , ℓ3 ) and ℓ′ = (ℓ′1 , ℓ′2 , ℓ′3 ). More explicitly, (
P {uℓa ia }ℓa |xijk , {G}ℓb ,ℓc , {muℓb ib }ℓb , {muℓc ic }ℓc =
e
− 21
PL1
u −muℓ i ℓ ,ℓ′ =1 ( ℓa ia a
a a
a
i h ) Sℓ−1,ℓ′ a
q
a ℓa ℓ′a
(2π)La Sℓa ,ℓ′a (24) where (ℓa , ℓb , ℓc ) = (ℓ1 , ℓ2 , ℓ3 ) and (ia , ib , ic ) = (i, j, k) should be cyclically changed and La should also be changed as (L1 , L2 , L3 ) accordingly.
3.5 Algorithm As mentioned above, uℓ1 i , uℓ2 j , and uℓ3 k cannot be fitted simultaneously but must be fitted one by one. To do this, we proposed the following procedure (Algorithm 1) where A† = (AT A)−1 AT is Moore–Penrose inverse. α is the multiplier of l2 -norm (do not confuse l2 of l2 -norm with index ℓ2 in G(ℓ1 ℓ2 ℓ3 ) or uℓ2 j ). The outline of Algorithm 1 is as follows. At first, uℓ1 i is tried to be computed with eq. (6). Orthogonalization as well as normalization are applied to uℓ1 i . G is tried to be computed with eq. (12). Then, uℓ2 j is tried to be computed with eq. (8). Orthogonalization as well as normalization are applied to uℓ2 j . G is tried to be computed with eq. (12). Finally, uℓ3 k is tried to be computed with eq. (10). Orthogonalization as well as normalization are applied to uℓ3 k . G is tried to be computed with eq. (12). These whole processes are repeated until all of uℓ1 i , uℓ2 j , uℓ3 k , and G are converged. The critical difference between HOOI and the Algorithm 1 is that the former employs SVD whereas the latter employs linear regression. In the sense that they both employ the scheme of alternative least squares, they are the same as others.
8
)
uℓ′ ia −mu ′ a
ℓa ia
Algorithm 1 Bayesian Tucker decomposition with linear regression Require: G, uℓ1 i , uℓ2 j , uℓ3 k computed by HOOI. 1: while not converged do 2: while ℓ1 ≤PL1 do [Solve eq. (6) ] Φjk,l1 ← ℓ2 ,ℓ3 G(ℓ1 ℓ2 ℓ3 )uℓ2 j uℓ3 k 3: while i ≤ N do P 4: if (α = 0) then uℓ1 i ← jk Φ†jk,ℓ1 xijk † P P P ′ ′ 5: else uℓ1 i ← ℓ′ jk Φjk,ℓ1 Φjk,ℓ1 + αI jk Φjk,ℓ1 xijk 1
6:
7:
8: 9:
10: 11: 12:
end if i←i+1 end while uℓ1 i ← uℓ1 i ⊥ uℓ′1 ̸=ℓ1 i [Orthogonalization] qP 2 uℓ1 i ← uℓ1 i / i uℓ1 i [Normalization] Φℓ1 ℓ2 ℓ3 ,ijk ← uℓ1 i uℓ2 j uℓ3 k P G(ℓ1 ℓ2 ℓ3 ) ← ijk Φ†ℓ1 ℓ2 ℓ3 ,ijk xijk [Solve Eq. (12)] ℓ1 ← ℓ1 + 1 end while while ℓ2 ≤PL2 do [Solve eq. (8) ] Φik,l2 ← ℓ1 ,ℓ3 G(ℓ1 ℓ2 ℓ3 )uℓ1 i uℓ3 k while j ≤ M do P if (α = 0) then uℓ2 j ← ik Φ†ik,ℓ2 xijk h † i P P ′ + αI else uℓ2 j ← ℓ′ Φ Φ ik,ℓ ik,ℓ 2 ik 2 2
13: 14:
15: 16:
17: 18: 19:
P ℓ2 ℓ′2
end if j ← j + 1 end while uℓ2 j ← uℓ2 j ⊥ uℓ′2 ̸=ℓ2 j [Orthogonalization] qP 2 uℓ2 j ← uℓ2 j / j uℓ2 j [Normalization] Φℓ1 ℓ2 ℓ3 ,ijk ← uℓ1 i uℓ2 j uℓ3 k P G(ℓ1 ℓ2 ℓ3 ) ← ijk Φ†ℓ1 ℓ2 ℓ3 ,ijk xijk [Solve Eq. (12)] ℓ2 ← ℓ2 + 1 end while while ℓ3 ≤ PL3 do [Solve Eq. (10)] Φij,l3 ← ℓ1 ,ℓ2 G(ℓ1 ℓ2 ℓ3 )uℓ1 i uℓ2 j while k ≤ K do P if (α = 0) then uℓ3 k ← ij Φ†ij,ℓ3 xijk † P P ′ else uℓ3 k ← ℓ′ Φ Φ + αI ij ij,ℓ3 ij,ℓ3 3
end if k ←k+1 21: end while uℓ3 k ← uℓ3 k ⊥ uℓ′3 ̸=ℓ3 k [Orthogonalization] qP 2 uℓ3 k ← uℓ3 k / k uℓ3 k [Normalization] Φℓ1 ℓ2 ℓ3 ,ijk ← uℓ1 i uℓ2 j uℓ3 k P G(ℓ1 ℓ2 ℓ3 ) ← ijk Φ†ℓ1 ℓ2 ℓ3 ,ijk xijk [Solve Eq. (12)] 9 ℓ3 ← ℓ3 + 1 22: end while 23: end while 20:
ℓ1 ℓ′1
ik Φik,ℓ2 xijk
P
′
ℓ3 ℓ′3
ij Φij,ℓ3 xijk ′
3.6 Replacement with HOOI Although we have proposed the algorithm, Algorithm 1, to perform BTuD, we occasionally found that HOOI can also give us the solution coincident with the proposed BTuD when α = 0. Empirically, if uℓ1 i is equal to muℓ1 i , uℓ2 j is equal to muℓ2 j , uℓ3 k is equal to muℓ3 k , and G(ℓ1 ℓ2 ℓ3 ) is equal to mG , they can be regarded as solutions for BTuD. Thus, when the solution computed by HOOI is coincident with muℓ1 i , muℓ2 j , muℓ3 k and mG , we do not execute Algorithm 1, but employ the solution by HOOI as it is. This means that, with starting initial values computed by HOOI, Algorithm 1 converges immediately. Also one should notice that this does not always mean the solution given by HOOI is identical to that given by Algorithm 1, since there are no uniqueness about the solutions coincident with the proposed BTM. In this study, we always employed the solution by HOOI, since we must anyway execute HOOI at the very first stage in Algorithm 1 and always found that it is the solution coincident with muℓ1 i , muℓ2 j , muℓ3 k and mG (in this study, we considered only α = 0). For HOOI, we employed tucker function in rTensor [14] package in R.
3.7 TD based unsupervised FE Before explaining how we can make use of the above BTuD for feature selection, we introduce the previously proposed TD based unsupervised FE [13] since how to select features is quite similar. Suppose that j and k are attributed to samples and i is attributed to features. For example, j represents the distinction between genders (e.g., M K j≤ M 2 is female and j > 2 is male) and k represents age (e.g., k ≤ 2 is young and K k > 2 is old) and is are properties of persons (e.g., income, weight, height, blood pressure and so on). Then xijk is the value of ith property of the person labeled by j, k. After getting TD, we can find some uℓ2 j which is distinct between genders and uℓ3 k which is distinct between ages. Next, by investigating |G(ℓ1 ℓ2 ℓ3 )| we can find which ℓ1 s have larger contribution. Assuming uℓ1 i selected by the ℓ2 and ℓ3 obey Gaussian (Null hypothesis), we can attribute P -values to is as
" Pi = Pχ2 >
X uℓ i 2 1
ℓ1
σ ℓ1
# (25)
where the summation is taken over only the selected ℓ1 s, Pχ2 [> x] is the cumulative χ2 distribution and σℓ1 is the standard deviation optimized such that the distribution of the resulting P -values is coincident with the Gaussian distribution as much as follows as shown in the below. P -values are corrected by BH (Benjamini-Hochberg) criterion and is associated with adjusted P -values less than the threshold value, say 0.05, are selected. In order to optimize SD, first we compute histogram of 1 − Pi as hs (1 − P ) that is the frequency of is in the sth bin excluding features to be selected, i.e., those associated with adjusted P -values less than the threshold value, say 0.01. Then we
10
compute standard deviation of hs as S
1X hs S s=1 v u S u1 X 2 σh = t (hs − ⟨hs ⟩) S s=1
⟨hs ⟩ =
(26)
(27)
where S is the total number of bins. Since P -values are computed as a function of σℓ1 , there should be the optimal σℓ1 that gives the minimum σh supposed to be most coincident with the null hypothesis that uℓ1 i obeys Gaussian since hs (1 − Pi ) should ′ be flat (i.e., hs (1 − Pi ) = NS where N ′ is the total number of the not selected features) if the null hypothesis is totally true. When we analyze not a tensor but a matrix (i.e., sinusoidal data and RCS-GCM in the below), we can use singular value decomposition (SVD) xij =
X
λℓ uℓi uℓj
(28)
ℓ
instead of TD. Since uℓi is always associated with uℓj in SVD, we can decide which uℓi should be used for the feature selection just after we decide which uℓj is of interest.
3.8 BTuD based unsupervised FE Up to the stage at which we decide which ℓ1 to use to select is, the processes are identical to those in TD based unsupervised FE. Then, P -values that uℓ1 i > (<)0 when muℓ1 i < (>)0 are attributed to i as
" Pi = Pχ2 >
X m2uℓ i 1
ℓ1
Sℓ1 ℓ1
# .
(29)
The summation is taken over only the selected ℓ1 s, P -values are corrected by BH criterion and is associated with adjusted P -values less then the threshold value, say 0.05, are selected.
3.9 Data sets 3.9.1 Synthetic data To test whether the results of BTuD can be used for feature selection, we prepared a tensor, xijk ∈ RN ×M ×K , as follows
xijk ∼
K N (µ, 1) i ≤ N1 ≤ N , j ≤ M 2 ,k ≤ 2 N (0, 1) otherwise
11
(30)
by which we can test whether BTuD based unsupervised FE can select i ≤ N1 correctly or not.
3.9.2 Sinusoidal data To test if the results of BTuD is used for feature selection, we also prepared a matrix, xij ∈ RN ×M , as the following
xij ∼
sin 2π 3j + ϵi , ϵi ∈ N (0, 1), i ≤ N1 ≤ N N (0, 1), N1 < i ≤ N .
(31)
by which we can test if BTuD based unsupervised FE can select i ≤ N1 correctly or not.
3.9.3 RCS-GCM GCM [15] is a globally coupled system of multiple chaotic systems xij+1 = f (xij , a), where j is regarded as time step, formulated as follows. N
xij+1 = (1 − g)f (xij , a) +
g X f (xi′ j , a) N ′
(32)
i =1
f (x, a) = 1 − ax2
(33)
where g is the coupling parameter and a is the non-linearity parameter, the increase of which results in the chaotic behaviour of xij ∈ RN ×M . By adjusting the parameters a, g of this system, we can generate a globally synchronized state, a partially synchronized state with correlations between parts of dimensions, and a chaotic state in which all dimensions create an independent turbulent state. In this study, we assume for the data-generating model of this world that all dimensions have different correlation lengths. The original GCM cannot achieve such a system. To generate variables that are a mixture of long correlated ordered states and short correlated random states, we introduce GCM with random parameters as N
xij+1 = gii f (xij , ai ) +
1 X gii′ f (xi′ j , ai′ ), N ′
(34)
i =1
gii′ = (1 − c)δii′ + cϵii′ ,
(35)
ai = a + (1 − a)ϵi ,
(36)
2
f (x, a) = 1 − ax ,
(37)
where ϵii′ and ϵi are uniform random numbers as ϵi , ϵii′ ∼ [0, 1]. This model extension has realized the behavior of the data generation model considered in this study, in which some dimensional groups behave synchronously with finite correlation lengths, whereas others behave chaotically with small correlation lengths. The specific parameters used in this study to generate data were a = 1.75, c = 0.04,
12
and N = 104 , such that a single f (x, a) falls in the chaotic region (a > 1.48). c, which expresses the strength of pairwise interactions between individual maps, is taken to be sufficiently small not to suppress the chaotic nature completely because of synchronization among individual maps and for f (x, a) to have the mixture of ordered and random states. js are taken to be 1 ≤ j ≤ 102 . Thus, the generated dataset is 4 2 xijk ∈ R10 ×10 . Initial values (xi0 ) are drawn from the same uniform distribution, [0, 1].
3.9.4 Gene expression The gene expression profiles used in this study were downloaded from the Gene Expression Omnibus (GEO) with GEO ID GSE142068. Twenty-four profiles named “GSE142068 count XXXXX.txt.gz” were downloaded, where “XXXXX” indicates one of the 24 tissues, i.e., AdrenalG, Aorta, BM (Bone marrow), Brain, Colon, Eye, Heart, Ileum, Jejunum, Kidney, Liver, Lung, Pancreas, ParotidG, PituitaryG, SkMuscle, Skin, Skull, Spleen, Stomach, Testis, Thymus, ThyroidG, and WAT (white adipose tissue), which were treated with 15 drugs: Alendronate, Acetaminophen, Aripiprazole, Asenapine, Cisplatin, Clozapine, Empagliflozin, Lenalidomide, Lurasidone, Olanzapine, Evolocumab, Risedronate, Sofosbuvir, and Teriparatide, and Wild type (WT). They were formatted as tensor, xijkm ∈ RN ×24×18×2 , for N genes, 24 tissues, 18 drug treatments, and two replicates. By applying HOSVD (higher order singular value decomposition) to xijkm , we get xijkm =
X
G(ℓ1 ℓ2 ℓ3 ℓ4 )uℓ1 j uℓ2 k uℓ3 m uℓ4 i
(38)
ℓ1 ℓ2 ℓ3 ℓ4
where G ∈ RN ×24×18×2 is the core tensor, uℓ1 j ∈ R24×24 , uℓ2 k ∈ R18×18 ,uℓ3 m ∈ R2×2 , and uℓ4 i ∈ RN ×N , represents singular valueP matrices that are also P orthogonal matrices. xijkm is considered to be standardized as i xijkm = 0 and i x2ijkm = N .
4 Results 4.1 Performance test of BTuD based unsupervised FE Since the proposal to interpret Tucker decomposition as sets of linear regression is new, to determine whether BTuD based unsupervised FE works correctly, we apply BTuD based unsupervised FE to the synthetic data sets with N = 1, 000, M = K = 20, N1 = 10, µ = 1. To get BTuD, we employed results given by HOOI (Higher Order Orthogonal Iteration of tensors) [1] with the setting “max iter=500,tol=1e-8” (we used tucker function in rTensor package [14]). We also restricted that ℓ1 ≤ 10 and ℓ2 , ℓ3 ≤ 5. HOOI successfully converged to muℓ1 i , muℓ2 j , and muℓ3 k (Here we did not use the proposed Algorithm 1). One might wonder why we can judge that HOOI converged to the solution required for BTuD and decided to use HOOI instead of proposed (time-consuming) algorithm, Algorithm 1. Since muℓ1 i , muℓ2 j , and muℓ3 k are the functions of uℓ1 i , uℓ2 j , uℓ3 k as shown 13
in eq. (18), they must be decided in the self-consistent manner; this means that we can check if solutions computed by HOOI can satisfy the conditions that the solution of BTuD must satisfy, based upon the self-consistency, muℓ1 i = uℓ1 i , muℓ2 j = uℓ2 j , and muℓ3 k = uℓ3 k . As a result we found that HOOI can converge to the solution that satisfies self-consistency and decided that we replace HOOI with the proposed time consuming Algorithm 1, also for the following examples other than the synthetic data. The results computed by HOOI are employed unless it is explicitly noticed that HOOI is not used. At first, we consider ℓ1 = 1 and compute Pi with eq. (29). The obtained 1,000 Pi values are corrected by BH criterion and is associated with adjusted P -values less than 0.05 are selected. As can be seen in Table 1, BTuD based unsupervised FE almost completely selected i ≤ 10 among 1,000 is for one hundred ensembles. One might wonder why considering ℓ1 = 1 enables us to select i ≤ N1 in spite of the unsupervised method since we did not provide any information about the properties of dataset at all. We define L2 X L3 X Yℓ1 jk = G(ℓ1 ℓ2 ℓ3 )uℓ2 j uℓ3 k . (39) ℓ2 =1 ℓ3 =1
Thus, if uℓ2 j and uℓ3 k are coincident with xijk , i ≤ N1 , Yℓ1 jk with larger absolute G(ℓ1 ℓ2 ℓ3 ) is also coincident with xijk , i ≤ N1 . Since xijk ≃
L1 X
Yℓ1 jk uℓ1 i ,
(40)
ℓ1 =1
uℓ1 i should take larger absolute values only for i ≤ N1 as well such that Yℓ1 jk is coincident with xijk , i ≤ N1 . This results in the smaller P -values attributed to i ≤ N1 by eq. (29) and i ≤ N1 will be selected. In fact, Y1jk as well as u1j and u1k are coincident with xijk , i ≤ N1 . For example, K M K we checked if u1j and u1k are distinct between j ≤ M 2 , k ≤ 2 and j > 2 , k > 2 , which means that u1j and u1k are coincident with xijk , i ≤ N1 . As can be seen in Table 2, in at least 85 out of 100 ensembles, they are significantly distinct and G(111) takes larger absolute values (not shown here), which means that Y1jk is supposed to be coincident with xijk , i ≤ N1 as well. In addition to this, we also checked if u1j u1k s K themselves are distinct between j ≤ M 2 , k ≤ 2 and others. Then as can be seen in Table 3, more than 90 out of 100 ensembles, they are distinct; in the remaining 10 out K of 100 ensembles, we can also find Y1jk s are distinct between j ≤ M 2 , k ≤ 2 and others with considering ℓ2 > 1 and ℓ3 > 1 since not only G(111) but also other G(1ℓ2 ℓ3 )s have negligible contributions as well. Thus the reason why absolute values of u1i s for i ≤ N1 are much larger than those for i > N1 is because u1j u1k s or Y1jk s are distinct K between j ≤ M 2 , k ≤ 2 and others (Fig. 3). K Why are (u1j , u1k )s or Y1jk s distinct between j ≤ M 2 , k ≤ 2 and others despite the unsupervised nature? When the number of uℓ1 i , uℓ2 j and uℓ3 k are restricted, it is better to select those coincident with the majority of xijk to minimize residuals. K Since xijk , i ≤ N1 share the distinction between j ≤ M 2 , k ≤ 2 and others whereas 14
0.1
0.2
u1k
−0.4
−0.3
−0.3
−0.2
−0.1
−0.1 0.0
u1j
0.0
0.1
0.05 −0.05
u1i −0.15
200
400
600
800
1000
−0.5
−0.5
−0.25 0
1
i
2
1
j
2
k
Fig. 3 The results of BuTD applied to synthetic data. Left: u1i . Red ones are i ≤ N1 . middle: and 2 for j > M , right: boxplot of u1k with 1 for k ≤ K and 2 for boxplot of u1j with 1 for j ≤ M 2 2 2 K k > 2 . As expected, since xijk , i ≤ N1 are more likely represented as a product of two features, u1j u1k , than others (i.e., i > N1 ) , xijk , i ≤ N1 s are associated with larger absolute u1i because of larger absolute G(111) and are accordingly associated with smaller P -values and selected.
xijk , i > N1 have nothing to share since they are composed of pure random numbers, it is better to select those coincident with xijk , i ≤ N1 and uℓi , i ≤ N1 have larger absolute values, which results in smaller P -values attributed to uℓ1 i , i ≤ N1 that are finally selected. Table 1 Confusion matrix of BTuD based unsupervised FE applied to synthetic data. i ≤ N1 are associated with distinction between two classes and are supposed to be selected. Those associated with adjusted P -values less than 0.05 are selected. i ≤ N1 and the selected ones are consistent well. BTuD based unsupervised FE i > N1 i ≤ N1 Adjusted P -value
> 0.05 ≤ 0.05
990 0
0.09 9.91
Table 2 The confirmation if u1j as well as u1k are consistent with the distinction between two classes. The number of ensembles where u1j or u1k is consistent with the distinction between j≤ M and j > M or k ≤ K and k > K , respectively. In the 2 2 2 2 majority of ensembles, u1j as well as u1k are consistent with the distinction between two classes. P -values are computed by t test.
Adjusted P -value
u1j > 0.05 ≤ 0.05
u1k > 0.05 ≤ 0.05
The number of ensembles
15
12
85
88
Although one might wonder whether BTuD based unsupervised FE might be useless unless there are majority features as in the above synthetic data, it is not the case. 15
Table 3 The confirmation if (u1j , u1k ) is consistent with the distinction between two classes. The number of ensembles where (u1j , u1k ) is coincident with the distinction between j ≤ M ,k≤ K and others. In the 2 2 majority of ensembles, (u1j , u1k ) is consistent with the distinction between two classes. P -values are computed by t test.
Adjusted P -value
(u1j , u1k ) > 0.05 ≤ 0.05
The number of ensembles
8
92
To show that BTuD based unsupervised FE is effective even if there are no majority features, we applied BTuD based unsupervised FE to xijk that includes sinusoidal functions with distinct phases as a part (N = 10, 000, M = 100, K = 1, N1 = 1, 000). We also restricted that ℓ1 ≤ 10 and ℓ2 ≤ 2, ℓ3 ≤ 1. Since sinusoidal functions with distinct phases are distinct with each other and cannot be composed of majority features, one might think BTuD based unsupervised FE cannot detect them. In spite of the concern, BTuD based unsupervised FE successfully selects a set of sinusoidal functions (i ≤ N1 ) among 10,000 features (Table 4) when we consider ℓ1 = 1, 2 and compute Pi with eq. (29); the obtained 10,000 Pi s are corrected by BH criterion and is associated with adjusted P -values less than 0.05 are selected. The reason why it could occur is as follows; since sinusoidal functions, A sin(x + δ), with distinct phases, δ, that correspond to xij , i ≤ N1 can be represented as a linear combination of sin function and cos function, A sin(x + δ) = A cos δ sin x + A sin δ cos x,
(41)
majority features that are sinusoidal and allow us to select xij , i ≤ N1 are represented as well with u1j and u2j if u1j and u2j correspond to sin and cos, respectively. As can be seen in Fig. 4, u1j and u2j are sinusoidal. Since we use SVD instead of TD and SVD makes uℓj to directly correspond to uℓi , we can select xij , i ≤ N1 by using u1i and u2i to select is. Table 4 Confusion matrix of BTuD applied to sinusoidal data averaged over 100 ensembles. i ≤ N1 are associated with sinusoidal function and are supposed to be selected. Those associated with adjusted P -values less than 0.05 are selected. i ≤ N1 and the selected ones are consistent well. BTuD based unsupervised FE i > N1 i ≤ N1 Adjusted P -value
> 0.05 ≤ 0.05
16
8999.75 0.25
0.00 1000
0.00
u1i
0.02
0.04
0.15 0.05
u2j
−0.10 −0.05
−0.15 −0.10 −0.05 −0.02
0.00
0.00
u1j
0.05
0.10
0.10
0.04 0.02 0.00
u2i
−0.02 −0.04 −0.04
0
20
40
60
j
80
100
0
20
40
60
80
100
j
Fig. 4 The results of BuTD applied to sinusoidal data. Left: u1i vs u2i . Red ones are i ≤ N1 . As expected, since xij , i ≤ N1 s that are sinusoidal are more likely P represented as a linear P combination of two sinusoidal features, u1j (middle) and u2j (left), u1i ∝ j u1j xij and u2i ∝ j u2j xij should have larger absolute values than others (i > N1 ). Then u1i and u2i with i ≤ N1 are associated with smaller P -values and selected.
In the above two examples, we have shown that BTuD based unsupervised FE can be successful when we know what the correct answer is. In the following, we applied BTuD based unsupervised FE to two cases, whose correct answers we do not know in advance; one is the result of numerical model and another is that of real data, gene expression profiles. The first (numerical) example is so-called Randomized Coupling Strength GCM (RCS-GCM) [16]. By tuning parameters, RCS-GCM can generate the mixture of ordered (periodic) states and random (chaotic) states. The performance test of BTuD based unsupervised FE applied to RCS-GCM is whether BTuD based unsupervised FE can, in a fully unsupervised manner, distinguish ordered states from random states in RCS-GCM. As can be seen in Table 5, using Pi computed with eq. (29) with setting ℓ1 = 1, we can identify many i values which are apparently composed of three states (Fig. 5) in spite of that we intentionally select neither ordered state nor three states at all. The reason why three-state variables can be selected in the unsupervised manner is because u1j represents the three state vector that enables us to select xij s associated with three state unintentionally (Fig. 6). Since here we again analysed not a tensor but a matrix, correspondence between uℓ1 i and uℓ2 j is again always that ℓ1 = ℓ2 , since G(ℓ1 ℓ2 ) is always diagonal. Here we would like to emphasize that the separation between ordered states and others in the present study is much better than the previous one [16] where TD based unsupervised FE was employed. Thus, we can conclude that BTuD based unsupervised FE can deal with the data sets, correct answers of which we do not know. The second real example is gene expression profiles [17]. With applying BTuD based unsupervised FE to gene expression profiles, we can identify not one but multiple sets of genes each of which is associated with distinct properties (tissue specificities) as identified previously [18, 19]. Although the following explanation is almost the replication of the previous study [18], since the employed method is not TD based unsupervised FE but BTuD based unsupervised FE in the present study, we briefly explain the procedure. At first, we need to identify which uℓ2 j s are associated with tissue specificity. Fig. 7 shows the uℓ2 j , 1 ≤ ℓ2 ≤ 6. Although u1j does not exhibit clear tissue-specificity, other five uℓ2 j s have clear tissue specific expression. Among them, we selected ℓ2 = 2, 4, 5, 6 that simultaneously exhibit specificities for multiple tissues 17
Table 5 The performance of BTuD based unsupervised FE applied to RCS-GCM. The 1707 is are associated with adjusted P -values less than 0.05 and are selecteed. consistent with ordered state (see Fig. 5)
Adjusted P -value
(u1j , u1k ) > 0.05 ≤ 0.05
The number of is
8293
1707
7129 758 4870 3873 3241 5889 2409 5893 5687 5337 6440 2974 7172 3952 1914 5502 7084 2489 5510 1990 8226 7040 6044 2777 6928 5676 160 2540 626 5515 7230 3857 6764 6775 6003 1093 4595 1916 4214 4911 6062 5534 53 5075 161 8195 7073 3178 2291 7495 1137 4523 6323 4566 6609 5437 569 7978 3248 5365 147 3319 4450 862 6782 3864 3820 290 1994 5024 2040 4091 1527 3005 1954 7346 1436
57 60 42 51 87 96 72 66 78 79 47 65 74 94 89 12 24 36 27 18 3 39 37 31 22 98 83 16 10 4 100 67 76 44 32 61 59 5 43 71 52 26 14 88 70 8 20 58 85 38
2 11 23 17 65 71 77 62 95 92 86 50 56 47 38 26 29 1 12 18 21 24 30 42 45 60 54 99 96 63 69 84 75 81 7 13 22 25 28 34 46 100 85 97 49 52 61 67 73 79
769 897 954 439 348 1639 1047 860 145 1363 17 1672 175 1682 590 1157 1704 22 1333 155 814 164 708 157 524 511 1542 785 295 1505 1254 1411 887 1598 1648 143 467 1464 549 648 1524 1692 443 573 947 1480 1377 1223 4 1217 1663 1597 1410 620 1071 837 80 1340 1568 1616 1288 1301 996 1136 1516 1046 181 1023 815 1511 7 1517 718 1012 1459
u1j
−0.10
−0.05
0.00
0.05
0.10
0.15
Fig. 5 Heatmap of xij for 1707 is selected (left) and 8293 is not selected (right) by BTuD based unsupervised FE, respectively (Table 5). Rows (i) and columns (j) are clustered by the hierarchical clustering. It is obvious that is with three states are selected and not selected is lack ordered state.
0
20
40
60
80
100
j
Fig. 6
The result of BTD applied to RCS-GCM. Vertical axis is u1j and horizontal axes is j.
18
(Please notice that ℓ2 = 5, 6 are associated with the combination of two common tissues). Next we found that u2k and u3k exhibit distinction between drug treatments and controls and u1m exhibits independence of two replicates (not shown here). Thus we can seek |G(ℓ1 , ℓ2 , ℓ3 , 1)|, ℓ2 = 2, 4, 5, 6, ℓ3 = 2, 3 to find which ℓ1 s are associated with tissue specificities and distinction between drug treatments and controls, simultaneously (Fig. 8). The selected ℓ1 s are used to attribute P -values to i using eq. (29). Table 6 shows the number of selected genes. To validate genes selected by BTuD Table 6 The number of selected genes associated with tissue specificities and distinction between drug treatments and controls Adjusted P -value
> 0.05
≤ 0.05
ℓ2 = 2 ℓ2 = 4 ℓ2 = 5 ℓ2 = 6
24052 24120 23834 24028
369 301 587 393
l_2=2
l_2=3
−0.4
l_2=6
−0.2
0.2
0.2
0.4
0.6
0.6
l_2=5
−0.6
0.0
−0.4
0.0
l_2=4
−0.8
−0.8
0.0
−0.20
0.2
0.4
−0.10
0.6
0.0
0.00
l_2=1
Fig. 7 The results of BuTD applied to gene expression: uℓ2 j . Horizontal axes are 24 tissues. Top left: ℓ2 = 1, top middle: ℓ2 = 2, top right: ℓ2 = 3, bottom left ℓ2 = 4, bottom middle: ℓ2 = 5, bottom right: ℓ2 = 6. Please notice that ℓ2 = 5, 6 are associated with the distinct combination of two common tissues.
based unsupervised FE, we compare the genes selected by BTuD based unsupervised FE with those selected by TD based unsupervised FE [19]. Since genes selected by TD based unsupervised FE were biologically evaluated and proven to be reasonable, if the 19
l2=4
20 0
1 2 3 4 5 6 7 8 9 10
1 2 3 4 5 6 7 8 9 10 l1
l2=5
l2=6 100
l1
0
20
20
40
40
60
60
80
80 100
0
50
40
60
100
80
150
200
120
l2=2
1 2 3 4 5 6 7 8 9 10 l1
1 2 3 4 5 6 7 8 9 10 l1
Fig. 8 The results of BuTD applied to gene expression:|G(ℓ1 , ℓ2 , ℓ3 , 1)|. Horizontal axes are ℓ1 . Top left: ℓ2 = 2, top right: ℓ2 = 4, bottom left ℓ2 = 5, bottom right: ℓ2 = 6. Vertical arrows show the selected ℓ1 s. Red: ℓ3 = 2, blue: ℓ3 = 3.
overlap between genes selected by BTuD based unsupervised FE and those selected by TD based unsupervised FE is significant, we can judge that BTuD based unsupervised FE selected reasonable genes. Figure 9 shows the Venn diagram between genes selected by BTuD based unsupervised FE and those selected by TD based unsupervised FE; since the former is almost always the subset of the latter, the overlap between them is definitely significant (Since ℓ2 = 5, 6 are associated with the combination of two common tissues, they are grouped together). Thus, BTuD based unsupervised FE is not only applicable to real data sets, but also can deal with data sets including sets of variables associated with multiple properties.
20
l2=2
l2=4
l2=5,6
BTD_b BTD
TD
BTD
TD
TD_a 563
6 BTD_a
370
103
302
17 33
22
362
302
72
122
TD_b 61
14
191
Fig. 9 Venn diagram of the selected genes between BTuD based unsupervised FE (present study) and TD based unsupervised FE (previous study [19]). Left: ℓ2 = 2, middle: ℓ2 = 4, right:ℓ2 = 5, 6. ℓ2 = 5 and ℓ2 = 6 correspond to BTuD/TD a and BTuD/TD b, respectively. Since ℓ2 = 5, 6 are associated with the combination of two common tissues, they are treated together.
In conclusion, BTuD based unsupervised FE can be applied to a wide range of data sets without pre-assigned labels or prior knowledges.
5 Discussion 5.1 Equivalence between BTuD based unsupervised FE and TD based unsupervised FE Since BTuD-based unsupervised FE is practically identical to TD-based unsupervised FE especially in the last example, gene expression, it is important to clarify the reason of this coincidence; TD-based unsupervised FE is known to be successfully applied to many problems [1]. P -values are attributed to features with eqs. (25) and (29) for TD based unsupervised FE and BTuD based unsupervised FE, respectively. In TD based unsupervised FE, mean and SD are estimated based upon the distribution of not selected features that are assumed to obey Gaussian distribution in BTuD based unsupervised FE. Thus Gaussian distribution attributed to not selected features are empirically identical and thus the same P -values are attributed to not selected features in both methods. The Gaussian distributions with common SD are attributed to the selected features in both methods, but they are centered distinctly. However, if we attribute P -values of zero to selected features as N (0|uℓ1 i , SD) in BTuD based unsupervised FE that are identical to N (uℓ1 i |0, SD) in TD based unsupervised FE, the attributed P -values are identical as well. Thus, both methods attribute the same P -values to selected and not selected features, respectively. This means that both methods are empirically identical (Fig. 10). One should also notice that the equivalence between TD and BTuD based unsupervised FE stands only when TD and BTuD give the common uℓ1 i by which we can select is. In Fig. 9, we were forced to treat ℓ2 = 5, 6 as one group, since uℓ1 i , ℓ = 5, 6 used to attribute P -values to is are not identical between TD and BTuD. Nevertheless, ℓ2 = 5, 6 as one group are coincident between TD and BTuD; it is biologically 21
BTuD based unsupervised FE
TD based unsupervised FE
7
7
6
6
5
5 4
i
i
4
−1
3
3
2
2
1
1
0
1
−1
2
uℓ1 i
0
1
2
uℓ1 i
Fig. 10 The schematic figure that explains the equivalence between BTuD based unsupervised FE (left) and TD based unsupervised FE (right). Blue and red dots corresponds to uℓ1 i attributed to not selected and selected features, respectively. For BTuD based unsupervised FE, Gaussian distributions having the common SD and distinct means are attributed to features to be selected (red) and those not to be selected (blue), respectively. Since the latter has zero within the 95 % confidence intervals, they are rejected whereas the former is not rejected since zero is outside the 95 % confidence interval. For TD based unsupervised FE, no Gaussian distribution is attributed to individual features, but two Gaussian distribution that have distinct means between features not to be selected (blue) and those to be selected (red) are assumed to exist. The latter is decided if they are regarded to be outliers based upon Gaussian distribution attributed to features not to be selected as a whole (blue). Since SD is common for all the Gaussian distributions, the same P -values are attributed to individual features in both methods and the same features are selected in both method.
reasonable since ℓ2 = 5, 6 are associated with the distinct combination of two common tissues specificity (Fig. 7).
6 Conclusions In this paper, we have proposed a new implementation of Bayesian Tucker decomposition (BTuD) to perform unsupervised feature selection. BTuD was practically was practically implemented by HOOI and we can attributed P -values to features with following Bayesian theory. BTuD based unsupervised FE works well for various datasets. BTuD is also expected to coincide with the previously proposed TD based unsupervised FE. 22
7 Acknowledgement This work was supported by JSPS KAKENHI 22K13979, 23K28150, 24K15168, 24K22309, 24H00247, 25K00986, 25H01470; JST; and PRESTO Grant numbers JPMJPR212A, JST CREST JPMJCR2431, and NEDO JPNP22100843-0.
References [1] Taguchi, Y.-h.: Unsupervised Feature Extraction Applied to Bioinformatics: A PCA Based and TD Based Approach, 2nd edn. Unsupervised and SemiSupervised Learning. Springer, Switzerland (2024) [2] Cheng, L., Chen, Z., Wu, Y.-C.: Bayesian Tensor Decomposition for Signal Processing and Machine Learning, 1st edn. Springer, Cham, Switzerland (2023) [3] Kanagawa, H., Suzuki, T., Kobayashi, H., Shimizu, N., Tagami, Y.: Gaussian process nonparametric tensor estimator and its minimax optimality. In: Balcan, M.F., Weinberger, K.Q. (eds.) Proceedings of The 33rd International Conference on Machine Learning. Proceedings of Machine Learning Research, vol. 48, pp. 1632–1641. PMLR, New York, New York, USA (2016). https://proceedings.mlr.press/v48/kanagawa16.html [4] Wimalawarne, K., Sugiyama, M., Tomioka, R.: Multitask learning meets tensor factorization: task imputation via convex optimization. In: Ghahramani, Z., Welling, M., Cortes, C., Lawrence, N., Weinberger, K.Q. (eds.) Advances in Neural Information Processing Systems, vol. 27. Curran Associates, Inc., ??? (2014). https://proceedings.neurips.cc/paper files/paper/2014/file/6d9cb7de5e8ac30bd5e8734bc96a35c1Paper.pdf [5] Takayama, H., Zhao, Q., Hontani, H., Yokota, T.: Bayesian Tensor Completion and Decomposition with Automatic CP Rank Determination Using MGP Shrinkage Prior. SN Computer Science 3(3), 225 (2022) https://doi.org/10.1007/ s42979-022-01119-8 [6] Zhao, Q., Zhang, L., Cichocki, A.: Bayesian Sparse Tucker Models for Dimension Reduction and Tensor Completion (2015). https://arxiv.org/abs/1505.02343 [7] Mørup, M., Hansen, L.K.: Automatic relevance determination for multi-way models. Journal of Chemometrics 23(7-8), 352–363 (2009) https://doi.org/10.1002/cem.1223 https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/pdf/10.1002/cem.1223
[8] Tong, X., Cheng, L., Wu, Y.-C.: Bayesian tensor tucker completion with a flexible core. IEEE Transactions on Signal Processing 71, 4077–4091 (2023) https://doi. org/10.1109/TSP.2023.3327845
23
[9] Tipping, M.E., Bishop, C.M.: Probabilistic principal component analysis. Journal of the Royal Statistical Society. Series B (Statistical Methodology) 61(3), 611–622 (1999). Accessed 2024-09-16 [10] Bishop, C.M.: Pattern Recognition and Machine Learning, 1st edn. Information Science and Statistics. Springer, New York, NY (2006) [11] Ho, T.K.: Random decision forests. In: Proceedings of 3rd International Conference on Document Analysis and Recognition, vol. 1, pp. 278–282 (1995). IEEE [12] Hartigan, J.A., Hartigan, P.M.: The dip test of unimodality. The Annals of Statistics 13(1), 70–84 (1985). Accessed 2025-01-07 [13] Barata, J.C.A., Hussein, M.S.: The Moore–Penrose Pseudoinverse: A Tutorial Review of the Theory. Brazilian Journal of Physics 42(1-2), 146–165 (2012) https: //doi.org/10.1007/s13538-011-0052-z [14] Li, J., Bien, J., Wells, M.T.: rTensor: An R package for multidimensional array (tensor) unfolding, multiplication, and decomposition. Journal of Statistical Software 87(10), 1–31 (2018) https://doi.org/10.18637/jss.v087.i10 [15] Kaneko, K.: Globally coupled chaos violates the law of large numbers but not the central-limit theorem. Phys. Rev. Lett. 65, 1391–1394 (1990) https://doi.org/10. 1103/PhysRevLett.65.1391 [16] Mototake, Y.-i., Taguchi, Y.-h.: Signal identification without signal formulation (2023). https://arxiv.org/abs/2304.06522 [17] Kozawa, S., Sagawa, F., Endo, S., De Almeida, G.M., Mitsuishi, Y., Sato, T.N.: Predicting human clinical outcomes using mouse multi-organ transcriptome. iScience 23(2), 100791 (2020) https://doi.org/10.1016/j.isci.2019.100791 [18] Taguchi, Y.-h., Turki, T.: Universal nature of drug treatment responses in drug-tissue-wide model-animal experiments using tensor decomposition-based unsupervised feature extraction. Frontiers in Genetics 11 (2020) https://doi.org/ 10.3389/fgene.2020.00695 [19] Taguchi, Y.-h., Turki, T.: Adapted tensor decomposition and PCA based unsupervised feature extraction select more biologically reasonable differentially expressed genes than conventional methods. Scientific Reports 12(1), 17438 (2022) https: //doi.org/10.1038/s41598-022-21474-z
24