Non-negative Factorization of the
Occurrence Tensor from Financial Contracts

arXiv:1612.03350v1 [] 10 Dec 2016

1

Zheng Xu1 , Furong Huang3 , Louiqa Raschid2 , Tom Goldstein1
Department of Computer Science, 2 Smith School of Business, University of Maryland
3
Microsoft Research, New York City

Abstract
We propose an algorithm for the non-negative factorization of an occurrence tensor built
from heterogeneous networks. We use `0 norm to model sparse errors over discrete values
(occurrences), and use decomposed factors to model the embedded groups of nodes. An
efficient splitting method is developed to optimize the nonconvex and nonsmooth objective.
We study both synthetic problems and a new dataset built from financial documents, resMBS.

1

Introduction

Tensor factorization is a powerful approach to a myriad of unsupervised learning problems [1, 14, 17].
We propose an efficient algorithm called non-negative occurrence tensor factorization (NOTF) for
analyzing heterogeneous networks. As an application of NOTF, we study the resMBS dataset[5, 19],
which consists of the relationships that financial institutions (FIs, e.g., Bank of America) play roles
(e.g., issuer) in financial contracts (FCs). ResMBS is automatically extracted from a collection of
public financial contracts . We can represent resMBS as a three-mode (FC, FI, Role) occurrence tensor
with non-negative discrete values corresponding to the confidence of the relationship extraction.
The occurrence tensor has several characteristics. First, the tensor has positive discrete values.
Second, the observed tensor contains sparse noise corresponding to errors caused by inaccurate
extraction. Most importantly, we expect that FIs will play specific roles across multiple FCs to create
a community. The observed occurrence tensor could then be decomposed into a low-rank tensor
together with sparse noise. This decomposition can be used to clean and complete the observation,
and also allow us to understand the underlying (FC, FI, Role) communities.
The proposed NOTF extends the CANDECOMP/PARAFAC (CP) decomposition [6, 12] of a tensor.
After decomposing the occurrence tensor, each rank-one tensor component represents a (FC, FI, Role)
community. Hence the decomposed factors are constrained to be non-negative. Instead of the `2 norm
for factorization of real valued tensors with Gaussian noise, the `0 norm is considered for the discrete
values and the sparse errors. We develop an algorithm based on the alternating direction method of
multipliers (ADMM) [4, 22] to optimize the nonconvex and nonsmooth objective of NOTF.
We briefly discuss some related research. While tensor decomposition dates back to the 1920s, robust
tensor decomposition solutions have been presented in some recent papers [14, 2, 8, 11]. Standard
CP decomposition has been applied to heterogenous networks [16]. Discrete valued tensors have
been studied in [18]. ADMM has been extensively used as a solver for robust tensor recovery [8, 11],
and also for standard non-negative tensor factorization[15]. Non-`2 norms have been used for CP
decomposition of real valued tensors in computer vision [13, 7]. ADMM for `0 norm was empirically
studied in [21]. Probabilistic communities were discussed for financial documents in [19].

2

Non-negative Occurrence Tensor Factorization

We use notations similar to [14]. Vectors and matrices are denoted by lowercase and capital letters,
respectively. Higher mode tensors are denoted by Euler script letters. Three mode tensors X ∈
RN1 ×N2 ×N3 are used as an example in this paper. Fibers are column vectors extracted from tensors
by fixing every index but one, e.g., x:jk , xi:k , xij: . Mode-d fibers are arranged to get the mode-d

unfolding matrix X(d) , e.g., X(1) = [x:jk ] ∈ RN1 ×N2 N3 . Denote the vector outer product by ◦, then
PR
the CP decomposition factorizes a tensor into a sum of rank-one tensors as X = r=1 ar ◦ br ◦ cr . If
we denote A = [ar ] ∈ RN1 ×R , B = [br ] ∈ RN2 ×R and C = [cr ] ∈ RN3 ×R , the Kronecker product
A ⊗ B ∈ RN1 N2 ×RR , and the Khatri-Rao product A B = [ar ⊗ br ] ∈ RN1 N2 ×R , then we can
PR
compactly represent CP as X = r=1 ar ◦ br ◦ cr = [[A, B, C]], which is equivalent to the unfolding
representation X(1) = A(C B)T , X(2) = B(C A)T , and X(3) = C(B A)T .
We are seeking to recover a tensor X = [[A, B, C]] that is at most rank R from the noisy observations
O. Each rank-one tensor ar ◦ br ◦ cr captures a community and ar ≥ 0, br ≥ 0, cr ≥ 0 represent the
weights of nodes (e.g., FI, FC, and roles) in the community. For discrete valued tensors, we minimize
the sparse error measured with the `0 norm rather than the `2 loss in standard decomposition,
min k[[A, B, C]] − Ok0 , subject to A ≥ 0, B ≥ 0, C ≥ 0,

A,B,C

(1)

P
where kXk0 = i,j,k 11{z: |z|>0} (xijk ) is the counts of nonzero values in a tensor, 11S is the indicator
function of the set S: 11S (v) = 1, if v ∈ S, and 11S (v) = 0, otherwise.
We minimize the NOTF objective (1) by introducing an intermediate variable U,
min kUk0 + ι{z: z≥0} (A, B, C), subject to U = [[A, B, C]] − O,

U,A,B,C

(2)

where ιS is the characteristic function of the set S; ιS (v) = 0, if v ∈ S, and ιS (v) = ∞, otherwise.
We then apply alternating direction method of multipliers (ADMM) [9, 4, 10, 22] by introducing dual
variables λ and alternatively solving subproblems of U and A, B, C, with p indexing iterations,
τ
2
U − [[Ap , B p , C p ]] + O + λp F
(3)
Up+1 = arg min kUk0 +
U
2
τ
2
Ap+1 , B p+1 , C p+1 = arg min ι{z: z≥0} (A, B, C) +
Up+1 − [[A, B, C]] + O + λp F (4)
A,B,C
2
λp+1 = λp + Up+1 − [[Ap+1 , B p+1 , C p+1 ]] + O,
(5)
qP
2
where kXkF =
i,j,k xijk , and τ is a hyperparameter called the penalty parameter.
Subproblem (3) can be solved by the proximal operator of the `0 norm, known as hard-thresholding,
Up+1 = hard([[Ap , B p , C p ]] − O − λp , 1/τ ),
arg minX kXk0 + 1/2tkX − Zk2F

