ACMTF-R: supervised multi-omics data integration uncovering shared and distinct outcome-associated variation

19 The rapid growth of high -dimensional biological data has necessitated advanced data fusion 20 techniques to integrate and interpret complex multi -omics and longitudinal datasets. Advanced 21 Coupled Matrix and Tensor Factorization (ACMTF) has emerged as a powerful framework for 22 uncovering common, local, and distinct sources of variation across datasets. However, ACMTF lacks 23 the ability to model variation linked to a dependent variable, limiting its applicability to studies 24 investigating biological phenotypes. N-way Partial Least Squares (NPLS) is a supervised method that 25 identifies variation in relation to a dependent variable but lacks the ability to identify common, 26 local and distinct sources of variation across multiple datasets. To bridge the gap between data 27 exploration and prediction, we introduce ACMTF-Regression (ACMTF-R), an extension of ACMTF 28 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 2 that incorporates a regression term, allowing for the simultaneous decomposition of multi-way data 29 while explicitly capturing variation associated with an outcome variable. 30 We present a detailed mathematical formulation of ACMTF-R, including its optimisation algorithm 31 and implementation. Through extensive simulations, we systematically evaluate its ability to 32 recover a small 𝒚-related component shared between multiple blocks, its robustness to noise, and 33 the impact of the tuning parameter ( 𝜋) which controls the balance between data exploration and 34 outcome prediction. Our results demonstrate that ACMTF -R can robustly identify the 𝒚-related 35 component, correctly identifying outcome-associated shared and distinct variation, distinguishing 36 it from existing approaches such as N-way Partial Least Squares and ACMTF. 37 To validate its applicability in a real -world setting, we apply ACMTF -R to a multi -omics dataset 38 integrating human milk microbiome, human milk metabolome, and infant faecal microbiome data, 39 investigating how maternal pre -pregnancy BMI affects microbial and metabolic signatures. 40 ACMTF-R successfully identifies novel mother-infant relationships associated with maternal pre-41 pregnancy BMI, underscoring its utility in multi-omics research. Our findings establish ACMTF-R 42 as a versatile tool for multi-way data fusion, offering new insights into complex biological systems 43 by integrating common, local, and distinct variation in the context of a dependent variable. 44 Key Words 45 Coupled matrix and tensor factorization, multi-omics data integration, tensor decomposition 46

47 The advent of high -throughput omics technologies has led to a rapid growth of complex, high -48 dimensional biological datasets, including genomics [1–3], transcriptomics [4], proteomics [5], 49 metabolomics [4], and microbiom e [6]. When measured on the same samples, t hese multi-omics 50 datasets give researchers t he opportunity to obtain a systems -level understanding of biological 51 processes [5]. However, their integration and analysis pose significant challenges due to data 52 heterogeneity [5, 7], high dimensionality [4, 6], and complex interactions between different omics 53 datasets [8, 9] . Advanced data analytical methods are required to extract meaningful biological 54 insights by distinguishing variation that is shared across datasets from variation that is dataset -55 specific [10–13]. 56 A critical need in the analysis of multi-omics data is identifying common, local and distinct sources 57 of variation (Figure 1; [8–10, 14]). Common variation is shared across all datasets, local variation is 58 shared between a subset of datasets, and d istinct variation is unique to a dataset. Identifying the 59 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 3 Common, Local and Distinct (CLD) structure of a multi -omics dataset elucidates the complex 60 biological processes that occur within the system. 61 Advanced Coupled Matrix and Tensor Factorization (ACMTF) has emerged as a powerful multi -62 way data integration approach capable of uncovering common, local, and distinct sources of 63 variation across datasets [10, 15–18]. However, ACMTF is unsupervised and therefore does not allow 64 the identification of shared and distinct variation that relates to a dependent variable. N-way Partial 65 Least Squares (N-PLS) [19] is a supervised method that identifies variation in relation to a dependent 66 variable, but is limited to analysing one dataset at a time and thus cannot simultaneously identify 67 common, local, and distinct variation across multiple datasets. 68 To bridge the gap between data exploration and prediction, we introduce Advanced Coupled Matrix 69 and Tensor Factorization Regression (ACMTF -R). ACMTF-R extends ACMTF by incorporating a 70 regression term into the objective function, enabling the extraction of common, local, and distinct 71 variation while simultaneously predicting an outcome . By introducing a tuning parameter ( 𝜋), 72 which has been used successfully in other approaches such as Principal Covariates Regression 73 (PCovR) [20, 21] and Multiway Covariates Regression ( MCovR) [22], the solution can be steered 74 towards explaining the data or towards predicting 𝒚. 75 In this study, we first provide a detailed overview of the ACMTF framework, and the 76 methodological development represented by ACMTF -R. We then present a simulation -based 77 approach to systematically evaluate the ability of ACMTF-R to capture a small, hidden, 𝒚-related 78 component under different noise conditions and tuning parameter values. Finally, we apply 79 ACMTF-R to a real -world dataset, integrating human milk (HM) microbiome, HM metabolome, 80 and infant gut microbiome data to study how maternal pre-pregnancy body mass index (ppBMI) is 81 associated with these data. Our findings highlight the potential of ACMTF-R as a versatile tool for 82 multi-omics data integration, facilitating the discovery of biologically meaningful common, local 83 and distinct variation, while accommodating both exploratory and predictive research objectives. 84 85 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 4 86 Figure 1: common, local and distinct variation. Common variation 𝐶123 is shared between all 87 datasets 𝑿(1), 𝑿(2) and 𝑿(3). Local variation is shared between some of the datasets. Hence 𝐿12 is 88 shared between 𝑿(1) and 𝑿(2), but not 𝑿(3). Similarly, 𝐿13 is shared between 𝑿(1) and 𝑿(3), but not 89 𝑿(2), and 𝐿23 is shared between 𝑿(2) and 𝑿(3), but not 𝑿(1). Distinct variation is unique to one 90 dataset: 𝐷1 for block 𝑿(1), 𝐷2 for block 𝑿(2), and 𝐷3 for block 𝑿(3). Together, they make up the 91 Common, Local and Distinct (CLD) structure of the data. 92 93 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 5

94 Notation and definitions 95 We briefly define the mathematical notation that will be used throughout this paper. The notation 96 proposed by Kiers [23] is followed with some minor extensions for multi-omics. Scalars are denoted 97 by lower-case letters (e.g., 𝑎, 𝑏, 𝑐). Column vectors are denoted by boldface lowercase letters (e.g., 98 𝒂, 𝒃, 𝒄). Matrices are denoted by boldface capital letters (e.g., 𝑿, 𝒀, 𝒁). Three-way arrays are denoted 99 by underlined boldface capital letters (e.g., 𝑿, 𝒀, 𝒁). While 𝑿 can more generally be used to indicate 100 four-way or higher-order arrays, this paper does not discuss such data. 101 The characters 𝐼, 𝐽, and 𝐾 are reserved to indicate the first, second and third mode of an array. A 102 two-way array 𝑿 will be assumed to be of size 𝐼 × 𝐽, while a three-way array 𝑿 will be assumed to 103 be of size 𝐼 × 𝐽 × 𝐾. Similarly, lowercase characters 𝑖, 𝑗, and 𝑘 will be used as indices for the first, 104 second, and third modes, respectively. Columns of a matrix are vectors with a subscript of their 105 index (e.g., 𝒙1, 𝒙2, … , 𝒙𝐽 for some matrix 𝑿). Similarly, the i th-row and jth-column entry is denoted 106 as a scalar with subscripts of their indices (e.g., 𝑥𝑖𝑗). For three-way arrays, the entries per mode are 107 given as subscripts (e.g., 𝑥𝑖𝑗𝑘). 108 The index 𝑝 = 1, … , 𝑃 is used to indicate the block number, which is denoted by a superscript (e.g., 109 𝑿(1), 𝑿(2), … , 𝑿(𝑃)) and the size is indicated likewise (e.g., the size of 𝑿(𝑝) is 𝐼 × 𝐽(𝑝) × 𝐾(𝑝)). In 110 ACMTF and ACMTF-R 𝑿(1), 𝑿(2), … , 𝑿(𝑃) can be any mixture of two-way or three-way arrays, but 111 this study focuses only on the case where all blocks are three-way arrays with a shared subject mode. 112 Hence the superscript (𝑝) is not used to describe the size of the first mode 𝐼. 113 Mathematically it is convenient to matricise a three-way array 𝑿, as this transformation facilitates 114 the description of multi -way models using matrix notation [23–25]. There are three different 115 variants of matricisation, depending on which mode of the three-way array is preserved as rows in 116 the resulting matrix. Hence, first mode matrici sation yields a matrix 𝑿𝑎 of size 𝐼 × 𝐽𝐾. Similarly, 117 second mode matricization produces a matrix 𝑿𝑏 of size 𝐽 × 𝐼𝐾, while third mode matricization 118 generates a matrix 𝑿𝑐 of size 𝐾 × 𝐼𝐽. 119 The vectorisation of an 𝐼 × 𝐽 matrix 𝑿 is denoted vec(𝑿), such that vec(𝑿) = [𝒙1 𝑇 𝒙2 𝑇 … 𝒙𝐽 𝑇] 𝑇 . The 120 Hadamard product [26] of two equally sized arrays 𝑿 and 𝒀 is denoted 𝑿 ∗ 𝒀, such that (𝑿 ∗ 𝒀)𝑖𝑗 =121 𝑥𝑖𝑗𝑦𝑖𝑗. The Kronecker product, denoted ⨂, of two matrices 𝑿 with size 𝐼 × 𝐽 and 𝒀 with size 𝐾 × 𝐿 122 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 6 produces a block matrix containing all pairwise products of their entries following (𝑿 ⨂ 𝒀)𝑖𝑘,𝑗𝑙 =123 𝑥𝑖𝑗𝑦𝑘𝑙. Thus, we have 124 𝑿⨂𝒀 = [ 𝑥11𝒀 ⋯ 𝑥1𝐽𝒀 ⋮ ⋱ ⋮ 𝑥𝐼1𝒀 ⋯ 𝑥𝐼𝐽𝒀 ] (1) 125 where the outcome of 𝑿⨂𝒀 is size 𝐼𝐾 × 𝐽𝐿. 126 The Khatri-Rao product, denoted ⨀, of two matrices 𝑿 of size 𝐼 × 𝑅 and 𝒀 of size 𝐽 × 𝑅 is defined 127 as the column-wise Kronecker product [27, 28] 128 𝑿⨀𝒀 = [ | | ⋯ | 𝒙1⨂𝒚1 𝒙2⨂𝒚2 ⋯ 𝒙𝑅⨂𝒚𝑅 | | ⋯ | ] = [ | | ⋯ | vec(𝒚1𝒙1 𝑇) vec(𝒚2𝒙2 𝑇) ⋯ vec(𝒚𝑅𝒙𝑅 𝑇) | | ⋯ | ] 129 where the outcome of 𝑿⨀𝒀 has size 𝐼𝐽 × 𝑅. 130 The Frobenius (or Euclidean) norm of a three-way array 𝑿 is defined as 131 ‖𝑿‖ = √∑ ∑ ∑ 𝑥𝑖𝑗𝑘 2 𝐾 𝑘=1 𝐽 𝑗=1 𝐼 𝑖=1 (2) 132 and the notation ‖ ∙ ‖ will also be used to refer to the Frobenius norm for matrices and vectors, 133 respectively [24]. The Moore-Penrose inverse [29] of 𝑿 is denoted 𝑿+. 134 In this paper we follow the canonical polyadic form of defining multi-way arrays as the sum of rank-135 one arrays with normalized factor matrices and a scalar 𝜆 to capture the magnitude per term [15, 136 30–32]. The notation ⟦𝝀𝑝 ; 𝑨, 𝑩(𝑝), 𝑪(𝑝)⟧ is used to describe the model of the three -way array 𝑿(𝑝), 137 defined elementwise as 138 𝑥𝑖𝑗𝑘 (𝑝) = ∑ 𝜆𝑝𝑓𝑎𝑖𝑓𝑏𝑗𝑓 (𝑝)𝑐𝑘𝑓 (𝑝) + 𝑒𝑖𝑗𝑘 (𝑝) 𝐹 𝑓=1 (3) 139 where columns of the loading matrices 𝑨, 𝑩(𝑝), 𝑪(𝑝) are norm 1 and the column vector 𝝀𝑝 is used to 140 modify the size of each identified component 𝑓 = 1, … , 𝐹 to norm 𝜆𝑝𝑓. For a two -way array, 141 equation (3) reduces to ⟦𝑨, 𝑩⟧ = 𝑨𝑩𝑇, assuming that 𝝀 = 1 with size 1 × 𝐹 and that the singular 142 values are multiplied into 𝑨 to stay in line with common practice. 143 144 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 7 Advanced Coupled Matrix and Tensor Factorization (ACMTF) 145 Advanced Coupled Matrix and Tensor Factorization is a joint factor isation approach that is 146 applicable to mixtures of matrices and tensors [10, 15, 16, 33]. The method seeks a decomposition 147 of the datasets 𝑿(1), 𝑿(2), … , 𝑿(𝑃) as 𝑿̂(𝑝) = ⟦𝝀𝑝; 𝑨, 𝑩(𝑝), 𝑪(𝑝)⟧ such that loading matrices 148 corresponding to shared mode s are equal between the blocks. In this paper, we will assume that 149 only the first mode is shared between all blocks. In that case, the matrix containing the first mode 150 loadings 𝑨 will be the same for every block 𝑝 and the superscript can be left out. 151 The loss function for ACMTF is defined as 152 𝑓(α, β, ε, 𝚲, 𝑨, 𝑩(1), 𝑪(1), … , 𝑩(𝑃), 𝑪(𝑃)) = ∑‖𝑿(𝑝) − 𝑿̂(𝑝)‖ 2 𝑃 𝑝=1 + 𝛽 ∑ ∑ √𝜆𝑝𝑓 2 + 𝜖 𝐹 𝑓=1 𝑃 𝑝=1 + 𝛼 ∑ ∑(‖𝒂𝑓‖− 1) 2 𝐹 𝑓=1 𝑃 𝑝=1 + 𝛼 ∑ ∑ (‖𝒃𝑓 (𝑝)‖ − 1) 2 𝐹 𝑓=1 𝑃 𝑝=1 + 𝛼 ∑ ∑ (‖𝒄𝑓 (𝑝)‖ − 1) 2 𝐹 𝑓=1 𝑃 𝑝=1 (4) 153 where the 𝑃 × 𝐹 matrix 𝚲 = [𝝀1 𝝀2 … 𝝀𝑃]𝑇 scales the factors per component and per data block, 154 𝛼 > 0 is a penalty setting for norm-1 components, 𝛽 > 0 is a sparsity penalty setting for 𝚲, 𝐹 is the 155 number of components, and 𝑝 is the block index. Here the 1-norm loss term for the elements in 𝚲 156 is replaced with a differentiable approximation 𝛽‖𝜆𝑝𝑓‖1 ≈ 𝛽√𝜆𝑝𝑓 2 + 𝜖 for a sufficiently small 𝜖 >157 0 [33, 34]. The originally suggested default settings are 𝛼 = 1, 𝛽 = 10−3 and 𝜖 = 10−8 [10]. The 158 gradient of the loss function is reported in Supplementary Methods. 159 The ACMTF framework of putting the norm of a component into 𝚲 and constraining the loading 160 vectors to become norm one, allows the 𝚲 matrix to encode the Common, Local, and Distinct (CLD) 161 structure of the data (Figure 1). For example, the CLD-structure of a three-block case 𝑿(1), 𝑿(2), 𝑿(3) 162 can be encoded in 𝚲 through 163 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 8 𝚲 = 𝑿(1) 𝑿(2) 𝑿(3) 𝐶123 𝐿12 𝐿13 𝐿23 𝐷1 𝐷2 𝐷3 [ 𝜆11 𝜆12 λ13 0 𝜆15 0 0 𝜆21 𝜆22 0 𝜆24 0 𝜆26 0 𝜆31 0 𝜆33 𝜆34 0 0 𝜆37 ] 164 where 𝜆𝑝𝑓 > 0 corresponds to data block 𝑝 contributing to component 𝑓, 0 corresponds to a data 165 block not contributing to a component, 𝐶123 indicates a common component shared between all 166 blocks, 𝐿12, 𝐿13 and 𝐿23 indicate local components shared between some blocks, and 𝐷1, 𝐷2 and 𝐷3 167 indicate distinct components unique to one block. Sparsity is then achieved by applying an 𝐿1 168 penalty on the 𝚲-matrix. While this example uses 7 components to showcase the entire CLD -169 structure, multiple components may be needed to fully describe a specific type of variation in real 170 data (e.g., two components for 𝐶123) and some types of variation may be absent altogether. 171 Advanced Coupled Matrix and Tensor Factorization Regression (ACMTF-R) 172 We present Advanced Coupled Matrix and Tensor Factorization Regression (ACMTF-R) to describe 173 common, local, and distinct variation of interest in the data blocks in the context of a dependent 174 variable 𝒚 by adding a regression term and a tuning parameter 𝜋 to the loss function. The 175

