## Empirical Bayes Linked Matrix Decomposition

Eric F. Lock 1*

1* Division of Biostatistics and Health Data Science, School of Public Health, University of Minnesota, Minneapolis, 55455, MN, USA.

Corresponding author(s). E-mail(s): elock@umn.edu;

## Abstract

Data for several applications in diverse fields can be represented as multiple matrices that are linked across rows or columns. This is particularly common in molecular biomedical research, in which multiple molecular 'omics' technologies may capture different feature sets (e.g., corresponding to rows in a matrix) and/or different sample populations (corresponding to columns). This has motivated a large body of work on integrative matrix factorization approaches that identify and decompose low-dimensional signal that is shared across multiple matrices or specific to a given matrix. We propose an empirical variational Bayesian approach to this problem that has several advantages over existing techniques, including the flexibility to accommodate shared signal over any number of row or column sets (i.e., bidimensional integration), an intuitive model-based objective function that yields appropriate shrinkage for the inferred signals, and a relatively efficient estimation algorithm with no tuning parameters. A general result establishes conditions for the uniqueness of the underlying decomposition for a broad family of methods that includes the proposed approach. For scenarios with missing data, we describe an associated iterative imputation approach that is novel for the single-matrix context and a powerful approach for 'blockwise' imputation (in which an entire row or column is missing) in various linked matrix contexts. Extensive simulations show that the method performs very well under different scenarios with respect to recovering underlying low-rank signal, accurately decomposing shared and specific signals, and accurately imputing missing data. The approach is applied to gene expression and miRNA data from breast cancer tissue and normal breast tissue, for which it gives an informative decomposition of variation and outperforms alternative strategies for missing data imputation.

Keywords: Data integration, dimension reduction, low-rank factorization, missing data imputation, variational Bayes

## 1 Introduction

Low-rank matrix factorization techniques are widely used for many machine learning tasks such as data compression and dimension reduction, denoising to approximate underlying signal, and matrix completion. Often, instead of a single matrix, the data for a given application takes the form of multiple matrices that are linked by rows and/or columns. For example, in Section 9 we consider data from The Cancer Genome Atlas, in which gene expression (mRNA) and microRNA (miRNA) data are available from different platforms for breast cancer tumor samples and normal breast tissue samples from unrelated individuals. Thus, the data can be represented as four matrices: (1) mRNA for cancer tissues, (2) mRNA for normal tissues, (3) miRNA for cancer tissues, and (4) miRNA for normal tissues; these matrices share both rows (here, molecular features) columns (here, samples), giving a bidimensionally linked structure. Furthermore, these data have column-wise missingness, in which molecular data for one of the platforms (mRNA or miRNA) is entirely missing for some samples in each of the cancer and control cohorts. We are interested in imputing the missing data, as well as identifying patterns of variation that are shared or specific to the different data sources and sample cohorts; for these tasks, a low-rank factorization approach that integrates all of the matrices together while accounting for their differences is attractive.