(6)

where hard(Z, t) =
= Z ∗ I{z:|z|> 2t} (Z), with ∗ representing the
element-wise Hadamard product, and IS (X) = [11S (xijk )] r the element-wise indicator function.
√

Subproblem (4) is nonnegative tensor factorization, which can be solved by alternatively optimizing
one of A, B, or C when the other two are fixed, with q indexing iterations,
τ
k(Up+1 + O + λp )(1) − A(C p,q B p,q )T k2F
(7)
2
p+1
p
p,q
p,q
p,q T
p,q
p,q T
p,q †
= max{(U
+ O + λ )(1) (C
B )((C ) (C ) ∗ (B ) (B )) , 0}
(8)
τ
p,q+1
p+1
p
p,q
p,q+1 T 2
B
= arg min ι{z: z≥0} (B) + k(U
+ O + λ )(2) − B(C
A
) kF
(9)
B
2
p+1
p
p,q
p,q+1
p,q T
p,q
p,q+1 T
p,q+1 †
= max{(U
+ O + λ )(2) (C
A
)((C ) (C ) ∗ (A
) (A
)) , 0}
(10)
τ
C p,q+1 = arg min ι{z: z≥0} (C) + k(Up+1 + O + λp )(3) − C(B p,q+1 Ap,q+1 )T k2F
(11)
C
2
p+1
p
p,q+1
p,q+1
p,q+1 T
p,q+1
p,q+1 T
p,q+1 †
= max{(U
+ O + λ )(3) (B
A
)((B
) (B
) ∗ (A
) (A
)) , 0}. (12)
Ap,q+1 = arg min ι{z: z≥0} (A) +
A

Each subproblem is a constrained least squares problem the recovers the mode-d unfolding matrix
(Up+1 + O + λp ), starting from Ap,0 = Ap , B p,0 = B p , C p,0 = C p , and updating A, B, C according
to (7)-(12) until convergence, then Ap+1 = Ap,end , B p+1 = B p,end , C p+1 = C p,end . The updates
(7)-(12) for the subproblems (4) usually converge in less than ten iterations when warm started from
the previous iteration. Relative “residuals” are used to monitor the convergence of (3)-(5), and are
defined by
res1 =