of a tuning parameter has been used successfully in other approaches such as Principal 176 Covariates Regression (PCovR) [20, 21] and Multiway Covariates Regression (MCovR) [22] to steer 177 the solution towards explaining the data blocks or towards predicting 𝒚. This gives the user more 178 information about how the addition of 𝒚 affects the model of the data. We only define ACMTF-R 179 in the case where the first mode is shared across all data blocks and 𝒚. Hence 𝒚 is a column vector 180 of size 𝐼 × 1. The loss function of ACMTF-R is then defined as 181 𝑓(α, β, ε, π, 𝚲, 𝑨, 𝑩(1), 𝑪(1), … , 𝑩(𝑃), 𝑪(𝑃)) = π ∑‖𝑿(𝑝) − 𝑿̂(𝑝)‖ 2 𝑃 𝑝=1 + (1 − 𝜋)‖𝒚 − 𝑨𝝆‖2 + 𝛽 ∑ ∑ √𝜆𝑝𝑓 2 + 𝜖 𝐹 𝑓=1 𝑃 𝑝=1 + 𝛼 ∑ ∑(‖𝒂𝑓‖− 1) 2 𝐹 𝑓=1 𝑃 𝑝=1 + 𝛼 ∑ ∑ (‖𝒃𝑓 (𝑝)‖ − 1) 2 𝐹 𝑓=1 𝑃 𝑝=1 + 𝛼 ∑ ∑ (‖𝒄𝑓 (𝑝)‖ − 1) 2 𝐹 𝑓=1 𝑃 𝑝=1 (6) (6) 182 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 9 where 𝜋 is a tuning parameter used to focus the model on explaining the data blocks 𝑿(1), 𝑿(2), …, 183 𝑿(𝑃) versus predicting 𝒚, and 𝝆 is the 𝐹 × 1 vector of regression coefficients of the factor matrix of 184 the first mode 𝑨 onto 𝒚. Here 𝝆 is given when 𝑨 is defined, through 185 𝝆 = 𝑨+𝒚. (7) 186 The gradient of the loss function is reported in Supplementary Methods. 187 Prediction of 𝒚 using a new sample 188 A fitted ACMTF-R model can be used to predict 𝑦𝑛𝑒𝑤 for a new sample using a similar procedure as 189 Latent Root Regression [35, 36]. First, this requires a brief description of obtaining the subject mode 190 loadings for a new sample in CANDECOMP/PARAFAC (CP) [30, 32, 37, 38] . Subsequently, this 191

