{'fraction_non_alphanumeric': 0.04617201451696526, 'fraction_numerical': 0.015496284348176738, 'mean_word_length': 4.236651583710407, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Goodness-of-fit Test for Latent Block Models

Latent block models are used for probabilistic biclustering, which is shown to be an effective method for analyzing various relational data sets. However, there has been no statistical test method for determining the row and column cluster numbers of latent block models. Recent studies have constructed statistical-test-based methods for stochastic block models, which assume that the observed matrix is a square symmetric matrix and that the cluster assignments are the same for rows and columns. In this study, we developed a new goodness-of-fit test for latent block models to test whether an observed data matrix fits a given set of row and column cluster numbers, or it consists of more clusters in at least one direction of the row and the column. To construct the test method, we used a result from the random matrix theory for a sample covariance matrix. We experimentally demonstrated the effectiveness of the proposed method by showing the asymptotic behavior of the test statistic and measuring the test accuracy.


Introduction
Block modeling [1,2] is known to be effective in representing various relational data sets, such as the data sets of movie ratings [3], customer-product transactions [3], congressional voting [4], document-word relationships [5], and gene expressions [6]. Latent block models or LBMs [7] are used for probabilistic biclustering of such relational data matrices, where rows and columns represent different objects. For instance, suppose that a matrix A = (A ij ) ij ∈ R n×p represents the relationship between users and movies, where entry A ij is the rating of the j-th movie by the i-th user. In LBMs, we assume a regular-grid block structure behind the observed matrix A; i.e., both rows (users) and columns (movies) of matrix A are simultaneously decomposed into latent clusters. A block is defined as a combination of row and column clusters, and entries of the same block in matrix A are supposed to be i.i.d. random variables.
An open problem in using LBMs is that there has been no statistical criterion for determining the numbers of row and column clusters. Recently, statistical-test-based approaches [8,9,10] have been proposed for estimating the cluster number of stochastic block models (SBMs) [11]. SBMs are similar to LBMs in the sense that they assume a block structure behind an observed matrix; however, they are based on different assumptions from LBMs that an observed matrix is a square symmetric matrix and that the cluster assignments are the same for rows and columns [12]. In regard to the LBM setting, no statistical method has been constructed to determine row and column cluster numbers.
Aside from the test-based methods, several model selection approaches have been proposed based on cross-validation [13] or an information criterion [14,15,4]. However, these approaches have several limitations. (1) First, they cannot provide knowledge about the reliability of the result besides the finally estimated cluster numbers. Rather than minimizing the generalization error, in some cases, it is more appropriate to provide a probabilistic guarantee in reliability for the purpose of knowledge discovery. (2) Second, both the cross-validation-based and information-criterion-based methods depend on the clustering algorithm used. For instance, we can employ the Bayesian information criterion (BIC) for estimating the marginal likelihood only if the Fisher information matrix of the model is regular, which is not the case for block models. Constructing an information criterion that estimates the expectation of the generalization error for a wider class of models is generally difficult. (3) Finally, the above methods require relatively large computational complexity. Computation of an information criterion requires the process of approximating the posterior distribution by the Markov chain Monte Carlo (MCMC) method, and cross-validation requires the iterative calculation of the test error with different sets of partitions of the training and test data sets.
In this study, we proposed a new statistical test method for LBMs. To construct a hypothesis test with a theoretical guarantee, we used a result from random matrix theory. Recent studies on random matrix theory have revealed the asymptotic behavior of singular values of an n × p random matrix [16,17,18,19,20,21,22,23]. We assumed that if an observed matrix A is represented by the sum of block-wise mean effect P and i.i.d. noise Z with a standard deviation σ and its distribution satisfies certain conditions (e.g., Gaussian noise), the normalized maximum eigenvalue of 1 σ 2 Z Z converges in law to the Tracy-Widom distribution with index 1 [20,21,22,23]. Based on this result, we constructed a goodness-of-fit test for a given set of row and column cluster numbers of an LBM, using the maximum singular value of a standardized version of the observed matrix A (detailed description of the proposed method is presented in section 2). We proved that under the null hypothesis (i.e., observed matrix A consists of a given set of row and column cluster numbers), the proposed test statistic T converges in law to the Tracy-Widom distribution with index 1 (Theorem 4.1). We also showed that under the alternative hypothesis, test statistic T increases in proportion to m The proposed method solves the limitations of other model selection approaches. (1) Our statistical test method enables us to obtain knowledge about the reliability of the test results. When testing a given set of row and column cluster numbers, we can explicitly set the probability of Type I error (or false positive) as a significance level α. (2) Unlike the other model selection methods, the proposed method does not depend on the clustering algorithm as long as it satisfies the consistency condition (section 2). It only uses the output of a clustering algorithm to test a given set of cluster numbers; there is no need to modify the test method according to the clustering algorithm. (3) The proposed test method requires relatively small computational complexity. It does not require the MCMC procedure or partitioning into the training and test data sets. For these reasons, the proposed test-based method can be widely used for the purpose of knowledge discovery.
The next sections consist of the detailed explanation of the proposed test method for LBMs. In section 2, we describe the proposed goodness-of-fit test and its theoretical guarantee with the assumptions required for the problem setting. Next, we briefly review the related works and their differences from the proposed method in section 3. The main results are presented in section 4, where we prove the asymptotic properties of the proposed test statistic. In section 5, we experimentally demonstrate the effectiveness of the proposed test method by showing the asymptotic behavior of the test statistic and calculating the test accuracy. We discuss the results and limitations of the proposed method in section 6 and conclude the paper in section 7.
2 Problem setting and statistical model for goodness-of-fit test for latent block models Let A ∈ R n×p be an observed matrix with n × p entries. In this study, we assume n ≥ p; however, if n < p, we can still use the proposed method by applying it to the transpose matrix A instead of the original matrix A. In LBM, we assume that each entry of matrix A is generated based on its block or the combination of row and column clusters. Let g (1) i be the cluster of the i-th row of matrix A, and let g (2) j be the cluster of the j-th column of matrix A. To apply the result in [23], we assume that each entry A ij of A is generated as the sum of block-wise mean P ij and i.i.d. noise Z ij with zero mean and standard deviation σ. We assume that the noise distribution is symmetric (thus, E[|Z ij | 2k+1 ] = 0 holds for all k ∈ N) and that its moments are sub-Gaussian (i.e., E[|Z ij | 2k ] ≤ (C 0 k) k holds for all k ∈ N with some constant C 0 > 0). For the following discussion, we explain the case where the noise is generated by a normal distribution with zero mean and standard deviation σ for simplicity. Such a model can be applied to relational data matrices with continuous entries, and it has been called the Gaussian latent block model [24,25]. Here, B kh is the mean of entries in the (k, h)-th block under the null hypothesis: Let K and H, respectively, be the true cluster numbers for rows and columns of an observed matrix A. Our purpose is to estimate these cluster numbers from observed matrix A. To achieve this goal, we propose a goodness-of-fit test for a given set of cluster numbers (K 0 , H 0 ), with which we test whether the true cluster number (K, H) is (K 0 , H 0 ), or at least one of the given cluster numbers (K 0 or H 0 ) is smaller than the true cluster numbers (K or H). Formally, the null (N) and alternative (A) hypotheses are expressed as By sequentially testing the cluster numbers in the following order, we can estimate the cluster numbers of a given matrix A.
Testing the sets of cluster numbers in the above order requires 3n tests at most, with the assumption of n ≥ p.
In this study, we derive the test statistics under the following setting: (i). Let n and p, respectively, be the number of rows and columns of matrix A. We assume that both n and p increase in proportion to some sufficiently large number m (n, p ∝ m).
(ii). Let (K, H) be the minimum cluster numbers to represent the block structure of matrix A under the null hypothesis. We assume that the minimum number of rows and columns, which we denote as n min and p min , respectively, satisfies the following condition: In other words, we assume that there is no "too small" block in matrix A.
(iii). If the set of cluster numbers (K 0 , H 0 ) for estimating the cluster structure of matrix A is equal to the true cluster numbers (K, H), then we call it a realizable case. Otherwise (if at least one of them is smaller than the true cluster numbers), we call it an unrealizable case. We only consider the cases where K 0 ≤ K and H 0 ≤ H.
(iv). In the realizable case, we assume that a clustering algorithm for estimating the block structure of matrix A is consistent; i.e., the probability that its correctly recovers the true block structure converges to 1 in the limit of m → ∞. For instance, clustering algorithms in [26,27,28] have been proven to be consistent. By this assumption, the proposed method is independent of a specific clustering algorithm.

Relation to existing works
In this section, we briefly review the related works and explain the differences between them and the proposed method.

Model selection for block models
Statistical-test-based methods (for SBM) Recently, several methods have been proposed for testing the properties of a given observed matrix in relation to SBMs [8,9,29,10,30]. Particularly, the methods proposed in [8,9,10] have enabled us to estimate the number of blocks for SBMs. However, these methods differ from ours in the problem setting; they can be applied only to an SBM setting, where an observed matrix is a square symmetric matrix, and the cluster assignments are the same for rows and columns. There has been no method for estimating the block number for LBMs, where rows and columns (not necessarily square) of an observed matrix are simultaneously decomposed into clusters.
Cross-validation-based methods Cross-validation is a widely used method for model selection, where a data set is first split into training and test data sets, and then the best model with the minimum test error is determined. Recently, cross-validation methods for matrix data have been proposed [31,32,33,13] to determine the number of clusters in network data. Although the purpose of these methods and our method is similar, these methods differ from ours in that their target is the network data, where the observed matrix is square and its rows and columns represent the same node sets. Thus, the block structure is symmetric regardless of whether the network itself is directed or undirected). Moreover, unlike a statistical test, these methods cannot provide quantitative knowledge about the reliability of the selected model. Furthermore, the computational cost of cross-validation is generally high because it requires the iterative calculation of the test error with different data set partitions.
Information-criterion-based methods Another approach for determining the number of blocks in a matrix is to estimate the generalization error or marginal likelihood by some information criterion for given sets of block numbers. By using such information criteria, we can select a model in a statistically meaningful (non-heuristic) way. In regard to block models, many variants of BIC have been proposed [14,15,4,34,35]. Unlike our test-based method, which only requires a clustering algorithm to satisfy the consistency condition (section 2), an information criterion for a theoretical guarantee should be carefully chosen according to the given clustering algorithm. For instance, BIC can be employed for estimating the marginal likelihood only if the Fisher information matrix of the model is regular, which is not the case for block models.
To solve this problem, as an alternative criterion to BIC, the integrated completed likelihood (ICL) criterion has been used in many studies for estimating the number of blocks in LBMs [25,36,37]. In ICL, we first derive a marginal likelihood for a given set of an observed matrix and block assignments and then substitute the set of estimated block assignments to approximate the marginal likelihood. However, since ICL is computed based on a single estimator of block assignments, there is no guarantee for the goodness of the approximation of marginal likelihood.
Similar to cross-validation-based methods, information-criterion-based methods cannot provide a probabilistic guarantee for the reliability of the selected model, which is a disadvantage for the purpose of knowledge discovery. The computational cost also becomes a problem because the computation of an information criterion requires the process of approximating the posterior distribution by MCMC.
Other model selection methods Aside from the information criteria, several studies have proposed to determine the number of blocks in LBMs based on the co-clustering adjusted rand index [38], the extended modularity for biclustering [39], or the expected posterior loss for a given loss function [40]. Another approach is to define the posterior distribution not only on cluster assignments of rows and columns but also on row and column cluster numbers [41,42]. Unlike the model selection approaches, such nonparametric Bayesian methods can estimate the distribution of the block numbers. The bestfitted number of the blocks can be determined based on the posterior distribution (e.g., we can choose a MAP estimator [42]). However, in this case, the computational cost of MCMC is higher than that of the information-criterion-based methods because it requires a large number of iterations to approximate the posterior distribution both on the block assignments and the number of blocks.