[[Ap,q , B p,q , C p,q ]] − [[Ap,q−1 , B p,q−1 , C p,q−1 ]]
[[Ap,q−1 , B p,q−1 , C p,q−1 ]]

2

F

F

(13)

4

10

0

0

0.04

0.06

Noise ratio

0.08

0

10

50

0

2

0.02

0.04

0.06

Noise ratio

0.08

4

4

6

Reconstruction CP rank

8

10

(a)

0.06

0

0.02

0.04

0.06

Noise ratio

0.08

0.04
0.02
0

0.1

4

2

10

0.3

103

0.25

10

FP wrt X (l2)
FP wrt X (l1)
FP wrt X (l0)
FP wrt O (l2)
FP wrt O (l1)
FP wrt O (l0)

10
1
0

2

4

6

Reconstruction CP rank

8

2

1
0

(b)

2

0

0.02

0.04

0.06

Noise ratio

0.08

0.1

8

10

MSE wrt X (l2)
MSE wrt X (l1)
MSE wrt X (l0)
MSE wrt O (l2)
MSE wrt O (l1)
MSE wrt O (l0)

0.2

0.15

FN wrt X (l2)
FN wrt X (l1)
FN wrt X (l0)
FN wrt O (l2)
FN wrt O (l1)
FN wrt O (l0)

10

0

10

0.1

0.08

FN wrt X (l2)
FN wrt X (l1)
FN wrt X (l0)
FN wrt O (l2)
FN wrt O (l1)
FN wrt O (l0)

1
0

MSE wrt X (l2)
MSE wrt X (l1)
MSE wrt X (l0)
MSE wrt O (l2)
MSE wrt O (l1)
MSE wrt O (l0)

0.12

10

0.1

103

100

0

FP wrt X (l2)
FP wrt X (l1)
FP wrt X (l0)
FP wrt O (l2)
FP wrt O (l1)
FP wrt O (l0)

1

10

NTF (l2 norm)
NTF (l1 norm)
NOTF (l0 norm)

150

2

10

0.1

False positive counts

Convergence iterations

200

0.02

0.14

Mean square error

50

False negative counts

10

100

3

102

False negative counts

150

10

3

10

Mean square error

NTF (l2 norm)
NTF (l1 norm)
NOTF (l0 norm)

False positive counts

Convergence iterations

200

0.1

0.05
4

6

Reconstruction CP rank

(c)

8

10

0

0

2

4

6

Reconstruction CP rank

(d)

Figure 1: (a) Convergence iteration, (b) false positive count, (c) false negative count, and (d) mean square
error when varying the noise ratio (top) and varying the tensor CP rank R for NOTF (bottom) for the synthetic
dataset. We set rank R = 3 when varying the noise ratio (top) and set a noise ratio of 10% when varying the CP
rank (bottom). Note that reconstruction errors with respect to both groundtruth X and noisy observation O are
presented in (b)-(d).

(
res2 = max

[[Ap , B p , C p ]] − [[Ap−1 , B p−1 , C p−1 ]]
[[Ap−1 , B p−1 , C p−1 ]]

F

F

kλp − λp−1 kF
,
kλp−1 kF

)
.

(14)

The algorithm converges when res1 <  and res2 < , with typical  = 10−3 . The relative residual
res2 is inspired by the primal and dual residuals in [4]. Note that we are not seeking a unique
decomposition of A, B, C, but rather a stable low rank construction [[A, B, C]] of observation O that
has minimum sparse error.

3

Experiments on Synthetic Data

