arXiv is now an independent nonprofit! Learn more
License: arXiv.org perpetual non-exclusive license
arXiv:1809.03776v2 [cs.LG] 26 Sep 2018

Solving Non-identifiable Latent Feature Models

Ryota Suzuki Thanks: r-suzuki@fa.jp.nec.com    Shingo Takahashi Thanks: s-takahashi@fh.jp.nec.com    Murtuza Petladwala Thanks: murtuza@cq.jp.nec.com    Shigeru Kohmoto Thanks: s-koumoto@bq.jp.nec.com Affiliation: NEC Data Science Research Laboratories Affiliation: Tokyo, Japan
Abstract

Latent feature models (LFM)s are widely employed for extracting latent structures of data. While offering high, parameter estimation is difficult with LFMs because of the combinational nature of latent features, and non-identifiability is a particularly difficult problem when parameter estimation is not unique and there exists equivalent solutions. In this paper, a necessary and sufficient condition for non-identifiability is shown. The condition is significantly related to dependency of features, and this implies that non-identifiability may often occur in real-world applications. A novel method for parameter estimation that solves the non-identifiability problem is also proposed. This method can be combined as a post-process with existing methods and can find an appropriate solution by hopping efficiently through equivalent solutions. We have evaluated the effectiveness of the method on both synthetic and real-world datasets.

1 Introduction

Latent variable models are widely used for obtaining hidden data-structures. A mixture model is a leading example of such a latent variable model, in which each instance of data is classified into a latent class. Latent feature model (LFM) [1] is an extension of a mixture model, in which data is characterized into not only one but into a combination of latent features. LFM is used for various applications to extract such hidden structures of data in medical [2, 3], facial images [4], and social, gene, and document networks [5].

While LFM offers wide application, parameter estimation (unsupervised learning) is difficult with it because the optimization algorithms encounters its combinational non-convex nature; possible combinations of KK features arising in NN data increase exponentially in 2NK2^{NK}, and multi-modal structure of cost-functions (such as log-likelihood and model evidence) results in numerous local optima, preventing optimization algorithms to obtain global optima. Several approaches have been proposed to deal with non-convexity. Reed and Zoubin [6] focused on the submodularity, a discrete analog of convexity, of a cost function in a non-negative assumption of features, in which efficient greedy algorithm is available. Another approach has been proposed by Yen et al. [7], in which a convex relaxation is employed with a Lasso regularizer, and it may be solved as a certain class of semi-definite programming. Even though these methods avoid the problem of non-convexity, they still have to face the other difficulty, which is non-identifiability.

Identifiability in latent variable models represents the uniqueness of parameter estimation [8]. The parameters are correctly estimated if the solution of an optimization problem is unique. However, if there are multiple parameters that result in the same cost, the solution of the optimization problem may not be unique. In such a non-identifiable situation, error in parameter estimation may significantly worsen [9]. A special case of non-identifiability in LFM is shown in [10], in which two features have the same value, and their method solves the problem by balancing the size of the features. Additionally, [7] have shown a sufficient condition for identifiability and also that the condition holds with high probability under an assumption that features appear in an independently and identically distributed (i.i.d.) Bernoulli process. However, i.i.d. Bernoulli, assumed in most existing works [10, 11] including an Indian buffet process (IBP) [1, 4, 6, 12, 13], is too stringent an assumption in many real-world applications; there might be hidden constraints for which some features do not appear at the same time, and/or features may have a hidden hierarchical structure in which one subsumes another. We show that, in such a case of features’ having a dependency, parameter estimation is non-identifiable, and optimization methods may face the difficulty. Another case of non-identifiability we show is that of the existence of bias, which can be represented as a feature commonly appearing in all the data [2, 3].

Since difficulty in non-identifiability is due to the existence of solutions having the same cost, which we refer to here as equivalent solutions, optimization methods may find a solution that is not the true parameter but equivalent one. We have developed a hopping algorithm that efficiently finds equivalent solutions from one to another in succession, which maximizes prior probabilities without degrading likelihood.

In this paper, we present a solution to the problem of non-identifiability, and our contributions are as follows: first, we have derived a necessary and sufficient condition for non-identifiability in LFM. Secondly, we have derived sufficient conditions for non-identifiability that is, significantly, related to the dependency of features. Thirdly, we have developed a novel method for parameter estimation that can be combined as a post-process with existing methods and that can find appropriate solutions by hopping through equivalent solutions. Finally, we have also shown the effectiveness of the new method on both synthetic and real-world datasets.

2 Latent Feature Models

In a latent feature model (LFM), observed data is assumed to be represented as a combination of KK latent features. Let XN×DX\in\mathbb{R}^{N\times D} be an observed data matrix, and its rows represent instances of DD-dimensional observations. We assume that XX can be generated as:

X=ZW+ε,ZPZ(Z),WPW(W)X=ZW+\varepsilon,\;\;\;Z\sim P_{Z}(Z),\;\;\;W\sim P_{W}(W) (1)

where kk-th row of WK×DW\in\mathbb{R}^{K\times D} is the kk-th latent feature 𝒘k\bm{w}_{k}, and the kk-th column 𝒛k\bm{z}_{k} of unknown binary matrix Z{0,1}N×KZ\in\{0,1\}^{N\times K} represents the incidence of the kk-th latent feature along NN data, and ε\varepsilon is an unknown noise matrix. A typical class of LFMs is linear-Gaussian LFMs [10, 11]:

PZ=Bernoulli(𝝅),PW=𝒩(0,σW2),ε𝒩(0,σX2).P_{Z}=\mathrm{Bernoulli}(\bm{\pi}),\;P_{W}=\mathcal{N}(0,\sigma_{W}^{2}),\;\varepsilon\sim\mathcal{N}(0,\sigma_{X}^{2}). (2)

where 𝝅=(π1,,πK)\bm{\pi}=(\pi_{1},\cdots,\pi_{K}) is probability for feature incidence. If πk=1/2\pi_{k}=1/2, a maximum a posteriori (MAP) inference is obtained from following optimization problem:

argminZ{0,1}N×K,WK×D{XZWF2+τWF2},\mathop{\mathrm{arg~min}}_{Z\in\{0,1\}^{N\times K},W\in\mathbb{R}^{K\times D}}\left\{\|X-ZW\|_{F}^{2}+\tau\|W\|_{F}^{2}\right\}, (3)

where τ=σX2/σW2\tau=\sigma_{X}^{2}/\sigma_{W}^{2}, and it becomes maximum likelihood (ML) estimation in τ0\tau\to 0. IBP is also used for PZP_{Z} in Bayesian nonparametric settings [1, 4, 6, 12, 13], and the optimization problem (3) can also be obtained from MAP asymptotics of IBP linear-Gaussian LFM with 2\ell_{2} regularization [4, 7].

From a general point of view, LFM (1) can be seen as a class of matrix factorizations of XX into ZZ and WW with a constraint such that ZZ is binary [14]. Our analyses and methods deal not only with typical LFMs but also with such a general class of matrix factorization. Factorial hidden Markov Models (FHMMs) [15, 16] and non-negative LFMs [6, 17] are an example of such a class of matrix factorization, where additional constraints and prior probabilities are assumed. We focus on a general characteristics of matrix factorization (1), and our analyses and methods are applicable to a wide-range of matrix factorization problems.

Most existing works on LFMs, explicitly or implicitly, make two strong assumptions related to identifiability. The first common but strong assumption is the statistical independence of features [4, 6, 10, 11, 12, 13]. In these case, one assumes that the incidence of an individual feature will be independently generated in a Bernoulli process. As we will see in Sec. 3.5, the independence of features is sometimes too strong an assumption in actual situations, and an absence of independence will cause the problem of non-identifiability.

The second assumption related to identifiability is that a model has zero bias 𝔼[ε|Z]=0\mathbb{E}[\varepsilon|Z]=0 [4, 6, 10, 11, 12, 13, 14]. There might, however, be a background feature common to all instances in actual situations. Introducing bias term 𝒘bias\bm{w}_{\mathrm{bias}} to (1) is equivalent to an additional feature that is always active [2], and the existence of such a bias results in non-identifiability, as we will see in Sec. 3.5.

3 Analyses of Non-identifiability

3.1 Non-identifiability in LFMs

Figure 1: Whole feasible solutions of Z{0,1}4×3Z\in\{0,1\}^{4\times 3} for (a) identifiable and (b) non-identifiable cases. Blue dots correspond to individual solutions that are connected to one another by gray lines of Hamming distance one. Ground truth ZZ^{*} is plotted as a red circle and its components are shown in insets. (c) shows solutions equivalent to the ground truth of (b). We use noiseless data with features WW^{*}, as three of four synthetic images originally used in [1].

In this paper, as is also seen in [14, 7], we consider identifiability in LFMs in terms of uniqueness of the solution. Before presenting our theoretical analyses, let us briefly overview (non-)identifiability in LFMs with an example.

