Learning the Sparse and Low Rank PARAFAC Decomposition via the Elastic Net

TL;DR

This paper introduces a Bayesian model with elastic net regularization for sparse, low-rank PARAFAC tensor decomposition, capable of handling missing data, noise, and large-scale problems, with an efficient solution path for rank and sparsity selection.

Contribution

It proposes a novel Bayesian framework and algorithms for sparse, low-rank tensor decomposition that automatically determine the true rank and sparsity, improving robustness and scalability.

Findings

Effective on simulation data

Successful application to real data

Robust to noise and missing values

Abstract

In this article, we derive a Bayesian model to learning the sparse and low rank PARAFAC decomposition for the observed tensor with missing values via the elastic net, with property to find the true rank and sparse factor matrix which is robust to the noise. We formulate efficient block coordinate descent algorithm and admax stochastic block coordinate descent algorithm to solve it, which can be used to solve the large scale problem. To choose the appropriate rank and sparsity in PARAFAC decomposition, we will give a solution path by gradually increasing the regularization to increase the sparsity and decrease the rank. When we find the sparse structure of the factor matrix, we can fixed the sparse structure, using a small to regularization to decreasing the recovery error, and one can choose the proper decomposition from the solution path with sufficient sparse factor matrix with low recovery error. We test the power of our algorithm on the simulation data and real data, which show it is powerful.

Tables

Table 1. Table 1: Solution path by Algorithm 3 on the simultated data when α = 0.2 𝛼 0.2 \alpha=0.2 . R: rank, NZS: the number of zeros in the current solution, NZT: the number of the zeros such that the zero appears in the current solution and the true solution; IS1: which algorithm is used? 1–algorithm 1 and 0– algorithm 2 ;

$\lambda$	R	NZS	NZT	IS1	rel_err	iters
1e-10	3	3	2	1	8.8e-06	38
1e-08	3	3	2	0	6.4e-06	2
1e-04	3	3	2	1	1.8e-05	56
1e-08	3	3	2	0	3.5e-06	3
1e-01	3	8	4	1	3.7e-03	200
1e-08	3	8	4	0	2.3e-05	100
1e+00	2	5	5	1	2.3e-03	200
1e-08	2	5	5	0	4.4e-06	12
1e+01	2	7	7	1	1.3e-02	200
1e-08	2	7	7	0	3.2e-06	15
2e+01	2	8	8	1	2.5e-02	200
1e-08	2	8	8	0	4.0e-06	16
2e+02	2	10	10	1	1.9e-01	108
1e-08	2	10	10	0	3.3e-06	17
4e+02	2	11	11	1	3.2e-01	63
1e-08	2	11	11	0	4.4e-06	8
5e+02	2	12	12	1	3.8e-01	53
1e-08	2	12	12	0	1.6e-06	9
7e+02	1	3	3	1	6.1e-01	56
1e-08	1	3	3	0	5.4e-01	10
2e+03	0	0	0	1	1.0e+00	6

Table 2. Table 2: 50 times random test behavior when ( I 1 , I 2 , I 3 ) = ( 20 , 20 , 20 ) subscript 𝐼 1 subscript 𝐼 2 subscript 𝐼 3 20 20 20 (I_{1},I_{2},I_{3})=(20,20,20) , R = 10 𝑅 10 R=10 where the data with noise ( SNR = 20dB) and without missing entries and each algorithm using random initialization (given in above) and nvecs initialization(given below). And SPML stands for the Solution Path Method using algorithm 1 (Look algorithm), and And SPMS stands for the Solution Path Method using algorithm 2 (Sparse constrained Algorithm). And SPML(0.8) standes for that α = 0.8 𝛼 0.8 \alpha=0.8 . And SCORES(STD) stand for the mean score and the its standard deviation. REL_ERR stands for the mean of the full data relative error ‖ 𝐙 ¯ − 𝐗 ¯ ‖ F ‖ 𝐙 ¯ ‖ F subscript norm ¯ 𝐙 ¯ 𝐗 𝐹 subscript norm ¯ 𝐙 𝐹 \frac{||\underline{\bf Z}-\underline{\bf X}||_{F}}{||\underline{\bf Z}||_{F}} subtract the true relative error ‖ 𝐄 ¯ ‖ F ‖ 𝐙 ¯ ‖ F subscript norm ¯ 𝐄 𝐹 subscript norm ¯ 𝐙 𝐹 \frac{||\underline{\bf E}||_{F}}{||\underline{\bf Z}||_{F}} . TFNZ stands for the total number of zeros found by the selected solution; TFTNZ stands for total number of the zeros such that the zero appears in both the selected solution and the true solution; TNZ stands for total number of the zeros in the true solution. TIME stands for the the mean cost cpu time in 50 random testes, note that for each α 𝛼 \alpha , SPML and SPMS are computed in together, so the time are both the SPML and SPMS in each test.

Method	Scores(std)	rel_err	TFNZ	TFTNZ	TNZ	$\frac{\text{TFNZ }}{\text{TNZ}}$	$\frac{\text{TTNZ }}{\text{TNZ}}$	time
random	initialization
CP-AlS	0.566(0.360)	-2.436e-03	0	0	5078	0.000	0.000	0.6
LRTI	0.413(0.278)	1.918e-02	0	0	5078	0.000	0.000	1.2
SPML(0)	0.579(0.400)	1.824e-02	0	0	5009	0.000	0.000	73.7
SPMS(0)	0.648(0.387)	1.435e-02	0	0	5009	0.000	0.000	73.7
SPML(0.2)	0.922(0.141)	1.178e-02	2225	2130	5009	0.444	0.425	74.5
SPMS(0.2)	0.967(0.121)	-3.351e-03	2225	2225	5009	0.444	0.444	74.5
SPML(0.8)	0.934(0.164)	2.771e-03	3003	2933	5009	0.600	0.586	67.3
SPMS(0.8)	0.952(0.160)	-2.646e-03	3003	2989	5009	0.600	0.597	67.3
SPML(0.98)	0.946(0.144)	1.956e-03	3168	3083	5009	0.632	0.615	67.1
SPMS(0.98)	0.960(0.140)	-3.324e-03	3168	3112	5009	0.632	0.621	67.1
nvecs	initialization
CP-AlS	0.739(0.289)	-9.452e-03	0	0	5078	0.000	0.000	0.4
LRTI	0.647(0.240)	1.632e-02	2	0	5078	0.000	0.000	2.4
SPML(0)	0.784(0.324)	7.721e-03	0	0	5009	0.000	0.000	109.0
SPMS(0)	0.832(0.285)	4.766e-03	0	0	5009	0.000	0.000	109.0
SPML(0.2)	0.883(0.226)	1.360e-02	2274	2146	5009	0.454	0.428	113.0
SPMS(0.2)	0.927(0.224)	-2.580e-03	2274	2241	5009	0.454	0.447	113.0
SPML(0.8)	0.888(0.251)	4.696e-03	3090	2905	5009	0.617	0.580	105.3
SPMS(0.8)	0.909(0.252)	-1.160e-03	3090	2967	5009	0.617	0.592	105.3
SPML(0.98)	0.897(0.234)	4.672e-03	3344	3037	5009	0.668	0.606	104.4
SPMS(0.98)	0.921(0.230)	-1.498e-03	3344	3084	5009	0.668	0.616	104.4

Table 3. Table 3: 50 times random test behavior when ( I 1 , I 2 , I 3 ) = ( 50 , 50 , 50 ) subscript 𝐼 1 subscript 𝐼 2 subscript 𝐼 3 50 50 50 (I_{1},I_{2},I_{3})=(50,50,50) , R = 9 𝑅 9 R=9 where the data with noise ( SNR = 20dB) and with missing entries(random drop out 25% entries) and each algorithm using random initialization (given in above) and nvecs initialization(given below). And SPML stands for the Solution Path Method using algorithm 1 (Look algorithm), and SPMS stands for the Solution Path Method using algorithm 2 (Sparse constrained Algorithm) where the solution is coming from SPML. And SPML(0.8) standes for that α = 0.8 𝛼 0.8 \alpha=0.8 . And SCORES(STD) stand for the mean score and the its standard deviation. SPMLR stands for the Solution Path Method using algorithm 4 , And SPMSR stands for the Solution Path Method using algorithm 2 (Sparse constrained Algorithm) where the solution is coming from SPMLR. REL_ERR stands for the mean of the full data relative error ‖ 𝐙 ¯ − 𝐗 ¯ ‖ F ‖ 𝐙 ¯ ‖ F subscript norm ¯ 𝐙 ¯ 𝐗 𝐹 subscript norm ¯ 𝐙 𝐹 \frac{||\underline{\bf Z}-\underline{\bf X}||_{F}}{||\underline{\bf Z}||_{F}} subtract the true relative error ‖ 𝐄 ¯ ‖ F ‖ 𝐙 ¯ ‖ F subscript norm ¯ 𝐄 𝐹 subscript norm ¯ 𝐙 𝐹 \frac{||\underline{\bf E}||_{F}}{||\underline{\bf Z}||_{F}} . TFNZ stands for the total number of zeros found by the selected solution; TFTNZ stands for total number of the zeros such that the zero appears in both the selected solution and the true solution; TNZ stands for total number of the zeros in the true solution. TIME stands for the the mean cost cpu time in 50 random testes, note that for each α 𝛼 \alpha , SPML and SPMS are computed together, so the time was spent in both the SPML and SPMS in each test. NFTR stands for number of find true rank in the 30 testes.