There is a growing and dynamic literature on multi-matrix low-rank factorization approaches. The majority of existing methods assume that just one dimension, rows or columns, is shared across matrices, i.e., unidimensional integration. For example, the Joint and Individual Variation Explained (JIVE) approach (Lock et al., 2013; O'Connell and Lock, 2016) decomposes low-rank structure that is shared across matrices from variation that is specific to each matrix. Several related approaches yield a similar decomposition using different strategies or assumptions (Feng et al., 2018; Schouteden et al., 2014; Yang and Michailidis, 2016; Tang and Allen, 2021), while other approaches such as structural learning and integrative decomposition (SLIDE) (Gaynanova and Li, 2019) allow for structure that is partially shared across any subset of matrices. Such an approach is efficient and useful for interpretation, as for many applications identifying shared, partially shared, and specific underlying structure is of interest.

The bidimensional linked matrix factorization (BIDIFAC) approach (Park and Lock, 2020; Lock et al., 2022) was developed to accommodate bidimensionally linked integration scenarios, by identifying structure that is shared or specific across both row sets and column sets. The objective function for BIDIFAC involves a structured nuclear norm penalty, which approximates the data as a decomposition of low-rank modules that are shared across row or column subsets. This approach is attractive in the bidimensionally linked context because it has favorable convergence properties and generally leads to a uniquely-defined decomposition (Lock et al., 2022); methods for unidimensional integration without shrinkage penalties typically require orthogonality (Lock et al., 2013; Feng et al., 2018; Gaynanova and Li, 2019; Yi et al., 2023) or other constraints to assure uniqueness, which do not extend to the bidimensionally linked context. Furthermore, the application of the structured nuclear norm penalty for

unidimensional integration (termed UNIFAC) has been shown to outperform alternative approaches with respect to accurately decomposing shared and specific structures (Lock et al., 2022) and imputing missing data (Park and Lock, 2020) under certain scenarios. However, a drawback is that it tends to over-penalize the estimated structure, and is thus biased toward an over-conservative solution with respect to the underlying low-rank structure (Yi et al., 2023).

For a single matrix the tendency of the nuclear norm to over-shrink low-rank signal has been well-documented, and in contrast low-rank approximations without any further shrinkage tend to overfit (Josse and Sardy, 2016). Empirical Bayes approaches to low-rank factorization are an attractive compromise, in which the optimal levels of shrinkage are determined empirically within a model-based framework (Wang and Stephens, 2021). Under a variational Bayes approach, the true posterior distribution for Bayesian model is approximated by minimizing the free energy (reducing to Kullback-Leibler divergence) under simpler distributional assumptions (Fox and Roberts, 2012). In particular, an empirical variational Bayes approach to low-rank matrix factorization under an assumption of normally distributed factor matrices has desirable theoretical properties, and the resulting approximation for a scenario without missing data can be efficiently computed via a closed-form shrinkage formula on the observed singular values of the matrix (Nakajima et al., 2012, 2013).

In this article, we describe an empirical variational Bayesian approach to identify and decompose shared, partially shared, or specific low-rank structure in the context of unidimensionally or bidimensionally linked matrices. Moreover, we describe an accompanying iterative imputation approach analogous to an expectation-maximization (EM) algorithm that can handle entry-wise missing or row/column-wise missing values. The proposed approach is free of any tuning parameters, and is shown to share desirable properties of the BIDIFAC approach while addressing the issue of over-shrinkage.

## 2 Notation and Setting

Throughout this article bold uppercase characters (' X ') correspond to matrices and bold lowercase characters (' x ') correspond to vectors. Square brackets index the entries within a vector or matrix, e.g., ' X [ m,n ]'. The singular value decomposition (SVD) of a matrix X : M × N is given by U D V X X T X , where the diagonal entries of D X , D X [ r, r ], are the singular values of X ordered from largest to smallest for r = 1 , . . . , min( M,N ). The Frobenius norm is given by || · || F , so that || X || 2 F is the sum of the squared entries in X ; the nuclear norm is given by || · || ∗ , which is the sum of the singular values in a matrix:

$$| | \mathbf X | | _ { * } = \sum _ { r = 1 } ^ { \min ( M, N ) } D _ { \mathbf X } [ r, r ].$$

We consider the setting of a collection of matrices X ij : M i × N j for i = 1 , . . . , I and j = 1 , . . . , J . The column sets N ,. . . , N 1 J and row sets M ,... , M 1 L are consistent across the matrices. The subscript ' ' defines concatenation across row or column sets, · and the full data X ·· can be represented as a single M × N matrix where M =

M 1 + . . . + M I and N = N 1 + . . . + N J :

$$X _ { \centerdot } = \begin{bmatrix} X _ { 1 1 } & X _ { 1 2 } & \dots & X _ { 1 J } \\ \vdots & \vdots & \vdots & \vdots \\ X _ { I 1 } & X _ { I 2 } & \dots & X _ { I J } \end{bmatrix} \text{ where } X _ { i j } \text{ are } M _ { i } \times N _ { j }.$$

As for our motivating application in Section 9, the I column sets may correspond to different sample cohorts and the J row sets may correspond to features from different high-dimensional sources.

## 3 Single-matrix results

Here we review some established results for a single matrix that are critical to what follows. Propositions 1 and 2 describe solutions to the least squares matrix approximation problem under a rank-constrained or nuclear norm penalized criterion, respectively.

Proposition 1. For X : M × N , the minimizer of the least squares objective || X - || S 2 F under the constraint rank ( S ) = R ≤ min ( M,N ) is given by S = U D V X S T X where D S is diagonal with D S [ r, r ] = D X [ r, r ] for r = 1 , . . . , R and D S [ r, r ] = 0 for r &gt; R .

Proposition 2. For X : M × N , the minimizer of the nuclear norm penalized objective 1 2 || X -S || 2 F + λ || S || ∗ is given by S = U D V X S T X where D S is diagonal with D S [ r, r ] = max( D X [ r, r ] -λ, 0) for all r .

Proposition 1 is well-known, and a proof of Proposition 2 can be found in Mazumder et al. (2010). Because of how they operate on the singular values X , Proposition 1 represents a hard-thresholding and Proposition 2 a soft-thresholding approach to low-rank approximation. Proposition 3 below establishes that the soft-thresholding operator also solves an objective with L 2 -penalties on the factor components of S , while Corollary 4 further establishes that it is the posterior mode of a Bayesian model with normal priors on the factor components.

Proposition 3. For X : M × N ,

$$\underset { \substack { m i n \\ \iota \colon M \times H, \mathbb { V } \colon N \times H } } { \min } \| \mathbf X - \text{UV} ^ { T } \| _ { F } ^ { 2 } + \lambda ( \| \mathbf U \| _ { F } ^ { 2 } + \| \mathbf V \| _ { F } ^ { 2 } ) = \underset { \mathbf S \colon r a n k ( \mathbf S ) \leq H } { \min } \| \mathbf X - \mathbf S \| _ { F } ^ { 2 } + 2 \lambda \| \mathbf S \| _ { * }$$

Moreover, ̂ S = ̂ ̂ UV T , where ̂ S solves the right-hand side of (2) and { ̂ U V , ̂ } solves the left-hand side of (2) .

Corollary 4. For X : M × N , if U : M × H and V : N × H have independent Normal (0 , σ 2 /λ ) entries, and X = UV T + E where E has independent Normal (0 , σ 2 ) entries, then the mode of the posterior distribution p ( U V , | X ) satisfies ̂ S = ̂ ̂ UV T , where ̂ S solves the right-hand side of (2) .

A proof of Proposition 3 can be found in Mazumder et al. (2010), and Corollary 4 follows by noting that the log of the posterior density for the specified model is proportional to the left-hand side of (2). While the results hold for general H , in what follows we set H as the largest possible rank H = min( M,N ); the actual estimated rank will tend to be smaller due to singular value thresholding (as in Proposition 2). The distribution of singular values for the error matrix E has been well-studied, including the following results on the first singular value under general conditions from Rudelson and Vershynin (2010):

Proposition 5. If E : M × N has independent entries with mean 0 , variance σ 2 and finite fourth moment, then D [1 , 1] ≈ σ ( √ ( M ) + √ N ) as M,N → ∞ ; further, if the entries are Gaussian, E ( D [1 , 1]) ≤ σ ( √ ( M ) + √ N ) for any M,N .

These results can motivate choice of the nuclear norm penalty λ . A conservative approach, if the error variance σ 2 is known, is to select λ = σ ( √ M + √ N ); per Propositions 5 and 2, this will tend to shrink to 0 any singular values in X that result from error with no true low-rank signal. Another approach is to adopt an empirical Bayes estimate for λ based on the model in Corollary 4. However, a general limitation of the nuclear norm approach to low-rank matrix reconstruction is the use of a universal threshold on all singular values; in practice larger singular values should be penalized less, as they capture a relatively lower ratio of error and higher signal (Shabalin and Nobel, 2013). An analogous interpretation is that the entries of U and V should not have identical variance across columns, and this can be relaxed by a model in which they have column-specific variances: U [ · , r ] ∼ Normal( 0 , τ 2 U ,r I ) and V [ · , r ] ∼ Normal( 0 , τ 2 V ,r I ) for r = 1 , . . . , H . 1

For given τ = { τ U ,r τ V ,r } H r =1 , the approximation resulting from the posterior mode for p ( U V , | X , τ ) is given by thresholding the singular values of X , analogous to Corollary 4. To determine τ , directly maximizing the joint likelihood p ( X U V , , | τ ) over U V , and τ degenerates to τ → 0 in the global solution, though it can have a reasonable local mode that is a singular value thresholding operator (Nakajima and Sugiyama, 2014). Further, a marginal empirical Bayes estimator that optimizes p ( X | τ ) (effectively integrating over U and V ) is difficult to obtain computationally or analytically. As an alternative, one can use an empirical variational Bayes (EVB) approach. In general, for a variational Bayes approach (Attias, 1999), the true posterior of a Bayesian model with parameters Θ and data y , p (Θ | y ) is approximated by another distribution q (Θ) under simplifying assumptions to minimize the free energy

$$F ( q ) = E _ { q } \log \frac { q ( \Theta ) } { p ( \Theta, \mathbf y ) } = E _ { q } \log \frac { q ( \Theta ) } { p ( \Theta \, | \, \mathbf y ) } - \log p ( \mathbf y ),$$

where E q ( ) · is the expected value with respect to the distribution defined by q ( ). · In the context of low-rank approximation, a useful simplifying assumption for q ( ) is that · the left and right factor matrices are independent, q ( U V , ) = q u ( U ) q v ( V ). Thus, for

1 For most scenarios τ 2 U ,r and τ 2 V ,r are not independently identifiable, as the level of signal is determined by their product.

hyperparameters τ and σ the free energy to minimize is

$$F ( q \, | \, \sigma, \tau ) = E _ { q } \log \frac { q _ { u } ( \mathbf U ) q _ { v } ( \mathbf V ) } { p ( \mathbf X \, | \, \mathbf U, \mathbf V, \sigma ) p ( \mathbf U \, | \, \tau _ { \mathbf U } ) p ( \mathbf V \, | \, \tau _ { \mathbf V } ) },$$

and a point estimate for low-rank structure ˆ can be obtained via its expected value S under q ( ): ·

$$\hat { \mathbf S } = E _ { q } ( \mathbf U \mathbf V ^ { T } ).$$

Under an empirical variational Bayes approach, τ and σ are also estimated by minimizing the free energy in (3). For fully-observed X , Nakajima et al. (2012) showed that the global solution to the empirical variation Bayes objective can be obtained in closed form as a singular value thresholding operator; their result for fixed σ is given in Theorem 6, and for estimation of σ in Theorem 7.

Theorem 6. The approximation ˆ S (4) obtained by minimizing (3) over q and τ is ˆ = S S = U D V X S T X where D S is diagonal with

$$\mathbf D _ { \mathbf S } [ r, r ] = 0 \text{ if } D _ { \mathbf X } [ r, r ] < \sigma \sqrt { M + N + \sqrt { M N } ( \kappa + 1 / \kappa ) }$$

and

$$\mathbf D s [ r, r ] = \left ( \mathbf D _ { X } [ r, r ] ^ { 2 } - ( M + N ) \sigma ^ { 2 } + \sqrt { ( \mathbf D _ { X } [ r, r ] ^ { 2 } - ( M + N ) \sigma ^ { 2 } ) ^ { 2 } - 4 M N \sigma ^ { 4 } } \right ) / ( 2 \mathbf D _ { X } [ r, r ] ) \ o t h e r w i s e.$$

The value κ is given as the zero-cross point of the function

$$f ( x ) = x \sqrt { N / M } \, l o g \left ( x \sqrt { M / N } + 1 \right ) + x \sqrt { M / N } \, l o g \left ( x \sqrt { N / M } + 1 \right ) - 1.$$

Theorem 7. Assume without loss of generality that M &gt; N and α = N/M . The estimate ˆ σ that minimizes (3) is the global minimizer of

$$\Psi ( \sigma ) = \sum _ { r = 1 } ^ { N } \psi _ { 1 } \left ( \frac { \mathbf D _ { X } [ r, r ] ^ { 2 } } { M \sigma ^ { 2 } } \right ),$$

where

$$\psi _ { 2 } ( x ) = x - \log x + 1 _ { \{ x > c \} } \psi _ { 2 } ( x ), \\ \psi _ { 2 } ( x ) = \log \left ( \sqrt { \alpha } \psi _ { 3 } ( x ) + 1 \right ) + \alpha \log \left ( \psi _ { 3 } ( x ) / \sqrt { \alpha } + 1 \right ) - \alpha \psi _ { 3 } ( x ),$$

and

$$\psi _ { 3 } ( x ) = \frac { 1 } { 2 \sqrt { \alpha } } \left ( x - 1 - \alpha + \sqrt { ( x - ( 1 + \alpha ) ) ^ { 2 } - 4 \alpha } \right ).$$

Together, Theorems 7 and 6 provide a powerful and efficient way to approximate low-rank signal without any tuning parameters. The only substantial computing

requirement is in obtaining the SVD of X , as other steps involve univariate functions that are trivial to solve.

## 4 Bidimensional linked matrix factorization

BIDIFAC+ (Lock et al., 2022) was developed to decompose bidimensionally linked data as in (1) into modules of low-rank structure that are shared across column sets and/or row sets. That is,

$$X _ { \bullet } \approx \sum _ { k = 1 } ^ { K } \mathbf S _ { \bullet } ^ { ( k ) },$$

where each S ( k ) ·· is a concatenation of blocks S ( k ) ij as in (1) and

$$\mathbf S _ { i j } ^ { ( k ) } = \begin{cases} 0 _ { M _ { i } \times N _ { j } } & \text{if } \mathbf R [ i, k ] = 0 \text{ or } C [ j, k ] = 0 \\ \mathbf U _ { i } ^ { ( k ) } \mathbf V _ { j } ^ { ( k ) T } & \text{if } \mathbf R [ i, k ] = 1 \text{ and } C [ j, k ] = 1 \end{cases}.$$

Here, R and C are binary matrices that indicate the presence of module k ( S ( k ) ·· ) across the row and column sets; that is, S ( k ) ·· is low-rank structure specific to the submatrix defined by the row sets identified in R [ · , k ] and the column sets identified in C [ · , k ]. This general framework subsumes several other integrative decompositions of structured variation. For example, in the unidimensionally linked case with shared columns J = 1, JIVE (Lock et al., 2013) identifies components that are either shared across all matrices ( R [ i, k ] = 1 for all i ) or specific to a single matrix ( R [ i, k ] = 1 for just one i ), while SLIDE (Gaynanova and Li, 2019) allows for components that are shared across any number of matrices. The BIDIFAC method (Park and Lock, 2020) is a special case in which modules are either globally shared ( R [ i, k ] = C [ j, k ] = 1 for all i, j ), shared across all rows for a given column set ( R [ i, k ] = 1 for all i and C [ j, k ] = 1 for one j ), shared across all columns for a given row set ( C [ j, k ] = 1 for all j and R [ i, k ] = 1 for one i ), or specific to one matrix ( R [ i, k ] = 1 and C [ j, k ] = 1 for just one i, j ).

The BIDIFAC and BIDIFAC+ approaches solve a structured nuclear norm objective:

$$\underset { \{ \mathbf S _ { \bullet \, } ^ { ( k ) } \} _ { k = 1 } ^ { K } } { \arg \min } \ \frac { 1 } { 2 } | | \mathbf X _ { \bullet \, } - \sum _ { k = 1 } ^ { K } \mathbf S _ { \bullet \, } ^ { ( k ) } | | _ { F } ^ { 2 } + \sum _ { k = 1 } ^ { K } \lambda _ { k } | | \mathbf S _ { \bullet \, } ^ { ( k ) } | | _ { * },$$

subject to S ( k ) ij = 0 M i × N j if R [ i, k ] = 0 or C [ j, k ] = 0. For fixed R and C , objective (6) may be solved with an iterative soft-singular value thresholding algorithm that cycles through the estimation of each S ( k ) ·· with applications of Proposition 3. For smaller I and J , R and C can enumerate all possible submatrices of the row and column sets; the solution may still be sparse, with S ( k ) ·· = 0 if the submatrix corresponding to module k has no particular shared structure. Alternatively, if the number of row and

column sets I and J are large, Lock et al. (2022) describe an algorithm to adaptively estimate R and C to give the submatrices with the most significant low-rank structure during estimation.

The choice of the module-specific penalty parameters λ k is critical. By default, λ k are selected by an approach motivated by Proposition 5; under the assumption that each matrix X ij has error variance 1, λ k depends on the total number of non-zero rows and columns in its submatrix,

$$\lambda _ { k } = \sqrt { \sum _ { i = 1 } ^ { I } { \mathbf R } [ i, k ] M _ { i } } + \sqrt { \sum _ { j = 1 } ^ { J } { \mathbf C } [ j, k ] N _ { j } },$$

for k = 1 , . . . , K . This approach will shrink to 0 any singular values in the submatrix that are due to error only. In practice each data matrix X ij is scaled to have error variance 1 beforehand; if the error variance is unknown, then it can be taken as the total variance in X ij (a conservative choice) or via other approaches such as the median absolute deviation estimator in Gavish and Donoho (2017).

## 5 Empirical variational BIDIFAC (EV-BIDIFAC)

While the BIDIFAC+ approach has several advantages, one weakness is the required standardization by the noise variance in each data matrix, which is estimated pre-hoc and not within a cohesive framework. More critically, the use of nuclear norm penalties generally results in over-shrinkage of the signal, and the default theoretically-motivated choice of the λ k 's is particularly conservative. Thus, we aim to develop a more flexible empirical variational Bayes approach to the decomposition of bidimensionally linked matrices, analogous to that for a single matrix in Section 3. Our psuedo-Bayesian model is as follows for i = 1 , . . . I and j = 1 , . . . , J :

$$\mathbf X _ { i j } = \sum _ { k = 1 } ^ { K } \mathbf U _ { i } ^ { ( k ) } \mathbf V _ { j } ^ { ( k ) T } + \mathbf E _ { i, j } \text{ where } \mathbf E _ { i, j } [ l, m ] \stackrel { i i d } { \sim } \text{Normal} ( 0, \sigma _ { i j } ^ { 2 } ),$$

(7)

$$\mathbf U _ { i } ^ { ( k ) } [ \cdot, r ] & = 0 \text{ if } \mathbf R [ i, k ] = 0 \text{ and } \mathbf U _ { i } ^ { ( k ) } [ \cdot, r ] \sim \text{Normal} ( 0, \sigma _ { i j } \tau _ { u } ^ { 2 } [ k, r ] \mathbf I ) \text{ if } \mathbf R [ i, k ] = 1, \\ \mathbf V _ { j } ^ { ( k ) } [ \cdot, r ] & = 0 \text{ if } \mathbf C [ j, k ] = 0 \text{ and } \mathbf V _ { j } ^ { ( k ) } [ \cdot, r ] \sim \text{Normal} ( 0, \sigma _ { i j } \tau _ { v } ^ { 2 } [ k, r ] \mathbf I ) \text{ if } \mathbf C [ j, k ] = 1,$$

where I is an identity matrix of appropriate dimension. Using repeated applications of Corollary 4, the BIDIFAC+ objective can be shown to give the posterior mode of the model specified in (7) for which σ = { σ ij : i = 1 , . . . , I, j = 1 , . . . , J } and τ = { τ u , τ v } are fixed, and τ 2 u [ k, r ] = τ 2 v [ k, r ] = 1 /λ k for all r . We extend this framework by allowing σ and τ to be estimated within a cohesive empirical variational Bayes framework. We approximate the true posterior with a distribution q ( ), for which ·

the factor matrices are mutually independent:

$$q ( \mathbf U, \mathbf V ) = \prod _ { k = 1 } ^ { K } q _ { u, k } ( \mathbf U ^ { ( k ) } ) q _ { v, k } ( \mathbf V ^ { ( k ) } ).$$

The corresponding free energy is given by

$$F ( q \, | \, \sigma, \tau ) = E _ { q } \log \frac { \prod _ { k = 1 } ^ { K } q _ { u, k } ( \mathbf U ^ { ( k ) } ) q _ { v, k } ( \mathbf V ^ { ( k ) } ) } { p ( \mathbf X _ { \bullet } \, | \, \mathbf U, \mathbf V, \sigma ) p ( \{ \mathbf U ^ { ( k ) } \} _ { k = 1 } ^ { K } \, | \, \tau _ { \mathbf U } ) p ( \{ \mathbf V ^ { ( k ) } \} _ { k = 1 } ^ { K } \, | \, \tau _ { \mathbf V } ) }, \quad ( 8 )$$

with the densities p ( ) · in the denominator as specified in model (7). The resulting low-rank fit and decomposition is then given by

$$\hat { \mathbf S } _ { \centerdot } = E _ { q } \left ( \sum _ { k = 1 } ^ { K } \mathbf U _ { \cdot } ^ { ( k ) } \mathbf V _ { \cdot } ^ { ( k ) T } \right ) = \sum _ { k = 1 } ^ { K } \hat { \mathbf S } _ { \centerdot } ^ { ( k ) }$$

where for k = 1 , . . . , K ,

$$\hat { \mathbf S } _ { \bullet } ^ { ( k ) } = E _ { q _ { u }, q _ { v } } \left ( \mathbf U _ { \cdot } ^ { ( k ) } \mathbf V _ { \cdot } ^ { ( k ) T } \right ).$$

Algorithm 1 describes an iterative approach to compute the empirical variational Bayes solution. This algorithm converges to a blockwise minimum of the free energy (8) over q and τ with σ fixed at the distinct minimizer of (8) for each X ij separately. Alternatively, one can update the variances σ over the iterative algorithm to minimize free energy under the joint model. In Algorithm 1, the number of modules K and the row- and column-set inclusions for the modules, R and C , are assumed to be fixed. For a smaller number of linked matrices, one can set R and C to enumerate all of the possible submatrices of shared, partially shared, or specific structure: K=(2 I -1)(2 J -1). The estimated decomposition may still be sparse in that ˆ S ( k ) ·· = 0 if module k has no distinct low-rank structure, as demonstrated in Section 8.3. If enumerating all submatrices of the row- and column-sets is not feasible due to a large number matrices, one can specify K , R and C based on prior knowledge. Alternatively, we may approximate the fully enumerated decomposition with a large upper bound on the number distinct modules, K , and the row- and column-set inclusions R and C estimated empirically as in Lock et al. (2022).

## 6 Uniqueness

The uniqueness of decompositions into shared, partially shared, and specific low-rank structure requires careful consideration, as the terms in the model are generally not uniquely defined in general without further constraints. For unidimensional decompositions of low-rank structure, conditions of orthogonality are commonly used to assure the different components are uniquely identified (Lock et al., 2013; Feng et al., 2018; Gaynanova and Li, 2019; Yi et al., 2023); however, such conditions are not straightforward to extend to decompositions for bidimensionally linked data. Alternatively, Lock et al. (2022) showed that the terms in the BIDIFAC+ solution are uniquely defined

## Algorithm 1 Estimation of EV-BIDIFAC model for fully observed data

- 1. Estimate σ ij for each X ij , as in Theorem 7.
- 2. Define ˜ X ·· as in (1), with ˜ X ij = X ij /σ ij for i = 1 , . . . , I and j = 1 , . . . , J .
- 3. Initialize ˜ S ( k ) ·· = 0 for k = 1 , . . . , K
- 4. For k = 1 , . . . , K :

̸

- (a) Let ˜ X ( k ) ·· be the submatrix of ˜ X ·· -∑ k ′ = k ˜ S ( k ′ ) ·· over row sets indicated by R [ · , k ] and column sets C [ · , k ].
- (b) Update the non-zero submatrix in ˜ S ( k ) ·· by the application of the singular value thresholding in Theorem 6 to ˜ X ( k ) ·· .
- 5. Repeat Step 3. until convergence of the { ˜ S ( k ) ·· } K k =1 .
- 6. Unscale the estimated structure: ˆ S ( k ) ij = σ ij ˜ S ( k ) ij for k = 1 , . . . , K , j = 1 , . . . , J , i = 1 , . . . , I .

under certain conditions, without orthogonality, as the minimizer of the structured nuclear norm objective in (6). Here we extend this result to a much broader class of penalties, that includes the EV-BIDIFAC decomposition.

For fixed row-set and column-set indicators for the modules, R and C , let S ˆ X be the set of possible decompositions that yield a given approximation ˆ X ·· :

$$\mathbb { S } _ { \hat { X } } = \left \{ \{ \mathbf S _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } \, | \, \hat { \mathbf X } _ { \bullet } = \sum _ { k = 1 } ^ { K } \mathbf S _ { \bullet } ^ { ( k ) } \right \}.$$

In most setting the cardinality of | S ˆ X | is infinite without further conditions. Define P λ ( ) for a vector · λ as a weighted sum of the singular values of a matrix:

$$P _ { \lambda } ( \mathbf S ) = \sum \lambda [ r ] { \mathbf D } _ { \mathbf S } [ r, r ],$$

and let f pen ( ) give a structured penalty on the different modules of a decomposition,: ·

$$f _ { \text{pen} } ( \{ \mathbf S _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } ) = \sum _ { k = 1 } ^ { K } P _ { \lambda _ { k } } ( \mathbf S _ { \bullet } ^ { ( k ) } ),$$

with non-zero entries for each λ k . Theorem 8 gives sufficient conditions for the minimizer of this penalty to be uniquely determined.

Theorem 8. Consider { ˆ S ( k ) ·· } K k =1 ∈ S ˆ X and let U ( k ) · ˆ D ( k ) V ( k T ) · give the SVD of ˆ S ( k ) ·· for k = 1 , . . . , K . The following three properties uniquely identify { ˆ S ( k ) ·· } K k =1 .

- 1. { ˆ S ( k ) ·· } K k =1 minimizes f pen ( ) · over S ˆ X ,
- 2. { ˆ U ( k ) i [ · , r ] : R [ i, k ] = 1 and ˆ D ( k ) [ r, r ] &gt; 0 } are linearly independent for i = 1 , . . . I ,
- 3. { ˆ V ( k ) j [ · , r ] : C [ j, k ] = 1 and ˆ D ( k ) [ r, r ] &gt; 0 } are linearly independent for j = 1 , . . . , J .

The proof of Theorem 8 is given in Appendix A.1. Note that the linear independence conditions are straightforward to verify, and are likely to hold in practice if the ranks of the terms in { ˆ S ( k ) ·· } K k =1 are small relative to the matrix dimensions. As noted in Proposition 9, the terms in any accumulation point of Algorithm 1 also give a coordinate-wise optimum of f pen ( ) for a certain set of penalties · { λ k } K k =1 , and thus the EV-BIDIFAC solution is unique in this sense.

Proposition 9. If { ˆ S ( k ) ·· } K k =1 ∈ S ˆ X is an accumulation point of Algorithm 1, then it is a coordinate-wise optimum of f pen ( ) · over S ˆ X for some { λ k } K k =1 .

## 7 Missing data imputation

Now we consider the context in which the data are not fully observed and missing data must be imputed. Let M index observations in the full dataset X ·· that are not observed: M = { ( m,n ) : X ·· [ m,n ] is missing } . We develop an ExpectationMaximization (EM) approach for imputation that is analogous to previous approaches described for hard-thresholding (Kurucz et al., 2007), soft-thresholding (Mazumder et al., 2010), and BIDIFAC+ (Lock et al., 2022). The procedure for fixed estimates of the noise variance for each matrix, σ ij , is described in Algorithm 2.

- Algorithm 2 Estimation of EV-BIDIFAC model for missing data 1. Initialize ˜ S ( k ) ·· = 0 for k = 1 , . . . , K 2. Define ˜ X ·· as in (1), with ˜ X ij = X ij /σ ij for i = 1 , . . . , I and j = 1 , . . . , J . 3. Set ˜ X ·· [ m,n ] = ˜ S ·· [ m,n ] for ( m,n ) ∈ M , where ˜ S ·· = ∑ K k =1 ˜ S ( k ) ·· . 4. Update { ˜ S ( k ) ·· } K k =1 as in Step 4. of Algorithm 1. 5. Repeat Steps 3. and 4. until convergence of the { ˜ S ( k ) ·· } K k =1 . 6. Unscale the estimated structure: ˆ S ( k ) ij = σ ij ˜ S ( k ) ij for k = 1 , . . . , K , j = 1 , . . . , J i = 1 , . . . , I . 7. Impute X ·· [ m,n ] = ˆ S ·· [ m,n ] for ( m,n ) ∈ M .
- ,

Note that Step 4. in Algorithm 2 is a partial maximization step, as the free energy in (8) is conditionally minimized for each module given the current imputed values, while Step 5. is an expectation step replacing the missing data with their expected value (9) under the current model fit. Further, the algorithm can be shown to minimize the free energy over the missing data { X ·· [ m,n ] : ( m,n ) ∈ M} ; this result is given as Theorem 10, and is appealing because all unknowns (model hyperpameters, parameters and missing entries) are then estimated via optimizing the same unified free energy objective.

Theorem 10. Algorithm 2 iteratively minimizes the free energy (9) over q ( ) · , τ , and { X ·· [ m,n ] : ( m,n ) ∈ M} .

A proof of Theorem 10 is provided in Appendix A.2. For the choice of σ , we suggest two different approaches, depending on the type of missingness. If entries are missing from an entire row or columns of a given data matrix X ij (i.e., row- or columnwise missing), then σ ij can be estimated by the applying Theorem 2 to the complete submatrix { X ·· [ m,n ] : ( m,n ) / ∈ M} . If the matrix is missing entries without entire rows or columns, then we update σ ij iteratively during Algorithm 2 as follows at the end of each cycle: 1. apply Theorem 2 to the matrix X ij with currently imputed values to obtain σ ′ ij , 2. update σ ij = σ ′ ij ( M N i j M N i j -|M ij | ) where |M | ij is the number of missing entries in X ij .

## 8 Simulations

Here we present a collection of simulations to illustrate and compare the proposed approach with respect to various criteria. We explore several conditions covering different signal-to-noise ratios, missing data scenarios, and linked matrix patterns with shared and specific structure. Further simulation results (e.g., exploring different matrix dimensions) are available at https://ericfrazerlock.com/ev-bidifac-sims.html.

## 8.1 Single matrix

We first describe a simple simulation on a single matrix to illustrate the empirical variational Bayes approach to matrix factorization as a flexible compromise between a hard-thresholding and a conservative soft-thresholding approach. We consider a fully observed matrix X : 1000 × 100, with X = c UV T + E where the entries of E U , : M × 10 and V : N × 10 are each independent standard normal, and c is a manipulated value controlling the signal-to-noise ratio. We consider 10 values of c distributed on a log-uniform scale between c = 0 05 and . c = 1, and generate 10 replications of the simulation for each value of c . For each generated dataset, we apply the following approaches:

- · OPT: The oracle operator on the singular values when the structure S = c UV T is known: ˆ S OPT = U DV X ˆ T X where ˆ D is diagonal with entries that solve the OLS problem to minimize || S -U DV X ˆ T X || 2 F .
- · HT: The hard-thresholding approach described in Proposition 1, using the true rank R = 10.
- · EVB: the empirical variational Bayes approach applying Theorems 6 and 7.
- · NN: The nuclear norm penalized soft-thresholding approach described in Proposition 2, with λ = √ M + √ N (this correspond to a threshold motivated by Proposition 5 with true noise variance σ = 1).

For each setting and for each method we compute the relative squared error (RSE) and oracle-normalized squared error (ONSE)

$$R S E = \frac { | | \mathbf S - \hat { \mathbf S } | | _ { F } ^ { 2 } } { | | \mathbf S | | _ { F } ^ { 2 } } \text{ and O N S E } = \frac { | | \mathbf S - \hat { \mathbf S } | | _ { F } ^ { 2 } } { | | \mathbf S - \hat { \mathbf S } _ { O P T } | | _ { F } ^ { 2 } }. \quad \quad ( 1 0 )$$

Figure 1 shows the mean values for both RSE and ONSE under the different methods as a function of the signal-to-noise c . The hard-thresholding approach tends to overfit the signal, which can lead to particularly poor recovery when the signal is smaller relative to the noise. In contrast, the nuclear norm penalty is appropriately conservative in that it never overfits and performs better for lower signal; however, it tends to overshrink for higher signal and can have several fold higher error than hard-thresholding. The EVB approach is a flexible compromise between the two methods, performing comparably or better than both methods for all scenarios and substantially outperforming them for moderate levels of signal. Moreover, performance of the EVB method is impressively close to the optimal in most scenarios.

Fig. 1 Error in estimating the underlying low-rank signal for a single matrix under different methods, and under different signal-to-noise (s2n) ratios. The left-panel gives relative squared error, and the right-panel gives oracle normalized standard error. All axes are on a log-scale.

We consider another scenario with low-rank structure in which the underlying components have heterogeneous levels of signal. First, we simulate S ′ = UV ′ ′ T where the entries of U : M × 10 and V : 10 × N are independent standard normal. Then, we generate the signal as

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-eb2572fe-f9fe-437e-aee4-4af3620a960c.jpg)