In LFMs, non-identifiability difficulties are mainly due to discrete nature of ZZ. Since the optimization problem (3) for WW with fixed ZZ is convex, the MAP solution W^\hat{W} corresponding to Z^\hat{Z} (which is also an expectation 𝔼[W|Z^,X]\mathbb{E}[W|\hat{Z},X]) is uniquely calculated in closed form [1, 7]. The uniqueness of ZZ is more complicated. Figure 1 (a, b) shows every possible combinations of ZZ in two different XX’s. The axes of the figures correspond to the terms of (3). In the identifiable case (a), only the ground truth has zero-residual XZWF=0\|X-ZW\|_{F}=0. However, in the non-identifiable case (b), there are many solutions having zero-residual, and inference methods minimizing the residual, such as [7, 11], may fall into incorrect zero-residual solution rather than ground truth. Moreover, zero-residual solutions (Fig. 1c) are distant from one another at Hamming distances, and this results in a multi-modal nature in a cost function preventing iterative algorithms, such as Markov-chain Monte Carlo (MCMC)[1, 12] and variational Bayesian (VB)[18] methods, from converging to a global optimum.

3.2 General Matrix Factorizations

Let us next consider conditions for (non-)identifiability and properties of equivalent solutions. Since LFMs are a class of matrix factorizations, we start from identifiability in general matrix factorization [19, 7, 14] to get an overall picture.

Definition 1 (Identifiability)

Let 𝒵N×K,𝒲K×D\mathcal{Z}\subset\mathbb{R}^{N\times K},\mathcal{W}\subset\mathbb{R}^{K\times D} be sets of matrices. We say that a pair of matrices (Z,W)𝒵×𝒲(Z,W)\in\mathcal{Z}\times\mathcal{W} is identifiable if for all (Z,W)𝒵×𝒲(Z^{\prime},W^{\prime})\in\mathcal{Z}\times\mathcal{W}, ZW=ZWZW=Z^{\prime}W^{\prime} implies {Z:,k}k=1K={Z:,k}k=1K\{Z_{:,k}\}_{k=1}^{K}=\{Z^{\prime}_{:,k}\}_{k=1}^{K} and {Wk,:}k=1K={Wk,:}k=1K\{W_{k,:}\}_{k=1}^{K}=\{W^{\prime}_{k,:}\}_{k=1}^{K}.

We can consider wide-ranging classes of matrix factorization by choosing 𝒵\mathcal{Z} and 𝒲\mathcal{W}. The LFM is the case of 𝒵={0,1}N×K\mathcal{Z}=\{0,1\}^{N\times K} and 𝒲=K×D\mathcal{W}=\mathbb{R}^{K\times D}. We can further say that a matrix factorization of XN×DX\in\mathbb{R}^{N\times D} is identifiable if one of the minimizers of a residual XZWF2\|X-ZW\|_{F}^{2} (i.e., one of the ML solutions) is identifiable. In noiseless settings (ε=0\varepsilon\!=\!0), the identifiability of XX is consistent with that of the ground truth (Z,W)(Z^{*},W^{*}).

It is known that most equivalent solutions of matrix factorization ZW=ZWZ^{\prime}W^{\prime}=ZW have a specific form. For a context of non-negative matrix factorization (NMF), where 𝒵=+N×K\mathcal{Z}=\mathbb{R}_{+}^{N\times K} and 𝒲=+K×D\mathcal{W}=\mathbb{R}_{+}^{K\times D}, Laurberg, et al. [19] showed that all equivalent NMF solutions ZW=ZWZ^{\prime}W^{\prime}=ZW have a form Z=ZU,W=U1WZ^{\prime}=ZU,W^{\prime}=U^{-1}W if rank(ZW)=K\mathrm{rank}(ZW)=K. This assertion is also true in the general matrix factorization in Definition 1 since their proof only used properties of linear space.

Theorem 1

Laurberg, et al. [19] Let Z𝒵N×KZ\in\mathcal{Z}\subset\mathbb{R}^{N\times K} and W𝒲K×DW\in\mathcal{W}\subset\mathbb{R}^{K\times D}. Assume rank(ZW)=K\mathrm{rank}(ZW)=K. Then for any (Z,W)𝒵×𝒲(Z^{\prime},W^{\prime})\in\mathcal{Z}\times\mathcal{W}, ZW=ZWZW=Z^{\prime}W^{\prime} holds if and only if there exists a regular matrix UK×KU\in\mathbb{R}^{K\times K} such that Z=ZU,W=U1WZ^{\prime}=ZU,W^{\prime}=U^{-1}W.

3.3 Equivalence Classes

Let us consider a set of equivalent solutions of an LFM to quantify (non-)identifiability. Let [(Z,W)][(Z,W)] be an equivalence class of an equivalence relation "\sim" defined as (Z,W)(Z,W)ZW=ZW(Z,W)\sim(Z^{\prime},W^{\prime})\,\Leftrightarrow\,ZW=Z^{\prime}W^{\prime}. Namely, [(Z,W)]={(Z,W)𝒵×𝒲|ZW=ZW}[(Z,W)]=\left\{(Z^{\prime},W^{\prime})\in\mathcal{Z}\times\mathcal{W}\middle|ZW=Z^{\prime}W^{\prime}\right\}.

Assuming rank(ZW)=K\mathrm{rank}(ZW)=K, Theorem 1 guarantees that all elements of [(Z,W)][(Z,W)] will be represented as (ZU,U1W)({Z}U,U^{-1}{W}). Then we can consider LFM (𝒵={0,1}N×K,𝒲=K×D\mathcal{Z}\!=\!\{0,1\}^{N\times K},\mathcal{W}\!=\!\mathbb{R}^{K\times D}), the equivalence class to be expressible as [(Z,W)]={(ZU,U1W)|UH(Z)}[(Z,W)]=\left\{({Z}U,U^{-1}{W})\middle|U\in H({Z})\right\}, where

H(Z)={UK×K|detU0,ZU{0,1}N×K}.H({Z})=\left\{U\in\mathbb{R}^{K\times K}\middle|\det U\neq 0,{Z}U\in\{0,1\}^{N\times K}\right\}.

Non-identifiability in the LFM is characterized by distinct elements UH(Z)U\in H(Z) other than permutation matrices U𝒮KU\not\in\mathcal{S}_{K}, where 𝒮K\mathcal{S}_{K} is symmetric group of degree KK. To reduce the degree of freedom in permutation, we assume a quotient set11 1 In a precise sense, H(Z)/𝒮KH(Z)/\mathcal{S}_{K} stand s for H(Z)/𝒮KH(Z)/\mathop{\sim_{\mathcal{S}_{K}}}, the quotient of H(Z)H(Z) by an equivalence relation “𝒮K\sim_{\mathcal{S}_{K}}” such that U𝒮KUσ𝒮K,U=UσU\mathop{\sim}_{\mathcal{S}_{K}}U^{\prime}\Leftrightarrow{}^{\exists}\sigma\in\mathcal{S}_{K},U^{\prime}=U\sigma. H(Z)/𝒮KH(Z)/\mathcal{S}_{K}, which is a set of sets, in which each set consists of transform matrices UH(Z)U\in H(Z) having the same column entries but in different orders.

3.4 Equivalent Condition for Identifiability in LFMs

By using the notation above, non-identifiability can be quantified by a cardinality |H(Z)/𝒮K||H(Z)/\mathcal{S}_{K}|, which represents the distinct number of equivalent solutions avoiding duplication in permutations. Additionally, we get |H(Z)/𝒮K|=1|H(Z)/\mathcal{S}_{K}|=1 as a necessary and sufficient condition for the identifiability under the assumption of rank(ZW)=K\mathrm{rank}(ZW)\!=\!K. We consider H(Z)/𝒮KH(Z)/\mathcal{S}_{K} to be trivial if H(Z)/𝒮KH(Z)/\mathcal{S}_{K} has only one element . While the assumption rank(ZW)=K\mathrm{rank}(ZW)\!=\!K denotes both rank(Z)=K\mathrm{rank}(Z)\!=\!K and rank(W)=K\mathrm{rank}(W)\!=\!K, we further prove that the former is unnecessary, and we get a following theorem providing a necessary and sufficient condition for the identifiability:

Theorem 2

Let Z{0,1}N×K,WK×DZ\in\{0,1\}^{N\times K},W\in\mathbb{R}^{K\times D}. Assume rank(W)=K\mathrm{rank}(W)=K. Then (Z,W)(Z,W) is identifiable if and only if H(Z)/𝒮KH({Z})/\mathcal{S}_{K} is trivial.

Note that Theorem 2 is a stronger result than the “identifiability condition” mentioned in [7], which supplies only a sufficient condition for identifiability in the case of rank(ZW)=K\mathrm{rank}(ZW)\!=\!K. While the assumption rank(W)=K\mathrm{rank}(W)\!=\!K in Theorem 2 is not necessarily required for the identifiability22 2 We can show an example in which identifiability holds despite rank(W)<K\mathrm{rank}(W)\!<\!K. See Supplemental Materials. , we assume rank(W)=K\mathrm{rank}(W)=K in the following discussion in consideration for DKD\gg K in many applications.

Finally, in Figure 2, we summarize the relationship between the conditions we have derived: (Full-Z) rank(Z)=K\mathrm{rank}(Z)=K, (ID) the identifiability, and (Trivial) |H(Z)/𝒮K|=1|H(Z)/\mathcal{S}_{K}|=1.