Test statistic for determining the set of cluster numbers
Under the null hypothesis, if the block-wise mean effect P and the standard deviation σ in (1) are given, by subtracting P from the original matrix A and dividing it by σ, we obtain a normalized matrixÃ * , which is given byÃ Based on this normalization, each entry ofÃ * in (4) independently follows the standard normal distribution N (0, 1). Therefore, according to the results in [23] and the assumption of n = n(p) ≥ p and n/p → γ ≥ 1 in the limit of p → ∞, the distribution of the normalized maximum eigenvalue T * of matrix (Ã * ) Ã * converges to the Tracy-Widom distribution with index 1 (T W 1 ) in the limit of p → ∞: where λ 1 is the maximum eigenvalue of matrix (Ã * ) Ã * and Typically, we do not know the true cluster numbers (K, H) or the true cluster assignments g (1) and g (2) , and we can only estimate their values based on the observed matrix A. Given a set of cluster numbers (K 0 , H 0 ), suppose that we have obtained estimated cluster assignmentsĝ (1) (for rows) and g (2) (for columns). Based on the clustering result (ĝ (1) ,ĝ (2) ), we estimate the block-wise mean bŷ where Here, I k is the set of indices of rows that are assigned to the k-th cluster, and J h is the set of indices of columns that are assigned to the h-th cluster. The consistency assumption (iv) guarantees that if the cluster numbers (K 0 , H 0 ) are set to the true numbers (K, H), the probability that above cluster assignments (I k ) k and (J h ) h are correct converges to 1 in the limit of m → ∞. By using the estimated block-wise meanP , we estimate the standard deviation of each entry of matrix A bŷ With the estimatorsP andσ, we define a matrixÃ as follows: We define the test statistic T for estimating the row and column cluster numbers as the normalized maximum eigenvalue of matrixÃ Ã : whereλ 1 is the maximum eigenvalue of matrixÃ Ã , and a and b are defined as in (6).
Theorem 4.1 (Realizable case). We assume that the matrix size satisfies the following condition: n = n(p) ≥ p and n/p → γ ≥ 1 in the limit of p → ∞. Suppose (K 0 , H 0 ) = (K, H). Under the assumption (iv) on the consistency of the clustering algorithm, in the limit of p → ∞, where T is defined as in (12).
Proof. We denote the operator norm as · op , In the subsequent discussion, we use the following notation: • q: the constant element in the same (k, h)-th block in matrix P : q ≡ B kh .
•q: the constant element in the same (k, h)-th block in matrixP , which is calculated based on the observed matrix A and the correct cluster assignments:q ≡B kh .
•q: the constant element in the same (k, h)-th block in matrixP (orB kh ), where the block (9) is given by an output of the clustering algorithm that satisfies the consistency assumption (iv), and P andB are given by the sample mean of the (k, h)-th block of the observed matrix A according to (7):q ≡B kh .
We first derive the joint probability of the event F m that the clustering algorithm outputs the correct cluster assignments (i.e.,q =q) and event G m,C that |q −q| ≤ C/m holds. The joint probability is expressed by Since the number of entries in the block is proportional to m 2 (assuming (ii)), from the central limit theorem, |q −q| √ m 2 converges to a normal distribution N (0, σ 2 ). From Prokhorov's theorem [43], Moreover, according to the consistency assumption (iv), if (K 0 , H 0 ) = (K, H), the probability of event F m converges to 1 in the limit of m → ∞. Combining these results with (15), we obtain In other words, with a high probability, the clustering algorithm detects the block structure correctly and |q −q| ≤ C/m holds in the limit of m → ∞. Furthermore, Combining the above inequality and (17), we obtain where 1 ≡ KH . By using the fact that Frobenius norm upper bounds operator norm, we obtain From (19) and (20), From (5), the maximum eigenvalue λ 1 of matrix 1 σ 2 Z Z is in the order of O p (m). In other words, the largest singular value of matrix Z is in the order of O p ( √ m). By the subadditivity of the operator norm and A −P op = Z + (P −P ) op , we obtain By (21) and (22), In the second equation of (23), we usedσ = σ + O p  • Based on the consistency assumption (iv), the probability of the eventσ MLE =σ MLE converges to 1 in the limit of m → ∞.
Here, Z op is equal to the largest singular value of Z, which is the square root of the maximum eigenvalue of Z Z. Therefore, the maximum eigenvalue of 1 σ 2 Z Z is equal to 1 σ 2 Z 2 op . By combining these results and (5), from Slutsky's theorem (note that both a and b are constants determined only by the matrix size), This is equivalent to the statement of Theorem 4.1.
where T is defined as in (12).
Proof. In the subsequent discussion, we use the following notation (see also Figure 1): • P : A matrix that consists of the true block structure and whose entries are the population block-wise means.
•P : A matrix that consists of the estimated block structure and whose entries are the population block-wise means, which can be calculated using P and the estimated block structure.
•P : A matrix that consists of the estimated block structure and whose entries are the block-wise sample means, which can be calculated using the observed matrix A and the estimated block structure.
We first consider the relationship between matrices P andP under the assumption of unrealizability (K 0 < K or H 0 < H). We assume the case of K 0 < K without loss of generality. In this case, for all k ∈ {1, · · · , K}, at least one estimated block contains n k K0 or more rows whose true block is k (Figure  1 left). Since K 0 < K, at least one estimated block contains two or more row sets that belong to mutually different blocks in the true block structure, and each of them contains nmin K0 or more rows. Therefore, in the estimated block structure, there exists at least one block, which consists of two or more submatrices (defined by mutually different blocks in the true block structure) with a size of at least nmin K0 × pmin H0 (Figure 1 right). Let X 1 and X 2 be the submatrices and let q 1 and q 2 be the entries of the blocks defined by the row and column indices of X 1 and X 2 , respectively. Here, we assume q 1 > q 2 without the loss of generality. In matrixP , both of these submatrices have the same entriesq.
. Therefore, for anyq, there exists at least one submatrix (X in Figure 1) in P with a size of at least nmin K0 × pmin H0 , where all the entries are q and Next, in regard to the relationship between matricesP andP , where Z is a submatrix of Z with the size of n 1 × p 1 , and u 1 = (1, 1, · · · , 1) ∈ R n1 and u 2 = (1, 1, · · · , 1) ∈ R p1 . Here,q is the constant entry of the block in matrixP , which includes the same row and column indices as Z . To derive the final equation in (27), we used the fact that Z 2 op is equal to the maximum eigenvalue of Z Z, which is O p (m) from (5). Since Z does not depend on the cluster structure, (27) holds uniformly over all clustering algorithms and clustered blocks.
Here, for all q,q, andq, the following inequality holds: By combining (27) and (28), we obtain The only non-zero (and thus, the largest) singular value of matrixX is |q −q| √np , wheren andp are the row and column sizes ofX, respectively. Since the largest singular value of a matrix is equal to its operator norm and the operator norm of a submatrix is not larger than that of the original matrix, By combining (29) and (30), we have Therefore, from (26), we obtain By the subadditivity of the operator norm and A −P op = Z + (P −P ) op , we obtain From (32) and the fact that the order of Z op is O p ( √ m) (based on the same discussion as the proof of Theorem 4.1), we obtain Next, we check the order of the estimated standard deviationσ, which is given bŷ By the subadditivity of Frobenius norm, we obtain Here, Z F = √ npσ MLE , whereσ MLE is the maximum likelihood estimator of standard deviation σ in the realizable case. Based on the Cramér's theorem ( Here, to derive the second to last equation, we used the assumption that the number of row and column blocks are finite constants and (27) holds for the block with the maximum difference betweenP andP .
In the final equation, we used the fact that max k,h,k ,h |B kh − B k h | is bounded by a finite constant. Based on (36) and (37), we obtain Z + (P −P ) F = O p (m). By using (35), we finally obtain the order of the estimated standard deviation:σ = O p (1). By combining this fact and (34), we obtain Here, 1 σ A −P op is equal to the largest singular value of matrixÃ in (11). Therefore, the maximum eigenvalueλ 1 of matrixÃ Ã , which is equal to 1 σ A −P op 2 ), satisfies The test statistic is T =λ 1−a b . Using the definition (6), we obtain a = O p (m) and where we used the assumption n ≥ p and definitions C 1 ≡ n/m and C 2 ≡ p/m. By combining these results and (40), we obtain , which concludes the proof.
, where T is defined as in (12).
Proof. We define P ,P , andP as in Theorem 4.2. By using the subadditivity of operator norm and A −P op = Z + (P −P ) op , we have A −P op ≤ Z op + P −P op . Combining this result with the fact that Frobenius norm upper bounds operator norm, we obtain Based on the same discussion as in Theorem 4.2, we obtain Z F = O p (m), P −P F = O p (m) and thus the following equation holds: Next, we check the order of 1/σ, which is given by By combining (34) and the fact that Frobenius norm upper bounds operator norm, we obtain Therefore, the following equation holds: By using (45) and (48), we obtain 1 σ A −P op = O p (m). Here, 1 σ A −P op is equal to the largest singular value of matrixÃ in (11). Therefore, the maximum eigenvalueλ 1 of matrixÃ Ã , which is equal to Using the definition (6), we obtain a = O p (m) and where we used the assumption n ≥ p and definitions C 1 ≡ n/m and C 2 ≡ p/m. According to Theorems 4.1, 4.2 and 4.3, we perform a one-sided test on the test statistic T . We define the goodness-of-fit test for a set of cluster numbers (K 0 , H 0 ) at the significance level of α by where t(α) is α upper quantile of the Tracy-Widom distribution with index 1. By sequentially testing the cluster numbers based on the above rejection rule (50) as in (2), we obtain a set of estimated cluster numbers (K,Ĥ) for a given matrix A.