$$\mathbf S = \mathbf U _ { S ^ { \prime } } \mathbf D _ { S } \mathbf V _ { S ^ { \prime } } ^ { T }$$

where D S is diagonal with entries sampled randomly from a log-uniform distribution between 0 05 . √ MN and √ MN , and X = S + E with entries of E independent standard normal. Thus, some singular values of the true low-rank structure S may correspond to high signal, others may be of moderate signal, and others may be of weak or undetectable signal. We generate 100 replications in this case, and for each we apply the EVB, HT, and NN methods described above, as well as the following additional approaches:

- · HT-OPT: the oracle hard-thresholding approach, in which the rank is selected that results in minimal squared error between the estimated signal and the true signal.
- · RMT: the low-rank reconstruction approach proposed in Shabalin and Nobel (2013), based on minimizing squared error loss using asymptotic random matrix theory, with the true error variance σ 2 = 1.
- · NN-OPT: the oracle soft-thresholding approach, in which the nuclear norm penalty is selected that results in minimal squared error between the estimated signal and the true signal.

Boxplots showing distribution of RSE and ONSE for recovering the underlying structure across replications are shown in Figure 2. The EVB and RMT methods have performance close to the oracle, and vastly outperform hard-thresholding and nuclear norm penalized approaches. The relatively poor performance of HT-OPT suggests that singular values benefit from some shrinkage, while the poor performance of NN-OPT indicates that the level of shrinkage should not be uniform across singular values.