Method	Scores(std)	rel_err	TFNZ	TFTNZ	TNZ	$\frac{\text{TFNZ }}{\text{TNZ}}$	$\frac{\text{TTNZ }}{\text{TNZ}}$	time	NFTR
random	initialization
LRTI	0.726(0.256)	6.641e-03	0	0	6985	0.000	0.000	222.6	10
SPMLR(0.2)	0.851(0.154)	1.008e-01	567	507	7064	0.080	0.072	458.9	30
SPMSR(0.2)	1.000(0.000)	-3.376e-04	567	567	7064	0.080	0.080	458.9	30
SPMLR(0.8)	0.910(0.114)	3.316e-02	2846	2619	7064	0.403	0.371	494.6	30
SPMSR(0.8)	1.000(0.000)	-3.200e-04	2846	2846	7064	0.403	0.403	494.6	30
SPMLR(0.98)	0.957(0.093)	1.216e-02	4658	4479	7064	0.659	0.634	541.8	30
SPMSR(0.98)	1.000(0.000)	-3.037e-04	4658	4658	7064	0.659	0.659	541.8	30
SPML(0.2)	0.943(0.104)	6.954e-02	6722	6454	7064	0.952	0.914	322.1	30
SPMS(0.2)	1.000(0.000)	-2.863e-04	6722	6722	7064	0.952	0.952	322.1	30
SPML(0.8)	0.990(0.040)	1.612e-02	7039	6992	7064	0.996	0.990	322.2	30
SPMS(0.8)	1.000(0.000)	-2.826e-04	7039	7039	7064	0.996	0.996	322.2	30
SPML(0.98)	0.983(0.056)	1.427e-02	7032	6948	7064	0.995	0.984	323.1	30
SPMS(0.98)	1.000(0.000)	-2.827e-04	7032	7032	7064	0.995	0.995	323.1	30
nvecs	initialization
LRTI	0.670(0.236)	7.102e-02	2	0	6985	0.000	0.000	260.4	6
SPMLR(0.2)	0.836(0.158)	1.067e-01	575	511	7064	0.081	0.072	444.3	30
SPMSR(0.2)	1.000(0.000)	-3.382e-04	575	575	7064	0.081	0.081	444.3	30
SPMLR(0.8)	0.936(0.110)	3.260e-02	2799	2649	7064	0.396	0.375	487.3	30
SPMSR(0.8)	1.000(0.000)	-3.205e-04	2799	2799	7064	0.396	0.396	487.3	30
SPMLR(0.98)	0.954(0.106)	1.160e-02	4307	4134	7064	0.610	0.585	543.6	30
SPMSR(0.98)	1.000(0.000)	-3.068e-04	4307	4307	7064	0.610	0.610	543.6	30
SPML(0.2)	0.943(0.104)	6.954e-02	6722	6454	7064	0.952	0.914	657.3	30
SPMS(0.2)	1.000(0.000)	-2.863e-04	6722	6722	7064	0.952	0.952	657.3	30
SPML(0.8)	0.876(0.296)	1.541e-02	6781	6196	7064	0.960	0.877	610.1	26
SPMS(0.8)	0.886(0.297)	-2.858e-04	6781	6243	7064	0.960	0.884	610.1	26
SPML(0.98)	0.844(0.333)	1.339e-02	7000	6138	7064	0.991	0.869	590.4	26
SPMS(0.98)	0.868(0.336)	-2.888e-04	7000	6274	7064	0.991	0.888	590.4	26

Table 4. Table 4: Performance comparison of the algorithms in COIL-20 image clustering, where Rank is the initial rank, where R and N in the first row in the parenthesis are stand for the random initialization and nvecs initialization respectively.

Algorithms	CP-ALS(R)	CP-ALS(N)	SPMLR(0.8R)	SPMSR(0.8R) )	SPMLR(0.8N)	SPMSR(0.8N)
Rank=20
RelErr	0.294	0.291	0.316	0.301	0.309	0.299
NumZeros	0	0	31	31	26	26
accuracy(%)	0.637	0.602	0.700	0.724	0.760	0.762
time(s)	514.8	478.1	40882.5	40882.5	45547.4	45547.4
Rank=30
RelErr	0.265	0.263	0.296	0.278	0.284	0.272
NumZeros	0	0	30	30	23	23
accuracy(%)	0.561	0.670	0.717	0.752	0.676	0.751
time(s)	550.7	566.8	43123.7	43123.7	47380.1	47380.1

Table 5. Table 5: Solution path by Algorithm 3 on the simultated data when α = 0.2 𝛼 0.2 \alpha=0.2 . R: rank, NZS: the number of zeros in the current solution, NZT: the number of the zeros such that the zero appears in the current solution and the true solution; IS1: which algorithm is used? 1–algorithm 1 and 0– algorithm 2 ;

$\lambda$	R	NZS	NZT	IS1	rel_err	iters
1e-10	3	3	2	1	8.8e-06	38
1e-08	3	3	2	0	6.4e-06	2
1e-09	3	3	2	1	6.4e-06	2
1e-08	3	3	2	1	5.3e-06	1
1e-07	3	3	2	1	4.4e-06	1
1e-06	3	3	2	1	2.4e-06	3
1e-05	3	4	2	1	2.1e-06	3
1e-08	3	4	2	0	1.6e-06	1
1e-04	3	3	2	1	1.8e-05	56
1e-08	3	3	2	0	3.5e-06	3
1e-03	3	3	2	1	1.7e-04	200
1e-02	3	3	2	1	1.2e-03	200
1e-01	3	8	4	1	3.7e-03	200
1e-08	3	8	4	0	2.3e-05	100
1e+00	2	5	5	1	2.3e-03	200
1e-08	2	5	5	0	4.4e-06	12
1e+01	2	7	7	1	1.3e-02	200
1e-08	2	7	7	0	3.2e-06	15
2e+01	2	8	8	1	2.5e-02	200
1e-08	2	8	8	0	4.0e-06	16
3e+01	2	8	8	1	3.6e-02	200
4e+01	2	8	8	1	4.8e-02	200
5e+01	2	8	8	1	5.8e-02	200
6e+01	2	8	8	1	6.9e-02	200
7e+01	2	8	8	1	7.9e-02	178
8e+01	2	8	8	1	8.8e-02	160
9e+01	2	8	8	1	9.8e-02	146
1e+02	2	8	8	1	1.1e-01	134
2e+02	2	10	10	1	1.9e-01	108
1e-08	2	10	10	0	3.3e-06	17
3e+02	2	10	10	1	2.6e-01	63
4e+02	2	11	11	1	3.2e-01	63
4e+02	2	11	11	0	4.4e-06	8
5e+02	2	12	12	1	3.8e-01	53
1e-08	2	12	12	0	1.6e-06	9
6e+02	2	12	12	1	4.4e-01	45
7e+02	1	3	3	1	6.1e-01	56
1e-08	1	3	3	0	5.4e-01	10
8e+02	1	3	3	1	6.3e-01	34
9e+02	1	3	3	1	6.4e-01	31
1e+03	1	3	3	1	6.5e-01	28
2e+03	0	0	0	1	1.0e+00	6

Equations

$${\underline{\mathbf{Z}}}_{i_1,i_2,\dots ,i_N}={\underline{\mathbf{X}}}_{i_1,i_2,\dots ,i_N}+{\underline{\mathbf{E}}}_{i_1,i_2,\dots ,i_N}$$

$${\underline{\mathbf{X}}}_{i_1,i_2,\dots ,i_N}=\sum _{r=1}^R{\mathbf{A}}^{(1)}(i_1,r){\mathbf{A}}^{(2)}(i_2,r)\cdots {\mathbf{A}}^{(N)}(i_N,r)$$

$$\underline{\mathbf{X}}$$

$$=\sum _{r=1}^R{\gamma }_r({\mathbf{u}}_r^{(1)}\circ {\mathbf{u}}_r^{(2)}\circ \cdots \circ {\mathbf{u}}_r^{(N)})$$

$$p_n({\mathbf{a}}_r^{(n)})={\beta }_{n}exp(-1/2{\mathbf{a}}_r^{(n)}{\mathbf{R}}_n^{-1}{\mathbf{a}}_r^{(n)}-{\mu }_n||{\mathbf{a}}_r^{(n)}|{|}_1)$$

$$\theta :=\text{Tr}({\mathbf{R}}_1)=\text{Tr}({\mathbf{R}}_2)=\cdots =\text{Tr}({\mathbf{R}}_N)$$

$$\mu :={\mu }_1={\mu }_2=\dots ={\mu }_N$$

$$L(\underline{\mathbf{X}})=\frac{1}{2{\sigma }^2}||(\underline{\mathbf{Z}}-\sum _{r=1}^R{\mathbf{a}}_r^{(1)}\circ {\mathbf{a}}_r^{(2)}\circ \cdots \circ {\mathbf{a}}_r^{(N)})\text{varoast}\underline{\mathbf{\Delta }}|{|}_F^2$$