Figure 2: Relationship between conditions under the assumption rank(W)=K\mathrm{rank}(W)=K (Full-W). Full-Z: rank(Z)=K\mathrm{rank}(Z)=K, ID: Identifiability in LFMs (Definition 1), Trivial: |H(Z)/𝒮K|=1|H(Z)/\mathcal{S}_{K}|=1 (Theorem 2), i.i.d.: Zi.i.d.Bernoulli(p)Z\sim\mathrm{i.i.d.}\,Bernoulli(p) with 0<p<1,N0<p<1,N\to\infty (Theorem 3), PDC: pairwise dependency conditions (Theorem 4).

3.5 Sufficient Conditions for (Non-)Identifiability

Non-identifiability in LFMs is significantly related to dependency between features. We illustrate here two sufficient conditions, both for identifiability and non-identifiability (i.i.d. and PDC in Figure 2), related to independence and dependence of features, respectively.

The first condition, the sufficient condition for identifiability, is the statistical independence of the features. [7] have shown that an LFM is identifiable with high probability if entries of ZZ are i.i.d. Bernoulli p=0.5p\!=\!0.5. More generally, identifiability holds with probability one in NN\!\to\!\infty for any <p<10\!<\!p\!<\!1. We show this in the following theorem:

Theorem 3

H(Z)/𝒮KH({Z})/\mathcal{S}_{K} is trivial if ZZ is a binary matrix s.t. {Zn,:}n=1N={0,1}K\{Z_{n,:}\}_{n=1}^{N}=\{0,1\}^{K}.

If ZZ is i.i.d. Bernoulli with 0<p<10<p<1, probabilities for every combination of each rowP(Zn,:)P(Z_{n,:}) will be non-zero. Therefore, in the limit of NN\to\infty, all the combinations {0,1}K\{0,1\}^{K} may appear in rows of ZZ with probability one, and this results in identifiability via Theorem 3. In other words, from a contraposition of Theorem 3, it can be said that possible non-identifiability is due to a lack of observed combinations. In real-world applications observing all 2K2^{K} combinations is rarely possible, and, as we will see later, some combinations may never appear in rows of ZZ even if NN\to\infty because of hidden dependency between features.

The second condition is in regard to non-identifiability. We propose three pairwise dependency conditions (PDCs) sufficient for non-identifiability,

Theorem 4

Let Z0,1N×K{Z}\in{0,1}^{N\times K} be a binary matrix. |H(Z)/𝒮K|3|H({Z})/\mathcal{S}_{K}|\geq 3 holds if there exists a distinct pair i,j(ij)i,j\;(i\neq j) of features satisfying one of following conditions for all n=1,,Nn=1,\cdots,N:

  • PDC1: (zn,i=1)(zn,j=0)(z_{n,i}=1)\Rightarrow(z_{n,j}=0).

  • PDC2: (zn,i=1)(zn,j=1)(z_{n,i}=1)\Rightarrow(z_{n,j}=1).

  • PDC3: (zn,i=0)(zn,j=0)(z_{n,i}=0)\Rightarrow(z_{n,j}=0).

The PDCs in Theorem 4 often appear in real-world applications33 3 We can show the commonality of PDCs from our survey on actual datasets. See Supplemental Materials. , and, unfortunately, they are unknown in most cases. Inference methods may then suffer from non-identifiability whenever there exists at least a pair of features PDC holds. For instance, a typical case of PDC1 would be disjoint features. Assuming, for example, that these features correspond to characteristics of cats, then a pair of features (i,j)=("male","black")(i,j)=("\textrm{male}","\textrm{black}") may appear at the same time, but (i,j)=("male","female")(i,j)=("\textrm{male}","\textrm{female}") will not appear concurrently. In the case of PDC2 and 3, typical cases would be latent hierarchical structures. For example, in considering features (i,j)=("cats","mammals")(i,j)=("\textrm{cats}","\textrm{mammals}"), the feature "mammals" is always active whenever "cats" is active since cats are mammals.

Another example of non-identifiability is the existence of a bias term. An LFM with bias term is equivalent to an unbiased LFM with an extra feature that is always active [2, 3]. The existence of such a bias feature is followed by PDC2 because the feature is always active regardless of other features. Let (Z,W)(Z^{\prime},W^{\prime}) be an equivalent solution corresponding to the transform matrix U=I+𝒆i𝒆biasT2𝒆i𝒆iTU=I+\bm{e}_{i}\bm{e}_{\mathrm{bias}}^{\mathrm{T}}-2\bm{e}_{i}\bm{e}_{i}^{\mathrm{T}}, where the index "bias" refer to the bias feature, absence and the presence of ii-th feature 𝒛i\bm{z}_{i}^{\prime} is inverted from 𝒛i\bm{z}_{i} and a sign of the ii-th feature is flipped 𝒘i=𝒘i\bm{w}_{i}^{\prime}=-\bm{w}_{i} in swap of raising the level of the bias 𝒘bias=𝒘bias+𝒘i\bm{w}_{\mathrm{bias}}^{\prime}=\bm{w}_{\mathrm{bias}}+\bm{w}_{i}. Such an “inverted” solution has been obtained by some algorithms, including Gibbs sampler [20], K-features [4], and possibly other LFM algorithms. In the case of a biased LFM, we can prove that there is a lower bound to the number of equivalent solutions |H(Z)/𝒮K|(K+1)×2K1|H(Z)/\mathcal{S}_{K}|\geq(K+1)\times 2^{K-1} (see Supplemental Materials).

4 Hopping through Equivalent Solutions

We now consider finding a superior solution among equivalent solutions. Once an arbitrary estimator finds a (not necessarily optimal) MAP solution (Z^,W^)(\hat{Z},\hat{W}) for the model (1), there might be some equivalent solutions (Z,W)[(Z^,W^)](Z^{\prime},W^{\prime})\in[(\hat{Z},\hat{W})]. Although they have the same residual, some of those might be close to the ground truth (Z,W)(Z^{*},W^{*}) but others might be far from it. Our method obtains a superior one among them having a maximal prior probability. In it, for efficiency, we assume transform matrix UU to be an integer matrix, i.e., we sample UU from a subset of H(Z^)H(\hat{Z}):

H~(Z^)=H(Z^)K×K.\tilde{H}(\hat{Z})=H(\hat{Z})\cap\mathbb{Z}^{K\times K}. (4)

Although there might be UH(Z)U\in H(Z) that is not an integer matrix in some cases, we can show that H~(Z)=H(Z)\tilde{H}(Z)=H(Z) holds in many cases in consideration of NKN\gg K (see Supplemental Materials).

Let us next introduce a quadratic form f(Z)=12n,kzn,k(zn,k1)f(Z)=\frac{1}{2}\sum_{n,k}z_{n,k}(z_{n,k}-1), then (4) is denoted as H~(Z)={UK×K|detU0,f(ZU)=0}\tilde{H}(Z)=\left\{U\in\mathbb{Z}^{K\times K}\middle|\det U\neq 0,f(ZU)=0\right\}. We can assume f(ZU)f(ZU) to be a cost function that measures how different ZUZU is from an binary matrix since it is non-negative integer and becomes zero if and only if ZUZU is a binary. Further, f(ZU)f(ZU) is an upper bound to the number of ZUZU components other than 0 or 1, and is equal to it if 1ZU2-1\preceq ZU\preceq 2.

4.1 Sampling Equivalent Solutions

Algorithm 1 Sampling Columns of  UU
Input: binary matrix ZZ, sample size NsN_{s}
 Initialize 𝒰=\mathcal{U}=\emptyset.
Φ,Σ,ΨT=svd(Z)\Phi,\Sigma,\Psi^{\mathrm{T}}=svd(Z);   Λ=ΨΣ1\Lambda=\Psi\Sigma^{-1};   𝝁=12ΛTZT𝟏\bm{\mu}=\frac{1}{2}\Lambda^{\mathrm{T}}Z^{\mathrm{T}}\bm{1}.
for i=1i=1 to NsN_{s} do
  Sample f(i)f^{(i)} from arbitrary distribution.
  Sample 𝒔(i)\bm{s}^{(i)} uniformly from a (K1)(K\!-\!1)-sphere of center 𝝁\bm{\mu} and radius 𝝁2+2{f(i)}2\sqrt{\|\bm{\mu}\|^{2}+2\{f^{(i)}\}^{2}}.
  𝒖(i)=round(Λ𝒔(i))\bm{u}^{(i)}=\mathrm{round}(\Lambda\bm{s}^{(i)})
  𝒰=𝒰{𝒖(i)}\mathcal{U}=\mathcal{U}\cup\{\bm{u}^{(i)}\}   if   𝒖(i)𝟎\bm{u}^{(i)}\neq\bm{0}.
end for
Algorithm 2 Equivalence Hopper
Input: initial matrices ZZ, WW,   𝒰\mathcal{U} sampled by Algorithm 1
 Initialize U0=IU_{0}=I
for m=1,,Mm=1,\cdots,M do
  J={U:=Um1+(𝒖(i)𝒖k)𝒆kT|𝒖(i)𝒰;k=1,,K;rank(U)=K}J=\left\{U^{\prime}:=U_{m-1}+(\bm{u}^{(i)}-\bm{u}_{k})\bm{e}_{k}^{\mathrm{T}}\;\middle|\;\bm{u}^{(i)}\in\mathcal{U};\;k=1,\!\cdots,\!K;\;\mathrm{rank}(U^{\prime})=K\right\}
  Sample UmJU_{m}\in J according to the probability (7).