Fig. 2 Error in estimating underlying low-rank structure in which the rank-1 components have heterogenous signal sizes. The left-panel gives relative squared error (RSE), and the right-panel gives oracle normalized standard error (ONSE). All axes are on a log-scale.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-293174f7-0692-41ee-a0c0-54b33819297e.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-bae4a68b-4cc4-4089-b1cc-39cd0407e3df.jpg)

We also compare missing data imputation performance of the EVB approach using Algorithm 2. For each of the 100 simulated datasets generated, we consider scenarios in which 20%, 50%, or 80% of the entries are randomly held out as missing and imputed. In addition to EVB imputation, we apply similar EM-imputation methods for the soft-thresholding (corresponding to softImpute of Mazumder et al. (2010)) and hard thresholding (Kurucz et al., 2007) approaches. The RMT approach is not considered here, as it does not have an associated missing data imputation procedure. For the NN-OPT and HT-OPT approaches, we select the value of λ or R , respectively, that

Fig. 3 Missing data imputation accuracy for different levels of missingness. The left column gives RSE miss and the right gives ONSE miss .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-5c53615d-b9c2-425c-83c7-5a27f4578350.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-b51d2637-f176-4a00-9713-a087715603b5.jpg)

RSE 502 missing

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-88c86ab6-059e-478c-8389-df8cf931cd16.jpg)

RSE 8016 missing

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-d83da0a7-4d84-4672-b3e7-e3626150fcbd.jpg)

Normalized RSE 501 missing

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-39ef5770-9fcc-4933-90b7-133d2c41d977.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-98b3f43a-52f0-46d9-a140-d45696c03a95.jpg)

15

minimizes squared imputation error for the signal of the missing data

$$R S E _ { \text{miss} } = \frac { \sum _ { ( m, n ) \in \mathcal { M } } ( \mathbf S [ m, n ] - \hat { \mathbf S } [ m, n ] ) ^ { 2 } } { \sum _ { ( m, n ) \in \mathcal { M } } \mathbf S [ m, n ] ^ { 2 } }.$$

For oracle normalized imputation error ONSE miss , we divide RSE miss by the RSE for the OPT method in (10). The results are shown in Figure 3, and demonstrate that the EVB imputation approach is stable for different levels of missingness and substantially outperforms the other approaches.

## 8.2 Two linked matrices

Here we generate two matrices, X 1 : 500 × 100 and X 2 : 500 × 100, that are linked by columns: X · = [ X 1 ; X 2 ] : 1000 × 100. We consider this context primarily to allow comparison of the EV-BIDIFAC approach with several other methods that identify shared and unshared low-rank structure in this setting.

First, we perform an illustrative simulation akin to that in Figure 1 for a single matrix. We simulate shared structure S 0 = c U V 0 T 0 and individual structures S 1 = c U V 1 T 1 S 2 = c U V 2 T 2 , where U 0 : 1000 × 5, V 0 : 100 × 5, U 1 : 500 × 5, U 2 : 500 × 5, V 1 : 100 × 5, and V 2 : 100 × 5 each have independent entries from a standard normal distribution. The observed data are then generated as X · = S 0 +[ S 1 ; S 2 ] + E where entries of E are independent standard normal. Thus c again controls the signal to noise ratio and we consider 10 values of c distributed on a log-uniform scale between c = 0 05 and . c = 1, and generate 10 replications for each value of c . Here we consider estimates for the following methods:

- · EVB: the EV-BIDIFAC approach described in Algorithm 2
- · UNIFAC: the BIDIFAC method (Park and Lock, 2020) for vertical integration using the true noise variance σ 2 = 1.
- · JIVE: the JIVE method implemented in O'Connell and Lock (2016) for true shared and unshared ranks R = 5 , R 1 = 5 , R 2 = 5.

Here JIVE is analogous to HT and UNIFAC to NN for a single matrix, as they involve iterative applications HT or NN, respectively, to the shared and unshared components. For each method we compute RSE for the overall structure ˆ S = ˆ S 0 + [ ˆ S 1 ; ˆ S 2 ] as in (10), and also compute the relative decomposition squared error (RDSE) as

$$R D S E = \frac { | | \mathbf S _ { 0 } - \hat { \mathbf S } _ { 0 } | | _ { F } ^ { 2 } + | | \mathbf S _ { 1 } - \hat { \mathbf S } _ { 1 } | | _ { F } ^ { 2 } + | | \mathbf S _ { 2 } - \hat { \mathbf S } _ { 2 } | | _ { F } ^ { 2 } } { | | \mathbf S _ { 0 } | | _ { F } ^ { 2 } + | | \mathbf S _ { 1 } | | _ { F } ^ { 2 } + | | \mathbf S _ { 2 } | | _ { F } ^ { 2 } }.$$

The results are shown in Figure 4. The results for overall RSE are similar to that for a single matrix, with UNIFAC appropriately conservative for lower signal but substantially underperforming JIVE for higher signal due to over-shrinkage. However, per the RDSE results, UNIFAC can still identify the shared and unshared components more accurately with higher signal - this phenomenon was noted previously (Lock et al., 2022) and is in part due to orthogonality restrictions imposed by JIVE and related methods that are required for uniqueness of the decomposition without shrinkage. The

EVB approach again acts as a flexible compromise between the two methods, and substantially outperforms both for many s2n scenarios, particularly with respect to RDSE.

Fig. 4 Error in estimating the underlying low-rank signal for two linked matrices under different signal-to-noise (s2n) ratios. The left-panel gives the RSE, and the right-panel gives RDSE. All axes are on a log-scale.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-a1c2f0b0-21f3-4385-8d91-ab3295768c6d.jpg)

Now, we compare RSE and RDSE under a scenario with heterogenous signal levels for the factors in each of the shared and unshared structures. The rank-5 signal matrices S 0 , S 1 , and S 2 are each generated as in (11) with singular values again generated uniformly from a log-uniform distribution between 0 05 . √ MN ∗ 0 05 and . √ MN for shared structure and between 0 05 . √ M N i and √ M N i for specific structure, and 100 replications are generated for this scenario. In addition to the EVB, JIVE, UNIFAC methods we also consider the following approaches:

- · SLIDE: the structural learning and integrative decomposition (SLIDE) method (Gaynanova and Li, 2019), with the true number of shared and unshared components (i.e., the ranks) specified.
- · HNN.OPT: the hierarchical nuclear norm (HNN) approach (Yi et al., 2023), with tuning parameters selected on a grid to minimize actual RSE for the true signal.

Figure 5 gives distribution of RSE and RDSE over the replications for the different methods. The EVB approach universally outperforms other methods in terms of both RSE and RDSE, and the benefit is large (by a factor of 2 or more) for both criteria in several cases. While the HNN method is motivated by a model with shared and unshared low-rank structure, it does not provide an explicit decomposition of the structure and thus it is not shown for RDSE.

Fig. 5 RSE (left) and RDSE (right) for low-rank structure with heterogenous signal levels for two linked matrices.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-442ec61a-1d67-4798-bb14-ebafe7c2c7f8.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-3dc86c86-6de6-43ee-a4d9-83ac678dd3f7.jpg)

## 8.3 Bidimensionally linked matrices

Now, we consider the full bidimensionally linked framework in (1), with I = J = 2, N 1 = N 2 = 50 and D 1 = D 2 = 500. Data were generated as

$$\mathbf X _ { \bullet \bullet } = \sum _ { k = 1 } ^ { K } \mathbf S _ { \bullet \bullet } ^ { ( k ) } + \mathbb { E },$$

where the entries of E are independent standard normal. There are K = 9 possible distinct modules of low-rank structure in this setting: globally shared, shared across rows for either of the of the J = 2 column sets, shared across columns for either of the I = 2 row sets, or structure specific to any of the 4 matrices X ij . We generate 100 datasets, where for each replication 5 of the 9 possible modules are randomly selected to have non-zero low-rank structure and the others are set to zero, S ( k ) ·· = 0. All nonzero modules are generated to have rank 2 structure on the given submatrix as in (11). For each simulated dataset, we apply the following methods:

- · EB-BIDI: the approach in Algorithm 1, with R and C enumerating all K = 9 possible modules.
- · EB-SEP: the EVB approach applied to each matrix X ij separately
- · BIDIFAC: the BIDIFAC algorithm described in (Lock et al., 2022) with R and C enumerating all possible modules and using the true noise variance σ 2 = 1.
- · EB-JOINT: the EVB approach applied to the overall matrix X ·· .

For each method, we compute the overall RSE for S ·· = ∑ K k =1 S ( k ) ·· as in (10), and for EB-BIDI and BIDIFAC we compute the RDSE as

$$R D S E = \frac { \sum _ { k = 1 } ^ { 9 } | | \mathbf S _ { \bullet } ^ { ( k ) } - \hat { \mathbf S } _ { \bullet } ^ { ( k ) } | | _ { F } ^ { 2 } } { \sum _ { k = 1 } ^ { 9 } | | \mathbf S _ { \bullet } ^ { ( k ) } | | _ { F } ^ { 2 } }.$$

̸

The results are summarized in Figure 6. Overall RSE is substantially better for EBBIDI vs. other methods; in comparison to EB-SEP and EB-JOINT, this illustrates the advantages of a parsimonious decomposition of the low-rank signal into the different shared and unshared components with respect recovering the overall signals. Moreover, EB-BIDI has substantially better RDSE compared to BIDIFAC, demonstrating the advantages of the empirical variational Bayes objective with respect to accurately decomposing the structure into its different components. Moreover, EB-BIDI correctly estimated modules with no true structure ( S ( k ) ·· = 0 ) as ˆ S ( k ) ·· = 0 for all of the 400 total instances (100 total datasets with 4 possible modules set to 0 in each case); on the other hand, it estimated ˆ S ( k ) ·· = 0 in only 19 / 500 = 3 8% of instances when . S ( k ) ·· = 0, which occured when S ( k ) ·· was difficult to detect due to a small signal size. This demonstrates how the approach can correctly recover the underlying rank-sparsity structure, only estimating non-zero structure for a module if it truly has distinct low-rank structure on its given rows and columns.

Fig. 6 RSE (left) and RDSE (right) for the scenario with 2 × 2 bidimensionally linked matrices.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-9509383e-fdea-41df-a836-2122ce0bf832.jpg)

We extend this simulation design to consider missing data imputation accuracy for the generated datasets. We consider two different missingness scenarios: an entrywise missing scenario in which 20% of the entries are randomly set to missing in each matrix, and a blockwise missing scenario in which 10% of the columns and 10% of the rows are randomly selected to have all of their entries missing in each matrix. In addition to imputation using the previous methods, we also consider the NN.OPT and

HT.OPT methods as described in Section 8.1. For entrywise imputation the RSE miss is computed as in (12). For blockwise imputation, the true underlying structure is taken to be the sum of column-shared structure for missing rows and the sum of row-shared structure for missing columns, as this is the only estimable signal when an entire row or column is missing. That is, letting S C ·· = ∑ k : C [ · ,k ]= 1 S ( k ) ·· and S R ·· = ∑ k : R [ · ,k ]= 1 S ( k ) ·· ,

$$R S E _ { m i s s } = \frac { \sum _ { ( m, n ) \in \mathcal { M } ^ { R } } ( \mathbf S _ { \bullet } ^ { C } [ m, n ] - \hat { \mathbf S } _ { \bullet } [ m, n ] ) ^ { 2 } + \sum _ { ( m, n ) \in \mathcal { M } ^ { C } } ( \mathbf S _ { \bullet } ^ { R } [ m, n ] - \hat { \mathbf S } _ { \bullet } [ m, n ] ) ^ { 2 } } { \sum _ { ( m, n ) \in \mathcal { M } ^ { R } } \mathbf S _ { \bullet } ^ { C } [ m, n ] ^ { 2 } + \sum _ { ( m, n ) \in \mathcal { M } ^ { C } } \mathbf S _ { \bullet } ^ { R } [ m, n ] ^ { 2 } },$$

where M R indexes row-missing entries and M C indexes column missing entries. The results are shown in Figure 7. The EB-BIDI method substantially and consistently outperforms the other approaches with respect to missing data imputation, and the advantages are particularly large for blockwise missing imputation. Note that the EBSEP approach does not estimate any structure for blockwise missing data, as it is not estimable when considering each matrix separately.

Fig. 7 RSE for entrywise missing data imputation (left) and blockwise missing data imputation (right), for the scenario with 2 × 2 bidimensionally linked matrices.

## Entry Missing imputation

## Block Missing imputation

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-635b7ed4-f364-43a6-90ce-4bea54429e62.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-4c9ac488-1a58-47ed-b085-5435ea3fa243.jpg)

## 9 Application to BRCA Data

Here we describe an application to gene expression (mRNA) and microRNA (miRNA) data for breast cancer (BRCA) from the Cancer Genome Atlas (TCGA) (TCGA Research Network et al., 2012). The data were processed as described in Park and Lock (2020), resulting in D 1 = 500 mRNAs and D 2 = 500 miRNAs for N 1 = 660 cancer tumor samples and N 2 = 86 non-cancerous breast tissue samples from unrelated individuals. Thus, the data can be represented as bidimensionally linked matrices as

in (1), where X 11 is mRNA data for cancer samples, X 12 is mRNA data for normal samples, X 21 is miRNA data for cancer samples and X 22 is miRNA data for normal samples. The decomposition of data into different modules of shared structured variation is very well-motivated in this context, as low-rank structure on any of the K = 9 possible modules is plausible and interpretable. For example, it is reasonable to expect shared structure across mRNA and miRNAs, as miRNAs can regulate mRNA. It is also reasonable to expect some shared structure across tumor and normal samples, as breast tumors arise from normal tissue and so they likely share some patterns of molecular variation. However, it is also reasonable to expect specific structure in each case, e.g., patterns of molecular variability that are 'cancer-specific' (i.e., not present in normal tissue) are of particular interest. The BIDIFAC approach decomposes such shared and unique structure across both row sets (mRNA/miRNA) and column sets (tumor/normal), while other methods do not. Moreover, as demonstrated in our simulations, the empirical variational (EV) Bayes approach captures underlying low-rank structure more efficiently and accurately than other low-rank approximation methods. Thus, we expect our EV-BIDIFAC approach to decompose the underlying shared and unshared signals in an informative way, and accurately recover the underlying signal in the full dataset X ·· better than alternative methods.