We test the NOTF algorithm on a synthetic dataset constructed as follows: (1) We create random
sparse matrices A ∈ R50×3 , B ∈ R20×3 , C ∈ R10×3 , where the sparse ratios (ratio of zero to nonzero values) are 70.67%, 55% and 30%, respectively. Each nonzero value is uniformly sampled from
range (0, 1). (2) Create the ground truth low rank matrix X = I{z:|z|>0} ([[A, B, C]]) ∈ R50×20×10 ;
it is the indicator tensor of nonzero values in the CP reconstruction from A, B, C; [[A, B, C]] has CP
rank 3. Note that X is a discrete (binary) valued tensor with sparsity ratio 75.75%. (3) Create the
observation tensor O by adding noise that flips a small portion of the binary values in X.
NOTF reconstructs a low rank matrix [[Â, B̂, Ĉ]] with at most CP rank R, where A ∈ R50×R , B ∈
R20×R , C ∈ R10×R . We vary the ratio of noise, i.e., the percentage of flips in O, up to 10%,
and the rank parameter R up to 10. We report on convergence iterations, false positive and false
negative counts when reconstructing the binary values in tensor X and O, and the mean square
error for reconstructing X and O. We compare NOTF with `0 norm in (1) with non-negative tensor
factorization (NTF) baselines using the `1 norm and `2 norms. Both NOTF and baseline methods
are initialized with a CP decomposition of the observation tensor O, λ0 = 0, with penalty parameter
τ = 10, and implementated in Matlab using the Tensor toolbox [3].
Fig. 1 presents the results of varying the noise ratio up to 10% (top) with rank parameter R = 3, and
varying the rank parameter R up to 10 (bottom) with noise of 10%, for the synthetic dataset. We
observe that the proposed NOTF solution with `0 norm performs well on the discrete measures (false
positive and false negative counts in (b) and (c)); we consider non-zeros as positives and zeros as
negatives in a tensor. NOTF achieves zero false positives and a relatively low false negative count
over all noise ratios. The `2 baseline achieves zero false negatives, but the false positive counts are
quite large. The `1 baseline achieves larger errors than NOTF on both positives and negatives, and is
slower to converge.
Fig. 1 (d) shows that NOTF with `0 norm does not outperform the baselines for mean square error; this
is not surprising. The discrete measurements are more important when reconstructing an occurrence
3

tensor. The error between the recovered low rank tensor [[Â, B̂, Ĉ]] and the observation O grows with
the noise ratio, while the error between [[Â, B̂, Ĉ]] and the groundtruth X is relatively stable. This
suggests that the recovered tensor [[Â, B̂, Ĉ]] can be used to de-noise the observation.
In Fig. 1 (bottom), we set a 10% noise ratio and vary the rank parameter R for the recovered tensor
[[Â, B̂, Ĉ]]. Note that R is an upper bound of the CP rank of [[Â, B̂, Ĉ]]. NOTF could achieve
zero for both false positives and false negatives when R = 6, which means that the ground truth
can be completely recovered from the noisy observation. However, NOTF becomes unstable with
larger R and leads to large false negative counts. A possible reason is that the initialization by CP
decomposition of O becomes less stable when large R is used, and the least squares in (7)-(12) are
often ill-posed and hard to solve. We finally observe that modeling the sparse error by the `0 norm
brings an additional benefit in that the recovered Â, B̂, and Ĉ are sparse; this leads to a clearer
interpretation of each rank-one tensor as a community.

4

Experiments on the resMBS dataset