end for

Let 𝒖k\bm{u}_{k} be the kk-th column of UU, and Λ\Lambda be a K×KK\times K matrix such that ΛΛT=(ZTZ)1\Lambda\Lambda^{\mathrm{T}}=(Z^{\mathrm{T}}Z)^{-1}, which is calculated by e.g., singular value decomposition (SVD) of ZZ, 𝒔k=Λ1𝒖k\bm{s}_{k}=\Lambda^{-1}\bm{u}_{k}, and 𝝁=12ΛTZT𝟏\bm{\mu}=\frac{1}{2}\Lambda^{\mathrm{T}}Z^{\mathrm{T}}\bm{1}, where 𝟏=(1,,1)T\bm{1}=(1,\cdots,1)^{\mathrm{T}}. then f(ZU)f(ZU) for kk-th column can be evaluated as:

f(Z𝒖k)=12(𝒖kTZTZ𝒖k𝒖kTZT𝟏)=12𝒔k𝝁212𝝁2.f(Z\bm{u}_{k})=\frac{1}{2}\left(\bm{u}_{k}^{\mathrm{T}}Z^{\mathrm{T}}Z\bm{u}_{k}-\bm{u}_{k}^{\mathrm{T}}Z^{\mathrm{T}}\bm{1}\right)=\frac{1}{2}\|\bm{s}_{k}-\bm{\mu}\|^{2}-\frac{1}{2}\|\bm{\mu}\|^{2}. (5)

If UH~(Z)U\in\tilde{H}(Z), (5) is evaluated to be zero, then 𝒔k\bm{s}_{k} will be on the (K1K\!-\!1)-sphere of center 𝝁\bm{\mu} and radius 𝝁\|\bm{\mu}\|. A possible 𝒖k\bm{u}_{k} can then be obtained by sampling 𝒔\bm{s} uniformly from the (K1K\!-\!1)-sphere, taking the nearest integer 𝒖=round(Λ𝒔)\bm{u}\!=\!\mathrm{round}(\Lambda\bm{s}), and accept 𝒖k=𝒖\bm{u}_{k}=\bm{u}   if   𝒖𝟎f(Z𝒖)=0\bm{u}\neq\bm{0}\wedge f(Z\bm{u})\!=\!0.

The strict method mentioned above (which samples UU strictly from H~(Z)\tilde{H}({Z})) may, unfortunately, fail in some cases because Z{Z}, estimated with an arbitrary algorithm, may have randomness and its flipped component may spoil equivalent solutions by breaking PDCs in Theorem 4. To handle such a randomness, we employ a tolerance to the equivalent conditions. If the transformed matrix Z=ZUZ^{\prime}=ZU includes an integer other than 0 or 1, it will be rejected from equivalent solutions. However, if there exists another binary matrix Z~\tilde{Z}^{\prime} which is close to ZZ^{\prime}, then Z~\tilde{Z}^{\prime} will be a nearly equivalent solution to ZZ. We measure this closeness by f(ZU)f(ZU) since it is an approximation for the number of non-binary components as mentioned above, and it can be calculated by (5) without scanning all NN rows of ZUZU. In the Algorithm 1 shows this tolerant method. We first sample ff^{*} from some distribution, and then sample 𝒖k\bm{u}_{k} such that f(Z𝒖k)=ff(Z\bm{u}_{k})=f^{*} by sampling 𝒔\bm{s} uniformly from the (K1K\!-\!1)-sphere of radius 𝝁2+2f2\sqrt{\|\bm{\mu}\|^{2}+{2f^{*}}^{2}} instead of 𝝁\|\bm{\mu}\|. By choosing the distribution of ff^{*}, we can tune the tolerance for the protrusion of ZU{Z}U from binary matrices. We employ discrete exponential distribution with parameter λ\lambda for this role. The strict case is a limit of λ\lambda\to\infty.

Figure 3: Number of equivalent solutions |H~(Z)/𝒮K||\tilde{H}(Z)/\mathcal{S}_{K}| found by Algorithm 1 for K=6K=6. ZZ is sampled from (a) i.i.d Bernoulli(0.5), (b) with a bias feature, and (c,d) with 1, 3 pair(s) of features with PDC, respectively. Boxes, whiskers and red lines represent, respectively, quartiles, extrema, and medians of 50 trials.

Once 𝒰\mathcal{U} is sampled by Algorithm 1, equivalent solutions are obtained by selecting KK columns of UU from 𝒰\mathcal{U} so that UU is a regular matrix. Figure 3 shows the number of equivalent solutions found from H~(Z)/𝒮K\tilde{H}(Z)/\mathcal{S}_{K} by using the strict version (λ\lambda\to\infty) of Algorithm 1. In the i.i.d. case (Fig. 3a), the number of equivalent solutions diminishes rapidly to one with an increasing NN. This is consistent with results mentioned in [7]. Figure 3 (b) to (d) show cases of ZZ sampled with a bias feature and PDCs (Theorem 4). In these three cases, solutions remains multiple even in N=200N=200, keeping the problem non-identifiable. In the case of a bias feature (Fig. 3b), the number of solutions converges to the theoretical lower bound (K+1)×2K2=112(K+1)\times 2^{K-2}=112, and when a pair of features holds a PDC (Fig. 3c), it converges to 3, the lower bound shown in Theorem 4. When three pairs of features have PDCs (Fig. 3d), the number of solutions fluctuate up to 10310^{3} while a combination of 3 independent PDCs will result in 33=273^{3}=27 solutions. This implies that multiple PDCs sharing a same feature give rise to another PDC, e.g., PDC2 for (i,j)=(1,2)(i,j)=(1,2) and (2,3)(2,3) implies PDC2 for (1,3)(1,3).

4.2 Optimizing over Equivalent Solutions

Now we consider obtaining appropriate solution among equivalent solutions. We select UU~(Z^)U\in\tilde{U}(\hat{Z}) so that obtained solution (Z,W)=(Z^U,U1W^)(Z^{\prime},W^{\prime})=(\hat{Z}U,U^{-1}\hat{W}) is more appropriate, i.e., having higher prior probability without degrading likelihood. So we employ the following cost function:

g(U,Z^,W^)=logPZ(Z^U)+logPW(U1W^)+γf(Z^U).g(U;\hat{Z},\hat{W})=\log P_{Z}(\hat{Z}U)+\log P_{W}(U^{-1}\hat{W})+\gamma f(\hat{Z}U). (6)

The first two terms of (6) correspond to log-priors in (1). The second term with parameter γ\gamma plays a role keeping (Z,W)({Z}^{\prime},{W}^{\prime}) nearly equivalent to (Z^,W^)(\hat{Z},\hat{W}) and maintaining little change in the likelihood.

We can solve the optimization problem by MCMC-based sampling method (Algorithm 2), sampling UU by updating its columns successively according to a Boltzmann distribution with parameter β\beta:

P(U|Z^,W^)exp{βg(U,Z^,W^)},P(U|\hat{Z},\hat{W})\propto\exp\{-\beta g(U;\hat{Z},\hat{W})\}, (7)

In the limit of β\beta\to\infty, Algorithm 2 becomes a greedy local-search algorithm, which has a lower-bound in its optimality for typical linear-Gaussian LFMs (see Supplemental Materials). However, sampling method (β<\beta<\infty) works well in many cases since the feasible space is small enough.

5 Experiments

In this section, we demonstrate utility of our method by applying it to both synthetic and actual data as a post-process combined with existing algorithms. Once an estimation (Z^,W^)(\hat{Z},\hat{W}) is obtained by such an algorithm, we execute the Equivalence Hopper (Algorithm 2) with (Z^,W^)(\hat{Z},\hat{W}) as an input. We employed a cost function (6) with linear-Gaussian priors (2) for consistency with baseline methods. We evaluate our method with state-of-the-art algorithms both for a parametric approach, LatentLasso [7], and a Bayesian nonparametric approach with non-negative constraints, MEIBP [6].

We examine our method with both synthetic and actual datasets. For the synthetic data, we use simulated images also used in [7], where each feature is ×3030\!\times\!30 image and its randomly selected ×77\!\times\!7 region is set as44 4 For MEIBP, we use absolute values instead so that non-negative constraint of MEIBP makes sense. 𝒩(0,1)\mathcal{N}(0,1). We employ PDC constraints for the synthetic ZZ to examine effectiveness of our method for non-identifiability. For actual-data experiments, we use the UK-DALE [21] a dataset for Non-intrusive Load Monitoring [22], and the Piano transcription dataset [23]. From the UK-DALE dataset, we extracted raw current waveforms for every 2 minutes from house-1 / 2015 / week-1 data. And we used Bach_850 from Piano dataset by taking a power spectrum. Since UK-DALE data contains negative values and MEIBP is not applicable, we only applied LatentLasso to it.