Fig. 8 Heatmaps of the BRCA data (left) and the full low-rank structure estimated by EV-BIDIFAC. Higher values are colored red, lower values are colored blue, and missing columns are colored black.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-2908bc78-cd1c-4ca8-b2c6-4214412f5275.jpg)

̸

We applied EV-BIDIFAC wherein R and C enumerate all K = 9 possible modules, and significant distinct low-rank structure was identified for each module ( ˆ S ( k ) ·· = 0 for all k ). Figure 8 shows heatmaps of the original data (including missing data for miRNA and mRNA) and the resulting overall fit ˆ S ·· = ∑ K k =1 ˆ S ( k ) ·· . While all estimated structures are low-rank, it captures the original data well and imputation of the missing columns appears reasonable. Considering the relative amount of variation present in the different modules reveals the extent of shared variation across mRNA and miRNA,

and across tumor and normal cells. Interestingly, only 1.47% of the low-rank signal is in the module shared globally across all submatrices, while 59.8% of the signal is shared across mRNA and miRNA but specific to tumors. This suggests that the substantial coordinated activity between mRNA and miRNA in breast tumors is cancer-specific and not present in the normal cells they derive from.

We further compare imputation accuracy for data under entry-wise, column-wise or row-wise missingness under a cross-validation scheme. The accuracy of entry-wise imputation is of interest because it is a good proxy for the accuracy of the the overall fit of the estimated structure. Column-wise imputation accuracy is of interest for two reasons: (1) the actual data have column-wise missingness, as not all samples have both mRNA and miRNA data in the full TCGA cohort, and (2) it provides a robust measure of the level of association between the two data sources, addressing how well can we predict gene expression with miRNA data, and vice versa. Row-wise imputation accuracy is also of interest, as it provides a measure of similarity in molecular structure between cancer and normal samples (e.g., how well can we predict levels of a missing gene in cancer samples, given its associations in normal samples).

We create 20 folds in the data, where for each fold 5% of the columns, 5% of the rows, and 5% of the remaining entries in each matrix are randomly withheld as missing. For each fold, the missing data are imputed using the EV-BIDIFAC, BIDIFAC, EB-SEP or EB-JOINT approaches described in Section 8.3. We further consider nuclear norm penalized imputation (NN) and hard-thresholding (HT) using softImpute (Hastie and Mazumder, 2021), jointly for the full matrix (NN-JOINT and HT-JOINT) or separately for each submatrix (NN-SEP and HT-SEP). In each case the nuclear norm penalty is selected via random matrix theory, and the rank for hard-thresholding is determined via cross-validation by optimizing MSE for an additional 5% of randomly withheld entries. As an alternative to low-rank methods we consider k-nearest neighbor (KNN) imputation using the package impute (Hastie et al., 2021) on the full matrix. Lastly, we consider a structured matrix completion approach (Cai et al., 2016) designed to impute a contiguous missing submatrix under a low-rank approximation using the StructureMC package (Liu and Zhang, 2019). As StructureMC imputes only one submatrix, it does not apply to randomly missing entries, and we hold out one set of rows or columns at a time to impute from the full matrix, with all other entries observed. In each case, the relative squared error for the true observed data is computed and averaged across the four datasets yielding

$$\text{MRSE} _ { \text{miss} } = \frac { 1 } { 4 } \sum _ { i = 1 } ^ { 2 } \sum _ { j = 1 } ^ { 2 } \frac { \sum _ { ( m, n ) \in \mathcal { M } _ { i j } } ( \mathbf X _ { i j } [ m, n ] - \hat { \mathbf S } _ { i j } [ m, n ] ) ^ { 2 } } { \sum _ { ( m, n ) \in \mathcal { M } _ { i j } } \mathbf X _ { i j } [ m, n ] ^ { 2 } }.$$

The MRSE miss under different imputation scenarios are summarized in Table 1. The EV-BIDIFAC approach performs relatively well across all scenarios. For entrywise missing data EV-BIDIFAC performs much better than BIDIFAC (likely due to over-shrinkage of the signal for BIDIFAC) and similarly to EB-SEP, indicating that structure can be well-estimated by considering each dataset separately; however, separate estimation does not allow for row- or column imputation and does not estimate

the shared variation that is of interest. The MRSE miss for row-missing data is relatively modest (0 90) compared to entry-wise missing (0 38), indicating that there is . . little shared information between tumor and normal cells. The joint approach EBJOINT on X ·· performs relatively poorly across all scenarios as it tends to overfit shared variation when there substantial specific variation across platforms (mRNA and miRNA) and cohorts (tumor and normal); this illustrates the benefits of a comprehensive decomposition of shared and specific variation across the constituent matrices.

.

Table 1 Imputation accuracy (MRSE miss ) for BRCA and normal tissue mRNA and miRNA data. The lowest values in each column are shown in bold font; if multiple values are bolded, then their means were not significantly different at a 0.05 level using a paired t-test over the 20 folds

| Method Name   | Entry-missing   |   Column-missing |   Row-missing | Overall   |
|---------------|-----------------|------------------|---------------|-----------|
| EV-BIDIFAC    | 0.389           |            0.756 |         0.903 | 0.682     |
| BIDIFAC       | 0.526           |            0.87  |         0.909 | 0.768     |
| EB-SEP        | 0.381           |            1     |         1     | 0.79      |
| EB-JOINT      | 0.510           |            1.01  |         1.67  | 1.06      |
| NN-SEP        | 0.536           |            1     |         1     | 0.845     |
| NN-JOINT      | 0.614           |            0.891 |         0.881 | 0.796     |
| HT-SEP        | 0.459           |            1     |         1.03  | 0.831     |
| HT-JOINT      | 0.520           |            4.5   |        13.1   | 6.05      |
| KNN           | 0.646           |            1.26  |         1.01  | 0.974     |
| StructureMC   | -               |            0.977 |         1.01  | -         |

## 10 Conclusion and discussion

The proposed EV-BIDIFAC approach can be considered as a flexible compromise between hard- and soft- thresholding techniques for the linked matrix decomposition problem, and has impressive performance with respect to 1. estimating underlying low-rank structure, 2. decomposing underlying structure into its shared, partiallyshared and specific components, and 3. missing data imputation under a variety of scenarios. An appealing aspect of the method is that it derives from a unified and intuitive model with a corresponding objective function that does not require any prehoc choices or parameters to tune. In simulation comparisons the proposed approach tended to outperform other methods even when their parameters are set to the true values (e.g., true rank(s) or true error variance) or ideal values (e.g., to minimize error with knowledge of the true structure).

While we have focused on obtaining a point estimate for underlying structure that minimizes free energy, this could also be used to empirically estimate hyperparameters and initialize a sampling algorithm for the associated fully Bayesian model. This would capture uncertainty in the underlying factorization, which can be used for multiple imputation of missing values or to propagate uncertainty if the lower-dimensional components are used in subsequent statistical modeling. A related Bayesian model

is described in Samorodnitsky et al. (2024), which also illustrates the importance of uncertainty propagation.

Often, data take the form of higher-order tensors instead of matrices, and lowrank tensor factorizations are useful to uncover their underlying structure (Kolda and Bader, 2009). The models considered here may be extended to the tensor setting, e.g., by allowing factors for each dimension of a tensor to have a distribution analogous to that for U and V in the matrix setting. However, the algorithms are not straightforward to extend, as there is likely not a closed-form solution for a tensor akin to that in Theorem 6 for a matrix. Extending this approach to the setting of a single tensor or linked tensors is an interesting direction for future work.

## Supplementary materials

Reproducible R scripts for all results presented in this manuscript are available in a supplementary zipped folder. Software to perform the EV-BIDIFAC method is available at https://github.com/lockEF/bidifac. Additional simulation study results are available at https://ericfrazerlock.com/ev-bidifac-sims.html.

## Acknowledgement

This work was supported by the NIH National Institute of General Medical Sciences (NIGMS) grant R01-GM130622.

## Appendix A Proofs

## A.1 Proof of Theorem 8

The proof for Theorem 8 follows a similar structure to that for Theorem 1 in Lock et al. (2022). We first establish four lemmas, Lemmas 11, 12, 13, and 14, which are useful to prove the main result.

Lemma 11. Take two decompositions { ˆ S ( k ) ·· } K k =1 ∈ S ˆ X and { ˜ S ( k ) ·· } K k =1 ∈ S ˆ X , and assume that both minimize the structured penalty for given { λ } K k =1 :

$$f _ { p e n } ( \{ \tilde { \mathbf S } _ { \bullet \, } ^ { ( k ) } \} _ { k = 1 } ^ { K } ) = f _ { p e n } \left ( \{ \tilde { \mathbf S } _ { \bullet \, } ^ { ( k ) } \} _ { k = 1 } ^ { K } \right ) = m i n \, f _ { p e n } ( \{ \mathbf S _ { \bullet \, } ^ { ( k ) } \} _ { k = 1 } ^ { K } ).$$

Then, for any α ∈ [0 , 1] ,

$$| | \alpha \hat { \mathbf S } _ { \bullet } ^ { ( k ) } + ( 1 - \alpha ) \tilde { \mathbf S } _ { \bullet } ^ { ( k ) } | | _ { * } = \alpha | | \hat { \mathbf S } _ { \bullet } ^ { ( k ) } | | _ { * } + ( 1 - \alpha ) | | \tilde { \mathbf S } _ { \bullet } ^ { ( k ) } | | _ { * }$$

for k = 1 , . . . , K .

Proof. Because S ˆ X is a convex space and f pen is a convex function, the set of minimizers of f pen over S ˆ X is also convex. Thus,

$$f _ { \text{pen} } \left ( \{ \alpha \hat { \mathbf S } _ { \bullet } ^ { ( k ) } + ( 1 - \alpha ) \tilde { \mathbf S } _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } \right ) = \min _ { \mathbb { S } _ { \mathfrak x } } f _ { \text{pen} } ( \{ \mathbf S _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } ).$$