We explore the roles played by financial institutions (FIs) across multiple contracts (FCs) using NOTF
with the `0 norm. ResMBS[5, 20, 19] contains extracted relationship of FI (e.g., Bank of America)
playing a role (e.g., issuer) for a specific financial contract. The discrete values of occurrence tensor
O ∈ R971×85×27 indicate the counts of extractions of the specific (FC, FI, Role) occurrence from
documents issued in 2005. O is sparse (1.02% non-zero values) and extraction noise is estimated to
be ≤ 0.2%. We describe some observations here and present the relevant figures in Section 6 due to
space limitations.
We vary the CP rank parameter R and reconstruct tensors for both the discrete observation O and
its binary version. Performance is similar for both while it is notably slower to reconstruct discrete
values (Fig. 2).We note that resMBS is challenging as the tensor is sparse. The false positives are
relatively stable while the false negatives decrease as R increases. With R = 20, the total error
count between the reconstructed tensor and the noisy observation is 3002; this roughly matches the
expected errors of the information extractor. The histogram (Fig. 3 (left)) shows that errors for each
FC is in a reasonable range (0, 20) with a mean of 3.
At last, we examine the discovered communities by NOTF for resMBS. Each rank-one tensor
ar ◦ br ◦ cr , r = 1, . . . , R represents a community. Fig. 3 (right) presents the nonzero ratio and Fig. 4
presents the distribution of the reconstructed tensor component Â = [ar ], B̂ = [br ], Ĉ = [cr ]. An
interesting observation is that the communities are “centered” around FIs, i.e., each community only
contains one or two FIs. Some FIs could play various roles and appear in various FCs, while some
FIs only play a limited number of roles in a limited number of FCs.

5

Discussion and future work

We present non-negative occurrence tensor factorization (NOTF) for analyzing heterogeneous networks. CP tensor decomposition is adapted to discover the embedded communities. The `0 norm is
used to model the discrete tensor values and sparse errors, and the objective is solved with an efficient
splitting optimization algorithm. NOTF is applied to both synthetic data and a new heterogeneous
bipartite graph, resMBS, representing financial role relationships extracted from financial contracts.
Preliminary results are promising and suggest that NOTF can be used to de-noise the occurrence
tensors and identify communities in resMBS.
There are several directions for future work. The `0 norm is known to be difficult to optimize. The
`p norm (0 < p < 1) satisfies the KL inequality, is often used as a surrogate, and may provide a
theoretical convergence guarantee. The penalty parameter τ is crucial for both convergence speed and
solution quality for nonconvex problems; adaptive ADMM[22, 21] which automates the selection
of τ, achieves promising practical performance. To deal with the high sparsity of the resMBS
tensor, domain-specific constraints (e.g., each FC should contain an FI play role “Issuer” ) may
boost performance. Finally, it is interesting to apply NOTF for analyzing some other heterogeneous
networks that could be represented with an occurrence tensor.
4

Acknowledgments
ZX and TG were supported by US NSF grant CCF-1535902 and by US ONR grant N00014-15-12676. ZX and LR were supported by NSF grants CNS1305368 and DBI1147144, and NIST award
70NANB15H194.

References
[1] A. Anandkumar, R. Ge, D. Hsu, S. M. Kakade, and M. Telgarsky. Tensor decompositions for learning
latent variable models. Journal of Machine Learning Research, 15(1):2773–2832, 2014.
[2] A. Anandkumar, P. Jain, Y. Shi, and U. Niranjan. Tensor vs matrix methods: Robust tensor decomposition
under block sparse perturbations. arXiv preprint, 2015.
[3] B. W. Bader and T. G. Kolda. Algorithm 862: MATLAB tensor classes for fast algorithm prototyping. ACM
Transactions on Mathematical Software, 32(4):635–653, December 2006. doi: 10.1145/1186785.1186794.
[4] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning
via the alternating direction method of multipliers. Found. and Trends in Mach. Learning, 3:1–122, 2011.
[5] D. Burdick, S. De, L. Raschid, M. Shao, Z. Xu, and E. Zotkina. resMBS: Constructing a financial supply
chain graph from financial prospecti. In SIGMOD DSMM workshop. ACM, 2016.
[6] J. D. Carroll and J.-J. Chang. Analysis of individual differences in multidimensional scaling via an n-way
generalization of “eckart-young” decomposition. Psychometrika, 35(3):283–319, 1970.
[7] X. Chen, Z. Han, Y. Wang, Q. Zhao, D. Meng, and Y. Tang. Robust tensor factorization with unknown noise.
In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 5213–5221,
2016.
[8] D. Goldfarb and Z. Qin. Robust low-rank tensor recovery: Models and algorithms. SIAM Journal on
Matrix Analysis and Applications, 35(1):225–253, 2014.
[9] T. Goldstein and S. Osher. The split Bregman method for L1-regularized problems. SIAM Journal on
Imaging Sciences, 2(2):323–343, 2009.
[10] T. Goldstein, B. O’Donoghue, S. Setzer, and R. Baraniuk. Fast alternating direction optimization methods.
SIAM Journal on Imaging Sciences, 7(3):1588–1623, 2014.
[11] Q. Gu, H. Gui, and J. Han. Robust tensor decomposition with gross corruption. In Advances in Neural
Information Processing Systems, pages 1422–1430, 2014.
[12] R. A. Harshman. Foundations of the parafac procedure: Models and conditions for an" explanatory"
multi-modal factor analysis. 1970.
[13] H. Huang and C. Ding. Robust tensor factorization using r 1 norm. In Computer Vision and Pattern
Recognition, 2008. CVPR 2008. IEEE Conference on, pages 1–8. IEEE, 2008.
[14] T. G. Kolda and B. W. Bader. Tensor decompositions and applications. SIAM review, 51(3):455–500, 2009.
[15] A. P. Liavas and N. D. Sidiropoulos. Parallel algorithms for constrained tensor factorization via alternating
direction method of multipliers. IEEE Transactions on Signal Processing, 63(20):5450–5463, 2015.
[16] K. Maruhashi, F. Guo, and C. Faloutsos. Multiaspectforensics: Pattern mining on large-scale heterogeneous
networks with tensor analysis. In Advances in Social Networks Analysis and Mining (ASONAM), 2011
International Conference on, pages 203–210. IEEE, 2011.
[17] E. E. Papalexakis, C. Faloutsos, and N. D. Sidiropoulos. Tensors for data mining and data fusion: Models,
applications, and scalable algorithms. ACM Transactions on Intelligent Systems and Technology (TIST), 8
(2):16, 2016.
[18] A. Schein, J. Paisley, D. M. Blei, and H. Wallach. Bayesian poisson tensor factorization for inferring multilateral relations from sparse dyadic event counts. In Proceedings of the 21th ACM SIGKDD International
Conference on Knowledge Discovery and Data Mining, pages 1045–1054. ACM, 2015.
[19] Z. Xu and L. Raschid. Probabilistic financial community models with latent dirichlet allocation for financial
supply chains. In SIGMOD DSMM workshop. ACM, 2016.