is extended to ACMTF-R, after which the prediction step is described. 192 In the CP case, finding the subject mode loadings 𝒂𝑛𝑒𝑤 with size 𝐹 × 1 for a sample 𝑿𝑛𝑒𝑤 with size 193 𝐽 × 𝐾 requires solving the problem 194 arg min 𝒂𝑛𝑒𝑤 ‖𝑿𝑛𝑒𝑤 − 𝑩diag(𝒂𝑛𝑒𝑤)𝑪𝑇 ‖𝐹 2 = arg min 𝒂 ‖vec(𝑿𝑛𝑒𝑤) − (𝑪⨀𝑩)𝒂𝑛𝑒𝑤‖𝐹 2 (9) 195 where 𝑩 and 𝑪 are the feature and time mode loadings of a previously fitted model [32]. This 196 problem statement comes down to the least squares solution 197 𝒂𝑛𝑒𝑤 = (𝑪⨀𝑩)+vec(𝑿𝑛𝑒𝑤) = 𝒁+vec(𝑿𝑛𝑒𝑤) (10) 198 where the block matrix 𝒁 = 𝑪⨀𝑩 with size 𝐽𝐾 × 𝐹 is used for convenience and where 𝒁+ is the 199 Moore-Penrose inverse of 𝒁. 200 In the case of ACMTF-R, the shared subject mode loadings across the blocks need to be identified. 201 This is done by finding the block matrix 𝒁(𝑝) with size 𝐾(𝑝)𝐽(𝑝) × 𝐹 through 202 𝒁(𝑝) = 𝝀𝑝𝑇⨀𝑪(𝑝)⨀𝑩(𝑝) (11) 203 which can be concatenated across all blocks into one matrix for all data blocks 204 1, … , 𝑃 simultaneously as 205 Z = ( 𝒁(1) 𝒁(2) ⋮ 𝒁(𝑃) ) (12) 206 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 10 where 𝒁 is of size (𝐾(1)𝐽(1) + 𝐾(2)𝐽(2) + ⋯ + 𝐾(𝑃)𝐽(𝑃)) × 𝐹. The shared subject mode loadings 𝒂𝑛𝑒𝑤 207 with size 𝐹 × 1 for a new sample 𝑿𝑛𝑒𝑤 can then be found by reusing Equation 10. 208 Finally, the prediction of 𝑦𝑛𝑒𝑤 is done using the ACMTF-R regression coefficients 𝝆 through 209 𝑦𝑛𝑒𝑤 = 𝝆𝑻𝒂𝑛𝑒𝑤. (12) 210 When multiple new samples are obtained, this procedure is performed per sample. The procedure 211 is performed by the npred() function in the supplied R package. 212 Implementation and stopping criteria 213 The decomposition of ACMTF and ACMTF -R is achieved through an all -at-once optimi sation 214 algorithm, originally developed in [15] for MATLAB, now available in the CMTFtoolbox package 215 for R on CRAN. The all-at-once optimisation is achieved by defining the loss and gradient function 216 and subsequently using the nonlinear conjugate gradient (NCG) method with Hestenes -Stiefel 217 updates and the Moré -Thuente line search as originally suggested [15, 18] . For the provided R 218 package, this was implemented using the mize package (v0.2.4, [39]). The provided package also 219 supports the Limited-memory Broyden-Fletcher-Goldfarb-Shanno algorithm (L-BFGS) to speed up 220 the line search as an experimental feature, but this setting is not used for any results in this paper 221 [40–42]. 222 ACMTF-R uses the same default values as ACMTF for the minimum value of the loss function 223 (default 10−10) and the minimum change in loss function between two function evaluations (default 224 10−10) as stopping criteria [33]. The all-at-once optimization implementation through the mize 225 package also defines three other stopping criteria: the maximum number of iterations to calculate 226 (default 10 000), the maximum number of function evaluation s that are allowed (default 10 000), 227 the minimum l2-norm of the gradient vector (default 10−10), and the absolute value of the size of 228 the parameter update (default 10−10). Fitting an ACMTF -R model for any of the datasets in this 229 paper takes less than a minute of computational time on an average computer (Microsoft Windows 230 11 Home v10.0.22631 with an 12 th Gen Intel® Core™ i5-12400F running 2500 MHz with 6 cores 231 or 12 logical processors and 16 Gb RAM). Multi -core parallelisation has been implemented as part 232 of the provided R package to allow many randomly initiali sed models to be fitted simultaneously. 233 This is needed to efficiently find the appropriate number of components through cross -validation, 234 as well as to find the global minimum for a given number of components. 235 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 11 Simulation approach 236 We demonstrate the effect of the ACMTF-R model parameters using various simulations (Figure 2). 237 We follow the guidelines provided in previous work [10]. We generate factor matrices 𝑨 ∈ ℝ𝐼×𝐹, 238 𝑩(1) ∈ ℝ𝐽(1)×𝐹, 𝑪(1) ∈ ℝ𝐾(1)×𝐹, 𝑩(2) ∈ ℝ𝐽(2)×𝐹, 𝑪(2) ∈ ℝ𝐾(2)×𝐹, 𝑩(3) ∈ ℝ𝐽(3)×𝐹 and 𝑪(3) ∈ ℝ𝐾(3)×𝐹 239 with random entries drawn from the standard normal distribution 𝑁~(0,1). The columns of the 240 factor matrices are normalized to unit norm. We set 𝐼 = 100, 𝐽(1) = 𝐽(2) = 𝐽(3) = 200, and 𝐾(1) =241 𝐾(2) = 𝐾(3) = 5 to resemble a typical pre-processed multi-omics dataset (for example, a mixture of 242 microbiome and metabolomics data). The factor matrices are used to create three third-order tensors 243 𝑿(1) ∈ ℝ100×200×5, 𝑿(2) ∈ ℝ100×200×5, 𝑿(3) ∈ ℝ100×200×5 through 244 𝑿(1) = ⟦𝝀1; 𝑨, 𝑩(1), 𝑪(1)⟧, 245 𝑿(2) = ⟦𝝀2; 𝑨, 𝑩(2), 𝑪(2)⟧, 246 and 247 𝑿(3) = ⟦𝝀3; 𝑨, 𝑩(3), 𝑪(3)⟧. 248 where we encode the required CLD -structure into 𝚲 consisting of eight components total: two 249 common components shared across all blocks, three local components, and a distinct component for 250 each block. This corresponds to the 𝚲-matrix 251 𝚲 = [ 1 1 1 𝛿 0 1 0 0 1 1 1 0 1 0 1 0 1 1 0 𝛿 1 0 0 1 ] 252 where all components are norm 1, except for the second local component whose size 𝛿 varies per 253 simulation. 254 Next, randomly distributed noise is added to the third-order tensors through 255 𝑿𝑛𝑜𝑖𝑠𝑦 (𝑝) = 𝑿(𝑝) + 𝜂𝑋 ‖𝑿(𝑝)‖ ‖𝑵(𝑝)‖𝑵(𝑝) (13) 256 where 𝜂𝑋 indicates the noise level, which is the same across all blocks. 257 We generate 𝒚 equal to the subject mode loadings of the second local component 258 𝒚 = 𝒂𝐿12 259 and some noise is added through 260 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 12 𝒚𝑛𝑜𝑖𝑠𝑦 = 𝒚 + 𝜂𝑦 ‖𝒚‖ ‖𝒏𝑦‖𝒏𝒚 (14) 261 where 𝜂𝑦 indicates the noise level on 𝒚. The outcome 𝒚𝑛𝑜𝑖𝑠𝑦 is then normalized to norm 1. 262 263 Figure 2: overview of the simulation setup of the study. All simulations contain a CLD-structure of 264 2 global components shared between blocks 𝑿(1), 𝑿(2), 𝑿(3), 3 local components , and a distinct 265 component unique to each block. This yields 8 components total for ACMTF -R to identify. All 266 components are norm 1, except for the second local component whose norm 𝛿 varies per simulation. 267 𝒚 is equal to the subject mode loadings of the second local component . Noise and terms 268 corresponding to components that have no contribution to a block have been removed for visual 269 clarity. Subject mode loadings (black lines) are equal between the data blocks per component but 270 are different between components. 271 272 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 13 Simulation 1: ACMTF-R can detect a small, hidden component by leveraging 𝒚 273 In the first simulation, the ability of ACMTF-R to detect a small, hidden component by leveraging 274 𝒚 was examined. This was done using the three -block simulation with 𝛿 = 0.03, 40% noise on 275 𝑿(1), 𝑿(2), 𝑿(3) (𝜂𝑋 = 0.40), and 5% noise on 𝒚 (𝜂𝑦 = 0.05). The data were then block scaled to 276 Frobenius norm 1. One hundred randomly initialised eight -component ACMTF models were 277 compared with one hundred randomly initialised eight -component ACMTF -R models using 278 different values of 𝜋 ∈ {0.1, 0.2, … ,0.9}, keeping all other parameter and convergence settings at 279 their default values. Subsequently, factor recovery was assessed using Factor Match Score and 280 Tucker Congruence Coefficient. The line search algorithm in most ACMTF -R models was 281 terminated due to the absolute change for the size of the parameter update stopping criterion, while 282 ACMTF was terminated only due to the minimum value of the loss function being reached 283 (Supplementary Figure 1). 284 Simulation 2: the benefit of supervision using 𝒚 versus the size of the hidden component 285 In the second simulation, the improvement of ACMTF -R compared to ACMTF was assessed by 286 changing the size of the small, hidden component. This was done using the three-block simulation 287 with 288 𝛿 ∈ {0.01, 0.02, 0.03, 0.04, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5}. 289 The noise on 𝑿(1), 𝑿(2), 𝑿(3) was set at 40% ( 𝜂𝑋 = 0.40), and the noise on 𝒚 was 5% (𝜂𝑦 = 0.05). 290 The data were then block scaled to Frobenius norm 1. One hundred randomly initialised eight -291 component ACMTF models were compared with one hundred randomly initialised eight -292 component ACMTF-R models using different values of 𝜋 ∈ {0.1, 0.2, … ,0.9} for each case, keeping 293 all other parameter and convergence settings at their default values. Subsequently, factor recovery 294 was assessed using Factor Match Score and Tucker Congruence Coefficient. The line search 295 algorithm in most ACMTF-R models was terminated due to the absolute change for the size of the 296 parameter update stopping criterion, while most ACMTF models were terminated only due to the 297 minimum value of the loss function being reached (Supplementary Figure 2). 298 299 300 301 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 14 Simulation 3: the benefit of supervision using 𝒚 versus noise on X 302 In the third simulation, the improvement of ACMTF -R compared to ACMTF was assessed by 303 changing the amount of noise on 𝑿(1), 𝑿(2), 𝑿(3). This was done using the three -block simulation 304 with 𝛿 = 0.03 and a variable amount of noise on 𝑿(1), 𝑿(2), 𝑿(3), 305 𝜂𝑋 ∈ {0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 10, 50, 100}, 306 and the noise on 𝒚 was 5% (𝜂𝑦 = 0.05). The data were then block scaled to Frobenius norm 1. One 307 hundred randomly initialised eight-component ACMTF models were compared with one hundred 308 randomly initialised eight -component ACMTF -R models using different values of 𝜋 ∈309 {0.1, 0.2, … ,0.9} for each case, keeping all other parameter and convergence settings at their default 310 values. Subsequently, factor recovery was assessed using Factor Match Score and Tucker 311 Congruence Coefficient. For most noise levels (𝜂𝑋 < 2), the line search algorithm in most ACMTF-312 R models was terminated due to the absolute change for the size of the parameter update stopping 313 criterion, while most ACMTF models were terminated only due to the minimum value of the loss 314 function being reached ( Supplementary Figure 3 ). At higher noise levels ( 𝜂𝑋 > 2), most ACMTF 315 models were terminated due to the relative change in the loss function being reached. 316 Simulation 4: the benefit of supervision using 𝒚 versus noise on 𝒚 317 In the fourth and final simulation, the improvement of ACMTF-R compared to ACMTF was assessed 318 by changing the amount of noise on 𝒚. This is done using the three-block simulation with 𝛿 = 0.03, 319 with 40% noise on 𝑿(1), 𝑿(2), 𝑿(3) (𝜂𝑋 = 0.40), and a variable amount of noise on 𝒚: 320 𝜂𝑦 ∈ {0.0, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4, 5, 10, 50, 100}. 321 The data were then block scaled to Frobenius norm 1. One hundred randomly initialised eight -322 component ACMTF models were compared with one hundred randomly initialised eight -323 component ACMTF-R models using different values of 𝜋 ∈ {0.1, 0.2, … ,0.9} for each case, keeping 324 all other parameter and convergence settings at their default values. Subsequently, factor recovery 325 was assessed using Factor Match Score and Tucker Congruence Coefficient. The line search 326 algorithm in most ACMTF-R models was terminated due to the absolute change for the size of the 327 parameter update stopping criterion, while ACMTF was terminated only due to the minimum value 328 of the loss function being reached (Supplementary Figure 4). 329 330 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 15 Simulation evaluation – Tucker Congruence Coefficient 331 In this paper, the recovered factors are compared against the input factors of the simulation using 332 the absolute Tucker Congruence Coefficient (TCC) [43, 44] through 333 |𝜙(𝒙, 𝒚)| = || ∑𝑥𝑖𝑖 𝑦𝑖 √∑𝑥𝑖 2𝑖 ∑𝑦𝑖 2𝑖 || (15) 334 where 𝒙 and 𝒚 are the recovered and true factors, respectively. The TCC can be interpreted as the 335 cosine of the angle between 𝒙 and 𝒚. A TCC value in the range of 0.85 -0.94 suggests that the 336 recovered and input factors are fairly similar, whereas a value of 0.95 -1 means that they are 337 essentially equal [45]. Although TCC is insensitive to scalar multiplication of the vectors (that is, 338 𝜙(𝒙, 𝒚) = 𝜙(𝑎𝒙, 𝑏𝒚) for any non-zero scalars 𝑎 and 𝑏), it is sensitive to sign flips between 𝒙 and 𝒚 339 [45]. For this reason, we use the absolute value instead. 340 Simulation evaluation – Factor Match Score 341 The Factor Match Score (FMS) is an adaptation of the Tucker Congruence Coefficient and indicates 342 the similarity of two multi-way models across all components. For model evaluation, FMS was used 343 to compare the decomposition against the input factors. This depends on the CLD -type of the 344 relevant factor. For the hidden local component 𝐿13 in the simulations, 𝐹𝑀𝑆𝐿13 is defined as 345 𝐹𝑀𝑆𝐿13 = |𝒂𝐿13 𝑇 𝒂̃𝐿13| ‖𝒂𝐿13‖‖𝒂̃𝐿13‖ |[𝒃𝐿13 (1) ] 𝑇 𝒃̃𝐿13 (1) | ‖𝒃𝐿13 (1) ‖ ‖𝒃̃𝐿13 (1) ‖ |[𝒄𝐿13 (1) ] 𝑇 𝒄̃𝐿13 (1) | ‖𝒄𝐿13 (1) ‖ ‖𝒄̃𝐿13 (1) ‖ |[𝒃𝐿13 (3) ] 𝑇 𝒃̃𝐿13 (3) | ‖𝒃𝐿13 (3) ‖ ‖𝒃̃𝐿13 (3) ‖ |[𝒄𝐿13 (3) ] 𝑇 𝒄̃𝐿13 (3) | ‖𝒄𝐿13 (3) ‖ ‖𝒄̃𝐿13 (3) ‖ (16) 346 where 𝒂𝐿13, 𝒃𝐿13, and 𝒄𝐿13 contain the true subject, feature and time loadings of component 𝐿13, 347 and 𝒂̃𝐿13, 𝒃̃𝐿13, and 𝒄̃𝐿13 are the best matching subject, feature and time loadings in the fitted model. 348 The FMS of the other components are defined equivalently. In this paper, the best matching 349 component is found using the Hungarian algorithm [46, 47] for all pairwise combinations of 350 components such that the FMS is maximised and the expected CLD structure is recovered 351 (Supplementary Methods). The FMS is a value in the range [0,1], where 𝐹𝑀𝑆 = 1 corresponds to 352 the model correctly recovering the input loadings. 353 The recovery of all input loadings can be quantified through 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙, defined as 354 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙 = 𝐹𝑀𝑆𝐶123 + 𝐹𝑀𝑆𝐶123 + 𝐹𝑀𝑆𝐿12 + 𝐹𝑀𝑆𝐿13 + 𝐹𝑀𝑆𝐿23 + 𝐹𝑀𝑆𝐷1 + 𝐹𝑀𝑆𝐷2 + 𝐹𝑀𝑆𝐷3 8 (17) 355 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 16 where 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙 = 1 corresponds to all input components being recovered correctly. Recovery of 356 the CLD structure is assessed using the Lambda Similarity Index, which is defined and reported in 357 the Supplementary Materials. 358 Example dataset: Jakobsen2025 359 We include an application study of ACMTF -R on longitudinally measured human milk (HM) 360 microbiome, HM metabolomics and infant faecal microbiome data. Full details on the preprocessing 361 of the data is outlined in the original paper [48]. 362 Briefly, the zOTU data and taxonomic information of the infant faecal microbiome and mother milk 363 microbiome were processed in R (v4.4.1, [49]). For the infant faecal microbiome data, features were 364 kept if they had ≤75% sparsity across the dataset. This step resulted in 93 out of 565 features being 365 selected. For the mother milk microbiome data, features were kept if they had ≤85% sparsity across 366 the dataset. This step resulted in 115 out of 707 features being selected. We performed a centred 367 log-ratio transformation with a pseudo -count of 1 to correct for compositionality [50, 51] . 368 Subsequently the data was converted to a three -way tensor, keeping missing samples as a row of 369 NAs. This resulted in three-way arrays of size 160 subjects × 93 microbial taxa × 3 time points and 370 169 subjects × 115 microbial taxa × 4 time points containing microbial abundances for the infant 371 faecal and mother milk microbiomes, respectively. 372 The human milk (HM) metabolomics data were processed using R (v4.4.1, [49]). Values below the 373 detection limit were imputed with a random value between 0 and the detection limit per metabolite 374 to preserve their distribution. Next, the dataset was (natural) log transformed to stabilise the 375 variance. The dataset was then converted to a three-way tensor of size 165 subjects × 70 metabolites 376 × 4 time points. 377 Since ACMTF -R assumes that the subject mode is shared between all data blocks, we took the 378 intersection of the subject identifiers across all three datasets and removed the other subjects. This 379 homogenised the subject mode size 𝐼 across all data blocks to 158 subjects. The data w ere then 380 centred across the subject mode such that variation per time point is only due to between -subject 381 variation, and scaled within the feature mode to make all features equally important for the 382 modelling procedure [52, 53]. As a result, the models will focus on inter-subject differences and the 383 features involved, while the time mode loadings indicate when these differences are largest. 384 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 17 After centring and scaling, the datasets were block scaled to Frobenius norm 1 to equali se the 385 variance across the blocks. The response variable 𝒚 of size 158 × 1 contained the maternal pre -386 pregnancy BMI for all mother-infant dyads, which was centred and subsequently scaled to norm 1 387 prior to modelling. ACMTF-R models were then created using the appropriate number of 388 components and regressing on maternal pre-pregnancy BMI. 389 Selecting the appropriate number of ACMTF components 390 This paper presents a two-armed approach for selecting the appropriate number of components for 391 ACMTF: a K-fold cross-validation scheme, and a random initialisation scheme. The full procedure 392 is available as the ACMTF_modelSelection() function in the supplied R package. 393 In the K-fold cross-validation scheme, one fold of the subjects is withheld for all blocks as test set 394 data, while randomly initialised ACMTF models are fitted on the training data using 𝑓 components, 395 keeping all other model parameters ( 𝛼, 𝛽, 𝜖) fixed. To avoid data leakage, the means, standard 396 deviations, and norms that are removed from the training data are also removed from the test set 397 data. Model stability is then reported by pairwise comparing the models with the lowest loss per 398 fold and calculating the 𝐹𝑀𝑆𝐶𝑉, defined as 399 𝐹𝑀𝑆𝐶𝑉 = 1 𝐹 ∑ |𝒃𝑓 𝑇𝒃̃𝑓| ‖𝒃𝑓‖‖𝒃̃𝑓‖ |𝒄𝑓 𝑇𝒄̃𝑓| ‖𝒄𝑓‖‖𝒄̃𝑓‖ 𝐹 𝑓=1 (18) 400 where 𝒃𝑓, and 𝒄𝑓 contain the feature and time loadings for component 𝑓 in one model, 𝒃̃𝑓, and 𝒄̃𝑓 401 are the best matching subject, feature and time loadings for component 𝑓 in a second model, and 402 the subject mode term is eliminated from the calculated due to not being comparable between CV 403 folds [15, 54]. The matching of components is achieved using the Hungarian algorithm [46, 47] for 404 all pairwise combinations of components such that the FMS is maximised. Hence, the distribution 405 of FMS values across all pairwise comparisons of CV (or jack -knife) folds give an indication of the 406 stability of the decomposition for a given number of components 𝐹. Previous work has established 407 that 95% of folds should have 𝐹𝑀𝑆𝐶𝑉 ≥ 0.95 for three-way data blocks and 𝐹𝑀𝑆𝐶𝑉 ≥ 0.9 for two-408 way data blocks in the ACMTF case for the decomposition to be considered replicable [54]. The 409 appropriate number of components is expected to have a stable 𝐹𝑀𝑆𝐶𝑉 near one across all 410 comparisons. 411 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 18 In the second arm of the approach, randomly initialised ACMTF models are fitted to the full data 412 using 𝑓 components. Model sensitivity to the starting values and to local minima is then assessed by 413 pairwise comparing the models and calculating the 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 through 414 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 = 1 𝐹 ∑ |𝒂𝑓 𝑇𝒂̃𝑓| ‖𝒂𝑓‖‖𝒂̃𝑓‖ 𝐹 𝑓=1 |𝒃𝑓 𝑇𝒃̃𝑓| ‖𝒃𝑓‖‖𝒃̃𝑓‖ |𝒄𝑓 𝑇𝒄̃𝑓| ‖𝒄𝑓‖‖𝒄̃𝑓‖ (19) 415 where 𝒂𝑓, 𝒃𝑓, and 𝒄𝑓 contain the subject, feature and time loadings for component 𝑓 in one model, 416 and 𝒂̃𝑓, 𝒃̃𝑓, and 𝒄̃𝑓 are the best matching subject, feature and time loadings for component 𝑓 in a 417 second model [15, 54] . The matching of components is again achieved using the Hungarian 418 algorithm [46, 47] for all pairwise combinations of components such that the FMS is maximised. 419 The distribution of FMS values across all pairwise comparisons of randomly initialised models for a 420 given number of components 𝐹 gives an indication of the stability of the decomposition. The 421 appropriate number of components is expected to have a stable 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 near one across all 422 randomly initialised models. 423 Additionally, it may be possible for ACMTF and ACMTF-R to split up common or local components 424 into distinct components with highly similar subject mode loadings when too many components 425 are chosen. This phenomenon can be detected by pairwise comparing all columns in the subject 426 mode per model and calculating the maximum absolute Tucker Congruence Coefficient 𝜑 [43–45]. 427 We define this as the Degeneracy Score, which is computed through 428 𝐷𝑒𝑔𝑒𝑛𝑒𝑟𝑎𝑐𝑦 𝑆𝑐𝑜𝑟𝑒= max 𝑝≠𝑞 |𝜑(𝒂𝑝, 𝒂𝑞)| (20) 429 where 𝑝 and 𝑞 indicate component numbers of the subject mode within the same model, and 𝑝 ≠430 𝑞. The Degeneracy Score is in [0, 1], will be close to one when a factor is split into two highly 431 collinear factors, and will have intermediate values otherwise. This procedure is only performed for 432 the randomly initialised ACMTF models since the complete subject mode is needed. The appropriate 433 number of components is expected to have a low Degeneracy Score across all randomly initialised 434 models. 435 Selecting the appropriate number of ACMTF-R components 436 This paper presents a similar two-armed approach for selecting the appropriate number of 437 components for ACMTF-R, reporting a few additional metrics related to the prediction and fit of 𝒚. 438 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 19 The full procedure is available as the ACMTFR_modelSelection() function in the supplied R 439 package. 440 In the cross-validation arm, the root mean squared error of cross-validation (RMSECV) is calculated 441 by using the models fitted on the training set to predict 𝒚 in the test set. The appropriate number of 442 components is expected to minimise the RMSECV while having a stable 𝐹𝑀𝑆𝐶𝑉 near one across all 443 comparisons. 444 In the random initialisation arm, the root mean squared error (RMSE) of 𝒚 and the variance 445 explained in 𝒚 are reported using the ACMTF -R model with the lowest loss. The appropriate 446 number of components is expected to minimise the RMSE while having a low Degeneracy Score 447 and a stable 𝐹𝑀𝑆𝑟𝑎𝑛𝑑𝑜𝑚 near one across all comparisons. 448 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 20