Since our method does not change likelihood except for a small change due to tolerance in Algorithm 1. We evaluate our method by means of following metrics:

  • Hamming Error: EHamm=minσ𝒮K1NKZσZ0E_{\mathrm{Hamm}}=\min_{\sigma\in\mathcal{S}_{K}}\frac{1}{NK}\|Z\sigma-Z^{*}\|_{0}

  • Regularizer: EReg=1KDWFE_{\mathrm{Reg}}=\frac{1}{KD}\|W\|_{F}

where the first one with ground truth ZZ^{*} is only available for synthetic data. The second corresponds to the logarithm of the prior PWP_{W}, which we minimize in Algorithm 2 via the cost function (6) and it is expected to offer sparser representation of the data, more representable and closer to true parameters.

Refer to caption

Figure 4: (a,b) Metrics for synthetic data (K=14,N=1000K=14,N=1000) with PDC constraints estimated by (a) LatentLasso and (b) MEIBP with/without Equivalence Hopper (Algorithm 2). Markers and error bars represent means and extrema of 5 trials. (c) Example WW for synthetic data (K=6K=6) of ground truth (upper), estimated by LatentLasso (middle), and with Equivalence Hopper (bottom).
Table 1: Metrics for actual-data experiments. The number of features K+K^{+} is constant for LatentLasso and estimated variables for MEIBP. Residuals XZWF\|X-ZW\|_{F} have not changed in 4 significant digits before and after Equivalence Hopper is applied.
Baseline w/ Equiv. Hopp.
Dataset Algorithm K+K^{+} XZWF\|X-ZW\|_{F} ERegE_{\mathrm{Reg}} ERegE_{\mathrm{Reg}}
UK-DALE LatentLasso 5 0.422 0.692 0.677 (-2.3%)
UK-DALE LatentLasso 10 0.247 1.574 1.480 (-6.0%)
Piano LatentLasso 20 0.830 0.699 0.698 (-0.2%)
Piano LatentLasso 40 0.698 0.492 0.491 (-0.3%)
Piano MEIBP 49 0.700 0.635 0.625 (-1.5%)

Figure 4 (a,b) shows the metrics for synthetic data with varying number of PDCs (NpdcN_{\mathrm{pdc}}). Accuracies in Hamming error of both two baseline methods (LatentLasso and MEIBP) degrade by increasing NpdcN_{\mathrm{pdc}}. By applying our method, Equivalence Hopper, accuracies are significantly improved, and especially in MEIBP, the error almost halves in average and the ground truth have been achieved within 5 trials at Npdc=3,12N_{\mathrm{pdc}}=3,12. In LatentLasso, while the best-case EHammE_{\mathrm{Hamm}} is mostly unchanged before and after applying Equivalence Hopper, a range of EHammE_{\mathrm{Hamm}} values significantly narrows keeping the average error lower and the estimation robust. A robustness of the estimation is quite important for unsupervised learning especially in non-identifiable case, in which we cannot evaluate the error without knowledge about ground truth. Further, Equivalence Hopper is worth applying in most cases rather than repeating the preceding methods until getting better results since our algorithm never worsen the result, and executes very fast, delivering results in less than 5 seconds (less than a single iteration of the preceding methods) for M=1000M=1000 iterations with Ns=1000N_{s}=1000. It is also remarkable that only few percentage improvement of ERegE_{\mathrm{Reg}} result in a drastic improvement in EHammE_{\mathrm{Hamm}}, which is consistent with the result in Figure 1(b,c) where completely different but equivalent solution has small difference in WF\|W\|_{F}.

Example features obtained by LatentLasso with/without Equivalence Hopper are shown in Figure 4 (c). In a raw estimation of LatentLasso (middle), there are some features, e.g. the second feature, having multiple features of the ground truth (top) in duplicate. With application of the Equivalence Hopper algorithm (bottom), such a duplication was suppressed, and we obtained a sparser solution. The fourth feature in Fig. 4(c) is an instance of “inverted” feature we discussed in Sec. 3, in which the sign of feature is flipped in the result of LatentLasso and is restored by Equivalence Hopper.

Finally, we show the actual-data experiments in Table 1. Our method obtained a better solution in a manner of a ERegE_{\mathrm{Reg}} than the preceding algorithm without any degradation in residual (i.e., likelihood) in all cases. While the change in ERegE_{\mathrm{Reg}} is in few percentage, it cannot be neglected because small change of ERegE_{\mathrm{Reg}} in synthetic data results in drastic improvement of Hamming error EHammE_{\mathrm{Hamm}}. We believe that small changes in ERegE_{\mathrm{Reg}} is important for comparison between equivalent solutions.

References

  • [1] T. Griffiths and Z. Ghahramani, “Infinite Latent Feature Models and the Indian Buffet Process,” in Gatsby Computational Neuroscience Unit Technical Report. Gatsby Computational Neuroscience Unit, May 2005.
  • [2] F. J. R. Ruiz, I. Valera, C. Blanco, and F. Perez-Cruz, “Bayesian Nonparametric Comorbidity Analysis of Psychiatric Disorders,” Journal of Machine Learning Research, vol. 15, pp. 1215–1247, 2014.
  • [3] I. Valera, M. F. Pradier, and Z. Ghahramani, “General Latent Feature Modeling for Data Exploration Tasks,” in 2017 ICML Workshop on Human Interpretability in Machine Learning (WHI 2017), Sydney, NSW, Australia, Jul. 2017, arXiv: 1707.08352. [Online]. Available: http://arxiv.org/abs/1707.08352
  • [4] T. Broderick, B. Kulis, and M. Jordan, “MAD-Bayes: MAP-based Asymptotic Derivations from Bayes,” in PMLR, Feb. 2013, pp. 226–234.
  • [5] J. Zhu and B. Chen, “Latent feature models for large-scale link prediction,” Big Data Analytics, vol. 2, p. 3, Feb. 2017.
  • [6] C. Reed and G. Zoubin, “Scaling the Indian Buffet Process via Submodular Maximization,” in PMLR, Feb. 2013, pp. 1013–1021.
  • [7] I. E.-H. Yen, W.-C. Lee, S.-E. Chang, A. S. Suggala, S.-D. Lin, and P. Ravikumar, “Latent Feature Lasso,” in PMLR, Jul. 2017, pp. 3949–3957.
  • [8] K. P. Murphy, Machine Learning: A Probabilistic Perspective. Cambridge, MA: The MIT Press, Aug. 2012.
  • [9] S. Watanabe, “Algebraic Analysis for Nonidentifiable Learning Machines,” Neural Computation, vol. 13, no. 4, pp. 899–933, Apr. 2001.
  • [10] K. Hayashi and R. Fujimaki, “Factorized Asymptotic Bayesian Inference for Latent Feature Models,” in Advances in Neural Information Processing Systems 26, C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, Eds., 2013, pp. 1214–1222.
  • [11] H.-Y. Tung and A. J. Smola, “Spectral Methods for Indian Buffet Process Inference,” in Advances in Neural Information Processing Systems 27, Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, Eds., 2014, pp. 1484–1492.
  • [12] F. Doshi-Velez and Z. Ghahramani, “Accelerated Sampling for the Indian Buffet Process,” in Proceedings of the 26th Annual International Conference on Machine Learning, ser. ICML ’09. New York, NY, USA: ACM, 2009, pp. 273–280.
  • [13] T. L. Griffiths and Z. Ghahramani, “The Indian Buffet Process: An Introduction and Review,” J. Mach. Learn. Res., vol. 12, pp. 1185–1224, Jul. 2011.
  • [14] M. Slawski, M. Hein, and P. Lutsik, “Matrix factorization with binary components,” in Advances in Neural Information Processing Systems 26, C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, and K. Q. Weinberger, Eds., 2013, pp. 3210–3218.
  • [15] Z. Ghahramani and M. I. Jordan, “Factorial Hidden Markov Models,” Machine Learning, vol. 29, no. 2-3, pp. 245–273, 1997.
  • [16] J. V. Gael, Y. W. Teh, and Z. J. C. Ghahramani, “The infinite factorial hidden markov model,” in NIPS, 2008.
  • [17] D. O’Malley, V. V. Vesselinov, B. S. Alexandrov, and L. B. Alexandrov, “Nonnegative/binary matrix factorization with a D-Wave quantum annealer,” arXiv:1704.01605 [quant-ph, stat], Apr. 2017.
  • [18] F. Doshi, K. T. Miller, J. Van Gael, and Y. Whye Teh, “Variational inference for the indian buffet process.” vol. 5, pp. 137–144.
  • [19] H. Laurberg, M. G. Christensen, M. D. Plumbley, L. K. Hansen, and S. H. Jensen, “Theorems on Positive Data: On the Uniqueness of NMF,” 2008.
  • [20] Z. Ghahramani and T. L. Griffiths, “Infinite latent feature models and the Indian buffet process,” in Advances in Neural Information Processing Systems 18, Y. Weiss, B. Schölkopf, and J. C. Platt, Eds. MIT Press, 2006, pp. 475–482.
  • [21] J. Kelly and W. Knottenbelt, “The UK-DALE dataset, domestic appliance-level electricity demand and whole-house demand from five UK homes,” vol. 2, no. 150007, 2015.
  • [22] G. W. Hart, “Nonintrusive appliance load monitoring,” Proceedings of the IEEE, vol. 80, no. 12, pp. 1870–1891, Dec. 1992.
  • [23] G. E. Poliner and D. P. W. Ellis, “A Discriminative Model for Polyphonic Piano Transcription,” EURASIP Journal on Advances in Signal Processing, vol. 2007, no. 1, p. 048317, Dec. 2006.
  • [24] J. S. Liu, F. Liang, and W. H. Wong, “The Multiple-Try Method and Local Optimization in Metropolis Sampling,” Journal of the American Statistical Association, vol. 95, no. 449, pp. 121–134, Mar. 2000.
  • [25] G. Storvik, “On the Flexibility of Metropolis Hastings Acceptance Probabilities in Auxiliary Variable Proposal Generation,” Scandinavian Journal of Statistics, vol. 38, no. 2, pp. 342–358, 2011.
  • [26] M. L. Fisher, G. L. Nemhauser, and L. A. Wolsey, “An analysis of approximations for maximizing submodular set functions II,” in Polyhedral Combinatorics, ser. Mathematical Programming Studies. Springer, Berlin, Heidelberg, 1978, pp. 73–87.