The result follows from the convexity of P λ , which implies that for any two matrices of equal size ˆ A and ˜ , A

$$P _ { \lambda } ( \alpha \hat { \mathbf A } + ( 1 - \alpha ) \tilde { \mathbf A } ) \leq \alpha P _ { \lambda } ( \hat { \mathbf A } ) + ( 1 - \alpha ) P _ { \lambda } ( \tilde { \mathbf A } ).$$

Applying (A1) to each additive term in f pen gives

$$( A 2 )$$

$$f _ { \text{pen} } \left ( \{ \alpha \hat { S } _ { \bullet } ^ { ( k ) } + ( 1 - \alpha ) \tilde { S } _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } \right ) & \leq \alpha f _ { \text{pen} } ( \{ \tilde { S } _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } ) + ( 1 - \alpha ) f _ { \text{pen} } ( \{ \tilde { S } _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } ) \\ & = \min _ { \tilde { S } _ { \tilde { X } } } f _ { \text{pen} } ( \{ \tilde { S } _ { \bullet } ^ { ( k ) } \} _ { k = 1 } ^ { K } ).$$

Because { α ˆ S ( k ) ·· +(1 -α ) ˜ S ( k ) ·· } K k =1 ∈ S ˆ X , the inequality in (A2) must be an equality, and it follows that the inequality (A1) must be an equality for each penalized term in the decomposition.

Lemma 12. Take two matrices ˆ A and ˜ A . If || ˆ + ˜ A A || ∗ = || ˆ A || ∗ + || ˜ A || ∗ , and UD V + T is the SVD of ˆ + ˜ A A , then ˆ A = ˆ ˆ ˆ UDV T where ˆ D is diagonal and || ˆ A || ∗ = || ˆ D || ∗ , and ˜ A = UDV ˜ T where ˜ D is diagonal and || ˜ A || ∗ = || ˜ D || ∗ .

Proof.

$$P r o o f. \text{ This result is proved in Lemma 2 of Lock et al. (2022).}$$

Lemma 13. Take two matrices ˆ A and ˜ A . If P λ ( ˆ A + ˜ ) = A P λ ( ˆ ) + A P λ ( ˜ ) A , and UD V + T is the SVD of ˆ + ˜ A A , then ˆ A = ˆ ˆ ˆ UDV T where ˆ D is diagonal and P λ ( ˆ ) = A P λ ( ˆ ) D , and ˜ A = UDV ˜ T where ˜ D is diagonal and P λ ( ˜ ) = A P λ ( ˜ ) D .

Proof. The result follows from repeated applications of Lemma 12 to the rank-1 terms of ˆ A and ˜ . A Note that P λ ( ˆ A + ˜ ) = A P λ ( ˆ ) + A P λ ( ˜ ) A implies

$$P _ { \lambda } \left ( \sum _ { r } \mathbf U _ { \mathbf A } [ \cdot, r ] \mathbf D _ { \mathbf A } [ r, r ] \mathbf V _ { \mathbf A } [ \cdot, r ] ^ { T } + \mathbf U _ { \mathbf A } [ \cdot, r ] \mathbf D _ { \mathbf A } [ r, r ] \mathbf V _ { \mathbf A } [ \cdot, r ] ^ { T } \right ) = \sum _ { r } \lambda [ r ] ( \mathbf D _ { \mathbf A } [ r, r ] + \mathbf D _ { \mathbf A } [ r, r ] ),$$

which requires that

$$| | \mathsf U _ { \hat { A } } [ \cdot, r ] \mathsf D _ { \hat { A } } [ r, r ] \mathsf V _ { \hat { A } } [ \cdot, r ] ^ { T } + \mathsf U _ { \hat { A } } [ \cdot, r ] \mathsf D _ { \hat { A } } [ r, r ] \mathsf V _ { \hat { A } } [ \cdot, r ] ^ { T } | | _ { * } \\ = | | \mathsf U _ { \hat { A } } [ \cdot, r ] \mathsf D _ { \hat { A } } [ r, r ] \mathsf V _ { \hat { A } } [ \cdot, r ] ^ { T } | | _ { * } + | | \mathsf U _ { \hat { A } } [ \cdot, r ] \mathsf D _ { \hat { A } } [ r, r ] \mathsf V _ { \hat { A } } [ \cdot, r ] ^ { T } | | _ { * }$$

for each r .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00237v1-figure-dddd1bf9-58ab-49ac-85e9-3762299dd0de.jpg)

Lemma 14. Take two decompositions { ˆ S ( k ) ·· } K k =1 ∈ S ˆ X and { ˜ S ( k ) ·· } K k =1 ∈ S ˆ X , and assume that both satisfy

$$f _ { p e n } ( \{ \tilde { \mathbf S } _ { \dots } ^ { ( k ) } \} _ { k = 1 } ^ { K } ) = f _ { p e n } \left ( \{ \tilde { \mathbf S } _ { \dots } ^ { ( k ) } \} _ { k = 1 } ^ { K } \right ) = m i n \, f _ { p e n } ( \{ \mathbf S _ { \dots } ^ { ( k ) } \} _ { k = 1 } ^ { K } ).$$

Then, ˆ S ( k ) ·· = U ( k ) · ˆ DV ( k T ) · and ˆ S ( k ) ·· = U ( k ) · ˜ D ( k ) V ( k T ) · where U ( k ) · : M × R k and V ( k ) · : N × R k have orthonormal columns, and ˆ D ( k ) and ˜ D ( k ) are diagonal.

Proof. The result follows as a direct corollary of Lemmas 11 and 13, as Lemma 11 implies P λ ( ˆ S ( k ) ·· + ˜ S ( k ) ·· ) = P λ ( ˆ S ( k ) ·· | )+ P λ ( ˜ S ( k ) ·· ) for each k , and then Lemma 13 implies the result.

Theorem 8 is then established as follows:

Proof. Take two decomposition { ˆ S ( k ) ·· } K k =1 and { ˜ S ( k ) ·· } K k =1 that satisfy properties 1., 2., and 3. of Theorem 1; we will show that { ˆ S ( k ) ·· } K k =1 = { ˜ S ( k ) ·· } K k =1 . For each k = 1 , . . . , K , write ˆ S ( k ) ·· = U ( k ) · ˆ DV ( k T ) · and ˆ S ( k ) ·· = U ( k ) · ˜ D ( k ) V ( k T ) · as in Lemma 14. Then, it suffices to show that ˆ D ( k ) [ r, r ] = ˜ D ( k ) [ r, r ] for all k, r .

First, consider module k = 1 with R [ · , 1] = [1 0 · · · 0] T and C [ · , 1] = [1 0 · · · 0] T . By way of contradiction, assume ˆ D (1) [1 , 1] &gt; 0 and ˜ D (1) [1 , 1] = 0. The linear independence of { V ( k ) j [ · , r ] : ˆ D ( k ) [ r, r ] &gt; 0 } and { V ( k ) j [ · , r ] : ˜ D ( k ) [ r, r ] &gt; 0 } implies that