449 Simulation 1: ACMTF-R can detect a small, hidden component by leveraging 𝒚 450 We evaluated how the tuning parameter 𝜋 influences the recovery of a local component that is 451 weakly present in some data blocks and is strongly reflected in 𝒚. This was done by applying 452 ACMTF-R to a simulated three-tensor dataset with a known structure containing 2 common, 3 local, 453 and 3 distinct components and 40% noise. All the components were norm 1, except for the second 454 local component (𝐿13), which was present across 𝑿(1) and 𝑿(3) with norm 𝛿 = 0.03. Since this 455 component was equal to 𝒚 with 5% noise, we expected that lowering 𝜋 would enhance its detection. 456 By systematically varying the tuning parameter 𝜋 from 1 (fully focused on explaining the data 457 blocks) to 0 (fully focused on predicting 𝒚), we examined whether outcome information could reveal 458 a latent structure that would otherwise remain hidden (Figure 3). One hundred randomly initialised 459 eight-component ACMTF-R and ACMTF models were created for each setting of 𝜋. 460 As assessed by the Factor Match Score, the lowest -loss ACMTF -R models using 0.6 ≤ 𝜋 ≤ 0.9 461 showed an improvement in factor recovery for all components compared to the lowest-loss ACMTF 462 model (Figure 3A, Supplementary Figure 5). This was also the case for the input factor recovery of 463 the hidden component 𝐿13 (Figure 3B, Supplementary Figure 6). Meanwhile, the fraction of models 464 that correctly identified the input CLD structure was comparable to ACMTF (Supplementary Figure 465 7-8). 466 Next, the ACMTF-R models with the lowest loss for each 𝜋 setting were selected for comparison 467 with the lowest -loss ACMTF model ( Figure 3C-H). The absolute Tucker Congruence Coefficient 468 (𝜙) was used to compare the best matching recovered factor loadings with the input loadings of 𝐿13 469 per mode (Supplementary Methods). This analysis revealed that the recovered shared subject mode 470 (Figure 3C) and time mode loadings of 𝑿(1) (Figure 3E) and 𝑿(3) (Figure 3G) were almost identical 471 to the input loadings in the ACMTF-R models for 0.6 ≤ 𝜋 ≤ 0.9, while ACMTF failed to identify 472 them. While the recovered feature mode loadings in 𝑿(1) (Figure 3D) and 𝑿(3) (Figure 3F) were not 473 exactly equal to the input factor of the hidden component, their orde r largely matched the input 474 (Kendall’s 𝜏 = 0.35 − 0.45 for 0.4 ≤ 𝜋 ≤ 0.9; Supplementary Table 1). 475 Taken together, ACMTF-R recovered the hidden component 𝐿13 across intermediate values of 𝜋 by 476 leveraging 𝒚, while ACMTF did not.477 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 21 478 Figure 3: recovery of a small, hidden component 𝑳𝟏𝟑 using ACMTF and ACMTF-R. One hundred 479 randomly initialised ACMTF -R models were fitted per setting of tuning parameter 𝜋. These were 480 compared with one hundred randomly initialised ACMTF models. (A) Overall factor recovery, 481 measured using 𝐹𝑀𝑆𝑜𝑣𝑒𝑟𝑎𝑙𝑙, was higher in ACMTF -R using 0.6 ≤ 𝜋 ≤ 0.9 compared to ACMTF. (B) 482 Likewise, the 𝐹𝑀𝑆𝐿13 maximum was higher for ACMTF-R than ACMTF. Red data points in (A,B) indicate 483 the model with the lowest loss, which were compared in (C -H). The absolute Tucker Congruence 484 Coefficient (𝜙) between the recovered loadings and the input (C) subject mode loadings, (D) feature 485 mode loadings of 𝑿(1), (E) time mode loadings of 𝑿(1), (F) feature mode loadings of 𝑿(3), and (G) time 486 mode loadings of 𝑿(3) for hidden component 𝐿13. For all modes, factor recovery using ACMTF-R with 487 intermediate values of 𝜋 was higher compared to ACMTF. 488 489 490 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 22 Simulation 2: the benefit of supervision using 𝒚 depends on the size of the hidden component 491 We next evaluated how the size of the hidden component 𝐿13 affected the improvement of ACMTF-492 R compared to ACMTF. This was done using the same three-tensor setup as in the first simulation, 493 except that the size of the second local component (𝐿13) was varied between 𝛿 = 0.01 to 𝛿 = 0.5. 494 We expected the benefit of supervision to decrease as the size of the hidden component 𝛿 increased. 495 As before, the tuning parameter 𝜋 was varied from 0.1 to 0.9 and one hundred randomly initialised 496 eight-component ACMTF-R and ACMTF models were created for each setting (Figure 4). 497 Input factor recovery across all components improved in ACMTF-R compared to ACMTF at 0.4 ≤498 𝜋 ≤ 0.9 when the hidden component was small ( 0.01 ≤ 𝛿 ≤ 0.05), while it deteriorated at lower 499 values of 𝜋 regardless of hidden component size due to overfitting (Figure 4A). Meanwhile, factor 500 recovery of the hidden component 𝐿13 improved in ACMTF -R at 0.4 ≤ 𝜋 ≤ 0.9 only at 501 intermediate component sizes (0.03 ≤ 𝛿 ≤ 0.05) (Figure 4B). At 𝛿 > 0.05, ACMTF and ACMTF-R 502 were both able to identify 𝐿13, while at 𝛿 < 0.03 neither method could recover it (Supplementary 503 Figure 9). 504 The recovered shared subject mode loadings of the hidden component showed a consistent recovery 505 at small component sizes (0.01 -0.05), regardless of 𝜋 (Figure 4C). Meanwhile, the recovery of the 506 feature and time mode loadings improved at similar component sizes only for intermediate values 507 of 𝜋 (0.4-0.9) (Figure 4D-G). Taken together, the simulation revealed that a hidden component that 508 is roughly 1/20th the size of the other components can be successfully recovered by supervising with 509 ACMTF-R using intermediate values of 𝜋. 510 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 23 511 Figure 4. Improvement in f actor recovery for varying sizes of the hidden component 𝑳𝟏𝟑 using 512 ACMTF-R. One hundred randomly initialised eight -component ACMTF -R and ACMTF models were 513 fitted for each combination of hidden component size ( 𝛿) and tuning (𝜋), after which the lowest-loss 514 models were compared. The improvement in factor recovery is computed by subtracting ACMTF 515 performance from ACMTF-R performance in all panels. Improvement per setting for (A) overall factor 516 recovery, (B) factor recovery of the hidden component, (C) shared subject mode recovery of the hidden 517 component, (D) feature mode recovery of 𝑿(1) for the hidden component, (E) time mode recovery of 518 𝑿(1) for the hidden component, (F) feature mode recovery of 𝑿(3) for the hidden component, (G) time 519 mode recovery of 𝑿(3) for the hidden component. Improvement compared to ACMTF is shown in 520 green, while deterioration is shown in purple. A version of this figure without subtraction of ACMTF 521 performance is shown in Supplementary Figure 9. 522 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 24 Simulation 3: the benefit of supervision using 𝒚 is consistent for biologically realistic noise levels 523 We evaluated how the amount of noise on the independent data blocks 𝑿(1), 𝑿(2), and 𝑿(3) affected 524 the improvement in recovering the hidden component 𝐿13 using ACMTF-R compared to ACMTF. 525 This was done using the same three-tensor setup as in the first simulation, except that the amount 526 of noise on the independent data blocks was varied from 𝜂𝑋 = 0.0 to 𝜂𝑋 = 100. We expected the 527 benefit of supervision to increase as the amount of noise increased. As before, the tuning parameter 528 𝜋 was varied from 0.1 to 0.9 and one hundred randomly initialised eight-component ACMTF-R and 529 ACMTF models were created for each setting (Figure 5). 530 Input factor recovery across all components improved in ACMTF -R compared to ACMTF at most 531 𝜋 values (0.3-0.9) at 0-50% noise (Figure 5A). However, at 10-30% noise, the improvement was low 532 due to both ACMTF and ACMTF-R recovering the hidden component (Supplementary Figure 10). 533 Likewise, factor recovery of the hidden component 𝐿13 improved only at intermediate noise levels 534 (30-50%) and at intermediate 𝜋 values (0.6 -0.9) ( Figure 5B ). When no noise was present, the 535 improvement in factor recovery of the hidden component was zero because both methods were 536 unable to identify it (Supplementary Figure 10). 537 The recovered shared subject mode loadings of the hidden component showed a consistent recovery 538 at all noise levels and values of 𝜋 (Figure 5C). The recovery of the feature mode loadings in block 539 𝑿(1) improved only at intermediate noise levels (20-50%) and intermediate 𝜋 values (0.6-0.9), while 540 the recovery of the feature mode loadings in block 𝑿(3) were largely independent of the noise level 541 (Figure 5D,F). Surprisingly, the recovery of the time mode loadings in block 𝑿(1) deteriorated at 542 high noise levels (50+%), while the recovery of the time mode loadings in block 𝑿(3) improved at 543 those levels (Figure 5E,G). Taken together, this simulation revealed that factor recovery by ACMTF-544 R improves at biologically realistic noise levels (>30%). 545 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 25 546 Figure 5. The improvement in factor recovery by ACMTF -R under different amounts of noise on X. 547 One hundred randomly initialised eight-component ACMTF-R and ACMTF models were fitted for each 548 combination of noise on X ( 𝜂𝑋) and tuning parameter (𝜋), after which the lowest -loss models were 549 compared. The improvement in factor recovery is computed by subtracting ACMTF performance from 550 ACMTF-R performance in all panels. Improvement per setting for (A) overall factor recovery, (B) factor 551 recovery of the hidden component, (C) shared subject mode recovery of the hidden component, (D) 552 feature mode recovery of 𝑿(1) for the hidden component, (E) time mode recovery of 𝑿(1) for the 553 hidden component, (F) feature mode recovery of 𝑿(3) for the hidden component, (G) time mode 554 recovery of 𝑿(3) for the hidden component. Improvement compared to ACMTF is shown in green, 555 while deterioration is shown in purple. A version of this figure without subtraction of ACMTF 556 performance is shown in Supplementary Figure 10. 557 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 26 Simulation 4: the benefit of supervision using 𝒚 is largely independent of the amount of noise on 𝒚 558 We evaluated how the amount of noise on 𝒚 affects the improvement in recovering the hidden 559 component 𝐿13 using ACMTF-R compared to ACMTF. This was done using the same three-tensor 560 setup as in the first simulation, except that the amount of noise on 𝒚 was varied from 𝜂𝑦 = 0.0 to 561 𝜂𝑦 = 100. We expected the benefit of supervision to decrease as the amount of noise on 𝒚 increased. 562 As before, the tuning parameter 𝜋 was varied from 0.1 to 0.9 and one hundred randomly initialised 563 eight-component ACMTF-R and ACMTF models were created for each setting (Figure 6). 564 Input factor recovery across all components improved in ACMTF -R compared to ACMTF at 565 intermediate 𝜋 values (0.4-0.9) for most noise levels (0-100%) (Figure 6A). Likewise, factor recovery 566 of the hidden component 𝐿13 improved at intermediate 𝜋 values (0.5-0.9) for most noise levels (0-567 100%) (Figure 6B). At lower 𝜋 values, factor recovery in ACMTF-R deteriorated due to the strong 568 focus on predicting 𝒚, due to ACMTF-R not finding the hidden component (Supplementary Figure 569 11). At very high noise levels (100+%), neither ACMTF nor ACMTF -R were able to recover the 570 hidden component. 571 The recovered shared subject mode loadings of the hidden component showed a consistent recovery 572 at most noise levels (0-100%) and intermediate values of 𝜋 (0.4-0.9) (Figure 6C). Similarly, recovery 573 of the feature and time mode loadings was consistent at most noise levels (0-100%) and intermediate 574 values of 𝜋 (0.4-0.9) (Figure 6D-G). Taken together, the simulation revealed that recovery of the 575 hidden component by ACMTF-R is insensitive to biologically realistic noise levels on 𝒚. 576 577 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 27 578 Figure 6. The improvement in factor recovery by ACMTF -R under different amounts of noise on 𝒚. 579 One hundred randomly initialised eight-component ACMTF-R and ACMTF models were fitted for each 580 combination noise on 𝒚 (𝜂𝑦) and tuning parameter (𝜋), after which the lowest -loss models were 581 compared. The improvement in factor recovery is computed by subtracting ACMTF performance from 582 ACMTF-R performance in all panels. Improvement per setting for (A) overall factor recovery, (B) factor 583 recovery of the hidden component, (C) shared subject mode recovery of the hidden component, (D) 584 feature mode recovery of 𝑿(1) for the hidden component, (E) time mode recovery of 𝑿(1) for the 585 hidden component, (F) feature mode recovery of 𝑿(3) for the hidden component, (G) time mode 586 recovery of 𝑿(3) for the hidden component. Improvement compared to ACMTF is shown in green, 587 while deterioration is shown in purple. A version of this figure without subtraction of ACMTF 588 performance is shown in Supplementary Figure 11. 589 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 28 Application: ACMTF-R enhances multi-omics integration compared to ACMTF and NPLS 590 To demonstrate the advantage of ACMTF -R for integrating multi -omics data, we analysed a 591 longitudinal dataset comprising human milk (HM) metabolomics, HM microbiome and infant faecal 592 microbiome data, using pre-pregnancy BMI (ppBMI) as the outcome (𝒚). To contextuali se 593 performance, we evaluated the results against ACMTF , which focuses only on the CLD structure, 594 and N-way Partial Least Squares (NPLS) models per block [48], which focus only on prediction. We 595 expected ACMTF-R to identify more strongly maternal ppBMI associated components compared to 596 ACMTF. 597 A combined random initialisation and cross-validation procedure was used to select the number of 598 components for ACMTF and ACMTF-R (Methods). The following tuning parameter values were 599 investigated: 𝜋 = 0.75, 𝜋 = 0.80, 𝜋 = 0.85, and 𝜋 = 0.90. For each 𝜋 value, 10 randomly initialised 600 models were fitted and a 10 -fold cross-validation was performed across models with one to five 601 components (Supplementary Figure 12). This analysis revealed that one component was optimal for 602 all ACMTF-R settings. For ACMTF, a similar procedure indicated that two components were 603 appropriate, each explaining a different source of variation (Supplementary Figure 13). For NPLS, 604 cross-validation indicated that one component minimised RMSECV for the infant faecal and HM 605 microbiome data, while two components were required for the HM metabolomics data 606 (Supplementary Figure 14). 607 After selecting the number of components, 100 randomly initialised one -component ACMTF-R 608 models were fitted for each 𝜋 value. The model with the lowest overall loss in each case was selected, 609 and the variance explained was computed ( Table 1). To contextualise performance, the ACMTF 610 model explained 14.34%, 9.24%, and 7.03% of the variation in the infant faecal microbiome, HM 611 microbiome, and HM metabolomics data blocks, respectively, while explaining 9.28% of the 612 variation in 𝒚. Additionally, the NPLS models explained 5.78%, 5.41%, and 8.25% of the variation 613 in the infant faecal microbiome, HM microbiome, and HM metabolomics datasets, respectively, 614 while explaining 11.97%, 11.22%, and 35.95% of the variation in 𝒚. 615 The ACMTF-R model with 𝜋 = 0.90 explained 2.08-8.31% of the variation in the independent data 616 blocks, and 24.08% of the variation in 𝒚. At 𝜋 = 0.85, this dropped to between 1.78 -5.38% of the 617 variation in the independent data blocks, and 62.06% in 𝒚. At 𝜋 = 0.80 and 𝜋 = 0.75, the models 618 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 29 were most heavily tuned towards predicting 𝒚, yielding 1 -2% of explained variation in the 619 independent data blocks and capturing 90+% of the variation in 𝒚 – a clear indication of overfitting. 620 Table 1. Variance explained per dataset for the fitted NPLS, ACMTF and ACMTF-R models. 𝐹: the 621 number of components. 622 Dataset Model 𝐹 Variation explained (%) Infant faecal microbiome Mother milk microbiome Mother milk metabolome ppBMI (𝒚) Infant faecal microbiome NPLS 1 5.78 - - 11.97 HM microbiome NPLS 1 - 5.41 - 11.22 HM metabolome NPLS 2 - - 8.25 35.95 All ACMTF-R (𝜋 = 0.75) 1 1.52 1.15 1.24 96.55 ACMTF-R (𝜋 = 0.80) 1 2.45 1.91 1.41 89.97 ACMTF-R (𝜋 = 0.85) 1 5.38 4.35 1.78 62.06 ACMTF-R (𝜋 = 0.90) 1 8.31 6.69 2.08 24.08 ACMTF 2 14.34 9.24 7.03 9.28 623 MLR was performed on the subject mode loadings to test for associations with subject metadata, 624 including maternal pre-pregnancy BMI, infant growth at six months, birth mode, maternal secretor 625 (𝑆𝑒), and maternal Lewis (𝐿𝑒) blood group status (Table 2). Maternal 𝑆𝑒 and 𝐿𝑒 status determine the 626 presence of specific fucosylated human milk oligosaccharides, which may influence the composition 627 of the infant faecal microbiome. 628 The subject mode loadings of the NPLS models all contained at least one component that was 629 associated with maternal pre-pregnancy BMI. While both components of the ACMTF model were 630 associated with both ppBMI and maternal secretor status , the first component was most strongly 631 associated with maternal secretor status and the second was most strongly associated with ppBMI . 632 This was not the case for the ACMTF-R models with 𝜋 = 0.90 and 𝜋 = 0.85, which both contained 633 subject mode loadings associated exclusively with ppBMI after correction . The subject mode 634 loadings of the ACMTF-R models with lower values of 𝜋 were strongly associated with ppBMI due 635 to overfitting. 636 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 30 Due to the extraction of a component exclusively related to maternal ppBMI, as well as the trade -637 off between explaining the independent data and predicting 𝒚, the ACMTF-R model using 𝜋 = 0.90 638 was selected for further interpretation. 639 Table 2. MLR results for the subject mode loadings of the created multi -way models . ppBMI: 640 maternal pre-pregnancy body mass index. 𝑆𝑒: maternal secretor status. 𝐿𝑒: maternal Lewis status. 641 WHZ: infant weight -of-height z-score at 6 months. Subject mode loadings were normalised to 642 norm 1 prior to MLR fitting for comparability of the results between the models. Values indicate 643 MLR coefficients times one thousand. Stars indicate the original p -values from the MLR model (*: 644 𝑝 < 0.05, **: 𝑝 < 0.01, ***: 𝑝 < 0.001), while hats indicate Benjamini-Hochberg corrected p-values 645 (^: 𝑞 < 0.05, ^^: 𝑞 < 0.01, ^^^: 𝑞 < 0.001). Coefficients with a significant p-value after correction 646 are made bold. The correlation matrix of the metadata is reported in Supplementary Table 2. 647 Dataset Model F ppBMI Birth mode Se Le WHZ Infant faecal microbiome NPLS 1 5.9***, ^^^ 26.2 18.7 -67.4 1.5 HM microbiome NPLS 1 5.3***, ^^^ 22.5 10.1 -25.4 -0.6 HM metabolome NPLS 1 8.7***, ^^^ 11.1 19.7 -78.3 -12.6*, ^ 2 -0.2 -12.8 -111***, ^^^ 43.0 7.0 All ACMTF-R (𝜋 = 0.75) 1 15.0***, ^^^ 6.7 3.0* -18.5 -0.9 ACMTF-R (𝜋 = 0.80) 1 -14.7***, ^^^ -10.8 -5.3* 31.7 1.4 ACMTF-R (𝜋 = 0.85) 1 -12.5***, ^^^ -18.1 -12.3* 59.6 2.4 ACMTF-R (𝜋 = 0.90) 1 -8.3***, ^^^ -21.5 -21.8 79.4 2.8 ACMTF 1 -2.9*, ^ -2.0 -74.8***, ^^^ 65.5 3.2 2 -4.4***, ^^ -40.1 -41.3*, ^ 72.3 1.0 648 The 𝚲-matrix of the ACMTF-R model using 𝜋 = 0.90 was examined to interpret the common, local, 649 and distinct (CLD) structure of the component. In this model, 𝜆 values of the only available 650 component were equal to 0. 34, 0.29, and 0.15, for the infant faecal microbiome, HM microbiome, 651 and HM metabolomics, respectively. Based on the relative magnitudes of the 𝜆 values, the 652 component appears to be primarily local to the infant faecal and HM microbiome blocks, with a 653 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 31 minor contribution of the HM metabolomics block. This suggests that ppBMI -associated variation 654 was mainly present in the microbiomes of mother and infant. 655 In the infant faecal microbiome, genera associated with higher ppBMI included Bifidobacterium 656 bifidum, Escherichia coli, Staphylococcus epidermidis, and B. breve (Figure 7A). Genera enriched 657 in infants of low-ppBMI mothers included Parabacteroides, Enterococcus and Streptococcus. In the 658 HM microbiome, higher ppBMI was associated with increased abundance of Staphylococcus and 659 Streptococcus spp., while genera associated with lower ppBMI included various taxa such as 660 Clostridium sensu stricto, Collinsella, and Lactobacillus (Figure 7C). In the HM metabolome, higher 661 ppBMI was associated with increased levels of many unknown HMOs, caffeine, and LDFT ( Figure 662 7E). Metabolites with increased levels in low -ppBMI mothers included several amino acids and 663 related compounds such as 2-aminobutyrates, aspartate, isoleucine, and methionine. 664 Due to the pre -processing of the data, the largest absolute time loading in the ACMTF -R model 665 indicates the time points at which inter -subject differences related to ppBMI are greatest ( Figure 666 7B,D,F). Overall, the time loadings were relatively flat, suggesting that ppBMI-associated variation 667 remained stable throughout the sampling period. In the HM microbiome, however, these 668 differences peaked at day 30 and declined thereafter. 669 Taken together, ACMTF -R was able to identify variation that was strongly and exclusively 670 associated with maternal ppBMI. Additionally, it revealed that this variation was primarily local to 671 the infant faecal and HM microbiome, with a limited contribution from HM metabolome. This 672 highlights the value of incorporating supervision to isolate outcome -relevant signals in complex 673 multi-omics data. 674 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 32 675 Figure 7. Feature and time mode loadings of the selected ACMTF-R model. Feature mode loadings 676 of the selected ACMTF -R model using 𝜋 = 0.90 describing (A) the infant faecal microbiome, (C) 677 the HM microbiome, and (E) the HM metabolomics datasets. Time loadings of the same model are 678 shown in the same order in (B,D,F). Signs of the feature mode loadings were adjusted post-hoc such 679 that positive values consistently corresponded to higher ppBMI and negative values to lower ppBMI. 680 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 33