Supplemental Materials

Appendix A Remarks on condition rank(W)=K\mathrm{rank}(W)=K

While a condition rank(Z)=K\mathrm{rank}(Z)\!=\!K is necessary for identifiability as shown in the proof of Theorem 2, a condition rank(W)=K\mathrm{rank}(W)\!=\!K is not necessarily required for identifiability. We can see this by assuming the following counterexample with K=2,D=1,N=4K=2,D=1,N=4:

Z=(00011011),W=(110),X=ZW=(011011),Z=\begin{pmatrix}0&0\\ 0&1\\ 1&0\\ 1&1\end{pmatrix},\;W=\begin{pmatrix}1\\ 10\end{pmatrix},\;X=ZW=\begin{pmatrix}0\\ 1\\ 10\\ 11\end{pmatrix},

then (Z,W)(Z,W) is identifiable despite rank(W)=1<K\mathrm{rank}(W)=1<K.

Appendix B Lower bound for the number of equivalent solutions in biased LFM

We show the number of equivalent solutions |H(Z)/𝒮K||H(Z)/\mathcal{S}_{K}| has a lower bound (K+1) 2K1(K+1)\,2^{K-1} if ZZ has a column that is always active.

Assume k=Kk=K be a bias feature (e.g., 𝒛K=𝟏\bm{z}_{K}=\bm{1}), and ff be as defined in Sec. 4. Then, f(Z𝒖)=0f(Z\bm{u})=0 holds for every 𝒖𝒰={𝒆1,,𝒆K,𝒆K𝒆1,,𝒆K𝒆K1\bm{u}\in\mathcal{U}=\{\bm{e}_{1},\cdots,\bm{e}_{K},\;\bm{e}_{K}\!-\!\bm{e}_{1},\cdots,\bm{e}_{K}\!-\!\bm{e}_{K-1}}, because 𝒛i,𝒛K𝒛i{0,1}N\bm{z}_{i},\,\bm{z}_{K}-\bm{z}_{i}\in\{0,1\}^{N}. Hence, we can construct UU by selecting its columns from 𝒰\mathcal{U} such that UU is regular matrix.

From the regularity of UU, it includes at least one non-zero element in each row. Therefore, let

𝒰i\displaystyle\mathcal{U}_{i} ={𝒆i,𝒆K𝒆i}fori=1,,K1,\displaystyle=\{\bm{e}_{i},\bm{e}_{K}-\bm{e}_{i}\}\;\;\;\mathrm{for}\;i=1,\cdots,K\!-\!1,
𝒰K\displaystyle\mathcal{U}_{K} ={𝒆K,𝒆K𝒆1,,𝒆K𝒆K1},\displaystyle=\{\bm{e}_{K},\bm{e}_{K}-\bm{e}_{1},\cdots,\bm{e}_{K}-\bm{e}_{K-1}\},

then, UU is regular iff all 𝒰1,,𝒰K\mathcal{U}_{1},\cdots,\mathcal{U}_{K} include at least one column of UU. We count up the number of UU by considering two cases:

  • The case of UU including 𝒆K\bm{e}_{K} in its columns,
    the other K1K-1 columns of UU are selected one from every 𝒰1,,𝒰K1\mathcal{U}_{1},\cdots,\mathcal{U}_{K-1}. Then, the number of combination is 2K12^{K-1}.

  • The case of UU not including 𝒆K\bm{e}_{K} in its columns,
    all KK columns of UU are selected at least one from every 𝒰1,,𝒰K1\mathcal{U}_{1},\cdots,\mathcal{U}_{K-1}. Then, two columns of UU are selected from a single 𝒰k\mathcal{U}_{k}, and the other K2K-2 columns are selected one from every 𝒰k(kk)\mathcal{U}_{k}^{\prime}\;(k^{\prime}\neq k). Then, the number of combination is (K1) 2K2(K-1)\,2^{K-2}.

Summing up the cases, we get a lower bound of equivalent solutions:

|H(Z)/𝒮K|2K1+(K1) 2K2=(K+1) 2K2.|H(Z)/\mathcal{S}_{K}|\geq 2^{K-1}+(K-1)\,2^{K-2}=(K+1)\,2^{K-2}.\square

Appendix C Sufficient condition for H(Z)=H~(Z)H(Z)=\tilde{H}(Z)

We show below some sufficient conditions for H~(Z)=H(Z)\tilde{H}(Z)=H(Z).

Theorem 5

Let Z{0,1}N×KZ\in\{0,1\}^{N\times K} be a binary matrix of rank KK. Assume UH(Z)U\in H(Z). UU is an interger matrix if one of following conditions holds:

  • a.

    min{n|zn,k=1}min{n|zn,k=1}\min\{n|z_{n,k}=1\}\neq\min\{n|z_{n,k^{\prime}}=1\} for kkk\neq k^{\prime},

  • b.

    ZZ has a K×KK\times K submatrix ζ\zeta s.t. |det(ζ)|=1|\det({\zeta})|=1.

The first condition in Theorem 5 assumes that the timing of the first appearance of each feature is different from that of each of the others. In other words, at most one new feature appears at the same time. The second condition is a more permissive condition since the first one follows it. The condition may hold when NKN\gg K since it holds if at least one of (NK)\binom{N}{K} submatrices in ZZ has an absolute determinant of one.

Appendix D Detail of Algorithm 2 for linear Gaussian LFMs

In linear-Gaussian LFMs, the cost function (6) is represented as:

g(U,Z^,W^)=τU1W^F2+γf(Z^U),g(U;\hat{Z},\hat{W})=\tau\|U^{-1}\hat{W}\|_{F}^{2}+\gamma f(\hat{Z}U), (8)

We use a MCMC-based method that samples UU according to (7) by updating each column in each iteration. Using Algorithm 1, we sample candidates 𝒰={𝒖(1),,𝒖(Ns)}\mathcal{U}=\{\bm{u}^{(1)},\cdots,\bm{u}^{(N_{s})}\} of the kk-th column of the next UU^{\prime}, and we select one according to distribution (7) such that UU^{\prime} is regular. Our method is similar to the multi-try method introduced in [24, 25], but in our method, the proposal of candidates 𝒰\mathcal{U} does not depend on the current state of UU. This means that we can use the same 𝒰\mathcal{U} both for sampling proposals and for calculating acceptance ratios, and this results in the acceptance ratio always being one if the current 𝒖k\bm{u}_{k} is among the candidates. While we could further reuse 𝒰\mathcal{U} over iterations, it is better to resample in several iterations for global convergence.

The most time-consuming step in Algorithm 2 is singularity determination of U=U+(𝒖(i)𝒖k)𝒆kT=U+Δ𝒖(i)𝒆kTU^{\prime}=U+(\bm{u}^{(i)}-\bm{u}_{k})\bm{e}_{k}^{\mathrm{T}}=U+\Delta\bm{u}^{(i)}\bm{e}_{k}^{\mathrm{T}} and probability calculation P(U|Z^,W^)P(U^{\prime}|\hat{Z},\hat{W}) in (7) for every 𝒖(i)𝒰\bm{u}^{(i)}\in\mathcal{U}. For singularity determination, we adapt the rank-1 update formula of determinants,

det(U)=det(U)(1+𝒆kTU1Δ𝒖(i)).\det(U^{\prime})=\det({U})\cdot\left(1+\bm{e}_{k}^{\mathrm{T}}U^{-1}\Delta\bm{u}^{(i)}\right). (9)

We can check the singularity of UU^{\prime} by (1+𝒆kTU1Δ𝒖(i))(1+\bm{e}_{k}^{\mathrm{T}}U^{-1}\Delta\bm{u}^{(i)}) to be zero. And for the probability calculation, we get

U1W^F2=U1W^F2+𝒆kTΩ(𝒗(i)2𝒆k2𝒗(i)),\|U^{\prime-1}\hat{W}\|_{F}^{2}=\|U^{-1}\hat{W}\|_{F}^{2}+\bm{e}_{k}^{\mathrm{T}}\,\Omega\,\left(\|\bm{v}^{(i)}\|^{2}\bm{e}_{k}-2\bm{v}^{(i)}\right),

where 𝒗(i)=U1Δ𝒖(i)1+𝒆kTU1Δ𝒖(i)\bm{v}^{(i)}=\frac{U^{-1}\Delta\bm{u}^{(i)}}{1+\bm{e}_{k}^{\mathrm{T}}U^{-1}\Delta\bm{u}^{(i)}}, and Ω=U1W^W^TUT\Omega=U^{-1}\hat{W}\hat{W}^{\mathrm{T}}U^{-\mathrm{T}}. Then we finally get

P\displaystyle P (U|Z^,W^)exp{β(τ𝒆kTΩ(𝒗(i)2𝒆k2𝒗(i))+γf(Z^𝒖(i)))}.\displaystyle(U^{\prime}|\hat{Z},\hat{W})\propto\exp\left\{-\beta\left(\tau\bm{e}_{k}^{\mathrm{T}}\Omega\left(\|\bm{v}^{(i)}\|^{2}\bm{e}_{k}-2\bm{v}^{(i)}\right)+\gamma f(\hat{Z}\bm{u}^{(i)})\right)\right\}. (10)

By reusing values of f(Z𝒖(i))f(Z\bm{u}^{(i)}) calculated in Algorithm 1, the calculation time of (9) and (10) is 𝒪(K)\mathcal{O}(K) for each candidate if U1U^{-1} and Ω\Omega are given. We can maintain these matrices incrementally in 𝒪(K2)\mathcal{O}(K^{2}) per iteration as U1=VU1U^{\prime-1}=VU^{-1} and Ω1=VΩVT\Omega^{\prime-1}=V\Omega V^{\mathrm{T}}, where V=I𝒗(i)𝒆kTV=I-\bm{v}^{(i)}\bm{e}_{k}^{\mathrm{T}}.

In total, calculation time with our method is 𝒪(NK2+DK2)\mathcal{O}(NK^{2}+DK^{2}) for initializing Ω\Omega and Λ\Lambda, 𝒪(NsK2)\mathcal{O}(N_{s}K^{2}) for resampling candidates (Algorithm 1), and 𝒪(NsK+K2)\mathcal{O}(N_{s}K+K^{2}) for each MCMC iteration (Algorithm 2). Our method is quite fast since there is no need to scan all NN data once initialized. The convergence of Algorithm 2 is proven by following theorem,

Theorem 6

Let q(U|U)q(U^{\prime}|U) be MCMC kernel defined as Algorithm 2. Then,

  • q(U|U)q(U^{\prime}|U) satisfies detailed balance condition.

  • q(U|U)q(U^{\prime}|U) is transitable over regular integer matrices.

After mm iterations of the MCMC step, the distribution of UU will converge to (7) in mm\to\infty.

Appendix E Lower-bound of optimality

Algorithm 2 has a lower bound on its optimality for typical linear-Gaussian LFMs in a limit of β\beta\to\infty because it becomes greedy local search algorithm for submodular maximization on a matroid constraint. Rewriting the cost function gg in (8) as:

g(U,Z^,W^)=τTr[(UUT)1W^W^T]+γk=1Kf(Z^𝒖k),g(U;\hat{Z},\hat{W})=\tau\mathrm{Tr}[(UU^{T})^{-1}\hat{W}\hat{W}^{T}]+\gamma\sum_{k=1}^{K}f(\hat{Z}\bm{u}_{k}),

its domain can be extended from regular matrices to rank-K ones with more than K columns. Assuming S𝒰S\subset\mathcal{U} be a set of 𝒖\bm{u} NOT included in columns of U, then possible SS’s form a matroid and the target function gg is supermodular on it. So, there exists linear set function hh so that ghg-h is monotone submodular, and the problem becomes a monotone submodular maximization on a matroid constraint. It is well known that local search yields 1/2-optimal solution for the problem [26].

Appendix F Survey on non-identifiability conditions

We can find out conditions for non-identifiability appears many real-world datasets. Table 2 shows the survey on datasets for multi-label classification datasets in LIBSVM library55 5 https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/multilabel.html . We can see that at least 6.1% (siam-competition2007) and in average 57.3% of 12K(K1)\frac{1}{2}K(K-1) pairs of features satisfy PDCs in Theorem 4. These results imply that LFMs will suffer from non-identifiability in most cases of actual applications.

Table 2: Statistics of datasets in LIBSVM.
Dataset NN KK # Pairs of PDC PDC Ratio
mediamill (exp1) 30,993 101 3,074 60.9 %
rcv1v2 (topics) 23,149 101 3,359 66.5 %
rcv1v2 (industries) 23,149 313 46,189 94.6 %
rcv1v2 (regions) 23,149 228 24,925 96.3 %
scene-classification 1,211 6 8 53.3 %
siam-competition2007 21,519 22 14 6.1 %
yeast 1,500 14 8 8.8 %

Appendix G Proofs for theorems

Theorem 1

Already shown in [19].

Theorem 2

In the case of rank(ZW)=K\mathrm{rank}(ZW)=K (i.e. rank(Z)=rank(W)=K\mathrm{rank}(Z)\!=\!\mathrm{rank}(W)\!=\!K) the assertion in Theorem 2 follows from Theorem 1 as discussed in the paper. Then we prove in the case of rank(Z)<rank(W)=K\mathrm{rank}(Z)\!<\!\mathrm{rank}(W)\!=\!K, where both the trivialness and the identifiability is false.

  • Non-trivialness of H(Z)/𝒮KH(Z)/\mathcal{S}_{K}

    In the case of ZZ having a column 𝒛l\bm{z}_{l} which is all zero, a matrix U0=I+a𝒆l𝒆lTU_{0}=I+a\bm{e}_{l}\bm{e}_{l}^{\mathrm{T}} for any a1,0a\neq-1,0 is in H(Z)H(Z) because detU0=a+10\det U_{0}=a+1\neq 0 and ZU0=Z{0,1}N×KZU_{0}=Z\in\{0,1\}^{N\times K}. And U0U_{0} is not a permutation matrix. Then |H(Z)/𝒮K|2|H(Z)/\mathcal{S}_{K}|\geq 2.

    In the case that no column in ZZ is zero vector, since rank(Z)<K\mathrm{rank}(Z)\!<K, there exists 𝒃K\{𝟎}\bm{b}\in\mathbb{R}^{K}\backslash\{\bm{0}\} s.t. Z𝒃=𝟎Z\bm{b}=\bm{0}. Let bmb_{m} be a non-zero component of 𝒃\bm{b}, then 𝒃𝒆m\bm{b}\neq\bm{e}_{m} because Z𝒃=𝟎Z\bm{b}=\bm{0} and Z𝒆m=𝒛m𝟎Z\bm{e}_{m}=\bm{z}_{m}\neq\bm{0}. A matrix U1=(𝒆1,,𝒆m1,𝒃,𝒆m+1,,𝒆K)U_{1}=(\bm{e}_{1},\cdots,\bm{e}_{m-1},\bm{b},\bm{e}_{m+1},\cdots,\bm{e}_{K}) is in H(Z)H(Z) because detU1=bm0\det U_{1}=b_{m}\neq 0 and ZU1=(𝒛1,,𝒛m1,𝟎,𝒛m+1,,𝒛K){0,1}N×KZU_{1}=(\bm{z}_{1},\cdots,\bm{z}_{m-1},\bm{0},\bm{z}_{m+1},\cdots,\bm{z}_{K})\in\{0,1\}^{N\times K}. And U1U_{1} is not a permutation matrix. Then |H(Z)/𝒮K|2|H(Z)/\mathcal{S}_{K}|\geq 2.

  • Non-identifiability of (Z,W)(Z,W)

    In the case of ZZ having a column 𝒛l\bm{z}_{l} which is all zero, (ZU0,U01W)(ZU_{0},U_{0}^{-1}W) is another solution since U01WWU_{0}^{-1}W\neq W (whereas ZU0=ZZU_{0}=Z).

    In the case that no column in ZZ is zero vector, (ZU1,U11W)(ZU_{1},U_{1}^{-1}W) is another solution because ZU1=(𝒛1,,𝒛m1,𝟎,𝒛m+1,,𝒛K)ZZU_{1}=(\bm{z}_{1},\cdots,\bm{z}_{m-1},\bm{0},\bm{z}_{m+1},\cdots,\bm{z}_{K})\neq Z. \square

Theorem 3

Assume {Zn,:}n=1N={0,1}K\{Z_{n,:}\}_{n=1}^{N}=\{0,1\}^{K} and UH(Z)U\in H(Z). Let 𝒖kT\bm{u}_{k}^{T} be a kk-th row or UU. Since 𝒆1T,,𝒆KT\bm{e}_{1}^{\mathrm{T}},\cdots,\bm{e}_{K}^{\mathrm{T}} and 𝟏T\bm{1}^{\mathrm{T}} is included in rows of ZZ, i.e.

n1,,nK,m,Znk,:=𝒆k,Zm,:=𝟏,{}^{\exists}n_{1},\cdots,n_{K},m,\;\;\;Z_{n_{k},:}=\bm{e}_{k},Z_{m,:}=\bm{1},

then

𝒖k=(𝒆kTU)T=(Znk,:U)T{0,1}K,\bm{u}_{k}=(\bm{e}_{k}^{\mathrm{T}}U)^{\mathrm{T}}=(Z_{n_{k},:}U)^{\mathrm{T}}\in\{0,1\}^{K},
k=1K𝒖k=(𝟏TU)T=(Zm,:U)T{0,1}K.\sum_{k=1}^{K}\bm{u}_{k}=(\bm{1}^{\mathrm{T}}U)^{\mathrm{T}}=(Z_{m,:}U)^{\mathrm{T}}\in\{0,1\}^{K}.

Therefore, UU is binary matrix and each row of UU sums to one. Also consider that UU is regular, UU is permutation matrix. Hence H(Z)=𝒮KH(Z)=\mathcal{S}_{K}. \square

Theorem 4

We provide examples of transfer matrices. Let 𝒆i\bm{e}_{i} be a vector of which the ii-th component is 1 and the others are 0, and

Rij=I+𝒆i𝒆jT,Qij=I+𝒆i𝒆jT2𝒆i𝒆iT.R_{ij}=I+\bm{e}_{i}\bm{e}_{j}^{\mathrm{T}},\;Q_{ij}=I+\bm{e}_{i}\bm{e}_{j}^{\mathrm{T}}-2\bm{e}_{i}\bm{e}_{i}^{\mathrm{T}}.

Then I,Rij,RjiH(Z)I,R_{ij},R_{ji}\in H(Z) in PDC1,  I,Rij1,QjiH(Z)I,R_{ij}^{-1},Q_{ji}\in H(Z) in PDC2,  and   I,Qij,Rji1H(Z)I,Q_{ij},R_{ji}^{-1}\in H(Z) in PDC3. \square

Note that PDC2 and PDC3 are essentially the same because one can be derived from the other by employing contraposition and exchanging ii and jj.

Theorem 5

We firstly prove the case of condition b., and then prove a. by using it.

  • b.

    ZZ has a K×KK\times K submatrix ζ\zeta s.t. |det(ζ)|=1|\det({\zeta})|=1.

    Assume UH(Z)U\in H(Z). Let ζ\zeta^{\prime} be a K×KK\times K submatrix of Z=ZUZ^{\prime}=ZU picking the same rows as ζ\zeta. Then ζ=ζU\zeta^{\prime}=\zeta U.

    As ζ\zeta is a binary matrix (hence, an integer matrix) and |det(ζ)|=1|\det({\zeta})|=1, ζ\zeta is a unimodular matrix, and it has an inverse ζ1\zeta^{-1} that is also a unimodular matrix. Hence U=ζ1ζU=\zeta^{-1}\zeta^{\prime} is integer matrix.

  • a.

    min{n|zn,k=1}min{n|zn,k=1}\min\{n|z_{n,k}=1\}\neq\min\{n|z_{n,k^{\prime}}=1\} for kkk\neq k^{\prime}.

    Let m(k)=min{n|zn,k=1}m(k)=\min\{n|z_{n,k}=1\}, and m(k)m(k) is well-defined because {n|zn,k=1}\{n|z_{n,k}=1\} is not empty since rankZ=K\mathrm{rank}\,Z=K. The condition assumed states that m(k)m(k) is injective.

    Consider a integer array (l1,,lK)(l_{1},\cdots,l_{K}) that is sorted from (1,,K)(1,\cdots,K) by ascending order of m(lk)m(l_{k}) values. Let K×KK\times K submatrix ζ\zeta of ZZ as ζk,:=Zm(lk),:\zeta_{k,:}=Z_{m(l_{k}),:} (i.e. kk-th row of ζ\zeta is a row of ZZ where the kk-th earliest appeared feature appears). Then ζ\zeta is a lower triangular matrix with diagonal components of one, having detζ=1\det\zeta=1. Hence UK×KU\in\mathbb{Z}^{K\times K} follows from c. \square

Theorem 6

Detailed Balance Condition

To prove the convergence of Algorithm 2, we first prove the detailed balance condition, a sufficient condition for MCMC kernel q(U|U)q(U^{\prime}|U) to keep the intended distribution P(U)P(U) invariant.

We factorize the MCMC kernel q(U|U)q(U^{\prime}|U) as

q(U|U)=q(U|U,k,𝒰)q(k)q(𝒰)q(U^{\prime}|U)=q(U^{\prime}|U,k,\mathcal{U})\,q(k)\,q(\mathcal{U})

where q(𝒰)q(\mathcal{U}) is a probability to obtain 𝒰\mathcal{U} by Algorithm 1, q(k)=1/Kq(k)=1/K, and

q(U|k,U,𝒰)=𝕀[UDk(U,𝒰)]P(U)UDk(U,𝒰)P(U),q(U^{\prime}|k,U,\mathcal{U})=\frac{\mathbb{I}\left[U^{\prime}\in D_{k}(U,\mathcal{U})\right]\cdot P(U^{\prime})}{\sum_{U^{\prime}\in D_{k}(U,\mathcal{U})}P(U^{\prime})},

where

Dk\displaystyle D_{k} (U,𝒰)\displaystyle(U,\mathcal{U})
={U+Δ𝒖|𝒖k+Δ𝒖𝒰,det(U+Δ𝒖)0}.\displaystyle=\left\{U+\Delta\bm{u}\middle|\bm{u}_{k}+\Delta\bm{u}\in\mathcal{U},\det(U+\Delta\bm{u})\neq 0\right\}.

We assume x=(k,𝒰)x=(k,\mathcal{U}) as auxiliary variables, and use Proposition 4 in [25] with

h(x|U,x,U)=δ(k,k)δ(𝒰,𝒰),h(x|U,x^{*},U^{\prime})=\delta(k,k^{*})\,\delta(\mathcal{U},\mathcal{U}^{*}),

and get an acceptance ratio

r(U,x,U,x)\displaystyle r(U;x^{*},U^{\prime},x) =P(U)q(U|U,k,𝒰)q(U|U,k,𝒰)P(U)\displaystyle=\frac{P(U^{\prime})\,q(U|U^{\prime},k^{*},\mathcal{U}^{*})}{q(U^{\prime}|U,k,\mathcal{U})\,P(U)}
={1,𝒖k𝒰,0,otherwise.\displaystyle=\begin{cases}1,&\bm{u}_{k}\in\mathcal{U},\\ 0,&otherwise.\end{cases}

Then Algorithm 2 satisfies a detailed balance condition. \square

Transitability

Since a value of q(U|U)q(U^{\prime}|U) becomes zero for some UU^{\prime}, we need to prove transitability for global convergence. By transitability we meant that: for arbitrary UU and U(0)U^{(0)}, a probability QT(U|U(0))Q_{T}(U|U^{(0)}) to obtain UU from initial state U(0)U^{(0)} after some finite TT MCMC iteration is non-zero.

Assume UU and U(0)U^{(0)} be integer regular matrices. Let V,V0V,V_{0} be a set of columns of U,U(0)U,U^{(0)}, respectively. If VV0V\neq V_{0}, we select 𝒗0V0\bm{v}_{0}\in V_{0} that is not in VV. There exists 𝒗0V\bm{v}_{0}^{\prime}\in V such that dimspan(V0\{𝒗0}{𝒗0})=K\dim\mathrm{span}(V_{0}\backslash\{\bm{v}_{0}\}\cup\{\bm{v}_{0}^{\prime}\})=K, because if there is no, 𝒗V,𝒗span(V0\{𝒗0}){}^{\forall}\bm{v}^{\prime}\in V,\bm{v}^{\prime}\in\mathrm{span}(V_{0}\backslash\{\bm{v}_{0}\}) yields dimspanV=dimspan(V0\{𝒗0})=K1\dim\mathrm{span}V=\dim\mathrm{span}(V_{0}\backslash\{\bm{v}_{0}\})=K-1 and it conflicts with dimspanV=rank(U)=K\dim\mathrm{span}V=\mathrm{rank}(U)=K. Then we set V1=V0\{𝒗0}{𝒗0}V_{1}=V_{0}\backslash\{\bm{v}_{0}\}\cup\{\bm{v}_{0}^{\prime}\}. Repeating above operation while ViVV_{i}\neq V, we obtain a sequence V0,V1,,VL=VV_{0},V_{1},\cdots,V_{L}=V, where LKL\leq K and Vi\{𝒗i}=Vi+1\{𝒗i}V_{i}\backslash\{\bm{v}_{i}\}=V_{i+1}\backslash\{\bm{v}_{i}^{\prime}\}.

If the distribution of ff^{*} in Algorithm 1 has a support covering {0,1,2,}\{0,1,2,\cdots\}, a probability of 𝒰\mathcal{U} to include both 𝒗i,𝒗i\bm{v}_{i},\bm{v}^{\prime}_{i} is non-zero. Therefore, in the ii-th MCMC step, a probability to obtain U(i)U^{(i)} (which have column entries ViV_{i}) from U(i1)U^{(i-1)} is non-zero. Hence after LL iteration, the probability QL(U,U(0))Q_{L}(U^{\prime},U^{(0)}) for UU^{\prime} having the same column entries as UU is non-zero.

Finally we consider a pseudo operation shuffling columns of U(i)U^{(i)} after each iteration (which is no effect at all, and even no need to be executed), we get QL(U,U(0))>0Q_{L}(U,U^{(0)})>0 for all regular integer matrices U,U(0)U,U^{(0)}. \square