Convergence of test statistic T to T W 1 distribution
Before applying our test method to relational data sets, we confirmed that test statistic T in (12) converges in law to the Tracy-Widom distribution with index 1, if we set the cluster numbers (K 0 , H 0 ) at the true cluster numbers of (K, H).
In the experiment, we applied the proposed statistical test in section 2 to synthetic data matrices with the ground truth block structure. We defined the ground truth cluster numbers (K, H), block-wise mean B, and the standard deviation σ as To generate data matrix A, we first defined the ground truth cluster assignments g (1) and g (2) . For each row of matrix A, we randomly chose its cluster index from the discrete uniform distribution on {1, 2, 3}. Similarly, we chose the cluster index of each column of matrix A from the discrete uniform distribution on {1, 2}. Then, we generated each entry of matrix A based on its block. Specifically, the (i, j)-th entry of matrix A was generated from a normal distribution N (B g (1) i ,g (2) j , σ). We estimated the cluster structure of each matrix A based on a hierarchical clustering algorithm (specifically, we used the Ward's method [44]), and based on the clustering result, we calculated the test statistic T in (12). Figure 2 shows the relationships between the normalized histograms of T for different matrix sizes (n, p) and the approximated Tracy-Widom distribution with index 1. To plot the T W 1 distribution, we used the following approximation [45]: where Γ(·) is Gamma function and k = 46.446, s = 0.186054, c = −9.84801.
(53) Figure 2 shows that in the limit of n → ∞, n/p → 2, test statistic T converges in law to the T W 1 distribution. We also checked the tail probability of test statistic T . Let t(α) be the α upper quantile of the T W 1 distribution. According to Table 2 in [46], t(0.01) ≈ 2.02345, t(0.05) ≈ 0.97931, and t(0.1) ≈ 0.45014. Figure 3 shows the ratios of the trials where T ≥ t(0.01), T ≥ t(0.05), and T ≥ t(0.1). In Figure 3, we see that the ratios, respectively, converge to 0.01, 0.05, and 0.1 as the matrix size (n, p) gets larger.

Behavior of test statistic T in unrealizable case
We also confirmed the behavior of test statistic T in the unrealizable case. Recall that in Theorems 4.2 and 4.3, if K 0 < K or H 0 < H, where T is defined as in (12).
In the experiment, we defined the ground truth cluster numbers (K, H), block-wise mean B, and standard deviation σ as To generate data matrix A, we first defined the ground truth cluster assignments g (1) and g (2) . For each row of matrix A, we randomly chose its cluster index from the discrete uniform distribution on {1, 2, 3, 4}. Similarly, we chose the cluster index of each column of matrix A from the discrete uniform distribution on {1, 2, 3}. Then, we generated each entry of matrix A based on its block. Specifically, the (i, j)-th entry of matrix A was generated from a normal distribution N (B g (1) i ,g (2) j , σ). For each matrix A, we applied a hierarchical clustering algorithm [44] to estimate its cluster structure. Figure 4 shows the behavior of the test statistic T for different matrix sizes (n, p). In Theorems 4.2 and 4.3, for large matrix sizes, test statistic T increases in proportion to n 5 3 for all settings of the block numbers, where K 0 < K or H 0 < H.

Accuracy of test for matrix data with ground truth block structures
To evaluate the accuracy of the proposed goodness-of-fit test for LBMs, we applied it to matrix data with ground truth block structures. We defined the ground truth cluster numbers as (K, H) = (4, 3) and the standard deviation as σ = 0.1. In regard to block-wise mean B, we tried 10 settings: Overall, for each of the above 100 settings (10 types for B and 10 types for (n, p)), we generated data matrix A for 1, 000 times, estimated their cluster structures, computed test statistic T , and performed the proposed goodness-of-fit test. The method for generating data matrices is the same as that used in section 5.2. Figure 5 shows the examples of the generated matrix A with the size of (n, p) = (40,30) for t = 1, · · · , 10, where the block-wise mean matrix B is defined as in (56) for a fixed t.
We estimated the cluster structure of each matrix A based on a hierarchical clustering algorithm [44]. Based on the clustering result, we calculated the test statistic T in (12) and applied our statistical test using a significance level of α = 0.01. Figure 6 shows the accuracy of tests with different matrix sizes. Each plotted line shows the ratio of trials where (K 0 , H 0 ) = (K, H), with a fixed block-wise mean B. For every setting of matrix B, the accuracy of the test increased with the matrix size. With the same matrix size (n, p), the accuracy decreased with the narrower range of the entries in matrix B because it gets more difficult to distinguish one cluster from another. Figure 6 shows that the test accuracy of 80% was attained with matrix size (n, p) ≥ (140, 105) and with the range of B equals to or broader than (0.34, 0.66) (that is, max k,h,k ,h |B kh − B k h | ≥ 0.32).

Discussion
In this section, we discuss the proposed test method in terms of the test statistic and the conditions for the generative model.
There are many test statistics for testing a given set of row and column cluster numbers. Although the most powerful test does not exist in our problem setting, where alternative hypothesis includes multiple models, it might be possible to construct a better test statistic in terms of the convergence rate. For instance, it might be better to combine the information of all the singular values, instead of using only the largest singular value. In this paper, Theorem 4.1 states the asymptotic behavior of the proposed test statistic. In practice, we typically obtain an observed matrix with a finite size. For a Gaussian case (i.e., each entry of a matrix independently follows N (0, 1)), the following statement holds [47]: Suppose n = n(p) > p and n/p → γ ∈ [1, ∞) in the limit of p → ∞. Then, for any s 0 , there exists N 0 ∈ N such that when max(n, p) ≥ N 0 and max(n, p) is even, for all s ≥ s 0 , where T * is defined as in (5) and C(·) is a continuous and non-increasing function. From the above inequality (57), if the clustering algorithm outputs the correct block assignments, the convergence rate of the normalized maximum eigenvalue T * of matrix (Ã * ) Ã * (whereÃ * is defined as in (4)) to the Tracy-Widom distribution with index 1 is O(m −2/3 ). However, since the distribution of T is unknown in the case where the correct block assignment is not obtained, the convergence rate of T is also unknown. Deriving the convergence rate of T by considering the above discussion is a future research topic.
In regard to the conditions for using the proposed test method, we assumed that an observed matrix is represented by the sum of a block-wise mean effect and a random noise matrix, whose entries independently and identically follow a symmetric sub-Gaussian distribution. The symmetric property is required as a sufficient condition for the test statistic convergence [23]; however, many distributions that are often used in practice do not satisfy this condition (e.g., binomial or multinomial distribution). Furthermore, there are proposed variants of latent block models with which we assume different block structures from a regular grid [48,49]. To construct test methods for the above settings is an important topic for future research. Another assumption in the proposed method is that for all the blocks, the size of rows and columns increases in proportion to m. In practice, there may be some cases where it is more appropriate to assume that the number of blocks increases with the matrix size. It will be useful to construct a test statistic, which is applicable to such cases and is dependent on the order of row and column cluster numbers.

Conclusion
Latent block models are effective tools for biclustering, where rows and columns of an observed matrix are simultaneously decomposed into clusters. Such a bicluster structure appears in various types of relational data, such as the customer-product transaction data or and the document-word relationship data. One open problem in using latent block models is that there has been no statistical test method for determining the number of blocks. In this study, we developed a goodness-of-fit test for latent block models based on a result from the random matrix theory. Since the maximum singular value of a symmetric sub-Gaussian random matrix with zero mean and unit variance converges in law to a Tracy-Widom distribution, we defined the test statistic as the normalized maximum singular value of an observed matrix standardized with the estimators of a block-wise mean and a standard deviation. We first proved that the proposed test statistic converges in law to the Tracy-Widom distribution under the consistency condition, if both row and column cluster numbers are set at the true numbers. We also showed that the test statistic increases in the order of O p (m asymptotic behavior of the test statistic and the test accuracy by using synthetic data matrices with ground truth block structures.