681 In this study, we presented ACMTF-Regression (ACMTF-R) as an extension of Advanced Coupled 682 Matrix and Tensor Factorization (ACMTF) that incorporates a regression term to model common, 683 local and distinct (CLD) variation related to a dependent variable. Through simulations, we 684 demonstrated its ability to capture a hidden component, its robustness to noise, and its flexibility in 685 tuning between data exploration and predicting. Our application demonstrated how ACMTF -R 686 finds novel relationships for the effect of maternal ppBMI on the HM metabolome, HM microbiome 687 and infant faecal microbiome. 688 ACMTF-R extends the capabilities of ACMTF by allowing explicit modelling of outcome-associated 689 variation, bridging the gap between exploratory data integration and supervised learning [10]. 690 Compared to NPLS, which is designed for predictive modelling but lacks the ability to separate 691 common, local, and distinct variation across multiple blocks, ACMTF -R provides a more nuanced 692 view of multi-way data structures [19]. Unlike classical ACMTF, which only uncovers underlying 693 data patterns, ACMTF -R allows for targeted analysis of biological variation of interest [33]. The 694 ability to tune between the two makes it particularly useful for cases where a balance between 695 exploration and prediction is necessary. 696 Our simulations revealed that ACMTF-R reliably recovers the underlying CLD-structure of multi-697 way datasets across a range of settings. We found that the tuning parameter (𝜋) governs the trade-698 off between reconstructing the input data and predicting the outcome variable. While 𝜋 = 1 699 corresponds to regular ACMTF, decreasing 𝜋 enhances the identification of factors associated with 700 the dependent variable, even revealing weak structures that would otherwise remain undetected. 701 However, tuning too far towards 𝒚 (𝜋 < 0.3) leads to overfitting and an inaccurate reconstruction. 702 The observation that intermediate values of 𝜋 (0.6-0.9) lead to recovery of the hidden local 703 component is largely in agreement with PCovR and MCovR [20, 22]. 704 Additionally, the simulations revealed that the hidden component needed to be sufficiently small 705 (𝛿 < 0.05) for ACMTF to fail in its recovery. However, this finding is largely dependent on the 706 relative sizes of the other components, which we have kept at Frobenius norm 1 in our simulations. 707 Further research is needed to investigate the recovery of a hidden component using ACMTF and 708 ACMTF-R in other scenarios. 709 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 34 Finally, the simulations showed that recovery of the hidden component using ACMTF-R was most 710 accurate when a biologically realistic amount of noise was present in the data ( 𝜂𝑋 ≤ 0.5, 𝜂𝑦 < 1). 711 The observation that neither ACMTF nor ACMTF-R can identify the hidden component in the no-712 noise case is surprising. This may occur when the noise has a smoothing effect on the loss space, 713 resulting in the global minimum being easier to identify, or reduces the collinearity of the loadings 714 in a statistical phenomenon comparable to ridge regression [55]. While the noise level for a real 715 dataset is often not known, it is reassuring that ACMTF -R can reliably capture the 𝒚-associated 716 subject mode loadings for most tested noise levels (Supplementary Figure 10-11). 717 Application to a longitudinal multi-omics dataset integrating human milk microbiome, human milk 718 metabolome, and infant gut microbiome data further highlighted the utility of ACMTF -R by 719 recovering a n exclusively ppBMI -associated component [48]. In the infant faecal microbiome, 720 elevated abundances of Enterococcus sp., E. coli, Clostridium sp., and Staphylococcus epidermidis 721 in infants born to high-ppBMI mothers are consistent with prior findings linking maternal obesity 722 to pro -inflammatory microbial profiles [56, 57] . In contrast, the increased abundance of 723 Bifidobacterium breve and Bifidobacterium bifidum in the same group contradicts earlier studies, 724 which generally associate these taxa with leaner maternal phenotypes [56, 57]. 725 Interestingly, the genera enriched in infants from low-ppBMI mothers represent a potential novel 726 finding, as comparable reports are currently lacking and warrant further investigation. In the HM 727 microbiome, enrichment of Streptococcus and Staphylococcus spp. in high -ppBMI mothers is in 728 agreement with previous research [58–60]. However, comparable data for taxa enriched in low -729 ppBMI mothers were not available for direct comparison. 730 In the HM metabolomics data, elevated levels of lactate and multiple unidentified HMOs were 731 observed in high-ppBMI mothers. While the relationship between ppBMI and HMO concentrations 732 has been investigated in several studies, the results remain inconsistent [61, 62]. One study did 733 report higher levels of a lactate derivative in the milk of obese mothers [61], suggesting a possible 734 link to altered maternal energy metabolism. Conversely, the observation of increased 735 concentrations of several amino acids in low-ppBMI mothers is supported by previous findings [61], 736 as well as by recent evidence showing enrichment of amino acid biosynthesis pathways in this group 737 [63]. These results show that ACMTF-R provides additional biological insights compared to ACMTF 738 and N -PLS, effectively disentangling common, local and distinct variation related to maternal 739 ppBMI. 740 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 35 While ACMTF-R demonstrates strong performance in identifying biologically meaningful patterns, 741 several challenges remain. The choice of the tuning parameter 𝜋 is dataset-dependent, requiring 742 careful selection through cross -validation or heuristic approaches. Additionally, in multi -omics 743 applications, different data blocks may contain varying amounts of biologically relevant 744 information. Allowing block -specific 𝜋 values could improve model flexibility, enabling certain 745 datasets to contribute more strongly to the dependent variable while ensuring others remain focused 746 on data structure. Future research should explore whether adaptive weighting strategies or block -747 specific hyperparameter tuning can further improve performance. 748 For ACMTF and ACMTF-R, the single guiding metric for selecting the number of components was 749 often the stability of the decomposition through the factor match score (FMS) [15, 54]. This paper 750 has expanded upon this approach by formalising FMS for a cross -validation and random 751 initialisation scheme, as well as developing the degeneracy score diagnostic to identify 752 overfactorisation. However, this approach is computationally intensive for the simulations. Future 753 work should investigate if this approach could be accelerated through the Limited -memory 754 Broyden-Fletcher-Goldfarb-Shanno [40–42] and AO-ADMM algorithms [64–66]. 755 While Factor Match Score (FMS) has been a useful diagnostic for model selection, it suffers from 756 the phenomenon of all blocks contributing to all components in ACMTF outlined above. Due to the 757 way the structural model of ACMTF is implemented, 𝐹 trilinear components are fitted per data 758 block, even when the rank of that block is lower than 𝐹 [33] This causes the creation of superfluous 759 feature and time mode loading vectors that contribute to the model due to having non-zero 𝜆 values. 760 As a result, the FMS is negatively impacted by considering these components in the calculation. 761 This problem may be avoided by more strongly penalising the 𝚲-matrix through replacement of the 762 L1 penalty term by a smooth L0 penalty, as presented in Relaxed ACMTF (RACMTF) [67, 68]. 763 Future research should investigate if this approach yields more stable ACMTF and ACMTF -R 764 models for multi-omics data. 765 The multi-omics application presented in this study highlights a critical shortcoming of ACMTF in 766 the identification of common, local, and distinct (CLD) sources of variation by revealing that every 767 data block contributes to all identified components. However, the relative magnitude of the 768 contributing blocks captured by 𝜆 still offer a useful interpretation. For example, a 𝚲-matrix 769 threshold could be applied post-hoc to identify CLD-structures [33]. Another approach could be to 770 pre-define the expected number of shared and distinct components, as suggested by Common and 771 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 36 Discriminative subspace Non-negative Tensor Factorization (CDNTF) [69]. Future research should 772 investigate if either approach or a combination of both might yield more interpretable ACMTF 773 models for multi-omics data. 774 Finally, treating 𝒚 as an additional data block for ACMTF could offer a different approach to 775 incorporating outcome variation into ACMTF -R. By embedding 𝒚 within the multi -way 776 decomposition framework, rather than optimizing it separately, the model might be able to better 777 capture latent structures linking 𝒚 to the input datasets. However, this approach would lose the 778 predictive aspect of ACMTF-R. Future work should focus on the comparability between ACMTF-R 779 and adding 𝒚 as a block to ACMTF. 780