$$+\sum _{r=1}^R\sum _{n=1}^N\left[\frac{1}{2}(({\mathbf{a}}_r^{(n)}{)}^T{\mathbf{R}}_n^{-1}{\mathbf{a}}_r^{(n)})+\mu ||{\mathbf{a}}_r^{(n)}|{|}_1\right]$$

$$\underset{\mathbf{A}}{\min }\frac{1}{2}||(\underline{\mathbf{Z}}-\sum _{r=1}^R{\mathbf{a}}_r^{(1)}\circ {\mathbf{a}}_r^{(2)}\circ \cdots \circ {\mathbf{a}}_r^{(N)})\text{varoast}\underline{\mathbf{\Delta }}|{|}_F^2$$

$$+\sum _{r=1}^R\sum _{n=1}^N\lambda \left[\frac{1-\alpha }{2}(({\mathbf{a}}_r^{(n)}{)}^T{\mathbf{R}}_n^{-1}{\mathbf{a}}_r^{(n)})+\alpha ||{\mathbf{a}}_r^{(n)}|{|}_1\right]$$

$$f(\mathbf{A}):=$$

$$+\sum _{r=1}^R\sum _{n=1}^N\lambda \left[\frac{1-\alpha }{2}(({\mathbf{a}}_r^{(n)}{)}^T{\mathbf{R}}_n^{-1}{\mathbf{a}}_r^{(n)})\right]$$

$$r(\mathbf{A})$$

$$F(\mathbf{A}):=f(\mathbf{A})+r(\mathbf{A})$$

$$F({\mathbf{a}}_1^{(n)}):=$$

$$+\lambda \left[\frac{1-\alpha }{2}(({\mathbf{a}}_1^{(n)}{)}^T{\mathbf{R}}_n^{-1}{\mathbf{a}}_1^{(n)})+\alpha ||{\mathbf{a}}_1^{(n)}|{|}_1\right]$$

$$F({\mathbf{a}}_1^{(n)}):=$$

$$+\lambda \left[\frac{1-\alpha }{2}(({\mathbf{a}}_1^{(n)}{)}^T{\mathbf{R}}_n^{-1}{\mathbf{a}}_1^{(n)})+\alpha ||{\mathbf{a}}_1^{(n)}|{|}_1\right]$$

$$F({\mathbf{a}}_1^{(n)}):=$$

$$+\lambda \left[\frac{1-\alpha }{2}(({\mathbf{a}}_1^{(n)}{)}^T{\mathbf{R}}_n^{-1}{\mathbf{a}}_1^{(n)})+\alpha ||{\mathbf{a}}_1^{(n)}|{|}_1\right]$$

$$F(\mathbf{x}):=$$

$$+\lambda \left[\frac{1-\alpha }{2}{\mathbf{x}}^T\mathbf{T}\mathbf{x}+\alpha ||\mathbf{x}|{|}_1\right]$$

$$F(\mathbf{x})$$

$$+\lambda \left[\frac{1-\alpha }{2}{\mathbf{x}}^T\mathbf{T}\mathbf{x}+\alpha ||\mathbf{x}|{|}_1\right]$$

$${\mathbf{x}}^{\ast }=\arg \underset{\mathbf{x}}{\min }F(\mathbf{x})$$

$$0\in \partial F({\mathbf{x}}^{\ast })$$

$$\partial F(\mathbf{x})$$

$${\bf H}=\left[\begin{array}[]{cccc}||{\bf h}\varoast{\boldsymbol{\delta}}_{1}||_{2}^{2}&&\\ &||{\bf h}\varoast{\boldsymbol{\delta}}_{2}||_{2}^{2}&&\\ &&\ddots&\\ &&&||{\bf h}\varoast{\boldsymbol{\delta}}_{I_{n}}||_{2}^{2}\\ \end{array}\right]$$

$$\mathbf{u}={\left[(\mathbf{h}\text{varoast}{\mathbf{\delta }}_1{)}^T({\mathbf{w}}_1\text{varoast}{\mathbf{\delta }}_1),\dots ,(\mathbf{h}\text{varoast}{\mathbf{\delta }}_{I_n}{)}^T({\mathbf{w}}_{I_n}\text{varoast}{\mathbf{\delta }}_{I_n})\right]}^T$$