$$\text{row} ( \mathbf X _ { \bullet } ) = \text{span} \{ \mathbf U _ { \bullet } ^ { ( k ) } [ \cdot, r ] \, \colon \hat { \mathbf D } ^ { ( k ) } [ r, r ] > 0 \} = \text{span} \{ \{ \mathbf U _ { \bullet } ^ { ( k ) } [ \cdot, r ] \, \colon \tilde { \mathbf D } ^ { ( k ) } [ r, r ] > 0 \}.$$

Thus, U (1) [ · , 1]] ∈ span {{ U ( k ) · [ · , r ] : ˜ D ( k ) [ r, r ] &gt; 0 } , and it follows from the orthogonality of U (1) [ · , 1] and { U (1) [ · , r ] , r &gt; 1 } that

$$\mathbf U _ { \cdot } ^ { ( 1 ) } [ \cdot, 1 ] \in \text{span} \{ \mathbf U _ { \cdot } ^ { ( k ) } [ \cdot, r ] \, \colon \tilde { \mathbf D } ^ { ( k ) } [ r, r ] > 0 \text{ and } k > 1 \}.$$

Moreover, because U (1) i = 0 for any i &gt; 1 and { U ( k ) i [ · , r ] : ˜ D ( k ) [ r, r ] &gt; 0 } are linearly independent it follows that

$$U _ { \cdot } ^ { ( 1 ) } [ \cdot, 1 ] \in \text{span} \{ U _ { \cdot } ^ { ( k ) } [ \cdot, r ] \, \colon \tilde { D } ^ { ( k ) } [ r, r ] > 0, \ k > 1, \text{ and } R [ i, k ] = 0 \text{ for any } i > 1 \}.$$

Note that (A3) implies U (1) 1 [ · , 1] ∈ row( X 12 + · · · +row( X 1 J ), however, this is contradicted by the linear independence of U (1) 1 [ · , 1] and { U ( k ) i [ · , r ] : ˆ D ( k ) [ r, r ] &gt; 0 , k &gt; 1 } . Thus, we conclude that ˜ D (1) [1 , 1] &gt; 0 implies ˜ D (1) [1 , 1] &gt; 0. Analogous arguments show that ˜ D ( k ) [ r, r ] &gt; 0 if and only if ˜ D ( k ) [ r, r ] &gt; 0 for any pair ( r, k ). It follows that { U ( k ) i [ · , r ] : ˆ D ( k ) [ r, r ] &gt; 0 or ˜ D ( k ) [ r, r ] &gt; 0 } are linearly independent for i = 1 , . . . I , and { V ( k ) j [ · , r ] : ˆ D ( k ) [ r, r ] &gt; 0 or ˜ D ( k ) [ r, r ] &gt; 0 } are linearly independent

for j = 1 , . . . , J . Thus,

$$\sum _ { k = 1 } ^ { K } \mathbf U _ { \cdot } ^ { ( k ) } ( \hat { \mathbf D } ^ { ( k ) } - \tilde { \mathbf D } ^ { ( k ) } ) \mathbf V _ { \cdot } ^ { ( k ) T } = \sum _ { k = 1 } ^ { K } \hat { \mathbf S } _ { \cdot \cdot } ^ { ( k ) } - \tilde { \mathbf S } _ { \cdot \cdot \cdot } ^ { ( k ) } = \mathbf X _ { \cdot \cdot } - \mathbf X _ { \cdot \cdot \cdot } = 0$$

implies that ˆ D ( k ) [ r, r ] = ˜ D ( k ) [ r, r ] for all k, r .

## A.2 Proof of Theorem 10

It is clear that the partial 'M-steps' in Algorithm 2 each partially optimize the free energy. It remains to be shown that the missing entries minimize free energy when they are replaced with their conditional expectation under the model. Note that the free energy (8) can be expressed as

$$E _ { q } \log \left ( \prod _ { k = 1 } ^ { K } q _ { u, k } ( \mathbf U _ { k } ) q _ { v, k } ( \mathbf V _ { k } ) \right ) - E _ { q } \log \left ( p ( \{ \mathbf U \} _ { k = 1 } ^ { K } \, | \, \tau _ { \mathbf U } ) p ( \{ \mathbf V \} _ { k = 1 } ^ { K } \, | \, \tau _ { \mathbf V } ) \right ) - E _ { q } \log \left ( p ( \{ \mathbf X.. \, | \, \mathbf U, \mathbf V, \sigma ) \,.$$

Consider the last term in the expression, decomposing the missing and non-missing entries of X ·· :

$$E _ { q } \log \left ( p ( \mathbf X _ { \bullet } \left | \, \mathbf U, \mathbf V, \sigma ) \right ) & = E _ { q } \log \left ( p ( \left \{ \mathbf X _ { \bullet } [ m, n ] \, \colon \left ( m, n \right ) \notin \mathcal { M } \right \} \, \mathbf U \,, \mathbf V, \sigma ) p \left \{ \mathbf X _ { \bullet } [ m, n ] \, \colon \left ( m, n \right ) \in \mathcal { M } \right \} \, \mathbf U \,, \mathbf V, \sigma ) \right ) \\ & = E _ { q } \log \left ( p ( \left \{ \mathbf X _ { \bullet } [ m, n ] \, \colon \left ( m, n \right ) \notin \mathcal { M } \right \} \, \mathbf U \,, \mathbf V, \sigma ) \right ) + \quad ( A 4 ) \\ & = E _ { q } \log \left ( p ( \left \{ \mathbf X _ { \bullet } [ m, n ] \, \colon \left ( m, n \right ) \in \mathcal { M } \right \} \, \mathbf U \,, \mathbf V, \sigma ) ).$$

Missing data { X ·· [ m,n ] : ( m,n ) ∈ M} are random variables that are independent and normally distributed given { U V , } as in (7). Thus, the last term in (A4) can be expressed as follows:

$$1$$

$$E _ { q } \log \left ( p ( \{ \mathbf X _ { \cdot \left [ m, n \right ] } \colon ( m, n ) \in \mathcal { M } \} \ | \ U, \mathbf V, \sigma ) \right ) & = c _ { 1 } + \sum _ { ( m, n ) \in \mathcal { M } } \frac { 1 } { 2 \sigma ^ { 2 } [ m, n ] } E _ { q } ( \mathbf X _ { \cdot \left [ m, n \right ] - \mathbf S _ { \cdot \left [ m, n \right ] } ) ^ { 2 } } \\ & \geq c _ { 1 } + \sum _ { ( m, n ) \in \mathcal { M } } \frac { 1 } { 2 \sigma ^ { 2 } [ m, n ] } E _ { q } ( \mathbf S _ { \cdot \left [ m, n \right ] ) - \mathbf S _ { \cdot \left [ m, n \right ] } ) ^ { 2 }.$$

Thus, for any given q , the free energy is minimized at X ·· [ m,n ] = E S q ·· [ m,n ] = ˆ S ·· [ m,n ] for ( m,n ) ∈ M .

.

## References

Attias, H.: Inferring parameters and structure of latent variable models by variational Bayes. In: Proceedings of the Fifteenth Conference on Uncertainty in Artificial Intelligence, pp. 21-30. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA (1999)

$$\cdot$$

| Cai, T., Cai, T.T., Zhang, A.: Structured matrix completion with applications to genomic data integration. Journal of the American Statistical Association 111 (514), 621-633 (2016)                   |
|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Feng, Q., Jiang, M., Hannig, J., Marron, J.: Angle-based joint and individual variation explained. Journal of multivariate analysis 166 , 241-265 (2018)                                               |
| Fox, C.W., Roberts, S.J.: A tutorial on variational Bayesian inference. Artificial intelligence review 38 , 85-95 (2012)                                                                               |
| Gavish, M., Donoho, D.L.: Optimal shrinkage of singular values. IEEE Transactions on Information Theory 63 (4), 2137-2152 (2017)                                                                       |
| Gaynanova, I., Li, G.: Structural learning and integrative decomposition of multi-view data. Biometrics 75 (4), 1121-1132 (2019)                                                                       |
| Hastie, T., Mazumder, R.: softImpute: Matrix Completion Via Iterative Soft- Thresholded SVD. (2021). R package version 1.4-1. https://CRAN.R-project.org/ package=softImpute                           |
| Hastie, T., Tibshirani, R., Narasimhan, B., Chu, G.: Impute: Impute: Imputation for Microarray Data. (2021). R package version 1.68.0                                                                  |
| Josse, J., Sardy, S.: Adaptive shrinkage of singular values. Statistics and Computing 26 , 715-724 (2016)                                                                                              |
| Kolda, T.G., Bader, B.W.: Tensor decompositions and applications. SIAM review 51 (3), 455-500 (2009)                                                                                                   |
| Kurucz, M., Bencz´r, u A.A., Csalog´ny, a K.: Methods for large scale SVD with missing values. In: Proceedings of KDD Cup and Workshop, vol. 12, pp. 31-38 (2007)                                      |
| Lock, E.F., Hoadley, K., Marron, J.S., Nobel, A.B.: Joint and Individual Variation Explained (JIVE) for integrated analysis of multiple data types. The Annals of Applied Statistics 7 (1), 523 (2013) |
| Lock, E.F., Park, J.Y., Hoadley, K.: Bidimensional linked matrix factorization for pan-omics pan-cancer analysis. Annals of Applied Statistics 16 (1), 193-215 (2022)                                  |
| Liu, Y., Zhang, A.: StructureMC: Structured Matrix Completion. (2019). R package version 1.0. https://CRAN.R-project.org/package=StructureMC                                                           |
| Mazumder, R., Hastie, T., Tibshirani, R.: Spectral regularization algorithms for learn- ing large incomplete matrices. The Journal of Machine Learning Research 11 , 2287-2322 (2010)                  |
| Nakajima, S., Sugiyama, M.: Analysis of empirical map and empirical partially bayes: Can they be alternatives to variational Bayes? In: Artificial Intelligence and                                    |

| Statistics, pp. 20-28 (2014). PMLR                                                                                                                                                                                                                                                                           |
|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Nakajima, S., Sugiyama, M., Babacan, S.D., Tomioka, R.: Global analytic solution of fully-observed variational Bayesian matrix factorization. Journal of Machine Learning Research 14 (1), 1-37 (2013)                                                                                                       |
| Nakajima, S., Tomioka, R., Sugiyama, M., Babacan, S.: Perfect dimensionality recovery by variational Bayesian pca. In: Pereira, F., Burges, C.J., Bottou, L., Wein- berger, K.Q. (eds.) Advances in Neural Information Processing Systems, vol. 25 (2012)                                                    |
| O'Connell, M.J., Lock, E.F.: R.JIVE for exploration of multi-source molecular data. Bioinformatics 32 (18), 2877-2879 (2016)                                                                                                                                                                                 |
| Park, J.Y., Lock, E.F.: Integrative factorization of bidimensionally linked matrices. Biometrics 76 (1), 61-74 (2020)                                                                                                                                                                                        |
| Rudelson, M., Vershynin, R.: Non-asymptotic theory of random matrices: extreme singular values. In: Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II-IV: Invited Lectures, pp. 1576-1602 (2010). World Scientific |
| Shabalin, A.A., Nobel, A.B.: Reconstruction of a low-rank matrix in the presence of Gaussian noise. Journal of Multivariate Analysis 118 , 67-76 (2013)                                                                                                                                                      |
| Schouteden, M., Van Deun, K., Wilderjans, T.F., Van Mechelen, I.: Performing disco- sca to search for distinctive and common information in linked data. Behavior research methods 46 (2), 576-587 (2014)                                                                                                    |
| Samorodnitsky, S., Wendt, C.H., Lock, E.F.: Bayesian simultaneous factorization and prediction using multi-omic data. Computational Statistics & Data Analysis 197 , 107974 (2024)                                                                                                                           |
| TCGA Research Network, et al. : Comprehensive molecular portraits of human breast tumors. Nature 490 (7418), 61 (2012)                                                                                                                                                                                       |
| Tang, T.M., Allen, G.I.: Integrated principal components analysis. Journal of Machine Learning Research 22 (198), 1-71 (2021)                                                                                                                                                                                |
| Wang, W., Stephens, M.: Empirical Bayes matrix factorization. Journal of Machine Learning Research 22 (120), 1-40 (2021)                                                                                                                                                                                     |
| Yang, Z., Michailidis, G.: A non-negative matrix factorization method for detecting modules in heterogeneous omics multi-modal data. Bioinformatics 32 (1), 1-8 (2016)                                                                                                                                       |
| Yi, S., Wong, R.K.W., Gaynanova, I.: Hierarchical nuclear norm penalization for multi- view data integration. Biometrics (2023)                                                                                                                                                                              |
Eric F. Lock
Empirical Bayes Linked Matrix Decomposition
Data for several applications in diverse fields can be represented as
multiple matrices that are linked across rows or columns. This is particularly
common in molecular biomedical research, in which multiple molecular "omics"
technologies may capture different feature sets (e.g., corresponding to rows in
a matrix) and/or different sample populations (corresponding to columns). This
has motivated a large body of work on integrative matrix factorization
approaches that identify and decompose low-dimensional signal that is shared
across multiple matrices or specific to a given matrix. We propose an empirical
variational Bayesian approach to this problem that has several advantages over
existing techniques, including the flexibility to accommodate shared signal
over any number of row or column sets (i.e., bidimensional integration), an
intuitive model-based objective function that yields appropriate shrinkage for
the inferred signals, and a relatively efficient estimation algorithm with no
tuning parameters. A general result establishes conditions for the uniqueness
of the underlying decomposition for a broad family of methods that includes the
proposed approach. For scenarios with missing data, we describe an associated
iterative imputation approach that is novel for the single-matrix context and a
powerful approach for "blockwise" imputation (in which an entire row or column
is missing) in various linked matrix contexts. Extensive simulations show that
the method performs very well under different scenarios with respect to
recovering underlying low-rank signal, accurately decomposing shared and
specific signals, and accurately imputing missing data. The approach is applied
to gene expression and miRNA data from breast cancer tissue and normal breast
tissue, for which it gives an informative decomposition of variation and
outperforms alternative strategies for missing data imputation.