781 Our results demonstrate that ACMTF -R is a powerful framework for multi -way data integration 782 that extends ACMTF by incorporating outcome-associated variation. Through simulations, we 783 showed that ACMTF -R can recover a small hidden component, is robust to noise, and provides 784 flexible tuning between data exploration and outcome prediction. Its application to multi -omics 785 data highlights how ACMTF-R finds novel relationships for the effect of maternal ppBMI on the 786 HM metabolome, HM microbiome and infant faecal microbiome. 787 Ethics approval and consent to participate 788 Not applicable. 789 Consent for publication 790 Not applicable. 791 Availability of the data and materials 792 The package used during the current study is available in the GitHub repository 793 https://github.com/GRvanderPloeg/CMTFtoolbox. The datasets generated and/or analysed during the 794 current study are available in the GitHub repository 795 https://github.com/GRvanderPloeg/CMTF_toolbox_project. 796 Competing interests 797 The authors declare that they have no competing interests. 798 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 37 Funding 799 GRvdP was funded by a grant from the University of Amsterdam, Research Priority Area on 800 Personal Microbiome Health. 801 Fred White was funded by a grant from the University of Amsterdam, Data Science Centre. 802 Authors’ contributions 803 G.R. van der Ploeg: conceptualization, visualization, formal analysis, methodology, programming, 804 writing – original draft; F.T.G. White: conceptualization, methodology, programming, writing – 805 original draft ; R.R. Jakobsen: programming, writing – review & editing ; J.A. Westerhuis: 806 conceptualization, methodology, supervision, formal analysis, writing – original draft; A. Heintz -807 Buschart: conceptualization, formal analysis, supervision, project administration, writing – original 808 draft; A.K. Smilde: conceptualization, methodology, supervision, funding acquisition, writing – 809 review & editing. 810