$$\text{Sign}(x{)}_i=\{\begin{array}{rlr}& 1& \text{ if }{x}_i>0\\ & [-1,1]& \text{ if }{x}_i=0\\ & -1& \text{ if }{x}_i<0\end{array}$$

$$\mathbf{d}:=\text{diag}(\mathbf{H}+\lambda (1-\alpha )\mathbf{T}):=(d_1,d_2,\dots ,d_{I_n}{)}^T$$

$$0\in \mathbf{x}-\mathbf{u}\text{varoslash}\mathbf{d}+\lambda \alpha \text{Sign}(\mathbf{x})\text{varoslash}\mathbf{d}$$

$$\mathbf{x}={\mathcal{T}}_{\lambda \alpha }(\mathbf{u})\text{varoslash}\mathbf{d}$$

$$\underset{\mathbf{A},\underline{\mathbf{X}}}{\min }$$

Full text

Keywords: PARAFAC decomposition, tensor imputation, elastic net

1 Introduction

Tensor decomposition has been used as a powerful tool on many fields (Kolda and Bader (2009)) for multiway data analysis. Among the representation models for tensor decompositions, PARAFAC decomposition is one of the simplest form. It aims to extract the data-dependent basis functions (i.e. factors), so that the observed tensor can be expressed by the multi-linear combination of these factors, which can reveal meaningful underlying structure of the data. The ideal situation for the PARAFAC decomposition is that we know (1) the true rank of the tensor, and (2) the sparse structure (i.e. locations of zeros) of the factors. Given a wrong estimation of the tensor’s rank or sparse structure of the factors, even the recovery error of the decomposition results are very low, there still exits a high generalization error as the learning procedure tries to fit the noise information. But in most of the real-case scenarios, we won’t be able to obtain the prior knowledge of the true rank of the tensor and the true sparse structure of the factors. Various literatures have proposed solutions for either of the above problems. For example, Liu et al. proposed the sparse non-negative tensor factorization using columnwise coordinate descent to get the sparse factors (Liu et al. (2012)). To get the proper rank estimation, a novel rank regularization PARAFAC decomposition was proposed in (Bazerque et al. (2013)). However, to the author’s best knowledge, there is no work that can handle both of the two problems simultaneously in an integrated framework. And in many works, they try to emphasize on the lower recovery error, few authors focus on finding the true underline factors of the observed data with noise. As we can see from our simulation examples, there are algorithms( e.g. CP-ALS), although they find a lower relative recovery error, but they deviate from the true underline factors. The reason is that these algorithm learn to fit the noise information to reduce the recovery error. Our algorithm is effective to de-noise from both our simulation data and real dataset. Moreover, there seems few researches to use stochastic scheme to calculate the large scale PARAFAC decomposition.

In response to the above problems, we develop a PARAFAC decomposition model capturing the true rank and producing the sparse factors by minimizing the least square lost with elastic net regularization.

The elastic net regularization is a convex combination of the $l_{2}$ regularization and $l_{1}$ regularization, where the $l_{2}$ regularization helps to capture the true rank, while the $l_{1}$ regularization trends to yield the sparse factors. For solving the minimization problem, we use the block coordinate descent algorithm and stochastic block coordinate descent algorithm, where the block coordinate descent algorithm have been shown to be an efficient scheme (Xu and Yin (2013)) especially on large scale problems. To identify the true rank and sparse factors, we generate a solution path by gradually increasing the magnitude of the regularization. And the initial dense decomposition is gradually transformed to a low rank and sparse decomposition. But when we find the rank and the appropriate sparse factors, the recovery error may be large because large regularization will shrink the solution closer to origin. One can fix the sparse structure of the factors and perform a constrained optimization with small regularization to decrease the recovery error. At last, we can choose the appropriate low rank and sparse PARAFAC decomposition from the refined solution of the solution path which features a low recovery error. We test our algorithm firstly on the synthetic data. The results shows that the algorithm captures the true rank with high probability and the sparse structure of the factors found by the algorithm is sufficiently close to the sparse structure of the true factors, and the sparse constraint algorithm with the sparse structure of the factors trends to find the true underline factors. And it finds out the meaningful structure of features when applied on the coil-20 dataset and the COPD data.

The rest of this paper is organized as follows: Our model based on a Bayesian framework is derived in Section 2. The coordinate descent algorithm to solve the optimization is proposed in Section 3. Theoretical analysis is given in Section 4. The solution path method is given in Section 5. The stochastic block coordinate descent algorithm to solve the optimization is proposed in Section 6. Results from applying the algorithm on the synthetic data and real data is summarized in Section 7, and conclude this paper in Section 8.

2 Bayesian PARAFAC model

The notations used through out the article are as follows: We use bold lowercase and capital letters for vectors ${\bf a}$, and matrices ${\bf A}$, respectively. Tensors are underlined, e.g., $\underline{\bf X}$. Both the matrix and tensor Frobenius norms are represented by $||\cdot||_{F}$ and use $\text{vec}(\cdot)$ denote the vectorization of a matrix. Symbols $\varotimes$, $\varodot$, $\varoast$, and $\circ$, denote the Kroneker, Kathri-Rao, Hadamard (entry-wise), and the outer product, respectively. $\varoslash$ denote the entry-wise divide operation, and $[N]:=\{1,2,\ldots,N\}$. $\text{Tr}(\cdot)$ stands for the trace operation.

The observation $\underline{\bf Z}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ is composed of the true solution $\underline{\bf X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ and the noise tensor $\underline{\bf E}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$.

[TABLE]

where $\underline{\bf E}_{i_{1},i_{2},\ldots,i_{N}}\sim\mathcal{N}(0,\sigma^{2}),\;\text{i.i.d.}$. We only observed a subset entries of $\underline{\bf Z}$, which is given by a binary tensor $\underline{\boldsymbol{\Delta}}\in\{0,1\}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$.

Assume that $\underline{\bf X}$ comes from the PARAFAC decomposition $\underline{\bf X}={\bf A}^{(1)}\circ{\bf A}^{(2)}\circ\cdots\circ{\bf A}^{(N)}$, which means that

[TABLE]

where the factor matrix ${\bf A}^{(n)}\in\mathbb{R}^{I_{n}\times R}$, for $n=1,2,\ldots,N$. Let ${\bf a}^{(n)}_{r},\;r=1,\cdots R$ be the columns of ${\bf A}^{(n)}$. Then we may express the decomposition as

[TABLE]

where ${\bf u}_{r}^{(n)}:={\bf a}_{r}^{(n)}/||{\bf a}_{r}^{(n)}||_{2}$, and character value $\gamma_{r}:=\prod_{n=1}^{N}||{\bf a}_{r}^{(n)}||_{2},r=1,\cdots R$. Denote ${\bf U}^{(n)}:=[{\bf u}_{1}^{(n)},{\bf u}_{1}^{(n)},\ldots,{\bf u}_{R}^{(n)}]$, ${\bf U}:=({\bf U}^{(1)},{\bf U}^{(2)},\ldots,{\bf U}^{(N)})$, $\boldsymbol{\gamma}:=(\gamma_{1},\gamma_{2},\ldots,\gamma_{R})$, ${\bf A}:=({\bf A}^{(1)},{\bf A}^{(2)},\ldots,{\bf A}^{(N)})$

For fixed $n$, since vectors ${\bf a}^{(n)}_{r}$ ,$n=1,\cdots N$ are interchangeable, identical distributions, they can be modeled as independent for each other, zero-mean with the density function,

[TABLE]

where $\beta_{n}$ is the normalization constant, and the density is the product of the density of Gaussian (mean [math], covariance ${\bf R}_{n}$) $1/\sqrt{\det(2\pi{\bf R}_{n})}\exp(-1/2{\bf a}^{(n)}_{r}{\bf R}_{n}^{-1}{\bf a}^{(n)}_{r})$ and the double exponential density $(\mu_{n}/2)^{I_{n}}\exp(-\mu_{n}||{\bf a}^{(n)}_{r}||_{1})$.

In addition ${\bf a}^{(n)}_{r},$ for $n=1,2,\ldots,N$ are assumed mutually independent. And since scale ambiguity is inherently present in the PARAFAC model, vectors ${\bf a}^{(n)}_{r},$ for $n=1,2,\ldots,N$ are set to have equal power, that is,

[TABLE]

And for the computation tractable and simplicity, we assume that ${\bf{R}}_{N}$ are diagonal. Under these assumptions, the negative logarithm of the posterior distribution (up to a constant) is

[TABLE]

Correspondingly, the MAP estimator is $\underline{\bf X}={\bf A}^{(1)}\circ{\bf A}^{(2)}\circ\cdots\circ{\bf A}^{(N)}$, where ${\bf A}$ is the solution of

[TABLE]

where $\lambda=\sigma^{2}(1+\mu)$, $\alpha=\frac{\mu}{1+\mu}$. Note that this regularization is elastic net as introduced in (Zou and Hastie (2005)).

3 Block Coordinate Descent for Optimization

In order to solve (7), we derive a coordinate descent-based algorithm, by firstly separating the cost function in (7) to the smooth part and the non-smooth part.

[TABLE]

[TABLE]

[TABLE]

The optimization problem in equation (7) is solved iteratively by updating one part at a time with all other parts fixed. In detail, we cyclically minimize the columns ${\bf a}^{(n)}_{r}$ for $r=1,\ldots,R$ and $n=1,\ldots N$. For example, if we consider do minimization about ${\bf a}_{1}^{(n)}$ and fixed all the other columns of ${\bf A}$, the following subproblem is then obtained:

[TABLE]

Setting $\underline{\bf W}=\underline{\bf Z}-\sum_{r=2}^{R}{\bf a}_{r}^{(1)}\circ{\bf a}_{r}^{(2)}\circ\cdots\circ{\bf a}_{r}^{(N)}$

[TABLE]

To make the calculation easier, we unfold the tensor to the matrix form. Let ${\bf W}_{(n)}\in\mathbb{R}^{I_{n}\times I_{1}I_{2}\cdots I_{n-1}I_{n+1}\cdots I_{N}}$ denote the matrix of unfolding the tensor $\underline{\bf W}$ along its mode $n$. And using the fact $({\bf a}_{1}^{(1)}\circ{\bf a}_{1}^{(2)}\circ\cdots\circ{\bf a}_{1}^{(N)})_{(n)}={\bf a}_{1}^{(n)}{\bf h}^{T}$, where ${\bf h}={\bf a}_{r}^{\varoast_{-n}}:={\bf a}^{(N)}_{r}\varotimes\cdots\varotimes{\bf a}^{(n+1)}_{r}\varotimes{\bf a}^{(n-1)}_{r}\varotimes\cdots\varotimes{\bf a}^{(1)}_{r}$. (12) becomes

[TABLE]

Denote ${\bf x}:={\bf a}^{(n)}_{1}$ , ${\bf T}:={\bf R}_{n}^{-1}$and drop out all subscript and superscript, we can rewrite (13) as

[TABLE]

which can be decomposed as

[TABLE]

where ${\bf w}_{i_{n}}^{T}$, ${\boldsymbol{\delta}}_{i_{n}}^{T}$, represent the ${i_{n}}$-th row of matrices ${\bf W}$, ${\boldsymbol{\Delta}}$, respectively. Note that $F$ is strongly convex at ${\bf x}$, the optimal solution

[TABLE]

satisfies the first order condition

[TABLE]

The subgradient of $F$ at ${\bf x}$ is

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Let

[TABLE]

The optimal condition (17) is actually equivalent to:

[TABLE]

So the solution is

[TABLE]

where $\mathcal{T}_{v}(t):=\text{sign}(t)\max\{|t|-v,0\}$ is the soft thresholding operator.

We formally give the algorithm 1 to solve (7).

4 Theory Analysis

4.1 Convergence results

Since our algorithm is a instance of the unified algorithm 1 in (Xu and Yin (2013)) , the convergence of algorithm 1 follows from Theory 2.8 and Theory 2.9 in (Xu and Yin (2013)) and the proof is given in the Appendix section 9.1. We has the following convergence theorem

Theorem 1

Let $\{{\bf A}^{k}\}$ be the sequence generated by Algorithm 1, where $\{{\bf A}^{k}\}$ is the solution $({\bf A}^{(1)},{\bf A}^{(2)},\ldots,{\bf A}^{(N)})$ in the k-th iteration in the repeat loop. Assume that $\{{\bf A}^{k}\}$ is bound. Then $\{{\bf A}^{k}\}$ converges to a critical point $\bar{\bf A}$, and the asymptotic convergence rates in Theory 2.9 in (Xu and Yin (2013)) apply.

4.2 Property of Elastic Regularization

4.2.1 $l_{2}$ regularization trend to find the true rank

First, the l2 regularization has the property to find the true rank of tensor $\underline{\bf X}$ (Bazerque et al. (2013)). To see this, note that when $\alpha=0$( This can be achieved $\mu=0$), the problem (7) becomes

[TABLE]

And this problem is equivalent to the following problem by the following Proposition 2(its proof is provided in the Appendix section 9.2)

[TABLE]

Proposition 2

The solution of (25) and (26) coincide, i.e. the optimal factors related by ${\bf a}_{r}^{(N)}=(\tilde{\gamma}_{r})^{\frac{1}{N}}{\bf R}_{n}^{1/2}\tilde{\bf u}_{r}^{(n)}$

To further stress the capability of (25) to produce a low-rank approximate tensor, consider transform (26) once more by rewritten it in the constrained-error form

[TABLE]

where $\tilde{\boldsymbol{\gamma}}:=(\tilde{\gamma}_{1},\tilde{\gamma}_{2},\ldots,\tilde{\gamma}_{R})$, $||\tilde{\boldsymbol{\gamma}}||_{\frac{2}{N}}:=(\sum_{r=1}^{R}|\tilde{\gamma}_{r}|^{\frac{2}{N}})^{\frac{N}{2}}$. For any vlue of $\eta$ there exists a corresponding Lagrange multiplier $\lambda$ such that (26) and (27) yield the same solution. And the $l_{\frac{N}{2}}$-norm $||\tilde{\boldsymbol{\gamma}}||_{\frac{2}{N}}$ in (27) produces a sparse vector $\tilde{\boldsymbol{\gamma}}$ when minimized (Chartrand (2007)). Note that the sparsity in vector $\tilde{\boldsymbol{\gamma}}$ implies the low rank of ${\underline{\bf X}}$.

The above arguments imply that when properly choose the parameter $\lambda$, the $l_{2}$ regularization problem (25) can find the true rank of the tensor $\underline{\bf X}$. i.e. when we give a overestimated rank $R>R^{*}$(where $R^{*}$ is the true rank of of tensor $\underline{\bf X}$) of tensor $\underline{\bf X}$, the perfect solution of (25) will shrink the redundant $R-R^{*}$ columns of the factor matrix ${\bf A}$ to [math].

4.2.2 $l_{1}$ regularization trend to find the true sparse structure of the factor matrix

Note that when there is only $l_{1}$ regularization

[TABLE]

For each mode factor matrix ${\bf A}^{(n)}$, it is a standard lasso problem, which implies that the solution $({\bf A}^{(n)})^{*}$ is sparse. Note that if when the true factor matrix $({\bf A}^{(n)})^{*}$ is sparse, with the properly choosed $\lambda$, we can reveal the true sparsity structure in $({\bf A}^{(n)})^{*}$. For application, the sparsity structure of $({\bf A}^{(n)})^{*}$ can help us to make meaningful explanation on its mode $n$, which standards for some attributes of the considered problem.

4.2.3 The elastic net give us a flexible model to find the true data structure in tensor

Combine the $l_{1}$ regularization and $l_{2}$ regularization we get the elastic net regularization, which can combine the advantages of both $l_{1}$ regularization and $l_{2}$ regularization, i.e. to find the true (low) rank and closed to the true sparse factor matrix of PARAFAC decomposition of an observed tensor data. It is helpful to understand the structure of the data and reveal the faces of the objects.

4.2.4 Estimate of the covariance matrix

To run the algorithm 1, the covariance must be postulated as a priori, or replace by their sample estimates. And often we don’t know the priori, so the proper sample estimates are very important, since it provides reasonable scaling in each dimension such that the algorithm performs well. Similarly in section C Covariance estimation in article (Bazerque et al. (2013)), we can bridge the covariance matrix with its kernel counterparts. Since that when the model is give in (4) the covariance of the the factor ${\bf a}^{(n)}_{r}$ is hard to evaluate, do not has an analytical expression. We can jump out the obstacles just assume our model is Gaussian now,

[TABLE]

Define the kernel similarity matrix in mode $n$ as

[TABLE]

i.e. the expectation of the inner product of the $i$-th slice and $j$-th slice of ${\bf X}$ in mode $n$. With some calculation(see in the supplementary material), we can get

[TABLE]

and

[TABLE]

From this we can get the covariance matrix estimate(just drop out the expectation)

[TABLE]

5 Solution Path by Warm Start in Practical Application

In this section, we will give a more practical algorithm which help us how to chose the proper sparse and low rank PARAFAC approximation of the observed tensor with missing values. From the Lagrange theory, the problem (7) is equivalent to

[TABLE]

For any value of $\lambda$ there exists a corresponding Lagrange multiplier $\eta$ such that (7) and (34) yield the same solution. Although there is no explicit equation to relate $\lambda$ and $\eta$, we know that small $\lambda$ implies large $\eta$, large ball-like solution space; large $\lambda$ implies small $\eta$, small ball-like solution space. So we can fixed the relation $\alpha$ of $l_{1}$ and $l_{2}$ regularization. And start from a small $\lambda$ toward to large $\lambda$, which can shrink the dense solution to the sparse solution(when $\lambda$ is very large, the solution is [math]). And this forms a solution path, the solutions are close to each other when $\lambda$ are close to each other, this means when we start from the previous solution, the algorithm can find quickly the close $\lambda$ solution. And for the lucky $\alpha$, we can meet the true data structure, i.e. the right sparse structure in the factor matrix. From our numerical test, the situation is just as we imagined. However, on the time when we find the right data structure, i.e. the right sparse structure in the factor matrix, the recovery error is not so low, since when $\lambda$ is large, the regularization will shrink the solution close to origin, and the best recovery error is achieved when $\lambda$ is small. This phenomenon inspired us a efficient method: using algorithm 1 to generate a solution path, and for a special solution in the solution path, fixed it sparse structure, solve a constrained problem with small $\lambda$ to get small recovery error — for fixed $\alpha$, we use Algorithm 1 to calculate a solution path, i.e. the solution generated by Algorithm 1 from a sequence increasing $\lambda$ (e.g $(10^{-10},10^{-9},\ldots,10^{10})$.

For a special solution ${\bf A}:=({\bf A}^{(1)},{\bf A}^{(2)},\ldots,{\bf A}^{(N)})$ , where ${\bf A}^{(n)}\in\mathbb{R}^{I_{n}\times R}$, calculate ${\bf U}$ and $\boldsymbol{\gamma}$ as in (3), reorder $\boldsymbol{\gamma}$ such that $\gamma_{1}\geq\gamma_{2}\geq\gamma_{R_{1}}>\gamma_{R_{1}+1}=\cdots=\gamma_{R}=0$ , reorder the columns of matrix ${\bf A}^{(n)}$ and ${\bf U}^{(n)}$ such that they coincide with the $\boldsymbol{\gamma}$. Let ${\bf B}:=({\bf B}^{(1)},{\bf B}^{(2)},\ldots,{\bf B}^{(N)})$ , where ${\bf B}^{(n)}\in\mathbb{R}^{I_{n}\times R_{1}}$ is sub-matrix of ${\bf A}^{(n)}$, which is formed by its first $R_{1}$ columns. Definde the sparse structure matrix ${\bf S}:=({\bf S}^{(1)},{\bf S}^{(2)},\ldots,{\bf S}^{(N)})$ , where ${\bf S}^{(n)}\in\mathbb{R}^{I_{n}\times R_{1}}$ of the factor matrix ${\bf B}$ as

[TABLE]

where $i\in[I_{n}]$, $j\in[R_{1}]$ and $\epsilon$ is a prespecified small number( e.g. $\epsilon=10^{-9}$).

For notation simplicity, we still use ${\bf A}$ denote ${\bf B}$ and $R$ denote $R_{1}$. If the true sparse structure matrix is ${\bf S}$, to solve (7), it is equivalent to solve the constrained problem,

[TABLE]

This problem can solved by algorithm 1 just renew ${\bf A}^{(n)}(i_{n},r)$ if ${\bf S}^{(n)}(i_{n},r)=1$, and set all ${\bf A}^{(n)}(i_{n},r)=0$ if ${\bf S}^{(n)}(i_{n},r)=0$. We formally give the Algorithm 2.

And the convergence theory still holds since we only update a subset elements of factor matrix ${\bf A}$.

To conclude, we give our finally algorithm 3 aim to find the true data structure and low recovery error. For convenience, we call the solution path method using algorithm 1 SPML, and the correspond using algorithm 2 SPMS where the solution comes from SPML.

6 Generalize the algorithm to large scale problem by using stochastic block coordinate descent algorithm by using admax method

For the real application, the data often have a large data dimnsion. the above algorithm may not compute fast, since for each iteration, we have $O(I_{1}*I_{2}*\cdots*I_{N})$ cost of operations. Fortunately, we can fixed this problem to some extent by using a stochastic scheme.

Note that for renew ${\bf a}^{(n)}_{r}$, we use all information to renew it, it it not necessary. To gain the insights. Now consider the problem ${\bf x}=a{\bf b}$, where ${\bf x}\in\mathbb{R}^{M\times 1}$ , $a\in R$ and ${\bf b}\in\mathbb{R}^{M\times 1}$. Now suppose we know the true ${\bf b}$, and we need to update $a$, the form the above algorithm, we renew $a$ by the formula(Now we drop out the regularity and suppose their is no noise and no missing entries) $a=\frac{\langle{\bf x},{\bf b}\rangle}{\langle{\bf b},{\bf b}\rangle}$. Note that we can also update $a$ by random choose a subset ${\bf s}\subseteq\{1,2,\ldots,M\}$, and update $a=\frac{\langle{\bf x_{\bf s}},{\bf b_{\bf s}}\rangle}{\langle{\bf b_{\bf s}},{\bf b_{\bf s}}\rangle}$. While when there are noise, such that we observed ${\bf z}={\bf x}+{\bf e}$, $\tilde{a}=\frac{\langle{\bf z},{\bf b}\rangle}{\langle{\bf b},{\bf b}\rangle}=\frac{\langle{\bf x}+{\bf e},{\bf b}\rangle}{\langle{\bf b},{\bf b}\rangle}=a+\frac{\langle{\bf e},{\bf b}\rangle}{\langle{\bf b},{\bf b}\rangle}$. Using a subset update $\hat{a}=\frac{\langle{\bf z}_{\bf s},{\bf b}_{\bf s}\rangle}{\langle{\bf b}_{\bf s},{\bf b}_{\bf s}\rangle}=\frac{\langle{\bf x}_{\bf s}+{\bf e}_{\bf s},{\bf b}_{\bf s}\rangle}{\langle{\bf b}_{\bf s},{\bf b}_{\bf s}\rangle}=a+\frac{\langle{\bf e}_{\bf s},{\bf b}_{\bf s}\rangle}{\langle{\bf b}_{\bf s},{\bf b}_{\bf s}\rangle}$. Suppose that $e\sim\mathcal{N}(0,\sigma^{2}{\bf I})$ is a white noise. then $\mathbb{E}\tilde{a}=\mathbb{E}\hat{a}$, and $\text{Var}(\tilde{a})=\sigma^{2}/||{\bf b}||_{2}^{2}$, $\text{Var}(\hat{a})=\sigma^{2}/||{\bf b}_{\bf s}||_{2}^{2}$. From this simple case, it inspire us that we can use a subset of ${\bf b}$ to update $a$, because it will lose some accuracy, higher variance, we should it a stochastic update scheme to renew $a$. Its key idea is keep the memory of the older information of the right direction and also take notice of the current of new information.

Now we begin to deduce our stochastic block coordinate descent method by using the Adamax scheme. Follow the same roads in the Section 3.

For example, if we consider do minimization about ${\bf a}_{1}^{(n)}$, fixed all the other columns of ${\bf A}$, and we choose subset ${\bf s}_{k}\subseteq\{1,2,\ldots,I_{k}\}$ for each mode $k=1,2,\ldots,N$. Denote

[TABLE]

the following subproblem is then obtained:

[TABLE]

Setting $\underline{\bf W}_{n}=\underline{\bf Z}_{n}-\sum_{r=2}^{R}{\bf a}_{r}^{(1)}({\bf s}_{1})\circ\cdots\circ{\bf a}_{r}^{(n-1)}({\bf s}_{n-1})\circ{\bf a}_{r}^{(n)}\circ{\bf a}_{r}^{(n+1)}({\bf s}_{n+1})\circ\cdots\circ{\bf a}_{r}^{(N)}({\bf s}_{N})$

[TABLE]

To make the calculation easier, we unfold the tensor to the matrix form. Let ${\bf W}_{(sn)}\in\mathbb{R}^{I_{n}\times|{\bf s_{1}}||{\bf s_{2}}|\cdots|{\bf s_{n-1}}||{\bf s_{n+1}}|\cdots|{\bf s_{N}}|}$ denote the matrix of unfolding the tensor $\underline{\bf W}_{n}$ along its mode $n$. And using the fact ${\bf a}_{1}^{(1)}({\bf s}_{1})\circ\cdots\circ{\bf a}_{1}^{(n-1)}({\bf s}_{n-1})\circ{\bf a}_{1}^{(n)}\circ{\bf a}_{1}^{(n+1)}({\bf s}_{n+1})\circ\cdots\circ{\bf a}_{1}^{(N)}({\bf s}_{N})={\bf a}_{1}^{(n)}{\bf h_{s}}^{T}$, where ${\bf h}_{s}:={\bf a}^{(N)}_{1}({\bf s}_{N})\varotimes\cdots\varotimes{\bf a}^{(n+1)}_{1}({\bf s}_{n+1})\varotimes{\bf a}^{(n-1)}_{1}({\bf s}_{n-1})\varotimes\cdots\varotimes{\bf a}^{(1)}_{1}({\bf s}_{1})$. (40) becomes

[TABLE]

Denote ${\bf x}_{s}:={\bf a}^{(n)}_{1}$ , ${\bf T}:={\bf R}_{n}^{-1}$and drop out all subscript and superscript, and to denote it is update by the subset information, we add a subscript $s$ for each symbol, we can rewrite (41) as

[TABLE]

which can be decomposed as

[TABLE]

where ${\bf w}_{s,i_{n}}^{T}$, ${\boldsymbol{\delta}}_{s,i_{n}}^{T}$, represent the ${i_{n}}$-th row of matrices ${\bf W}_{s}$, ${\boldsymbol{\Delta}}_{s}$, respectively. Note that $F$ is strongly convex at ${\bf x}_{s}$, the optimal solution

[TABLE]

satisfies the first order condition

[TABLE]

The subgradient of $F$ at ${\bf x}$ is

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Let

[TABLE]

The optimal condition (45) is actually equivalent to:

[TABLE]

So the solution is

[TABLE]

where $\mathcal{T}_{v}(t):=\text{sign}(t)\max\{|t|-v,0\}$ is the soft thresholding operator.

Note that although we solve for ${\bf a}_{1}^{(n)}$ using a subset information, however, we want it has the approximate the same effect as the original solution. Note that from equation (24) and (52), ${\bf u}_{s}\varoslash{\bf d}_{s}\approx{\bf u}\varoslash{\bf d}$ is no problem, however ${\lambda\alpha}\varoslash{\bf d}_{s}\not\approx{\lambda\alpha}\varoslash{\bf d}$. So the correct update should be like

[TABLE]

Since we do not have ${\bf d}$ in the subset update form, we now deduce an approximate formula for equation (53). Note that when $\underline{{\boldsymbol{\Delta}}}=\underline{\bf 1}$, ${\bf d}=||{\bf h}||_{2}^{2}{\bf 1}+\lambda(1-\alpha)\alpha\text{diag}({\bf T})$ and ${\bf d}_{s}=||{\bf h}_{s}||_{2}^{2}{\bf 1}+\lambda(1-\alpha)\alpha\text{diag}({\bf T})$, so ${\bf d}\approx\text{mean}(\frac{||{\bf h}||_{2}^{2}{\bf 1}+\lambda(1-\alpha)\alpha\text{diag}({\bf T})}{{\bf d}_{s}}){\bf d}_{s}$. Then the update can be approximated by

[TABLE]

Where the $||{\bf h}||_{2}^{2}=\prod_{i=1,i\neq n}^{N}||{\bf a}_{1}^{(i)}||_{2}^{2}$.

Note that for the Admax algorithm(Kingma and Ba (2014)), it calculate the gradient ${\bf g}_{1}^{(n)}$ of $F$ at ${\bf a}_{1}^{(n)(t-1)}$, then begin to update, but now we have the true solution of $\min_{{\bf a}_{1}^{(n)}}F({\bf a}_{1}^{(n)}$, we just set ${\bf g}_{1}^{(n)}={\bf a}_{1}^{(n)(t-1)}-{\bf x}_{s}$

We formally give the algorithm 4 to solve (7). For convenience, we call the solution path method using algorithm 4 SPMLR, and the correspond using algorithm 2 SPMSR where the solution comes from SPMLR.

The convergence theorem for Algorithm 4 is not clear. It belong to the class of the multiconvex stochastic scheme, we make the following conjecture.

Conjecture 3

Let $\{{\bf A}^{k}\}$ be the sequence generated by Algorithm 4, where $\{{\bf A}^{k}\}$ is the solution $({\bf A}^{(1)},{\bf A}^{(2)},\ldots,{\bf A}^{(N)})$ in the k-th iteration in the repeat loop. Assume that $\{{\bf A}^{k}\}$ is bound. Then $\{{\bf A}^{k}\}$ converges to a critical point $\bar{\bf A}$.

7 Numerical Test

7.1 Simulated Data

In this work, we generate the synthetic dataset to test the performance of the proposed algorithm. Without loss of generality, the simulation is done on the 3-way tensor data where $I:=I_{1}=I_{2}=I_{3}$ using the Bayesian model as described in section 2. For simplicity, we assume that ${\bf R}_{n}$ is a diagonal matrix, where its diagonal comes from uniform distribution. i.e.

[TABLE]

Note that when the ${\bf R}_{n}$ is diagonal,

[TABLE]

We can sample from the distribution by the rejection sample method. To produce the sparse factors, we set a gate $g$. When sampling out the factor matrix ${\bf A}$, we make it sparse by letting all its elements whose absolute value less than gate $g$ to [math].

7.1.1 Solution path

We firstly use a simple example to investigate the behavior of the solution path, by setting $\mu=0.1$ and $g=0.5$. $(I_{1},I_{2},I_{3}):=(6,6,6)$. $R=2$ and the true number of zeros in ${\bf A}$ is $12$. And we test our algorithm with initial estimate rank $R_{\text{est}}=3$ , $\alpha=0.2$, and the covariance ${\bf R}$ is estimated from equation (33) (where we only estimate the diagonal of ${\bf R}$, if one diagonal of ${\bf R}$ less than $1e-8$, we just set it to $1e-8$ to ensure stability). The initial estimate matrix ${\bf A}_{0}$ comes from standard Gaussian noise. The numbers of max iterations in Algorithm 1 and Algorithm 2 are $200$ and $100$, respectively.

In table 1 (the full table can be found in the Appendix Section 9.4), we see that as $\lambda$ increases, the algorithm will shrink the dense solution to [math], and decrease the rank $3$ to [math]. For all the the resulting factors, their sparsity increases and rank decreases along the solution path as more elements in the factors shrink to 0. When the true rank $2$ is found, the sparsity structure of the synthetic data is also found by the algorithm gradually, as indicated by the increasing number of zero elements found in the current solution and in the true solution, as shown in column NZS and NZT. When $\lambda=500$, the algorithm 1 finds the true sparse structure of the true solution ${\bf A}$. However the relative recovery error using $\lambda=500$ is quite high ($3.78e-01$). While Algorithm 2 can decrease the relative recovery error to $1.62e-6$ (find the true solution ${\bf A}$) in just 9 iterations.

We observed that when $\lambda$ is small, the algorithm converges very fast. But when $\lambda$ becomes large, the iterations will exceed the max number of iterations we set despite the fact that we used the warm start strategy. Yet it still contains the useful information for the true solution and makes it possible to decrease the relative recovery error from $3.78e-01$ to $1.62e-6$ only in 9 iterations.

7.1.2 Model Performance with Noise and No Missing Values

In this section, we test our algorithm with some baseline algorithms. The first algorithm is the CP-ALS method(which is provided by the Tensor Toolbox for MATLAB (Kolda and Bader (2006)) , and we reimplement it in python), and the second algorithm is the LRTI algorithm in (Bazerque et al. (2013)) (where the parameter $\mu$ in LRTI algorithm chosen by $\mu=0.01*\mu_{\text{max}}$ ).

To compare the numerical results, we use standard score metric to measure how well the ground truth is recovered by a CP decomposition in those case where the true factors are known (Tomasi and Bro (2006)). The score between two rank-one tensors $\underline{\bf X}={\bf a}\circ{\bf b}\circ{\bf c}$ and $\underline{\bf Y}={\bf p}\circ{\bf q}\circ{\bf r}$ is defined as:

[TABLE]

For $R>1$, we first sort the components such that $\gamma_{1}\geq\gamma_{2}\geq\cdots\geq\gamma_{R}$, then average the scores for all pairs of components.

After the above simple synthetic data experiment, we further design a more complex synthetic scheme with random noise added. Let $(I_{1},I_{2},I_{3})=(20,20,20)$, $R=10$, the number of test times is set to $50$, and the initial estimate of rank is $10$, and $\mu=0.1$, $g=0.5$. We add noise to the data with the SNR of 20dB, where the SNR is defined by $\text{SNR}:=10\log_{10}\frac{\text{Var}(\underline{\bf X})}{\text{Var}(\underline{\bf E})}$, and $\sigma^{2}=\text{Var}(\underline{\bf E})$. And the is no missing values $\underline{\boldsymbol{\Delta}}=\underline{\bf 1}$, since the CP-ALS algorithm cannot not hander the missing value entries.

For each test in our solution path method( SPML, SPMS. we don’t test on SPMLR, SPMSR algorithms, since the dimension is only $20$. but the solution are chosen by the same procedure in the next subsection), we select the solution in the solution according to the following steps: (1) Choose the estimate rank be the smallest found rank which not less than the true rank; (2) For those solutions whose rank is the estimate rank found by (1), choose the last solution which satisfies that the zeros in the current solution are also the zeros in the true solution, and the number of zeros in the current solution is not greater than the true number of zeros, if no such solution, go to step (3); (3) For those solutions whose rank is the estimate rank found by (1), choose the last solution which the number of zeros in the current solution less than the true number of zeros. (4) If not such a solution, we can choose a best solution from the path by ourself.

And for each algorithm, we use two ways to initialize the ${\bf A}$, the first one is to use the random standard Gaussian $\mathcal{N}(0,1)$, we call it random initialzation, and the second is to set ${\bf A}^{(n)}$ to be the leading $R$ left singular vectors of the mode-n unfolding, ${\bf Z}_{(n)}\varoast{\boldsymbol{\Delta}}_{(n)}$, and we call this nvecs initialization. Now we give our test result for both random initialization and nvecs initialization.

The results are summarized in the table 2.

And we also plot the scores for each algorithm with the random initialization in figure 1.

It can be seen that mean score value of our algorithm is better than the CP-ALS and LRTI algorithm. And we find our algorithm can find the sparse structure of the factor matrix. For the mean of the full data relative error $\frac{||\underline{\bf Z}-\underline{\bf X}||_{F}}{||\underline{\bf Z}||_{F}}$ subtract the true relative error $\frac{||\underline{\bf E}||_{F}}{||\underline{\bf Z}||_{F}}$, our method is as good as the CP-ALS method, and is better than the LRTI algorithm. and SPMS(0.2) with random initialization achieves the best mean scores 0.967. And the most spare solution is find by the algorithm SPMS(0.98) with random initialization. We find that the random initialization for our algorithm can achieve the better mean scores, while for the LRTI and CP-ALS algorithm, the nvecs initialization can achieve better mean scores. For the test time, that our algorithm may some costly, however, note that our algorithm calculate the whole solution path, the cost is relative small.

7.1.3 Model Performance with Noise and Missing Values

To gain more confidence for our mehthod, now we test on the synthetic data with noise and missing. Let $(I_{1},I_{2},I_{3})=(50,50,50)$, $R=9$, the number of test times is set to $30$, and the initial estimate of rank is $10$. We add noise to the data with the SNR of 20dB, where the SNR is defined by $\text{SNR}:=10\log_{10}\frac{\text{Var}(\underline{\bf X})}{\text{Var}(\underline{\bf E})}$, and $\sigma^{2}=\text{Var}(\underline{\bf E})$. And there are missing values(random drop out 25% entries). Since the CP-ALS algorithm cannot not hander the missing value entries, we only compare our method(both the solution path method(SPML, SPMS) and solution path method random(SPMLR , SPMSR)) with LRTI algorithm.

The results are summarized in the table 3.

And we also plot the scores for each algorithm in figure 2.

It can be seen that mean score value of our algorithm is better than the LRTI algorithm. And we find our algorithm can find the sparse structure of the factor matrix. For the mean of the full data relative error $\frac{||\underline{\bf Z}-\underline{\bf X}||_{F}}{||\underline{\bf Z}||_{F}}$ subtract the true relative error $\frac{||\underline{\bf E}||_{F}}{||\underline{\bf Z}||_{F}}$, our method is also better than the LRTI algorithm. Many of our algorithm find out the true factors the score value are 1.0. And the random initialization is better than the nvecs initialization for both our algorithm and the LRTI algorithm. Most choices of our algorithm find out the true rank perfectly, while LRTI only find 10 times of true rank in 30 testes. For the test time, that our algorithm may some costly, however, note that our algorithm calculate the whole solution path, the cost is relative small.

7.2 Application on Real Data

7.2.1 Coil-20 Data

In this section, we test the proposed algorithm on the coil-20 data (S. A. Nene and H. (1996)). This database contains 1440 images of 20 objects, each image is $128\times 128$ pixels, and each object is captured from varying angles with a 5-degree interval. We generated a tensor $\underline{\bf Z}$ with $128\times 128\times 1440$ using these images. And we compare our algorithms (SPMLR, SPMSR) with the CP-ALS and LRTI algorithm(Not given the result of the LRTI algorithm, since it behaves poor, the relative error is more than $0.9$). We don’t use the (SPML, SPMS) because it a little slow and the performance of SPMLR, SPMSR are as good as SPML, SPMS. And we use the two t-SNE (Hinton (2008)) components of ${\bf A}^{(3)}$ for visualization and clustering. The K-means algorithm was adopted for clustering. As K-means is prone to be affected by initial cluster centers, in each run we repeated clustering 20 times, each with the init ’k-means++’(see for the help for K-means method in scikit-learn(Pedregosa et al. (2011)) ). The performance averaged over 20 Monte-Carlo runs is detailed in Table 4. We give the t-SNE visualize of CP-ALS(R) and SPMSR(0.8N) when initial rank is 20 in Figure 3, and the t-SNE visualize of CP-ALS(N) and SPMSR(0.8R) when initial rank is 20 in Figure 4. From these results, we find that algorithm is good for find the underline true factors, although the relative recovery error is a little high than the the CP-ALS method, however, it can achieve a better cluster result (from both the table and the figure) than the CP-ALS method. Which that algorithm is trend to find the true factors.

7.2.2 COPD Data

In this work, we test the proposed algorithm on the baseline data from The Genetic Epidemiology of COPD (COPDGene*®*) cohort (Regan et al. (2010)). The data record consists of 10,300 subjects and 361 features for each subject, including self-administered questionnaires of demographic data and medical history, symptoms, medical record review, etc. All features are independently reviewed by certified professionals (the full data collection forms can be found www.COPDGene.org).

For the data preprocessing, we firstly remove the features which cannot be quantified, then discard the features with missing data covering more than 20 percent of the whole record, thus finally form a table of $213$ features. We then normalize the data by $\frac{{\bf x}-\min({\bf x})}{\max({\bf x})-\min({\bf x})}$ for each feature ${\bf x}$ except the ”GOLD value” feature which grades the COPD stages. As there are $7$ GOLD values( $-2,-1,\ldots,4$) for $7$ stages and the number of subjects within each stage is different (minimum is $108$), we perform subsampling on each stage for balancing the data.

To reveal the features which are closely related to the severity of COPD, we transform the 2D (subject by feature) table into a 3-way tensor $\underline{\bf Z}$, whose dimension is $(7,108,212)$, where the first dimension is the GOLD values, the second dimension is the subjects, and the third dimension is the features. We also keep records for the location of the missing values in tensor $\underline{\boldsymbol{\Delta}}$.

By applying Algorithm 3 on the 3-way tensor above, we generate a solution path and select the best solution from it. The rank of the best solution is $328$ and the relative recovery error is $0.137$. The sparseness (the number of zeros elements found / total number of elements in the factor matrices) is $0.26$. As shown in the visualization of the three factor matrices in Figure 5, 6, 7, all of them have sparse structure. We further calculate the sum of the top $6$ columns (which correspond the largest $6$ character values in ${\boldsymbol{\gamma}}$) of the third features’s factor matrix to get a vector of $212$ elements and sort them by their absolute value (see in Figure 8). The top $10$ features are EverSmokedCig, Blood_Other_Use, FEV1_FVC_utah, pre_FEV1_FVC, Resting_SaO2, ATS_ERS, FEV1pp_utah, SmokCigNow, HealthStatus, HighBloodPres. According to domain experts’ experience, most of these features are closely related the severity of COPD.

8 Conclusion

In this work, we derive a tensor decomposition model from the Bayesian framework to learn the sparse and low rank PARAFAC decomposition via elastic net regularization. An efficient block coordinate descent algorithm and stochastic block coordinate descent algorithm are proposed to solve the optimization problem, which can be applied to solve large scale problems. To increase the robustness of the algorithm for practical applications, we develop a solution path strategy with warm start to generate a gradually sparse and low rank solution, and we reduce the relative recovery error by using the sparse constrained algorithm 2. Evaluation on synthetic data shows our algorithm can capture the true rank and find a close approximation of the true true sparse structure with high probability and with the help of sparse structure, the constrained algorithm trends to find the true underline factors. Evaluation on the coil-20 data shows that our algorithm can extract the meaningful underline factors that can have better cluster performance. Evaluation on the real data shows that it can extract the sparse factors that reveal meaningful relationships between the features in the data and the severity of COPD.

In the future, we plan to expand the current algorithm to the N-way nonnegative PARAFAC decomposition in a similar structure.

9 Appendix

9.1 Proof of Theorem 1

**Proof ** From equation (8), (9), (10), our notations coincide with the notation in (Xu and Yin (2013)), (where in our algorithm, the block $i$ is just a column ${\bf a}^{(n)}_{r}$, for $r=1,2,\ldots,R,n=1,2,\ldots,N)$. Since our algorithm 1 is a special case of Algorithm 1 in (Xu and Yin (2013)) in which we only use the update rule (1.3a) in (Xu and Yin (2013)) , we only need to verify the assumptions in Theory 2.8 and Theory 2.9 in (Xu and Yin (2013)).

For the assumption 1 in (Xu and Yin (2013)): Obviously, $F$ is continuous in $\text{dom}(F)$ and $\inf_{\bf{A}}F(\bf{A})\geq 0$, and $F(\bf{A})$ has a Nash point (See (2.3) in (Xu and Yin (2013)) for definition).

For the assumption 2 in (Xu and Yin (2013)) :

[TABLE]

From the equation (18), it is easy to see that $\nabla^{2}_{{\bf a}^{(n)}_{r}}f({\bf a}^{(n)}_{r})={\bf H}+\lambda(1-\alpha){\bf R}_{n}^{-1}$, where H is defined by (19) with ${\bf h}={\bf a}_{r}^{{\varoast}_{-n}}:={\bf a}^{(N)}_{r}\varotimes\cdots\varotimes{\bf a}^{(n+1)}_{r}\varotimes{\bf a}^{(n-1)}_{r}\varotimes\cdots\varotimes{\bf a}^{(1)}_{r}$, and ${\boldsymbol{\delta}}_{i_{n}}^{T}$ is the $i_{n}$-th row of ${\boldsymbol{\Delta}}_{(n)}$. So $f({\bf a}^{(n)}_{r})$ is strongly convex with modulus $l{\bf I}\leq\nabla^{2}_{{\bf A}^{(n)}_{i_{n},r}}f({\bf A}^{(n)}_{i_{n},r})\leq L{\bf I}$, where $l:=\lambda(1-\alpha)\underset{n=1,2,\ldots,N,i_{n}=1,2,\ldots,I_{N}}{\min}{\bf R}_{n}^{-1}(i_{n},i_{n})>0$, since $\alpha=\frac{\mu}{1+\mu}<1$. And $L$ exist since we assume that $\{{\bf A}^{k}\}$ is bounded.

For the conditions in Lemma 2.6 in (Xu and Yin (2013)):

From the verification of the assumption 2 in (Xu and Yin (2013)) above, it is easy to see that $\nabla f$ is Lipschitz continuous on any bounded set 2. 2.

$F$ satisfies the KL inequility (2.14) in (Xu and Yin (2013)) at $\bar{\bf A}$ since $f$ satisfies the KL inequility (2.14) and $r$ satisfies KL inequility (2.14), so does their sum. 3. 3.

We may choose the initial estimate ${\bf A}$ sufficiently close to $\bar{\bf A}$, and $F({\bf A}^{k})>F(\bar{\bf A})$ for $k>0$, since $F({\bf a}^{(n)}_{r})$ is strongly convex in ${\bf a}^{(n)}_{r}$, so $F({\bf A}^{k})$ is strictly decreasing.

All the conditions in Theory 2.8 and Theory 2.9 in (Xu and Yin (2013)) are satisfied, so the conclusions follow out.

9.2 Proof of Proposition 2

**Proof ** Let $\tilde{\bf A}:=(\tilde{\bf A}^{(1)},\tilde{\bf A}^{(2)},\ldots,\tilde{\bf A}^{(N)})$, and $\tilde{\bf A}^{(n)}={\bf R}_{n}^{-1/2}{\bf A}^{(n)}$, $\tilde{\bf a}_{r}^{(n)}:={\bf R}_{n}^{-1/2}{\bf a}_{r}^{(n)}$, then we can rewrite (25) as

[TABLE]

We can standardize $\tilde{\bf a}_{r}^{(n)}$ to unit length, $\tilde{\bf u}_{r}^{(n)}:=\frac{\tilde{\bf a}_{r}^{(n)}}{||\tilde{\bf a}_{r}^{(n)}||_{2}}$, $\tilde{\bf U}^{(n)}:=(\tilde{\bf u}_{r}^{(1)},\tilde{\bf u}_{r}^{(2)},\ldots,\tilde{\bf u}_{r}^{(N)})$, $\tilde{\gamma}_{r}:=\prod_{n=1}^{N}||\tilde{\bf a}_{r}^{(n)}||_{2},\;r=1,2,\ldots R$, and let ${a}_{r}^{(n)}:=||\tilde{\bf a}_{r}^{(n)}||_{2}$. Then (59) is equivalent to

[TABLE]

Focus on the inner minimization w.r.t. norms $\{\tilde{a}^{(n)}_{r}\}$ for arbitrary directions $\{\tilde{\bf a}^{(n)}_{r}\}$ and fixed products $\tilde{\gamma}_{r}$ and $\underline{\bf X}$ , the (60) is equivalent to

[TABLE]

The arithmetic geometric-mean inequality gives the solution to (61), as it states that for scalars $\{\tilde{a}^{(n)}_{r}\}$, it holds that

[TABLE]

with equality when $\tilde{a}^{(1)}_{r}=\tilde{a}^{(2)}_{r}=\cdots=\tilde{a}^{(N)}_{r}$ so that the minimum of (61) is attained at $\tilde{a}^{(1)}_{r}=\tilde{a}^{(2)}_{r}=\cdots=\tilde{a}^{(N)}_{r}=({\tilde{\gamma}_{r}})^{\frac{1}{N}}$.

Then (60) is equivalent to

[TABLE]

So the solution of (25) and (63) coincide, i.e. the optimal factors related by ${\bf a}_{r}^{(N)}=(\tilde{\gamma}_{r})^{\frac{1}{N}}{\bf R}_{n}^{1/2}\tilde{\bf u}_{r}^{(n)}$

9.3 The relation of kernel similarity matrix and covariance matrix

Note that ${\bf X}_{(n)}={\bf A}_{n}({\bf A}_{N}\varodot\cdots\varodot{\bf A}_{n+1}\varodot{\bf A}_{n-1}\varodot\cdots\varodot{\bf A}_{1})^{T}$, let ${\bf B}_{n}:={\bf A}_{N}\varodot\cdots\varodot{\bf A}_{n+1}\varodot{\bf A}_{n-1}\varodot\cdots\varodot{\bf A}_{1}$, then

[TABLE]

Let $\theta_{n}:=\text{Tr}({\bf R}_{n})$. Then

[TABLE]

[TABLE]

Combine equation (30 ), (64 ), (65 ), (66 ), we get

[TABLE]

and

[TABLE]

Since from our assumption (5), $\theta_{n}=\theta,\;n=1,2,\ldots,N$, we get

[TABLE]

and

[TABLE]

9.4 Full Solution Path

We give the full solution path in table 5.