5

[20] Z. Xu, D. Burdick, and L. Raschid. Exploiting lists of names for named entity identification of financial
institutions from unstructured documents. arXiv preprint arXiv:1602.04427, 2016.
[21] Z. Xu, S. De, M. A. T. Figueiredo, C. Studer, and T. Goldstein. An empirical study of admm for nonconvex
problems. In NIPS workshop on nonconvex optimization, 2016.
[22] Z. Xu, M. A. Figueiredo, and T. Goldstein. Adaptive ADMM with spectral penalty parameter selection.
arXiv preprint arXiv:1605.07246, 2016.

6

Appendix: experimental results for resMBS
200

4

False positive counts

50

0

1.2

10

5

10

Reconstruction CP rank

15

0

20

1

0.8
0.6

1

0

Binary
Discrete

1.4

102

100

#104

1.6

103

150

Convergence iterations

1.8

10

Binary
Discrete

False negative counts

6

0.4

Binary
Discrete

0

5

10

Reconstruction CP rank

(a)

15

0.2

20

0

5

10

Reconstruction CP rank

(b)

15

20

(c)

Figure 2: (a) Convergence iteration, (b) false positive counts, and (c) false negative counts when vary the tensor
CP rank R for the reconstruction of resMBS dataset. Both binary and discrete tensor of resMBS are tested. Note
that reconstruction errors presented in (b)(c) are based on noisy observation O as the groundtruth is unknown.
0.5

300
Sparseness ratio (nonzeros)

250

#FC

200
150
100
50
0

0

5

10

15

0.4

0.3

0.2

0.1

0

20

A: FCs
B: FIs
C: roles

0

5

10

Error per FC

15

20

Community index

Figure 3: (left) Histogram of errors when constructed resMBS tensor with R = 20. (right) Nonzero ratio
for the constructed tensor component Â, B̂, Ĉ. Each rank-one tensor ar ◦ br ◦ cr , r = 1, . . . , R represents a
community.

200
300

0.15

FC index

400

10
20

0.1

600
700

0.05

800
900
5

10

Community index

(a)

15

20

0

5

18

10

14

15

10

20

6

16

1

30

0.8

40

500

20

1.2

50

0.6

60

0.4

70

0.2

80
5

10

Community index

15

20

0

Role index

0.2

FI index

100

12
8
4

25

2

5

Community index

(b)

Figure 4: The distribution of FC, FI, and roles in each community.

7

10

(c)

15

20

0