811 We want to thank Daniël Rademaker (RU/UvA) for their useful discussions and suggestions. 812

813 1. Capobianco E. Ten challenges for systems medicine. Front Genet. 2012;3:193. 814 2. Davis-Turak J, Courtney SM, Hazard ES, Glen WB, Da Silveira WA, Wesselman T, et al. Genomics 815 pipelines and data integration: challenges and opportunities in the research setting. Expert Rev Mol 816 Diagn. 2017;17:225–37. https://doi.org/10.1080/14737159.2017.1282822. 817 3. Zhao Z, Jin VX, Huang Y, Guda C, Ruan J. Frontiers in integrative genomics and translational 818 bioinformatics. BioMed Res Int. 2015;2015. 819 4. Bordbar A. Interpreting the deluge of omics data: new approaches offer new possibilities. Blood 820 Transfus. 2017;15:189. 821 5. Chen C, McGarvey PB, Huang H, Wu CH. Protein Bioinformatics Infrastructure for the Integration 822 and Analysis of Multiple High -Throughput “omics” Data. Adv Bioinforma. 2010;2010:1 –19. 823 https://doi.org/10.1155/2010/423589. 824 6. Misra BB, Langefeld C, Olivier M, Cox LA. Integrated omics: tools, advances and future approaches. 825 J Mol Endocrinol. 2019;62:R21–45. 826 7. Vitorino R. Transforming Clinical Research: The Power of High -Throughput Omics Integration. 827 Proteomes. 2024;12:25. https://doi.org/10.3390/proteomes12030025. 828 8. van der Kloet FM, Sebastián-León P, Conesa A, Smilde AK, Westerhuis JA. Separating common from 829 distinctive variation. BMC Bioinformatics. 2016;17:S195. https://doi.org/10.1186/s12859-016-1037-2. 830 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 38 9. Måge I, Smilde AK, van der Kloet FM. Performance of methods that separate common and distinct 831 variation in multiple data blocks. J Chemom. 2019;33:e3085. https://doi.org/10.1002/cem.3085. 832 10. Acar E, Papalexakis EE, Gürdeniz G, Rasmussen MA, Lawaetz AJ, Nilsson M, et al. Structure -833 revealing data fusion. BMC Bioinformatics. 2014;15:239. https://doi.org/10.1186/1471-2105-15-239. 834 11. Tan CS, Salim A, Ploner A, Lehtiö J, Chia KS, Pawitan Y. Correlating gene and protein expression data 835 using Correlated Factor Analysis. BMC Bioinformatics. 2009;10:272. https://doi.org/10.1186/1471 -836 2105-10-272. 837 12. Xiao X, Moreno -Moral A, Rotival M, Bottolo L, Petretto E. Multi -tissue analysis of co -expression 838 networks by higher -order generalized singular value decomposition identifies functionally coherent 839 transcriptional modules. PLoS Genet. 2014;10:e1004006. 840 13. Berger JA, Hautaniemi S, Mitra SK, Astola J. Jointly analyzing gene expression and copy number 841 data in breast cancer using data reduction models. IEEE/ACM Trans Comput Biol Bioinform. 2006;3:2–842 16. 843 14. Song Y, Westerhuis JA, Smilde AK. Separating common (global and local) and distinct variation in 844 multiple mixed types data sets. J Chemom. 2020;34:e3197. https://doi.org/10.1002/cem.3197. 845 15. Acar E, Kolda TG, Dunlavy DM. All -at-once Optimization for Coupled Matrix and Tensor 846 Factorizations. 2011. 847 16. Acar E, Bro R, Smilde A. Data Fusion in Metabolomics Using Coupled Matrix and Tensor 848 Factorizations. Proc IEEE. 2015;103:1602. https://doi.org/10.1109/JPROC.2015.2438719. 849 17. Acar E, Bro R, Smilde AK. Data Fusion in Metabolomics Using Coupled Matrix and Tensor 850 Factorizations. Proc IEEE. 2015;103:1602–20. https://doi.org/10.1109/JPROC.2015.2438719. 851 18. Acar E, Dunlavy DM, Kolda TG, Mørup M. Scalable tensor factorizations for incomplete data. 852 Chemom Intell Lab Syst. 2011;106:41–56. https://doi.org/10.1016/j.chemolab.2010.08.004. 853 19. Bro R. Multiway calibration. Multilinear PLS. J Chemom. 1996;10:47 –61. 854 https://doi.org/10.1002/(SICI)1099-128X(199601)10:13.0.CO;2-C. 855 20. Vervloet M, Kiers HAL, Noortgate WV den, Ceulemans E. PCovR: An R Package for Principal 856 Covariates Regression. J Stat Softw. 2015;65:1–14. https://doi.org/10.18637/jss.v065.i08. 857 21. De Jong S, Kiers HA. Principal covariates regression: part I. Theory. Chemom Intell Lab Syst. 858 1992;14:155–64. 859 22. Smilde AK, Kiers HAL. Multiway covariates regression models. J Chemom. 1999;13:31 –48. 860 https://doi.org/10.1002/(SICI)1099-128X(199901/02)13:13.0.CO;2-P. 861 23. Kiers HAL. Towards a standardized notation and terminology in multiway analysis. J Chemom. 862 2000;14:105–22. https://doi.org/10.1002/1099-128X(200005/06)14:33.0.CO;2-I. 863 24. Kolda TG. Multilinear operators for higher -order decompositions. Sandia National Laboratories 864 (SNL), Albuquerque, NM, and Livermore, CA …; 2006. 865 25. De Lathauwer L, De Moor B, Vandewalle J. A Multilinear Singular Value Decomposition. SIAM J 866 Matrix Anal Appl. 2000;21:1253–78. https://doi.org/10.1137/S0895479896305696. 867 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 39 26. Styan GP. Hadamard products and multivariate statistical analysis. Linear Algebra Its Appl. 868 1973;6:217–40. 869 27. McDonald RP. A simple comprehensive model for the analysis of covariance structures: Some 870 remarks on applications. Br J Math Stat Psychol. 1980;33:161 –83. https://doi.org/10.1111/j.2044 -871 8317.1980.tb00606.x. 872 28. Rao CR. GENERALIZED INVERSE OF A MATRIX AND ITS APPLICATIONS. In: Le Cam LM, Neyman J, 873 Scott EL, editors. Theory of Statistics. University of California Press; 1972. p. 601 –20. 874 https://doi.org/10.1525/9780520325883-032. 875 29. Ben -Israel A, Greville TN. Generalized inverses: theory and applications. Springer Science & 876 Business Media; 2006. 877 30. Carroll JD, Chang J -J. Analysis of individual differences in multidimensional scaling via an N -way 878 generalization of “Eckart-Young” decomposition. Psychometrika. 1970;35:283–319. 879 31. Harshman RA. FOUNDATIONS OF THE PARAFAC PROCEDURE: MODELS AND CONDITIONS FOR AN 880 “EXPLANATORY” MULTIMODAL FACTOR ANALYSIS. :84. 881 32. Bro R. MULTI-WAY ANALYSIS IN THE FOOD INDUSTRY. 882 33. Acar E, Papalexakis EE, Gürdeniz G, Rasmussen MA, Lawaetz AJ, Nilsson M, et al. Structure -883 revealing data fusion. BMC Bioinformatics. 2014;15:239. https://doi.org/10.1186/1471-2105-15-239. 884 34. Lee S-I, Lee H, Abbeel P, Ng AY. Efficient l1 regularized logistic regression. In: Aaai. 2006. p. 401–8. 885 35. Bertrand D, Qannari EM, Vigneau E. Latent root regression analysis: an alternative method to PLS. 886 Chemom Intell Lab Syst. 2001;58:227–34. https://doi.org/10.1016/S0169-7439(01)00161-7. 887 36. Webster JT, Gunst RF, Mason RL. Latent Root Regression Analysis. Technometrics. 1974;16:513 –888 22. https://doi.org/10.1080/00401706.1974.10489232. 889 37. Harshman RA. Foundations of the PARAFAC procedure: Models and conditions for an" explanatory" 890 multimodal factor analysis. 1970. 891 38. Bro R. PARAFAC. Tutorial and applications. Chemom Intell Lab Syst. 1997;:23. 892 39. Melville J. mize: Unconstrained Numerical Optimization Algorithms. 2019. 893 40. Liu DC, Nocedal J. On the limited memory BFGS method for large scale optimization. Math Program. 894 1989;45:503–28. https://doi.org/10.1007/BF01589116. 895 41. Malouf R. A comparison of algorithms for maximum entropy parameter estimation. In: proceedings 896 of the 6th conference on Natural language learning - Volume 20. USA: Association for Computational 897 Linguistics; 2002. p. 1–7. https://doi.org/10.3115/1118853.1118871. 898 42. Andrew G, Gao J. Scalable training of L1-regularized log-linear models. In: Proceedings of the 24th 899 international conference on Machine learning. New York, NY, USA: Association for Computing 900 Machinery; 2007. p. 33–40. https://doi.org/10.1145/1273496.1273501. 901 43. Burt C. The factorial study of temperamental traits. Br J Psychol. 1948. 902 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 40 44. Tucker LR. A method for synthesis of factor analysis studies. Educational Testing Service Princeton, 903 NJ; 1951. 904 45. Lorenzo-Seva U, ten Berge JMF. Tucker’s congruence coefficient as a meaningful index of factor 905 similarity. Methodol Eur J Res Methods Behav Soc Sci. 2006;2:57 –64. https://doi.org/10.1027/1614-906 2241.2.2.57. 907 46. Kuhn HW. The Hungarian method for the assignment problem. Nav Res Logist Q. 1955;2:83 –97. 908 https://doi.org/10.1002/nav.3800020109. 909 47. Kuhn HW. Variants of the hungarian method for assignment problems. Nav Res Logist Q. 910 1956;3:253–8. https://doi.org/10.1002/nav.3800030404. 911 48. Jakobsen R, Ploeg GR, Sundekilde U, Astono J, Poulsen K, Fuglsang J, et al. Supervised Modelling of 912 Longitudinal Human Milk and Infant Gut Microbiome Reveal Maternal Pre -Pregnancy BMI and Early 913 Life Growth Interactions. 2025. https://doi.org/10.21203/rs.3.rs-6244750/v1. 914 49. R Core Team R. R: A language and environment for statistical computing. 2013. 915 50. Gloor GB, Macklaim JM, Pawlowsky-Glahn V, Egozcue JJ. Microbiome Datasets Are Compositional: 916 And This Is Not Optional. Front Microbiol. 2017;8:2224. https://doi.org/10.3389/fmicb.2017.02224. 917 51. Aitchison J. The Statistical Analysis of Compositional Data. J R Stat Soc Ser B Methodol. 918 1982;44:139–60. https://doi.org/10.1111/j.2517-6161.1982.tb01195.x. 919 52. Bro R, Smilde AK. Centering and scaling in component analysis. J Chemom. 2003;17:16 –33. 920 https://doi.org/10.1002/cem.773. 921 53. van der Ploeg GR, Westerhuis JA, Heintz-Buschart A, Smilde AK. parafac4microbiome: Exploratory 922 analysis of longitudinal microbiome data using Parallel Factor Analysis. 2024. 923 https://doi.org/10.1101/2024.05.02.592191. 924 54. Li L, Yan S, Horner D, Rasmussen MA, Smilde AK, Acar E. Revealing static and dynamic biomarkers 925 from postprandial metabolomics data through coupled matrix and tensor factorizations. 2024. 926 55. Schott JR. Matrix analysis for statistics. John Wiley & Sons; 2016. 927 56. Singh SB, Madan J, Coker M, Hoen A, Baker ER, Karagas MR, et al. Associations of maternal pre -928 pregnancy BMI and gestational weight gain with the infant gut microbiome differ according to delivery 929 mode. Int J Obes 2005. 2019;44:23. 930 57. Collado MC, Isolauri E, Laitinen K, Salminen S. Effect of mother’s weight on infant’s microbiota 931 acquisition, composition, and activity during early infancy: a prospective follow -up study initiated in 932 early pregnancy. Am J Clin Nutr. 2010;92:1023–30. 933 58. Lundgren SN, Madan JC, Karagas MR, Morrison HG, Hoen AG, Christensen BC. Microbial 934 communities in human milk relate to measures of maternal weight. Front Microbiol. 2019;10:2886. 935 59. Cabrera -Rubio R, Collado MC, Laitinen K, Salminen S, Isolauri E, Mira A. The human milk 936 microbiome changes over lactation and is shaped by maternal weight and mode of delivery1234. Am 937 J Clin Nutr. 2012;96:544–51. https://doi.org/10.3945/ajcn.112.037382. 938 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint 41 60. Cortés-Macías E, Selma-Royo M, Rio-Aige K, Bäuerl C, Rodríguez-Lagunas MJ, Martínez-Costa C, et 939 al. Distinct breast milk microbiota, cytokine, and adipokine profiles are associated with infant growth 940 at 12 months: an in vitro host –microbe interaction mechanistic approach. Food Funct. 2023;14:148 –941 59. 942 61. Isganaitis E, Venditti S, Matthews TJ, Lerin C, Demerath EW, Fields DA. Maternal obesity and the 943 human milk metabolome: associations with infant body composition and postnatal weight gain. Am J 944 Clin Nutr. 2019;110:111–20. 945 62. Bardanzellu F, Puddu M, Peroni DG, Fanos V. The Human Breast Milk Metabolome in Overweight 946 and Obese Mothers. Front Immunol. 2020;11:1533. https://doi.org/10.3389/fimmu.2020.01533. 947 63. Sivalogan K, Liang D, Accardi C, Diaz-Artiga A, Hu X, Mollinedo E, et al. Human Milk Composition Is 948 Associated with Maternal Body Mass Index in a Cross-Sectional, Untargeted Metabolomics Analysis of 949 Human Milk from Guatemalan Mothers. Curr Dev Nutr. 2024;8:102144. 950 https://doi.org/10.1016/j.cdnut.2024.102144. 951 64. Schenker C, Cohen JE, Acar E. A Flexible Optimization Framework for Regularized Matrix -Tensor 952 Factorizations With Linear Couplings. IEEE J Sel Top Signal Process. 2021;15:506 –21. 953 https://doi.org/10.1109/jstsp.2020.3045848. 954 65. Schenker C, Cohen JE, Acar E. An Optimization Framework for Regularized Linearly Coupled Matrix-955 Tensor Factorization. In: 2020 28th European Signal Processing Conference (EUSIPCO). Amsterdam, 956 Netherlands: IEEE; 2021. p. 985–9. https://doi.org/10.23919/eusipco47968.2020.9287459. 957 66. Roald M. MatCoupLy: Learning coupled matrix factorizations with Python. SoftwareX. 958 2023;21:101292. https://doi.org/10.1016/j.softx.2022.101292. 959 67. Mohimani H, Babaie-Zadeh M, Jutten C. A fast approach for overcomplete sparse decomposition 960 based on smoothed $\backslashell^{0} $ norm. IEEE Trans Signal Process. 2008;57:289–301. 961 68. Rivet B, Duda M, Guérin-Dugué A, Jutten C, Comon P. Multimodal approach to estimate the ocular 962 movements during EEG recordings: A coupled tensor factorization method. In: 2015 37th Annual 963 International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC). 2015. p. 964 6983–6. https://doi.org/10.1109/EMBC.2015.7319999. 965 69. Liu W, Chan J, Bailey J, Leckie C, Ramamohanarao K. Mining Labelled Tensors by Discovering both 966 their Common and Discriminative Subspaces. In: Proceedings of the 2013 SIAM International 967 Conference on Data Mining (SDM). Society for Industrial and Applied Mathematics; 2013. p. 614 –22. 968 https://doi.org/10.1137/1.9781611972832.68. 969 970 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted July 31, 2025. ; https://doi.org/10.1101/2025.07.28.667162doi: bioRxiv preprint