\section{Introduction}\label{sec:intro}

Modern experiments are increasingly likely to produce complex data derived from multiple sources. One common example is a single group of $n$ \textit{objects} (i.e. cases, observations, patients) being observed across $K$ different \textit{views} or \textit{data blocks}, each with their own sets of $d_k,\ k=1,\ldots,K,$ \textit{traits} (i.e. variables, features, descriptors) and measurement methodologies. We call data collected and organized this way \textit{multi-block data}. Simple data analysis approaches would concatenate multiple data blocks into a single data matrix of $n$ data objects and $d_1+\cdots+d_K$ traits. However, these approaches ignore the often important relationships between the views and obscure insights on which information comes from which data block. We specifically aim to address this challenge by considering the data blocks as separate units and searching for \textit{shared structure} between them. Shared structure can be defined in several ways. Here and in some other contemporary approaches seen in Section~\ref{sec:litreview} it means common modes of variation modeled by low-rank matrices. More mathematical details of this modeling choice can be found in Section~\ref{sec:methodology}. Our proposed method, \textit{Data Integration Via Analysis of Subspaces} (DIVAS), incorporates state-of-the-art advances in matrix perturbation theory and optimization to provide insights about both shared and partially-shared joint structure between several data blocks.

DIVAS gives a novel approach for finding structure in a multi-block data set based on searching for shared subspaces between different collections of data blocks. The data blocks are intrinsically linked by having the same \textit{trait space} $\mathbb{R}^n$, so we primarily search for shared subspaces within the trait space. In contrast to other approaches, angles form the foundation of our analysis of the relationships between these subspaces. In particular, \textit{principal angles} are the measure of choice of proximity between subspaces. These angles undergird a rigorous framework of inference that provides relevant statistical significance determinations for the chosen subspaces.

Subspaces of trait space have corresponding induced subspaces of \textit{object space} $\mathbb{R}^{d_k}$ for each data block. Combined together, these subspaces can be decomposed into \textit{modes of variation}. A mode of variation is a rank 1 matrix formed from the outer product of two vectors: one in object space and one in trait space. In the terminology of Principal Components Analysis (PCA) these would be a loadings (direction) vector and a scores vector, respectively. Considering subspaces in terms of modes of variation is particularly useful for visualization. The scores vectors demonstrate relationships between the data objects, and the loadings vectors provide information about which traits are driving the variation. The ultimate result of a DIVAS exploration of a data set is a set of modes of variation for each data block associated with each block collection.

The rest of the paper proceeds as follows. The remainder of Section~\ref{sec:intro} continues with a motivating data example for DIVAS as well as background information on data integration in general. Section~\ref{sec:methodology} details DIVAS methodology and demonstrates the method's performance on a synthetic data set. Section~\ref{sec:data} contains a case study on twentieth century mortality using DIVAS. Section~\ref{sec:conclusion} summarizes some brief conclusions. The Appendices~\ref{sec:mp}, \ref{sec:paa}, and \ref{sec:noise} provide reference materials on random matrix theory, principal angle analysis, and residual matrix estimation respectively. MATLAB code is available for download at \url{https://github.com/jbprothero/DIVAS2021}.

\subsection{Motivating Example: Cancer Genomics}

One of the examples that motivates the development of DIVAS is a four-block data set containing different views of omics data from $n=616$ breast cancer patients from The Cancer Genome Atlas (TCGA) \citep{tcga}. We have a gene expression (GE) data block containing 16615 gene traits, a gene copy number (CN) data block with 24174 traits, a protein expression (RPPA) data block containing 187 protein traits, and a 0-1 mutation detection (Mut) block containing traits for 128 genes. Each patient is labeled as one of four breast cancer subtypes: Luminal A, Luminal B, Basal, or Her2-enriched.

We wish to obtain the entire hierarchy of joint structure among the four data blocks. Of particular interest are the four-way joint structure, three-way partially-shared joint structures, and partially-shared structure involving proteins or mutations. The biological understanding of the protein production pathway suggests that after accounting for the gene expression and gene copy number there should be no additional variation shared between mutations and proteins. Once all joint structure is cataloged, further conclusions can be drawn from loadings of the joint modes of variation. For example, the loadings from the gene expression data block would reveal which genes are involved in a certain cancer subtype if one of the joint modes of variation discriminates that subtype from the others. Note that DIVAS is an unsupervised method that does not make use of the class labels. Development of a supervised version of DIVAS remains an interesting open problem.  

Figure~\ref{fig:tcgascoresdiag} shows the DIVAS decomposition and angle diagnostics for the joint scores vectors for this data set. Detailed descriptions of the information plotted in the figure can be found in Section~\ref{sec:graphics}. For the purposes of this discussion we are only concerned with the ranks of partially shared joint structure displayed as numbers at the top of each colored box. DIVAS finds a single shared component between all four data blocks, a seven-dimensional subspace shared between all data blocks besides mutation, and lots of shared structure between pairs of data blocks not involving mutation. This result aligns well with biological expectations, particularly the large amount of structure shared uniquely between gene expression and copy number. 


\begin{figure}[h]
\includegraphics[width=\textwidth]{tcgaJointDiagnosticScores-11-18.png}
\caption{Joint structure breakdown and diagnostics for the breast cancer omics score vectors.}\label{fig:tcgascoresdiag}
\end{figure}


Figure~\ref{fig:tcga4way3wayscores} displays a scores scatter plot matrix of directions from the four-way joint and three-way joint components of the data. In these plots, each point corresponds to a single data object. Each cancer subtype is shown with a different color and point symbol: basal with red triangles, luminal A with blue asterisks, luminal B with cyan x'es, and Her2-enriched with magenta pluses. On-diagonal plots show the coefficients of projection of data objects onto the score vector indicated. The vertical axis provides a jitter for ease of visual interpretation. Solid curves in on-diagonal plots are kernel density estimates, with the black curve including all data objects and the colored curves corresponding with each respective subtype. Off-diagonal plots show scatter plots of coefficients of projection of the data objects onto the score vectors in that plot's respective axes labels. More information about these plots may be found in Chapter~1 of \citet{marron2021object}. We chose to include the first two directions in the basis for the three-way joint subspace along with the four-way joint direction for ease of visual interpretation. The four-way joint component separates basal cases from other cases. This is typically the first component found in any analysis of the modes of variation in breast cancer patients, as basal cell cancers have a very different gene expression profile than the other subtypes. The two-dimensional plot of the scores along the first two directions in the three-way joint subspace separates Her2 cases from Basal cases primarily, and from Luminal A cases secondarily. This indicates the potential for identifying useful genes, proteins, and copy number regions that drive the variation in this three-way joint subspace that distinguish Her2 cases from other breast cancer subtypes.

\begin{figure}[h]
\includegraphics[width=\textwidth]{tcga4way3wayscores-9-22.png}
\caption{Scores scatterplot matrix of four-way joint direction and first two three-way joint directions. The four-way joint subspace (top left) distinguishes basal (red triangles) from other subtypes and the three-way joint subspace (bottom right) distinguishes Her2 (magenta plusses) from other subtypes.}\label{fig:tcga4way3wayscores}
\end{figure}

\subsection{Data Integration Literature}\label{sec:litreview}

A time-honored multi-block data analysis method is Canonical Correlations Analysis (CCA), proposed by \citep{cca}. Given two blocks of data $\mathbf{X}$ and $\mathbf{Y}$, CCA seeks to maximize the Pearson correlation between vectors from the span of each data block in trait space.
The fundamental ideas of CCA have been thoroughly extended to more general settings. Several authors propose different generalizations of CCA for locating highly correlated structure between three or more data blocks, including \citep{horst1961mcca, gcca,nielsen2002mcca}. Some others have experimented with kernel methods for CCA, as in \citep{kcca} and \citep{cai2017kcca}.

In the context of machine learning, methods like CCA are termed \textit{multi-view methods}. Multi-view learning broadly includes methods like co-training for semi-supervised learning \citep{cotraining}, SVM-2K \citep{svm2k}, subspace learning \citep{mvsubspace}, and other multi-view extensions of paradigms like active learning and ensemble learning. See \citet{mvsumshort}, \citet{mvsumlong}, and \citet{li2019survey} for more details on these extensions.



Any CCA-based method is ultimately focused on finding jointly shared structure between each available block or view of the data. Oftentimes, given low-rank approximations of each data block we are also interested in a full factorization of the signal present in each block into a joint component shared between all blocks and an individual component unique to each block. One algorithm that produces such a factorization is Joint and Individual Variation Explained (JIVE) \citep{jive}. After initially choosing a signal rank for each data matrix using a permutation testing approach, the algorithm seeks to minimize residual energy by alternating between determining joint structure and individual structure. Broadly the algorithm accomplishes its goal, but the optimization problem can proceed slowly and there is no underlying inferential justification for the chosen boundary between joint and individual structure.


A later generation of JIVE, dubbed Angle-based JIVE (AJIVE) \citep{ajive}, was proposed to address the above shortcomings. Selection of joint structure happens in a quick, single step based on principal angle analysis and the delineation between joint and individual structure is based on a bound on the angles between original and perturbed subspaces found in \citet{wedin}. Initial rank selection, however, is performed ad hoc in a separate initial step, and as described in \citet{ajive} the perturbation angle bounds used can become extremely conservative under rank mis-specification. Additionally, neither JIVE nor AJIVE consider partially-shared joint structure between subsets of blocks.





Decomposition of data blocks into partially-shared joint structure components is one of the primary frontiers in contemporary data integration research. Two approaches to the problem, \citep{gaynapartial} and \citep{bass}, both model partially-shared information via structured sparsity in a basis matrix for the concatenated data blocks. \citet{gaynapartial} determine sparse structure via bi-cross-validation, while \citet{bass} determine sparse structure via a collection of Bayesian priors. A third approach is found in the recent unpublished manuscript \citep{yi2022pss}. Their method models partially-shared information via subspace intersections and a hierarchical matrix nuclear norm regularization scheme. This formulation eschews factorization into scores and loadings subspaces entirely in order to work with a convex optimization problem, to capture potentially non-orthogonal partially-shared structure across block collections, and to ensure identifiability. DIVAS maintains identifiability even among potentially non-orthogonal partially-shared structures via a sequential search through each collection of data blocks. DIVAS also involves a more challenging non-convex optimization problem based on factorized matrices, but in doing so achieves relevant statistical significance measurements for the resulting shared and partially-shared subspaces.






Another recent pursuit in data integration research is complete incorporation of all the information from the data blocks. In most cases the data blocks share the same number of data objects, so integrative analysis often takes place primarily in trait space, with corresponding information about the contributions of certain traits to shared structure being determined subsequently. In many cases, information about the contributing traits is just as pertinent as the shared structure itself, including if data blocks are bi-dimensionally linked as in \citep{lock2020bidimensional} or bi-dimensionally matched as in \citep{yuan2021doublematched}. This is often the case in bioinformatics where the traits represent measurements of particular genes and the primary goal is to identify genes or other biological factors that contribute to patterns observed across the data blocks. The above papers each propose their own method for incorporating trait information in the analysis in the situations where the data blocks are appropriately linked. Another methodology found in \citep{shu2021cdpa} attempts to incorporate trait information for more general multi-block situations, but it relies on a computationally-taxing row-matching algorithm as part of its procedure and the method cannot parse partially-shared joint structure. Our method utilizes subspace perturbation theory that applies along either dimension of the data blocks and in any scenario with mutli-block data. This makes trait information easy to incorporate throughout the algorithm that locates partially-shared structure.



\section{Methodology}\label{sec:methodology}

Let $\mathbf{X}_1,\dots,\mathbf{X}_K$ be data blocks each containing the same set of $n$ data objects and distinct sets of $d_1,\dots, d_K$ traits. In matrix calculations we follow the bioinformatics convention where matrix columns are data objects and matrix rows are traits (i.e. the matrix $\mathbf{X}_k$ has $d_k$ rows and $n$ columns). In our data model, shown in \eqref{eq:datamodel}, each data block is assumed to be a low-rank signal matrix $\mathbf{A}_k$ plus a full-rank noise matrix $\mathbf{E}_k$: 
\begin{equation}
\label{eq:datamodel}
\mathbf{X}_k = \mathbf{A}_k + \mathbf{E}_k.
\end{equation}
We assume each entry of $\mathbf{E}_k$ is independent with identical variance $\sigma^2$. For inferential purposes, we also assume \textit{rotational invariance} between signal and noise data matrices. The mathematical details of this assumption are provided in Section~\ref{sec:rotboot} where they are maximally relevant. 

Further discussion of the modeling assumptions for DIVAS requires notation for describing different collections of data blocks in detail. Consider the power set $2^{\{1,\dots,K\}}$ as a set of index sets, where each element $\mathbf{i}\in2^{\{1,\dots,K\}}$ represents a particular collection of data block indices. Each $\mathbf{i}$ indexes a shared structure among the data blocks $\mathbf{X}_k$ with $ k\in\mathbf{i}$. Denote by $|\mathbf{i}|$ the cardinality of $\mathbf{i}$.


We impose additional structure on each signal matrix $\mathbf{A}_k,\ k=1\ldots,K$ to model the partially-shared structure across the data blocks. We assume $\mathbf{A}_k$ is formed via a sum of low-rank signal matrices, each of which corresponds to the joint structure shared between the collection of blocks indicated by the index $\mathbf{i}$:
\begin{equation}
    \label{eq:signalmodel}
    \mathbf{A}_k = \sum_{\mathbf{i} | k\in\mathbf{i}} \mathbf{L}_{\mathbf{i},k}\mathbf{V}_{\mathbf{i}}^\top,
\end{equation}
where the sum extends over all index sets $\mathbf{i}\in 2^{\{1,\ldots,K\}}$ that satisfy $k\in\mathbf{i}$.
The $n\times r_\mathbf{i}$ \textit{scores matrices} $\mathbf{V}_{\mathbf{i}}$
model the shared structure in the trait space between the data blocks $\mathbf{X}_k$ with $k\in\mathbf{i}$. The $d_k\times r_\mathbf{i}$ \textit{loadings matrices} $\mathbf{L}_{\mathbf{i},k}$  contain the induced object space structure in each block $\mathbf{X}_k$ with $k\in\mathbf{i}$.
We require that the decomposition \eqref{eq:signalmodel} satisfies the following assumption. Here and for the rest of the manuscript the notation $\left[\mathbf{V}_{\mathbf{i}} \right]_{\mathbf{i}|\mathbf{i}\in S}$ denotes horizontal matrix concatenation $\left[ \mathbf{V}_{\mathbf{i}_1} \cdots \mathbf{V}_{\mathbf{i}_{|S|}} \right]$ of all matrices $\mathbf{V}_{\mathbf{i}}$ with $\mathbf{i}\in S$.
\begin{assumption}\label{asump:decomposition}
The factorized signal matrices $\mathbf{L}_{\mathbf{i}, k}$ and $\mathbf{V}_{\mathbf{i}}$ satisfy: 
\begin{enumerate}
    \item The columns of each $\mathbf{V}_{\mathbf{i}}$ are orthonormal.
    \item For two different block index sets $\mathbf{i}\neq\mathbf{j}$, if $\mathbf{i}\subset\mathbf{j}$ or $\mathbf{j}\subset\mathbf{i}$, then the subspaces spanned by the columns of $\mathbf{V}_{\mathbf{i}}$ and $\mathbf{V}_{\mathbf{j}}$ in the trait space are orthogonal.
    \item The matrix $\left[\mathbf{V}_{\mathbf{i}}\right]_{\mathbf{i}|\mathbf{i}\in 2^{\{1,\dots , K\}}}$ concatenated over all $\mathbf{i}\in 2^ {\{1,\ldots,K\}}$ has rank equal to its number of columns.
    \item  For all $k$, the matrix $\left[\mathbf{L}_{\mathbf{i},k}\right]_{\mathbf{i}\,|\,k\in\mathbf{i}}$ concatenated over all $\mathbf{i}$ so that $k\in\mathbf{i}\in 2^ {\{1,\ldots,K\}}$ has rank equal to its number of columns.
   
   
\end{enumerate}
\end{assumption}
Note that the columns of the loadings matrices $\mathbf{L}_{\mathbf{i}, k}$ are not required to be orthogonal and may have arbitrary magnitude in order to encode scale information. The $d_k\times n$ matrix $\mathbf{A}_{\mathbf{i},k}=\mathbf{L}_{\mathbf{i},k}\mathbf{V}_{\mathbf{i}}^\top$ has rank $r_{\mathbf{i}}$, the number of columns of $\mathbf{V}_\mathbf{i}$. The matrix $\mathbf{A}_{\mathbf{i},k}$  will be called the {\em partially shared joint structure} between blocks in $\mathbf{i}$. When $\mathbf{i}$ is a singleton we also call it {\em individual structure}, and when $\mathbf{i}=\{1,\ldots,K\}$ we also call it {\em fully joint structure}.

Under these assumptions, we prove a result about the existence and uniqueness of the joint structure decomposition in the absence of noise:
\begin{theorem}\label{th:construct}
For a set of signal matrices $\mathbf{A}_1,\dots ,\mathbf{A}_K$, there exists a set of matrices $\mathbf{L}_{\mathbf{i},k}, \mathbf{V}_{\mathbf{i}}$ satisfying \eqref{eq:signalmodel} and  Assumption~\ref{asump:decomposition}. The joint structure matrices $\mathbf{A}_{\mathbf{i},k}=\mathbf{L}_{\mathbf{i},k}\mathbf{V}_{\mathbf{i}}^\top$ are uniquely determined for all $\mathbf{i}\in 2^ {\{1,\ldots,K\}}$.
\end{theorem}
\begin{proof}
We proceed with a constructive proof by induction. Denote by $\mathcal{V}_{\mathbf{i}}$ the intersection of the subspaces spanned by the columns of signal matrices $\mathbf{A}_k$ with $k\in\mathbf{i}$ in trait space. 

Step 1: Consider the index set $\mathbf{i}=\{1,\ldots,K\}$. Choose $\mathbf{V}_{\{1,\dots , K\}}$ such that its columns form an orthonormal basis for $\mathcal{V}_{\{1,\dots , K\}}$.  Clearly the relevant parts of Assumption~\ref{asump:decomposition} are satisfied. Also notice that any vector orthogonal to $\mathbf{V}_{\{1,\dots , K\}}$ is in at most $K-1$ subspaces spanned by the columns of signal matrices $\mathbf{A}_k$.


Step 2: Let us assume that we have defined $\mathbf{V}_{\mathbf{i}}$ for all $K\geq |\mathbf{i}|\geq q$ that satisfy Assumption~\ref{asump:decomposition}, and,
for any $|\mathbf{i}|=q$, any vector orthogonal to $\left[\mathbf{V}_{\mathbf{j}}\right]_{\mathbf{j}\supset\mathbf{i}}$ is included in at most $q-1$ subspaces spanned by the columns of signal matrices $\mathbf{A}_k$, $k\in\mathbf{i}$.
For any $\mathbf{i}$ with $|\mathbf{i}|=q-1$ select  $\mathbf{V}_{\mathbf{i}}$ such that its columns form an orthonormal basis for $\mathcal{V}_{\mathbf{i}} \cap \left[\mathbf{V}_{\mathbf{j}}\right]_{\mathbf{j}\supset\mathbf{i}}^\perp$, the part of the space $\mathcal{V}_{\mathbf{i}}$ orthogonal to all $\mathbf{V}_{\mathbf{j}}$, $\mathbf{j}\supset\mathbf{i}$.
Each $\mathbf{V}_{\mathbf{i}}$ chosen this way satisfies parts 1 and 2 of Assumption~\ref{asump:decomposition} by construction.
Condition~3 is satisfied, because if there was a rank deficiency in a concatenated matrix $\left[\mathbf{V}_{\mathbf{i}}\right]_{\mathbf{i}\,|\, |\mathbf{i}|\geq q-1}$ there would be two indices $\mathbf{i}$, $\mathbf{j}$ such that $\mathbf{i}\neq\mathbf{j}$ and the spans of the matrices $\mathbf{V}_{\mathbf{i}}$ and $\mathbf{V}_{\mathbf{j}}$ share some vector in common. However this vector would already be included in $\mathbf{V}_{\mathbf{i}\cup\mathbf{j}}$ which is a contradiction.
Finally, for any $|\mathbf{i}|=q-1$, any vector orthogonal to $\left[\mathbf{V}_{\mathbf{j}}\right]_{\mathbf{j}\supset\mathbf{i}}$ is included in at most $q-2$ subspaces spanned by the columns of signal matrices $\mathbf{A}_k$, $k\in\mathbf{i}$.
This completes the inductive construction of the collection of $\mathbf{V}_{\mathbf{i}}$.

For each $k$, notice that the column span of $\left[\mathbf{V}_{\mathbf{i}}\right]_{\mathbf{i}|k\in\mathbf{i}}$ is $\mathcal V_{\{k\}}$, the space spanned by columns of $\mathbf{A}_k$. Thanks to Condition~3 of Assumption~\ref{asump:decomposition}, the matrices $\mathbf{L}_{\mathbf{i},k}$ are chosen as the unique solution of the equation formed using the concatenated matrices
$\left[\mathbf{L}_{\mathbf{i},k}\right]_{\mathbf{i}|k\in\mathbf{i}}\cdot \left[\mathbf{V}_{\mathbf{i}}\right]_{\mathbf{i}|k\in\mathbf{i}}^\top=\mathbf{A}_k.$ These $\mathbf{L}_{\mathbf{i},k}$ satisfy Condition~4 of Assumption~\ref{asump:decomposition} by construction.

Now assume there exists some other collection
$
    \mathbf{A}_k = \sum_{\mathbf{i} | k\in\mathbf{i}} \tilde{\mathbf{L}}_{\mathbf{i},k}\tilde{\mathbf{V}}_{\mathbf{i}}^\top
$
also satisfying Assumption~\ref{asump:decomposition}. Following  similar arguments as above we see that the columns spaces of $\tilde{\mathbf{V}}_{\mathbf{i}}$ and ${\mathbf{V}}_{\mathbf{i}}$ are the same, and therefore there exists an orthonormal $r_{\mathbf{i}}\times r_{\mathbf{i}}$  matrix $\mathbf{Q}_{\mathbf{i}}$, so that  $\tilde{\mathbf{V}}_{\mathbf{i}}=\mathbf{V}_{\mathbf{i}}\mathbf{Q}_{\mathbf{i}}$. Consequently
$\tilde{\mathbf{L}}_{\mathbf{i}, k}=\mathbf{L}_{\mathbf{i}, k}\mathbf{Q}_{\mathbf{i}}$ and
$\tilde{\mathbf{L}}_{\mathbf{i},k}\tilde{\mathbf{V}}_{\mathbf{i}}^\top=\mathbf{L}_{\mathbf{i},k}\mathbf{V}_{\mathbf{i}}^\top$.
\end{proof

Since  the matrices $\mathbf{L}_{\mathbf{i}, k}$ and $\mathbf{V}_{\mathbf{i}}$ are only determined up to basis rotation,
the most important information contained in $\mathbf{L}_{\mathbf{i}, k}$ and $\mathbf{V}_{\mathbf{i}}$ is the subspaces their columns span in object space and trait space respectively.   For interpretive purposes, it can be helpful to choose a particular informative basis for the shared subspaces and examine modes of variation of the data along those basis directions. These modes of variation of the data may be formed by outer-multiplying corresponding columns of suitably rotated $\mathbf{L}_{\mathbf{i},k}$ and $\mathbf{V}_{\mathbf{i}}$. In Section~\ref{sec:recon}, we discuss a procedure to choose such informative basis rotations $\mathbf{Q}_\mathbf{i}$ related to SVD.

Throughout the description of methodology we will use the following synthetic three-block data example to illustrate each step of DIVAS. Each block includes a different set of traits associated with 400 observations. To mimic challenging data situations with large disparities in trait set sizes, Block 1 has 200 traits, Block 2 has 400 traits, and Block 3 has $10000$ traits. Figure~\ref{3blocktruth} displays this synthetic data set using matrix \textit{heatmaps}. Heatmaps are a graphical display of matrix entry magnitude using color. Negative entries are shown in shades of blue and positive entries are shown in red, with color saturation indicating the magnitude of each entry. This means that entries close to zero are shown with a low-saturation white color. The color scaling ranges of each heatmap are shown in the color bar below each individual plot.

Each row demonstrates the formation of one of the data blocks via the data model in \eqref{eq:datamodel}. The left-most column of heatmaps shows the observed data blocks $\mathbf{X}_k$, and the right-most column of heatmaps shows the noise matrices $\mathbf{E}_k$, which in this example are i.i.d. Gaussian matrices. The middle columns display the various rank 1 components that sum to each block's signal matrix $\mathbf{A}_k$ as in \eqref{eq:signalmodel}. Each signal matrix is comprised of a fully-shared component (second column) and partially-shared components between two of the three data blocks (third, fourth, and fifth columns). As required by our model, these fully-shared and partially-shared components are constructed such that the trait space subspace of the fully-shared component is orthogonal to the corresponding trait space subspaces of each partially-shared component. However, the partially-shared subspaces are not mutually orthogonal. In fact, each pair of partially-shared trait subspaces each has a principal angle of 60 degrees between them in the trait space, $\mathbb{R}^{400}$. Adding the matrices in the middle columns across each row in the manner of \eqref{eq:signalmodel} combines the signal components into rank 3 signal matrices. Adding the noise matrices in the manner of \eqref{eq:datamodel} then produces the observed data blocks (first column).

\begin{figure}[h]
\centering
\includegraphics[scale=0.51]{3blocktruthAnnotated.png}
\caption{Heatmap view of $K=3$ synthetic example construction. Heatmaps share a common color scheme displayed in the color bars below each plot, with white representing 0 magnitude. The three blocks in the first column are the observed data, and are formed by adding up the other matrices in their respective row. The three blocks in the second column show the rank 1 fully-shared structure common to each block. The next three columns show the rank 1 partially-shared structure common to each subset of two blocks. The final column shows the noise matrices for each data block.}\label{3blocktruth}
\end{figure}


The procedure of DIVAS takes the observed data blocks as input and outputs a full breakdown of shared and partially-shared structure between them over three steps. The first step, described in Section~\ref{sec:sse}, extracts and estimates the dimension, magnitude, and direction of each block's signal subspace. The second step, described in Section~\ref{sec:jse}, combines the information from each block to locate shared directions between subspaces. The third step, described in Section~\ref{sec:recon}, uses those shared directions to form estimates $\hat{\mathbf{A}}_{\mathbf{i},k}$ of the partially shared joint structure matrices for each data block, e.g., estimates of the middle columns of Figure~\ref{3blocktruth}. Section~\ref{sec:graphics} describes the visual display of DIVAS output and corresponding diagnostic measurements. The properties and analyses described in Section \ref{sec:sse} apply to each data block independently. For notational simplicity, in that section we suppress the subscript $k$ indexing data blocks when referring to data-block-specific quantities for simplicity of notation. The indexing subscript will return in the remaining sections when information from different data blocks is combined together.

\subsection{Signal Subspace Extraction}\label{sec:sse}



The first step of DIVAS is to estimate the subspaces spanned in both object space and trait space by the signal matrix for each data block. When we observe $\mathbf{X}$, we are observing data that's been perturbed in both magnitude (singular values) and orientation (basis vectors) from the signal $\mathbf{A}$. The signal magnitude is readily recoverable from the data magnitude via the signal extraction procedure described subsequently in Section~\ref{sec:mpsignal}. The signal orientation is itself more challenging to estimate from the data orientation; there's no reason to favor one direction of rotation over another under the rotational invariance assumptions on the noise matrix $\mathbf{E}$. However, the key to DIVAS is to quantify a range of feasible signal orientations given the observed data orientation using bounds on principal angles. These techniques are described in Sections~\ref{sec:perturb} and \ref{sec:rotboot}. The results of this step for the synthetic data set presented in Figure~\ref{3blocktruth} are explored in Section~\ref{sec:syntheticsse}

\subsubsection{Signal Subspaces}\label{sec:mpsignal}


\citet{sn2013} demonstrated that, under the assumptions of \eqref{eq:datamodel}, any procedure for extracting signal from a data matrix $\mathbf{X}$ need only consider the singular values of the data matrix. We perform SVD on $\mathbf{X}$ to find $\mathbf{X} = \overbar{\mathbf{U}}\overbar{\mathbf{D}}\overbar{\mathbf{V}}^{\top}$. The columns of the matrices $\overbar{\mathbf{U}}$ and $\overbar{\mathbf{V}}$ are orthonormal bases for the subspaces spanned in object space and trait space of $\mathbf{X}$ respectively, and the diagonal entries of $\overbar{\mathbf{D}}$ are the singular values of $\mathbf{X}$. Denote these singular values as $\overbar{\nu}_1,\dots,\overbar{\nu}_{d\wedge n}$. Estimations of the signal matrix $\hat{\mathbf{A}}$ typically take the form of a decomposition in terms of rank 1 matrices/approximations that combine to form the estimated object space and trait space subspaces. The vectors $\overbar{\mathbf{u}}_i$ and $\overbar{\mathbf{v}}_i$ denote the $i$th columns of the matrices $\overbar{\mathbf{U}}$ and $\overbar{\mathbf{V}}$ respectively, and $\eta(\bullet)$ is a function from $\mathbb{R}^+$ to $\mathbb{R}^+$ for shrinking the singular values:
\begin{equation}
\label{eq:signalest}
\hat{\mathbf{A}} = \sum_{i=1}^{d\wedge n} \eta(\overbar{\nu}_i) \overbar{\mathbf{u}}_i \overbar{\mathbf{v}}_i^{\top}.
\end{equation}


Common choices for $\eta$ include \textit{soft thresholding}: $\eta_{soft}(\nu) = (\nu-c)\vee 0$, and \textit{hard thresholding}: $\eta_{hard}(\nu) = \nu \mathbb{I}_{\{\nu\geq c\}}$, for some constant $c$, and where $\mathbb{I}_{\{\bullet\}}$ represents an indicator function. In either case, any singular value smaller than $c$ is set to 0. This means both procedures have dimension-reducing effects on the estimated $\hat{\mathbf{A}}$, and only subspaces $\overbar{\mathbf{u}}_i \overbar{\mathbf{v}}_i^{\top}$ associated with nonzero transformed singular values contribute to the estimate. \citet{gd2014} outline optimal choices for $c$ for both soft and hard thresholding in terms of the \textit{aspect ratio} $\beta = \frac{d\wedge n}{d\vee n}$ and the standard deviation $\sigma$ of the noise matrix $\mathbf{E}$. 

In the additive noise data matrix model \eqref{eq:datamodel}, the presence of noise inflates the singular values associated with the signal component of the data matrix. Hard thresholding does not account for this phenomenon at all and soft thresholding often overcorrects by applying the same amount of shrinkage to each nonzero singular value. \citet{sn2013} and \citet{gd2017} each propose optimal thresholding functions based on the Marchenko-Pastur distribution (see Section~\ref{sec:mp}) under a variety of matrix norms. We use the operator-norm-optimal function $\eta^*$ from \citep{gd2017} shown in \eqref{eq:opShrink} for DIVAS. 
\begin{equation}
\label{eq:opShrink}
\eta^*(\nu) = \begin{cases}
\frac{1}{\sqrt{2}}\sqrt{\nu^2 - \beta - 1 + \sqrt{(\nu^2 - \beta - 1)^2 - 4\beta}}, & \nu \geq 1 + \sqrt{\beta} \\
0, & \nu < 1 + \sqrt{\beta}
\end{cases}
\end{equation}
Figure~\ref{threshes} demonstrates how this shrinkage function \eqref{eq:opShrink}, blue solid line, compromises between soft and hard thresholding for singular values with different magnitudes for a matrix with aspect ratio $\beta = 1$ and noise standard deviation $\sigma = 1$. Small values are thresholded according to optimal soft thresholding (magenta dot-dash line), but the shrinkage function approaches optimal hard thresholding (black dashed line) for larger values.

\begin{figure}[h]
\centering
\includegraphics[scale=0.4]{OpThresh.png}
\caption{Functions for hard thresholding, soft thresholding, and optimal shrinkage under operator norm loss for a square matrix ($\beta=1$). The optimal shrinkage function compromises between the other two approaches. Figure produced with code from \cite{gd2017}.}\label{threshes}
\end{figure}

Equation~\eqref{eq:opShrink} assumes noise standard deviation $\sigma = 1$. To use the shrinkage function in general settings, we must appropriately scale the singular values before and after shrinkage according to some estimate of the standard deviation of the noise $\hat{\sigma}$. \citet{sn2013} use a grid search over several candidate values for $\hat{\sigma}$ to find a value that minimizes the Kolmogorov-Smirnov distance between the non-signal singular values and the appropriate Marchenko-Pastur distribution. \citet{gd2017} opt for the simple, robust, closed-form estimate $\hat{\sigma} = \frac{\nu_{median}}{\sqrt{MP(\beta)_{0.5}}}$, where $\nu_{median}$ denotes the median singular value of $\mathbf{X}$ and $MP(\beta)_{0.5}$ denotes the median of the Marchenko-Pastur distribution with parameter $\beta$ (see Section~\ref{sec:mp}). We use the latter method for DIVAS noise standard deviation estimation.

Combining the previous equations, our estimate for the signal matrix $\hat{\mathbf{A}}$ for a given data block $\mathbf{X}$ is:
\begin{equation}
\label{eq:finala}
\hat{\mathbf{A}} =  \sum_{i=1}^{d\wedge n} \hat{\sigma}\eta^*(\overbar{\nu}_i/\hat{\sigma}) \overbar{\mathbf{u}}_i \overbar{\mathbf{v}}_i^{\top}.
\end{equation}
Let $\hat{\nu_i} = \hat{\sigma}\eta^*(\overbar{\nu}_i/\hat{\sigma})$ be the $i$th shrunken singular value of $\mathbf{X}$. Let $\hat{r}$ be the number of nonzero shrunken singular values, and therefore the estimated rank of $\mathbf{A}$. Let $\hat{\mathbf{U}}$ and $\hat{\mathbf{V}}$ be matrices containing the first $\hat{r}$ columns of $\overbar{\mathbf{U}}$ and $\overbar{\mathbf{V}}$, respectively. Using this notation and defining the matrix $\hat{\mathbf{D}}$ as the $\hat{r}\times\hat{r}$ diagonal matrix with diagonal entries equal to $\hat{\nu}_1,\dots,\hat{\nu}_{\hat{r}}$, we can also write $\hat{\mathbf{A}} = \hat{\mathbf{U}}\hat{\mathbf{D}}\hat{\mathbf{V}}^{\top}$. Note that $\hat{\mathbf{U}}$ is therefore an orthonormal basis for the subspace spanned in object space of $\hat{\mathbf{A}}$ and $\hat{\mathbf{V}}$ is an orthonormal basis for the subspace spanned in trait space of $\hat{\mathbf{A}}$.

\subsubsection{Angle Perturbation Theory} \label{sec:perturb}


The foundation of DIVAS is determining whether candidate directions $\mathbf{v}^{\star}\in \mathbb{R}^n$ lie in the trait space span $\textbf{TS}(\mathbf{A})$ of the signal matrix $\mathbf{A}$. If $\mathbf{A}$ was observable, this would simply amount to checking whether the angle $\theta$ between $\mathbf{v}^{\star}$ and $\textbf{TS}(\mathbf{A})$ was 0. Since $\mathbf{A}$ and $\theta$ are unobservable, we aim to estimate $\theta$ based on the observable estimated low-rank signal matrix $\hat{\mathbf{A}}$ and an estimate of the noise variation $\hat{\mathbf{E}}$ developed in Appendix~\ref{sec:noise}. Specifically, we want to choose a \textit{perturbation angle bound} $\hat{\phi}$ that defines a cone-shaped significance region around $\textbf{TS}(\hat{\mathbf{A}})$ which contains $\textbf{TS}(\mathbf{A})$ with high probability. Directions $\mathbf{v}^{\star}$ lying within that significance region would then be potential basis directions for $\textbf{TS}(\mathbf{A})$. Directions $\mathbf{v}^{\star}$ that lie within the significance regions of multiple data blocks would then be potential basis directions for partially-shared joint structure between those blocks. To arrive at such a data-block-wise uniform perturbation angle bound $\hat{\phi}$, we first look to bound the range of possible values for $\theta$ for one given candidate direction $\mathbf{v}^{\star}$.





The first step is construction of the range of values for $\theta$ based on the relationships between the projections of $\mathbf{v}^{\star}$ onto various subspaces. This is illustrated using a simple low-dimensional example in Figure~\ref{thetastars}. In particular, $\mathbf{v}^{\star}$ (green vector in both panels of Figure~\ref{thetastars}) is projected onto each of $\textbf{TS}(\mathbf{A})$ (translucent purple plane) and $\textbf{TS}(\hat{\mathbf{A}})$ (solid gold plane), with those projections denoted $\mathbf{v}^{\star}_{proj}$ (red solid lines) and $\hat{\mathbf{v}}^{\star}_{proj}$ (blue solid lines) respectively. Computations of these projections are based on the orthonormal basis matrices $\mathbf{V}$ for $\textbf{TS}(\mathbf{A})$ and $\hat{\mathbf{V}}$ for $\textbf{TS}(\hat{\mathbf{A}})$. Let $\hat{\theta}$ be the angle between $\mathbf{v}^{\star}$ and $\textbf{TS}(\hat{\mathbf{A}})$. We can write expressions for $\theta$ and $\hat{\theta}$ in terms of the above quantities as follows: $$\theta = \arccos\left(\frac{\langle\mathbf{v}^{\star},\mathbf{v}^{\star}_{proj}\rangle}{\parallel \mathbf{v}^{\star}\parallel\parallel\mathbf{v}^{\star}_{proj}\parallel}\right);\;\; \mathbf{v}^{\star}_{proj} = \mathbf{V}\Vb^{\top}\mathbf{v}^{\star}.$$
$$\hat{\theta} = \arccos\left(\frac{\langle\mathbf{v}^{\star},\hat{\mathbf{v}}^{\star}_{proj}\rangle}{\parallel \mathbf{v}^{\star}\parallel\parallel\hat{\mathbf{v}}^{\star}_{proj}\parallel}\right);\;\; \hat{\mathbf{v}}^{\star}_{proj} = \hat{\mathbf{V}}\hat{\mathbf{V}}^{\top}\mathbf{v}^{\star}.$$

We can construct bounds involving $\theta$ and $\hat{\theta}$ by considering further projections between $\textbf{TS}(\mathbf{A})$ and $\textbf{TS}(\hat{\mathbf{A}})$. The total angle traversed by projecting $\mathbf{v}^{\star}$ to $\textbf{TS}(\mathbf{A})$ and then projecting that result, $\mathbf{v}^{\star}_{proj}$, onto $\textbf{TS}(\hat{\mathbf{A}})$ (red dashed line in left panel of Figure~\ref{thetastars}) is at least as large as $\hat{\theta}$, the angle between $\mathbf{v}^{\star}$ and $\textbf{TS}(\hat{\mathbf{A}})$. Define $\theta^{\star}_1$ as the angle between $\mathbf{v}^{\star}_{proj}$ and $\textbf{TS}(\hat{\mathbf{A}})$. By the triangle inequality, $\hat{\theta} \leq \theta + \theta^{\star}_1$. Via an analogous projection of $\hat{\mathbf{v}}^{\star}_{proj}$ onto $\textbf{TS}(\mathbf{A})$ (blue dashed line in right panel of Figure~\ref{thetastars}), define $\theta^{\star}_2$ as the angle between $\hat{\mathbf{v}}^{\star}_{proj}$ and $\textbf{TS}(\mathbf{A})$. Then $\theta$, the angle between $\mathbf{v}^{\star}$ and $\textbf{TS}(\mathbf{A})$, is no larger than the total angle traversed by projecting $\mathbf{v}^{\star}$ onto $\textbf{TS}(\hat{\mathbf{A}})$ and then projecting that result, $\hat{\mathbf{v}}^{\star}_{proj}$, onto $\textbf{TS}(\mathbf{A})$: $\theta \leq \hat{\theta} + \theta^{\star}_2$. 

\begin{figure}[h]
\includegraphics[scale=0.3]{AnglePerturbationGraphic3.png}
\caption{Locations of $\theta$, $\hat{\theta}$, $\theta^{\star}_1$, and $\theta^{\star}_2$ in a low-dimensional example. Each panel demonstrates a different angle bound. Left: $\hat{\theta} \leq \theta + \theta^{\star}_1$. Right: $\theta \leq \hat{\theta} + \theta^{\star}_2$.}\label{thetastars}
\end{figure}

The above discussion of angles between subspaces summarizes the proof of the following theorem. More details can be found in the Ph.D. dissertation of \citet{mjdiss}.
\begin{theorem}
\label{thetastarMJ}
Let $\mathbf{X} = \mathbf{A} + \mathbf{E}$ be a $d\times n$ data matrix which is a sum of a signal matrix $\mathbf{A}$ and a noise matrix $\mathbf{E}$ under the assumptions of \eqref{eq:datamodel}. Given $\theta$, $\hat{\theta}$, $\theta_1^{\star}$, and $\theta_2^{\star}$ defined as above, we have:
\begin{equation}
\label{thetastarbounds}
(\hat{\theta} - \theta_1^{\star})_+ \leq \theta \leq \hat{\theta}+\theta_2^{\star}.
\end{equation}
\end{theorem}
Both inequalities in \eqref{thetastarbounds} will be used to rule out directions $\mathbf{v}^{\star}$ as a candidate for the joint spaces. These exclusions will be based on two distinct statistical arguments: a novel rotational bootstrap used for both bounds, and a distribution of angles between random directions used only for the upper bound.

If the angle to the estimated signal subspace $\hat{\theta}$ for a given candidate direction $\mathbf{v}^{\star}$ is less than $\theta_1^{\star}$, then the lower bound in \eqref{thetastarbounds} for the angle to the true signal subspace is 0, indicating that this direction can't be ruled out as lying in the true signal subspace. The angle $\theta^{\star}_1$ is behaving much like the desired perturbation angle bound, but as noted in \citet{mjdiss}, $\theta_1^{\star}$ is not directly estimable for a given direction. However, $\theta_1^{\star}$ is uniformly bounded from above for all $\mathbf{v}^{\star}$ by the maximum principal angle $\phi$ between $\textbf{TS}(\mathbf{A})$ and $\textbf{TS}(\hat{\mathbf{A}})$ (see Appendix~\ref{sec:paa}). Our chosen perturbation angle bound will therefore be a statistical estimate $\hat{\phi}$ of that maximum principal angle $\phi$. This estimation is performed via a rotational bootstrap as described in Section~\ref{sec:rotboot}.


Unlike $\theta_1^{\star}$, the angle $\theta_2^{\star}$ in the upper bound of \eqref{thetastarbounds} can be estimated for each $\mathbf{v}^{\star}$ using
\begin{align}
\begin{split}
\label{eq:thetatwostar}
\theta_2^{\star} 
&= \arccos\left(\frac{\parallel \mathbf{V}^{\top}\hat{\mathbf{V}}\hat{\mathbf{V}}^{\top}\mathbf{v}^{\star} \parallel}{\parallel \hat{\mathbf{V}}\hat{\mathbf{V}}^{\top}\mathbf{v}^{\star} \parallel}\right).
\end{split}
\end{align}
The only unknown quantity in this formula is the matrix $\mathbf{V}^{\top} \hat{\mathbf{V}}$, and we generate samples from the estimated distribution of this matrix as part of the rotational bootstrap in Section~\ref{sec:rotboot}. Therefore by recording those samples and using them in \eqref{eq:thetatwostar} we can generate from a bootstrap distribution of $\theta_2^{\star}$ and choose a high percentile denoted as $\hat{\theta}_2^{\star}$. 

The main use of this estimate is determining whether $\mathbf{v}^{\star}$ can be distinguished from an arbitrarily chosen direction based on the estimated upper bound $\hat{\theta}+\hat{\theta}^{\star}_2$. The rotational invariance property assumed of the signal and noise in \eqref{eq:datamodel} implies a natural null distribution for comparison. In particular, we choose the distribution of angles between a fixed arbitrary $\hat{r}$-dimensional subspace of $\mathbb{R}^n$ (recall that $\hat{r}$ is the estimated signal rank) and unit vectors chosen uniformly at random. We pick a \textit{random direction angle bound} $\theta_0$ as a low percentile of that null distribution. If $\hat{\theta}+\hat{\theta}^{\star}_2$ lies above $\theta_0$ for some direction $\mathbf{v}^{\star}$, then that direction cannot be distinguished from an arbitrarily chosen direction, which provides statistical evidence that $\mathbf{v}^{\star}$ is far from $\textbf{TS}(\mathbf{A})$.



The above derivations of a perturbation angle bound and other angle-based inference have taken place entirely in trait space. Analogous derivations can be carried out in object space, and the estimation of perturbation angle bounds in both spaces can take place simultaneously during the rotational bootstrap. When considering candidate directions $\mathbf{v}^{\star}$, we should additionally rule out directions whose corresponding basis directions in object space do not obey the object space perturbation angle bounds. Therefore both space's angle bounds play key roles in the optimization problem for locating joint structure between data blocks. This leads to more precise estimates of joint subspaces compared to methods like AJIVE \citep{ajive} that consider only trait space in their algorithms.

\subsubsection{Rotational Bootstrap}\label{sec:rotboot}

We estimate perturbation angle bounds for the object space and trait space of a data block using a novel \textit{rotational bootstrap}. This technique is designed to take advantage of the assumed rotational invariance property and  aims to estimate the distribution of principal angles between object space and trait space subspaces of $\mathbf{X}$ and $\mathbf{A}$ through random generation of replicate signal subspaces.

Recall $\mathbf{X}=\mathbf{A}+\mathbf{E}$ and as in Section~\ref{sec:perturb} the compact SVD of the rank $r$ signal matrix is $\mathbf{A}=\mathbf{U} \mathbf{D} \mathbf{V}^\top$ and the compact SVD of the rank $r$ approximation to the data $\mathbf{X}$ is 
$\hat\mathbf{A}=\hat\mathbf{U} \hat\mathbf{D} \hat\mathbf{V}^\top$. Next consider a random replication $\mathbf{X}^\diamond=\mathbf{U}^\diamond \mathbf{D} \mathbf{V}^{\diamond\top}+\mathbf{E}^\diamond$, where $\mathbf{E}^\diamond$ has the same distribution as $\mathbf{E}$, and $\mathbf{U}^\diamond$ and $\mathbf{V}^\diamond$ are random $d\times r$ and $n\times r$ orthonormal matrices, respectively. The corresponding compact SVD of the rank $r$ approximation to $\mathbf{X}^\diamond$ is $\hat\mathbf{A}^\diamond=\hat\mathbf{U}^\diamond \hat\mathbf{D}^\diamond \hat\mathbf{V}^{\diamond\top}$.
\begin{assumption}\label{as:rotation}
 Model \eqref{eq:datamodel} is rotationally invariant if the matrices  $\hat \mathbf{D}, \mathbf{V}^\top\hat\mathbf{V}, \hat\mathbf{V}^\top\mathbf{V}, \mathbf{U}^\top\hat\mathbf{U}, \hat\mathbf{U}^\top\mathbf{U}$ have the same distribution as the corresponding $\hat\mathbf{D}^\diamond, \mathbf{V}^{\diamond\top}\hat\mathbf{V}^\diamond, \hat\mathbf{V}^{\diamond\top}\mathbf{V}^\diamond, \mathbf{U}^{\diamond\top}\hat\mathbf{U}^\diamond, \hat\mathbf{U}^{\diamond\top}\mathbf{U}^{\diamond}$.
\end{assumption}
As discussed in Appendix~\ref{sec:paa} these matrices determine the principle angle structure between the spaced spanned by the low rank signal and its estimate in both trait and object spaces.
Theorem~7 of \citet{mjdiss}  shows that if the noise distribution is rotationally invariant, e.g., having i.i.d.~centered Gaussian entries, then the model is also rotationally invariant. Alternatively, Assumption~\ref{as:rotation} will be satisfied if the signal matrix $\mathbf{A}$ is considered random following a rotationally invariant prior distribution. 

The continuity of these distributions in the singular values $\mathbf{D}$ suggests use of a (parametric-like) bootstrap estimator of these quantities based on the estimated $\hat r\times\hat r$ singular value matrix $\hat\mathbf{D}$.
In particular, we form a bootstrap replication of a signal matrix $\mathbf{A}^{\circ} = \mathbf{U}^{\circ} \hat{\mathbf{D}} \mathbf{V}^{\circ\top}$, where $\mathbf{U}^{\circ}$ and $\mathbf{V}^{\circ}$  are random $d\times\hat{r}$ and $n\times\hat{r}$ orthonormal matrices, respectively. Using this randomly rotated estimated signal matrix along with an estimate $\hat{\mathbf{E}}$ of the noise matrix $\mathbf{E}$, we form a bootstrap replication of the data matrix $\mathbf{X}^{\circ} = \mathbf{A}^{\circ} + \hat{\mathbf{E}}$. If the rotational invariance assumption is satisfied, and the estimated $\hat\mathbf{D}$ is close to $\mathbf{D}$ this construction produces replicate signal and data matrices with principal angle structure drawn from a similar distribution as the unobserved principal angle structure between the true signal and data matrices.
An important, and perhaps surprising point is that the na\"ive noise matrix estimate $\hat{\mathbf{E}} = \mathbf{X} - \hat{\mathbf{A}}$ is not appropriate for use in this construction. This is because $\mathbf{X}-\hat{\mathbf{A}}$ has insufficient energy in the directions associated with $\hat{\mathbf{A}}$. Details of this including our proposed noise matrix estimator can be found in Appendix~\ref{sec:noise}.

As mentioned in Section~\ref{sec:perturb}, we will use estimates of the maximum principal angles between the subspaces spanned in object space and trait space by $\mathbf{X}$ and $\mathbf{A}$ as perturbation angle bounds. Through repeated replications of the randomly rotated signal and data matrices described in the previous paragraph, we generate bootstrap samples estimating the \textit{distribution} of principal angles between subspaces spanned by $\mathbf{X}$ and $\mathbf{A}$ in both object and trait space. With sufficiently many replications (we use $M = 400$) we can choose high quantiles (e.g. 0.95) of the empirical distributions of maximum principal angles as statistical perturbation angle bounds. Recall that $\hat{\phi}$ is the trait space perturbation angle bound, and denote the corresponding object space perturbation angle bound as $\hat{\psi}$. 

The procedure described in Section~\ref{sec:mpsignal} discriminates noise fairly well, but as our algorithm is based on angles we find that additional angle-based rank selection is necessary for good performance. Therefore as part of the rotational bootstrap we \textit{filter} the signal subspaces according to the random direction angle bound $\theta_0$ defined in Section~\ref{sec:perturb}. In particular we choose a filtered rank $\check{r}$ such that the estimated maximum principal angles between true and estimated signal don't exceed $\xi\theta_0$, where $\xi\in(0,0.5]$ is a tuning parameter. From experimentation with various data sets we've found a value based on the golden ratio $\xi=1-\frac{2}{1+\sqrt{5}}\approx 0.382$ to frequently be an effective choice, but this hyperparameter can be tuned up or down to include more or less information from the estimated signal matrix in the analysis. 

The limitation of $\xi\leq 0.5$ follows from the statistical inference framework laid out in Section~\ref{sec:perturb}. If the lower bound $(\hat{\theta} - \hat{\phi})_+$ is 0 and the upper bound $\hat{\theta} + \hat{\theta}_2^{\star}$ is simultaneously greater than $\theta_0$ for a given candidate direction $\mathbf{v}^{\star}$, the inference procedure says there is evidence that $\mathbf{v}^{\star}$ is both significantly close to the true signal subspace and indistinguishable from an arbitrary direction. This inference outcome is completely non-informative. In this case, both $\hat{\theta}$ and $\hat{\theta}_2^{\star}$ must be bounded from above by the maximum principal angle between estimated and true signal subspaces. Therefore this non-informative inference outcome is avoided by filtering the estimated signal subspace until the rotational-bootstrap-estimated maximum principal angle is at most $0.5\cdot\theta_0$.




The above description of the rotational bootstrap algorithm is formulated in Algorithm~\ref{alg:rotboot}. During each replication, random subspaces are generated from i.i.d. standard Gaussian matrices with the same centering operations used on the data. Note that orthogonalization of an i.i.d. random matrix in this fashion is identical to sampling from a rotationally uniform distribution of subspaces, according to Theorem 2.2.1 from \citet{chikuse2012statistics}. The inner \textbf{for} loop records maximum principal angles at each possible filtered rank from 1 to $\hat{r}$. After the outer \textbf{for} loop concludes, the algorithm filters the estimated signal rank to align with the chosen value of $\xi$ and chooses perturbation angle bounds for both object and trait space. Once the filtered rank $\check{r}$ is selected, 
we filter the columns of the estimated basis matrices for the signal object and trait space subspaces to correspond with the reduced rank. Let $\check{\mathbf{U}} = \overbar{\mathbf{U}}_{1:\check{r}}$ and $\check{\mathbf{V}} = \overbar{\mathbf{V}}_{1:\check{r}}$ be the final estimates of the signal object and trait space bases respectively.

\begin{algorithm}[h]
\caption{Rotational Bootstrap}\label{alg:rotboot}
\begin{algorithmic}
\State $objectAngles\gets 90*{\mathbf{1}}^{M\times \hat{r}}$; $traitAngles\gets 90*{\mathbf{1}}^{M\times \hat{r}}$
\ForAll{$m\in \{1,\dots,M\}$}
	\State $\mathbf{U}^{\circ}\gets \text{rand}^{d\times\hat{r}}$; $\mathbf{V}^{\circ}\gets \text{rand}^{n\times\hat{r}}$
	\State \Comment Replications must be orthogonal to constant function direction in appropriate spaces.
	\If{$\mathbf{X}$ is trait-centered}
		\State $\mathbf{U}^{\circ}\gets \left(\mathbf{I}^{d\times d} - \frac{1}{d}{\mathbf{1}}^{d\times d}\right) \mathbf{U}^{\circ}$
	\EndIf
	\If{$\mathbf{X}$ is object-centered}
		\State $\mathbf{V}^{\circ}\gets \left(\mathbf{I}^{n\times n} - \frac{1}{n}{\mathbf{1}}^{n\times n}\right) \mathbf{V}^{\circ}$
	\EndIf
	\State $\mathbf{U}^{\circ} \gets \text{orth}(\mathbf{U}^{\circ})$; $\mathbf{V}^{\circ} \gets \text{orth}(\mathbf{V}^{\circ})$; $\mathbf{A}^{\circ} \gets \mathbf{U}^{\circ} \hat{\mathbf{D}}\mathbf{V}^{\circ}$; $\mathbf{X}^{\circ} \gets \mathbf{A}^{\circ} + \hat{\mathbf{E}}$
	\State $[\overbar{\mathbf{U}}^{\circ},\sim,\overbar{\mathbf{V}}^{\circ}]\gets \text{SVD}(\mathbf{X}^{\circ})$.
	\ForAll{$j\in \{1,\dots,\hat{r}\}$}
		\State \Comment Smallest singular value equal to cosine of largest angle.
		\State $[\sim,\vec{\nu}_{object} ,\sim] \gets \text{SVD}\left(\mathbf{U}^{\circ T} \overbar{\mathbf{U}}^{\circ}_{1:j}\right)$; $[\sim,\vec{\nu}_{trait} ,\sim] \gets \text{SVD}\left(\mathbf{V}^{\circ T} \overbar{\mathbf{V}}^{\circ}_{1:j}\right)$
		\State $objectAngles[m,j] \gets\arccos(\min(\vec{\nu}_{object}))$; $traitAngles[m,j] \gets\arccos(\min(\vec{\nu}_{trait}))$
	\EndFor	
\EndFor
\State $objectAnglesSort\gets \text{sort}(objectAngles, col)$ $traitAnglesSort\gets \text{sort}(traitAngles, col)$ 
\State $\check{r} \gets \min\left(\sum_{j=1}^{\hat{r}} \mathbb{I}_{\{objectAnglesSort[0.95M,j] < \xi\theta_0\}}, \sum_{j=1}^{\hat{r}} \mathbb{I}_{\{traitAnglesSort[0.95M,j] < \xi\theta_0\}} \right)$
\State $\hat{\psi} \gets objectAnglesSort[0.95M,\check{r}]$; $\hat{\phi} \gets traitAnglesSort[0.95M,\check{r}]$ 
\end{algorithmic}
\end{algorithm}

\subsubsection{Signal Extraction for Synthetic Data}\label{sec:syntheticsse}

The results of signal space extraction on the synthetic data example from Figure~\ref{3blocktruth} are shown in Figure~\ref{toyexamplestep1}. Each heatmap shows the estimated signal matrix $\hat{\mathbf{A}}$ for the respective data block. The denoising of each data block appears to have been successful when comparing the visual impression of these heatmaps to the original data matrices shown in the first column of Figure~\ref{3blocktruth}. The signal rank is correctly chosen as 3 for all three blocks. In this case, angle-based rank filtering didn't further reduce the rank beyond the value determined from eigenvalue-based rank selection. 

\begin{figure}[h]
        \begin{minipage}[c]{0.3\linewidth}
            \centering
            \includegraphics[width=\textwidth]{A_1_hat_DIVAS.png}
        \end{minipage}
        \begin{minipage}[c]{0.3\linewidth}
            \centering
            \includegraphics[width=\textwidth]{A_2_hat_DIVAS.png}
        \end{minipage}
        \begin{minipage}[c]{0.3\linewidth}
            \centering
            \includegraphics[width=\textwidth]{A_3_hat_DIVAS.png}
        \end{minipage}
        \caption{Estimated signal matrices for each block in the toy example defined in Figure \ref{3blocktruth}. The heatmaps show good recovery of the original signal patterns. The trait and object spaces of each estimated matrix are rank 3.}\label{toyexamplestep1}
\end{figure}

Since we know the true signal object and trait subspaces for the synthetic data example, we can compare the estimated perturbation angle bounds to the actual angles between estimated and true signal subspaces to check the performance of the bounds. Tables~\ref{tbl:traittruthangles} and \ref{tbl:objecttruthangles} display the perturbation angle bounds and angles between the true and estimated signal subspaces in trait space and object space respectively. The calculated bounds exceed the true angles in all cases, so the true basis directions all lie within the cones of feasibility defined by the bounds. Since the bounds are calculated as uniform 95\% bounds, we'd expect to not cover the truth about 1 in 20 times, and the performance of the bounds in this case aligns with that expectation. 

Each synthetic data block has a similar signal-to-noise ratio, so the observed differences in perturbation angle bounds in the second column of each table are primarily explained by differences in matrix dimension. Recall that all four data blocks have 400 data objects, X1 has 200 traits, X2 has 400 traits, and X3 has 10000 traits. The trait space perturbation angle bounds decrease as the number of traits in the data blocks increase since we have a more precise idea of where the true signal subspace is with more trait vectors in the same trait space $\mathbb{R}^{400}$. The object space perturbation angle bounds increase as the number of traits in the data blocks increase since we have a more precise idea of where the true signal subspace is in object space with 400 object vectors in $\mathbb{R}^{200}$ than we do with 400 object vectors in $\mathbb{R}^{10000}$.

\begin{table}[h]
\centering
\begin{tabular}{cccccc}\hline
\thead{Data Block} & \thead{Trait Space \\ Angle Bound} & \thead{Angle to \\ 1,2,3 Truth} & \thead{Angle to \\ 1,2 Truth} & \thead{Angle to \{\mathbf{1}},3 Truth} & \thead{Angle to \\ 2,3 Truth} \\ \hline
1 & 11.7 & 9.2 & 8.5 & 6.1 & \\
2 & 8.6 & 6.9 & 5.6 & & 4.0 \\
3 & 2.8 & 2.5 & & 1.0 & 1.0 \\ \hline
\end{tabular}
\caption{Table of angles between estimated signal trait spaces and true signal trait spaces. All angles are within the calculated perturbation angle bounds.} \label{tbl:traittruthangles}
\end{table}

\begin{table}[h]
\centering
\begin{tabular}{cccccc}\hline
\thead{Data Block} & \thead{Object Space \\ Angle Bound} & \thead{Angle to \\ 1,2,3 Truth} & \thead{Angle to \\ 1,2 Truth} & \thead{Angle to \{\mathbf{1}},3 Truth} & \thead{Angle to \\ 2,3 Truth} \\ \hline
1 & 8.6 & 4.5 & 4.9 & 4.7 &\\
2 & 8.6 & 5.8 & 6.6 & & 4.0\\
3 & 13.1 & 7.9 & & 4.6 & 4.7\\ \hline
\end{tabular}
\caption{Table of angles between estimated signal object spaces and true signal object spaces. All angles are within the calculated perturbation angle bounds.} \label{tbl:objecttruthangles}
\end{table}

\subsection{Joint Subspace Estimation}\label{sec:jse}

Here we formally introduce the optimization problem for locating shared structure in DIVAS. The conceptual constraints and objective function are shown in \eqref{DivasOptHeur} and the full computational implementation is shown in \eqref{DivasOptCCSlackLoad}. 

For any given collection of blocks $\mathbf{i}$, the corresponding joint subspace should be near each of the included blocks in some sense. In DIVAS, proximity is evaluated in terms of angles between candidate directions and subspaces. In particular, during each phase of joint subspace estimation we minimize the angle between candidate directions $\mathbf{v}^{\star}$ and the estimated trait space subspaces of included blocks as our objective. This objective is expressed in terms of angle cosines in line 1 of \eqref{DivasOptHeur}. To ensure that a candidate direction lies in the true signal subspace of an included block $\mathbf{X}_k, k\in \mathbf{i}$ with high significance, the trait space angle between a candidate direction and the subspace spanned by the columns of $\check{\mathbf{V}}_k$ should be at most the trait space angle perturbation bound $\hat{\phi}_k$. Additionally, the object space angle between $\mathbf{X}_k\mathbf{v}^{\star}$ and the subspace spanned by the columns of $\check{\mathbf{U}}_k$ should be at most the object space angle perturbation bound $\hat{\psi}_k$. Finally, the angle between a candidate direction and an excluded block should be at least the trait space angle perturbation bound $\hat{\phi}_k$. These requirements are expressed as constraints for the optimization problem in lines 2-6 in \eqref{DivasOptHeur}, with subscripts $T$ and $O$ indicating angles in trait space and object space respectively. Crucially, our dimensionally flexible subspace-based angle perturbation approach to signal extraction allows object space information to be incorporated very naturally into the joint subspace estimation algorithm. This innovation enhances the significance and interpretability of loadings vectors found using DIVAS.

Following the proof of Theorem~\ref{th:construct}, we determine each block collection's potential joint structure in turn, starting with larger block collections and ending with singleton block collections. Within a joint structure search for a given block collection $\mathbf{i}$, joint subspace basis directions are found one at a time via successive solves of \eqref{DivasOptHeur}. If no new feasible direction is found, the search among the current block collection ends and the search among the next block collection begins. Candidate directions in trait space for a particular block collection must also obey orthogonality constraints expressed in conditions 1 and 2 of Assumption~\ref{asump:decomposition}. These conditions are concisely expressed in the constraint in line 7 of \eqref{DivasOptHeur}. The Gothic script symbol $\mathfrak{V}_{\mathbf{i}}$ is used to denote the current estimated trait space basis for the shared structure for block collection $\mathbf{i}$, i.e., an estimate of the $\mathbf{V}_\mathbf{i}$ from \eqref{eq:signalmodel}.
 
\begin{equation}\label{DivasOptHeur}
\begin{aligned} 
&\min_{\mathbf{v}^{\star}}& & - \sum_{k\in \mathbf{i}} \cos^2 \hat{\theta}_{Tk}  & & \\
&s.t.&  \hat{\theta}_{Tk} &= \angle(\mathbf{v}^{\star},\check{\mathbf{V}}_k) & \forall k & \\
&  & \hat{\theta}_{Ok} &= \angle(\mathbf{X}_k\mathbf{v}^{\star},\check{\mathbf{U}}_k) & \forall k & \\
&  &\hat{\theta}_{Tk} &\leq \hat{\phi}_k & \forall k &\in \mathbf{i} \\
&  &\hat{\theta}_{Tk} &> \hat{\phi}_k & \forall k &\in \mathbf{i}^c \\
&  &\hat{\theta}_{Ok} &\leq \hat{\psi}_k & \forall k &\in \mathbf{i} \\
&  &\mathbf{v}^{\star} &\perp \mathfrak{V}_{\mathbf{j}} & \forall \mathbf{j} &\supseteq \mathbf{i} \\
\end{aligned}
\end{equation}

Note that as basis directions are found for block collection $\mathbf{i}$, $\mathfrak{V}_{\mathbf{i}}$ grows larger, and so the constraint in line 7 becomes more stringent as more basis directions are located. Also as directions are located the angle constraints change. Each orthonormal basis for estimated trait space signal $\check{\mathbf{V}}_k$ is shrunk to only include directions in the null space of $[\mathfrak{V}_{\mathbf{j}}]_{\mathbf{j}\supseteq\mathbf{i}}$. This basis shrinking improves computation time and assists in the choice of basis directions satisfying our assumptions.

In practice this problem is solved via an iterative penalty convex-concave procedure as described in \citep{il2016} because the objective function and constraints can all be expressed as differences of convex functions: $f(\mathbf{v}^{\star}) - g(\mathbf{v}^{\star}) \leq 0$. These differences can then be appropriately linearized at a given initial condition point $\mathbf{v}_0$: $f(\mathbf{v}^{\star}) - \left(g(\mathbf{v}_{0}) + \nabla g(\mathbf{v}_{0})^{\top} (\mathbf{v}^{\star} - \mathbf{v}_{0})\right) \leq 0$. In any iteration besides the first while solving for the $i$th basis direction for block collection $\mathbf{i}$, $\mathbf{v}_0$ is set to the previous iteration's optimum. In the first iteration while solving for the $i$th basis direction for a block collection, we assign to $\mathbf{v}_0$ the $i$th flag mean of the estimated signal matrices of blocks in collection $\mathbf{i}$ from the SVD-based principal angle analysis step of the AJIVE algorithm proposed by \citet{ajive}. The initial condition is also projected in order to obey orthogonality constraints as necessary. 


The exact iterative procedure described in \citep{il2016} requires a feasible solution to execute. Our choice of initial condition will rarely be feasible, and finding a feasible initial condition is challenging in this setting. Therefore, we introduce positive slack variables on the right-hand-side of every constraint and penalize the slack variables in the objective. In particular the slack variables are penalized with a weight $\tau$ which changes on each iteration of the optimization problem, not unlike a learning rate in a machine learning algorithm. Notably, the values of the quadratic forms involved in the object space constraints are often much larger than those for the trait space constraints as the object space constraints include the full-energy data matrices $\mathbf{X}_k$. Therefore, we downweight the slack penalty on those constraints by the leading singular value $\overbar{\nu}_{1,k}$ of $\mathbf{X}_k$ so the optimization problem is not overly restricted by the object space constraints.

Due to the local linearization of angle proximity constraints, we must also enforce that candidate directions $\mathbf{v}^{\star}$ have unit length via the constraints on lines 5 and 6 of \eqref{DivasOptCCSlackLoad}. For further computational efficiency, if the algorithm reaches a point where all the angle-constraint slack variables are zero, it will stop early and add to the current basis a normalized version of the current iteration's intermediate solution. 

The full implementation with penalized slack variables, early stopping conditions, and object space information is displayed in \eqref{DivasOptCCSlackLoad}. The algorithm is implemented in Matlab and solved using CVX, a package for specifying and solving convex programs \cite{cvx, gb08}. 

\begin{equation}\label{DivasOptCCSlackLoad}
\begin{aligned} 
& \min_{\mathbf{v}^{\star}}& - 2\mathbf{v}_0^{\top}\left(\sum_{k\in \mathbf{i}} \check{\mathbf{V}}_k \check{\mathbf{V}}_k^{\top}\right)\mathbf{v}^{\star} + \mathbf{v}_0^{\top}\left(\sum_{k\in \mathbf{i}} \check{\mathbf{V}}_k \check{\mathbf{V}}_k^{\top}\right)\mathbf{v}_0 &+\tau_i\sum_{t=1}^{2K+2} s_t & & \\
& s.t. & \mathbf{v}^{\star T}\mathbf{v}^{\star} - 2 \frac{\mathbf{v}_0^{\top}\check{\mathbf{V}}_k \check{\mathbf{V}}_k^{\top}\mathbf{v}^{\star}}{\cos^2\left(\hat{\phi}_k\right)} + \frac{\mathbf{v}_0^{\top}\check{\mathbf{V}}_k \check{\mathbf{V}}_k^{\top}\mathbf{v}_0}{\cos^2\left(\hat{\phi}_k\right)} &\leq s_k & \forall k &\in \mathbf{i} \\
& & \frac{\mathbf{v}^{\star T}\check{\mathbf{V}}_k \check{\mathbf{V}}_k^{\top}\mathbf{v}^{\star}}{\cos^2\left(\hat{\phi}_k\right)} - 2\mathbf{v}_0^{\top}\mathbf{v}^{\star} + \mathbf{v}_0^{\top}\mathbf{v}_0 &\leq s_k  & \forall k &\in \mathbf{i}^c \\
& & \mathbf{v}^{\star T}\mathbf{X}_k^{\top}\mathbf{X}_k\mathbf{v}^{\star} - 2 \frac{ \mathbf{v}_0^{\top}\mathbf{X}_k^{\top}\check{\mathbf{U}}_k \check{\mathbf{U}}_k^{\top} \mathbf{X}_k \mathbf{v}^{\star}}{\cos^2\left(\hat{\psi}_k\right)} + \frac{\mathbf{v}_0^{\top}\mathbf{X}^{\top}_k\check{\mathbf{U}}_k \check{\mathbf{U}}_k^{\top}\mathbf{X}_k\mathbf{v}_0}{\cos^2\left(\hat{\psi}_k\right)} &\leq s_{K+k}/\overbar{\nu}_{1,k} & \forall k &\in \mathbf{i} \\
& & 1 - 2 \mathbf{v}_0^{\top} \mathbf{v}^{\star} +  \mathbf{v}_0^{\top}\mathbf{v}_0 &\leq s_{2K+1} & & \\
& & \mathbf{v}^{\star T} \mathbf{v}^{\star} - 1 &\leq s_{2K+2} & & \\
& & \mathfrak{V}_{\mathbf{j}}^{\top} \mathbf{v}^{\star} &= {\mathbf{0}} & \forall \mathbf{j} &\supseteq \mathbf{i} \\
\end{aligned}
\end{equation}

Figure~\ref{toyexamplestep2} shows the iterative process of the optimization problem for the synthetic data example displayed in Figure~\ref{3blocktruth}. On each panel, the horizontal axis represents the number of iterations and the vertical axis represents the angles in degrees to the panel's respective estimated signal subspace in trait space. Each horizontal green line represents the trait space perturbation angle bound. The blue paths represent the angles to the low rank approximations of trait spaces at each iteration. The red dashed paths represent the angle to the true trait subspaces at each iteration, which are known since this is a synthetic data set. Note that in the first row, the initial condition is infeasible for the $\mathbf{X}_3=$ X3 angle perturbation bound, demonstrating how the algorithm can use the flexibility afforded by the slack variables to explore the space before choosing a final solution.

\begin{figure}

\begin{center}
\begin{subfigure}[b]{1\textwidth}
	\centering
	\includegraphics[scale=0.3]{X1-X2-X3-joint-optV1.png}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
\centering
	\includegraphics[scale=0.3]{X1-X2-joint-optV1.png}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
\centering
	\includegraphics[scale=0.3]{X1-X3-joint-optV1.png}
\end{subfigure}
\begin{subfigure}[b]{1\textwidth}
\centering
	\includegraphics[scale=0.3]{X2-X3-joint-optV1.png}
\end{subfigure}
\caption{Iterative steps of the sequential optimizations locating joint structure between each possible combination of data blocks for the synthetic data from Figure~\ref{3blocktruth}. The horizontal axes represent iterations of the optimization problem and the vertical axes represent angles in degrees. Colored paths show progression of angles between the candidate direction and $\textbf{TS}(\hat{\mathbf{A}}_k)$ (blue) and $\textbf{TS}(\mathbf{A}_k)$ (red). Perturbation angle bounds $\hat{\phi}_k$ shown as green horizontal lines. From top to bottom: three-way joint, joint between blocks 1 and 2, joint between blocks 1 and 3, joint between blocks 2 and 3.}\label{toyexamplestep2}
\end{center}
\end{figure}

\subsection{Signal Reconstruction}\label{sec:recon}

 
Once we have located all possible joint structure, the remaining task is to reconstruct the signal matrix components for each data block. Recall that in Section~\ref{sec:jse} we denote the estimated orthonormal basis for the joint structure among blocks in collection $\mathbf{i}$ as $\mathfrak{V}_{\mathbf{i}}$. For a given data block $k$, we first horizontally concatenate all joint structure basis matrices found involving block $k$ into one matrix $\left[\mathfrak{V}_{\mathbf{i}}\right]_{\mathbf{i}|k\in \mathbf{i}}$. Then we form a linear regression problem to find the corresponding loadings vectors for block $k$ associated with the common scores vectors for block collection $\mathbf{i}$ in a similar fashion as the proof of Theorem~\ref{th:construct}.
In particular  $[\mathfrak{L}_{\mathbf{i},k}]_{\mathbf{i}|k\in\mathbf{i}}$ is chosen as the least-squares solution of the regression problem $\min_{\mathfrak{L}}\|\mathbf{X}_k - \mathfrak{L}\cdot [\mathfrak{V}_{\mathbf{i}}]^{\top}_{\mathbf{i}|k\in \mathbf{i}}\|_2^2$. This solution is unique when $[\mathfrak{V}_{\mathbf{i}}]_{\mathbf{i}|k\in \mathbf{i}}$  is full-rank. The columns of $[\mathfrak{L}_{\mathbf{i},k}]_{\mathbf{i}|k\in\mathbf{i}}$ can then be partitioned into loadings matrices $\mathfrak{L}_{\mathbf{i},k}$, each of which is associated with joint structure for one block collection $\mathbf{i}$ with $k\in\mathbf{i}$, i.e., $\mathfrak{L}_{\mathbf{i},k}$ is an estimate of $\mathbf{L}_{\mathbf{i},k}$ from \eqref{eq:signalmodel}.  The estimated partially shared joint structure between between data blocks in $\mathbf{i}$ is then simply $\hat{\mathbf{A}}_{\mathbf{i},k}=\mathfrak{L}_{\mathbf{i},k} \mathfrak{V}_{\mathbf{i}}^{\top}$.

As a brief remark, DIVAS may in certain cases select shared structure such that the rank of $\left[\mathfrak{V}_{\mathbf{i}}\right]_{\mathbf{i}|k\in \mathbf{i}}$ is larger than $\overbar{r}_k$. Since the subspaces spanned by $\mathfrak{V}_{\mathbf{i}}$ and $\mathfrak{V}_{\mathbf{j}}$ need not be orthogonal unless $\mathbf{i} \subseteq \mathbf{j}$ or $\mathbf{j} \subseteq \mathbf{i}$, more than $\overbar{r}_k$ total basis directions may be selected within the cone of feasibility for block $k$. 

Additional insight comes from further decomposition of $\hat{\mathbf{A}}_{\mathbf{i},k}$ into a sum of rank 1 modes of variation. Each mode comes from an outer product of corresponding columns of $\mathfrak{L}_{\mathbf{i},k}$ and $\mathfrak{V}_{\mathbf{i}}$. Since $\mathfrak{L}_{\mathbf{i},k}$ and $\mathfrak{V}_{\mathbf{i}}$ are only determined up to basis rotation, we first pick a maximally informative rotation matrix $\mathbf{Q}_{\mathbf{i}}$, and then examine modes of variation formed with the matrices $\mathfrak{L}_{\mathbf{i},k}\mathbf{Q}_{\mathbf{i}}$ and $\mathfrak{V}_{\mathbf{i}}\mathbf{Q}_{\mathbf{i}}$. We take an SVD of the projection of the stacked data matrix $[\mathbf{X}_k^{\top}]^{\top}_{k\in\mathbf{i}}$ onto the subspace spanned by $\mathfrak{V}_{\mathbf{i}}$ in trait space, and choose $\mathbf{Q}_{\mathbf{i}}$ as the matrix of right singular vectors from that calculation. This re-rotation can be thought of as sorting the modes of variation within the shared subspace in order of importance. 

Since DIVAS is based on angles, additional insight into the modes of variation is derived from angles between the loadings (columns of $\mathfrak{L}_{\mathbf{i},k}\mathbf{Q}_{\mathbf{i}})$ and scores (columns of $\mathfrak{V}_{\mathbf{i}}\mathbf{Q}_{\mathbf{i}}$) and the object and trait spaces spanned by the estimated low rank matrix $\hat\mathbf{A}_k$, respectively. In particular, for each estimated score vector and loadings vector, the angle to each estimated signal matrix $\hat{\theta}_k$ and the upper bound on the angle to the true signal $\hat{\theta}_k+\theta_2^{\star}$ for each direction in both trait space and object space are computed. See Section~\ref{sec:perturb} for a definition of $\theta_2^{\star}$. To calculate the upper bound for one of the vectors, we choose the 95th percentile of an empirical distribution of $\theta_2^{\star}$ generated using \eqref{eq:thetatwostar} and the cached matrices from the rotational bootstrap (See Section~\ref{sec:rotboot}). Scores for all blocks are in a shared trait space, and therefore in trait space we calculate these angles not only for included ($k\in\mathbf{i}$) but also for excluded ($k\notin\mathbf{i}$) blocks for each block collection $\mathbf{i}$. The angle to the included block is expected to be small and the angle to the excluded block is expected to be large, though not necessarily 90 degrees. If some score vector has an upper bound below the random direction bound $\theta_0$ for an excluded block, then the corresponding mode of variation is correlated with that excluded block even though it is not joint with that block. The object spaces are block specific, and therefore for loadings we calculate angles to the included ($k\in\mathbf{i}$) data blocks only.  These diagnostic angles form the crux of the overall diagnostic displays that we describe next.

\subsection{DIVAS Diagnostic Graphics}\label{sec:graphics}

We compile all angle-based diagnostics for DIVAS into comprehensive displays. These displays can be seen in Figure~\ref{fig:tcgascoresdiag} for breast cancer omics data, in Figures~\ref{fig:rbbToyExample} and \ref{fig:rbbToyExampleLoad} for the synthetic data set shown in Figure~\ref{3blocktruth}, and in Figures~\ref{fig:mortscoresdiag} and \ref{fig:mortloadingsdiag} for 20th century mortality data. We explain the interpretation of these displays in this section using the synthetic data example.

Figure~\ref{fig:rbbToyExample} shows the diagnostic angles for the joint scores vectors found for the synthetic data example from Figure~\ref{3blocktruth}, and Figure~\ref{fig:rbbToyExampleLoad} shows those same diagnostic angles for the joint loadings vectors. Each row of boxes corresponds to a data block and the various block collections appear in the columns. Boxes for included blocks in a given column are colored-in while boxes for excluded blocks are white. The number in each colored box specifies the rank of the estimated joint subspace between the blocks included in that column. Block collections where no partially shared joint structure was found are labeled with a 0 and grayed out. The last column labeled \textit{Ranks} contains the dimensions of key subspaces for each data block $k=1,\dots , K$. The \textit{final} rank is the dimension of the subspace spanned by all structure involving that data block, i.e. the rank of $\left[\mathfrak{V}_{\mathbf{i}}\right]_{\mathbf{i}|k\in \mathbf{i}}$. The \textit{filtered} rank is the dimension of the estimated signal subspace in both object and trait space for that data block, i.e. $\overbar{r}_k$. The \textit{maximum} rank is the largest possible dimension spanned by structure involving that data block, i.e. $d_k\wedge n$. These three ranks will usually appear in ascending order, but as discussed in Section~\ref{sec:recon}, the final rank is sometimes larger than the filtered rank.

To explain the interpretation of the comprehensive information in DIVAS diagnostic displays, we first focus on the top-left corner of Figure~\ref{fig:rbbToyExample}. Each box is a scatter plot, with the horizontal axis indicating basis direction index and the vertical axis indicating angle from $0^{\circ}$ at the bottom to $90^{\circ}$ at the top. Within a box of this figure, each candidate direction found for that column's joint structure is represented by two points: $\times$ and $\bullet$. The $\times$ represents the angle $\hat{\theta}_k$ between the candidate direction and the corresponding estimated signal matrix for data block $k$, and the $\bullet$ represents the upper bound $\hat{\theta}_k + \hat{\theta}_2^{\star}$ on $\theta_k$, the angle between the direction and the true subspace. The dashed line represents the perturbation angle bound $\hat{\phi}_k$ (for trait space) or $\hat{\psi}_k$ (for object space) and the dot-dash line represents the random direction angle bound $\theta_{0,k}$. The numerical values of those angle bounds are given to the right of each group of columns. As per the inferential framework laid out in Section~\ref{sec:perturb}, a $\times$ below the dashed line indicates strong evidence that the direction can't be ruled out as joint structure for that data block, and a $\bullet$ above the dot-dash line indicates strong statistical evidence that the direction can't be distinguished from an arbitrarily chosen direction with respect to that data block. Due to the rank filtering procedure that takes place during the rotational bootstrap, no direction has both a $\times$ below the dashed line and a $\bullet$ above the dot-dash line.

Based on the placements of $\times$ and $\bullet$ in each colored box in Figures~\ref{fig:rbbToyExample} and \ref{fig:rbbToyExampleLoad}, we have strong evidence that each piece of estimated joint structure located by DIVAS is statistically significant in both trait space and object space, respectively. Specifically, all $\times$ are below the perturbation angle bound dashed line within their respective boxes. We also gain additional insight about the angular relationships between the two-way joint subspaces via the angles to the excluded blocks in Figure~\ref{fig:rbbToyExample}. Each $\bullet$ lies below the random direction angle bound dot-dash line in columns 3-5 of the display, indicating strong evidence that the chosen joint subspaces are distinguishable from arbitrary directions with respect to the excluded block in each column. In fact, the true joint subspaces were constructed to have pairwise principal angles of $60^{\circ}$ between them, so this statistical rejection of arbitrariness is not surprising.

\begin{figure}[h]
\includegraphics[width=\textwidth]{Demo-Scores-Diagnostic-Graphic-11-18.png}
\caption{Summary of joint structure diagnostics for the trait spaces of the synthetic data example. All joint structure located is statistically significant, and angles to excluded blocks confirm underlying angular relationships between two-way shared subspaces. Effective Numbers of Cases (ENC) values in last row align well with the true score vectors.}\label{fig:rbbToyExample}
\end{figure}

\begin{figure}[h]
\includegraphics[width=\textwidth]{Demo-Loadings-Diagnostic-Graphic-11-18.png}
\caption{Summary of joint structure diagnostics for the object spaces of the synthetic example. All joint structure located is statistically significant. Effective Contributions of Traits (ECT) in last row align with expectations per the proportion of colored rows in each heatmap of Figure~\ref{3blocktruth}.}\label{fig:rbbToyExampleLoad}
\end{figure}

There are some situations where DIVAS basis directions appear statistically significant from an angular perspective but depend on a very small number of observations or traits. Useful insight comes from summarizing the contributions of each observation and/or trait to the shared structure. In the case of observations, we quantify their involvement with a summary statistic called the \textit{Effective Number of Cases} (ENC) based on ideas in importance sampling from \citep{kish}. Let $v_j$ for $j\in\{1,\dots, n\}$ be the entries of a chosen direction $\mathbf{v}^{\star}$. Note that the entries are scaled so $\mathbf{v}^{\star}$ has norm 1 (i.e. $\sum_{j=1}^n v_j^2 = 1$). In this case, the ENC is:
\begin{equation*}
ENC = \frac{1}{\sum_{j=1}^n v_j^4}
\end{equation*}
If one entry $v_j$ is $\pm 1$ while the rest are 0, meaning a single observation determines the direction, then the ENC evaluates to 1. If all entries $v_j$ have the same magnitude $\pm\frac{1}{\sqrt{n}}$, meaning all observations have equal influence on the direction, then the ENC evaluates to $n$. Any chosen direction will fall somewhere between those two extremes.

The last row of Figure~\ref{fig:rbbToyExample} shows the ENC for each joint scores direction found for the synthetic data example. Each box is again a scatter plot, with the horizontal axis indicating basis direction index and the vertical axis indicating the ENC from 1 to $n$ on a logarithmic scale. Each ENC value is shown with a $+$. All the values for this example are very close to $n=400$, indicating near equal contribution from each observation in all the scores vectors. This aligns with expectations since the entries of the true shared scores directions all have equal magnitude.

In the case of summarizing individual trait contributions, we use an analogous metric that takes into account the differing magnitudes of loadings vectors within data blocks and the different dimensions of loadings vectors between data blocks. To differentiate the two metrics we call this one \textit{Effective Contribution of Traits} (ECT). ECT performs the same operation as ENC, except it uses the entries $l_m$ for $m\in\{1,\dots, d_k\}$  of candidate loadings directions $\mathbf{l}^{\star}_k$. Furthermore, it scales the result by both the magnitude $||\mathbf{l}_k^{\star}|| = \sum_{m=1}^{d_k} l_m^2$ of the candidate direction and by the number of traits $d_k$ to allow for comparisons between data blocks:
\begin{equation*}
ECT = \frac 1{d_k} \frac{\left(\sum_{m=1}^{d_k} l_m^2\right)^2}{\sum_{m=1}^{d_k} l_m^4}.
\end{equation*}

The last row of Figure~\ref{fig:rbbToyExampleLoad} shows the ECT for each loadings direction found for the synthetic data example. Within each box, the horizontal axis indicates basis direction index and the vertical axis indicates ECT percentage ranging from $0\%$ to $100\%$. Each block has its own loadings direction for each basis direction, so each block's ECT is shown with a number corresponding with that block's index in the analysis. In all cases, the contribution percentages align quite well with the percentage of traits involved in each piece of true joint structure as per Figure~\ref{3blocktruth}. For example, half of the traits in X3 have the characteristic pinstripe pattern of the fully joint structure in Figure~\ref{3blocktruth}, and the ``3" in the bottom-left box of Figure~\ref{fig:rbbToyExampleLoad} sits right around $50\%$.


\section{Case Study: Twentieth Century Mortality}\label{sec:data}








\citet{ooda} consider a data matrix containing mortality rates (proportion of the population of a given age that died in a year) of Spanish males from 1908 to 2002. We expand on this initial analysis by incorporating three additional data blocks: one for Spanish women and two more for Swiss men and women, and by increasing the end of the time frame to 2018. We are interested in how mortality rates changed over the course of this time period as a function of age. Hence, we will treat each year as a data object and the mortality rates of each age as a trait. We consider ages 12 to 90 to avoid zero counts for particularly high and low ages. Data was downloaded on April 8, 2021 from the Human Mortality Database  \citep{mortalitydatabase}.





To appropriately handle the multiple orders of magnitude present in mortality proportions, we transform each entry of the data blocks with a \textit{logit} function $f(x) = \log\left(\frac{x}{1-x}\right)$. After the logit transformation, each data block was \textit{double-centered}: the mean vectors in both object space and trait space were removed from each data block. As per the discussion in \citet{centering}, this type of centering is effective when all the data blocks share a common mode of variation in the trait mean direction. The object and trait means for each data block are displayed as curves in Figure~\ref{fig:mortswissmeans}. Each curve is colored according to year using a rainbow color scheme starting at magenta and blue, through orange and red. In the object mean panels (top row), we see the overall mortality profile across ages for each country and gender. Males exhibit a slightly higher increase in mortality upon entering adulthood than females in both countries due to increased risk-taking behaviors at that age \citep{mortsex2000, mortsex2009}. We also observe systematic anomalies in the mortality rates for older Spanish individuals that are not present in the Swiss data. As discussed in \citet{ooda}, these anomalies are manifestations of an age-rounding effect, and reflect major early differences in demographic record keeping practices between the two countries. The distinct rainbow sequence in each trait mean panel (bottom row) shows steady overall decreases in average mortality rate over time. The worst year of the 20th century flu pandemic, 1918, appears prominently at the top of each trait mean panel in violet. Spanish data block panels (especially males) have out-of-sequence light blue lines in their plots due to a civil war in the late 1930s.

\begin{figure}[h]
\includegraphics[width=\textwidth]{BlockMeanCurves-10-10.png}
\caption{Object means (top row) and trait means (bottom row) for each data block. Both male data blocks have a more dramatic increase in mortality for young adults than the female data blocks. Both Spanish data blocks display effects of record-keeping round-offs absent from Swiss data blocks. All trait means capture the overall improvement in mortality rate over time across the population.}\label{fig:mortswissmeans}
\end{figure}






Figure~\ref{fig:mortscoresdiag} shows the scores diagnostic graphic for the DIVAS decomposition of the mortality data, and Figure~\ref{fig:mortloadingsdiag} shows the angle diagnostics for the corresponding loadings vectors (see Section~\ref{sec:recon}). Since all four data blocks have identical trait and object dimensions and relatively similar variation, all the perturbation angles are also similar to each other. DIVAS finds a two-dimensional four-way shared component, a one-dimensional three-way component shared between each block besides Swiss females, and a six-dimensional shared component between Spanish males and females. Intuitive reasons for these findings are explained below via discussion of the modes of variation.


\begin{figure}[h]
\includegraphics[width=\textwidth]{MortalityDiagnosticScores-11-18.png}
\caption{Joint structure breakdown and diagnostics for the mortality data score vectors. Perhaps surprising is some amount of three-way partially shared joint structure. Spanish men and women have complex two-way partially shared joint structure due to age rounding.}\label{fig:mortscoresdiag}
\end{figure}

\begin{figure}[h]
\includegraphics[width=\textwidth]{MortalityDiagnosticLoadings-11-18.png}
\caption{Joint structure breakdown and diagnostics for the mortality data loadings vectors.}\label{fig:mortloadingsdiag}
\end{figure}





We further investigate the joint structure by visualizing the joint modes of variation about the mean in curve plots throughout Figures~\ref{fig:mortfourwayload}-\ref{fig:mortspanishtwowayload}. Figure~\ref{fig:mortfourwayload} shows such a visualization for the four-way joint structure. In this and subsequent mode of variation figures, each row of panels corresponds to a different basis direction and each column of panels corresponds to a data block. The final plot in each row shows the entries of the common score vector corresponding to that mode of variation. Each row of panels in this figure contains one trend in mortality that was found to be common across both countries and genders. The first mode is a contrast between older and younger individuals that manifests as a change in slope over time. In particular, while mortality rates decreased for all ages over the 20th century as per the trend seen in the trait mean, this mode of variation shows that decrease was more pronounced for younger individuals. The second mode is primarily a contrast between younger adults and middle-age adults that takes place between the 1970s and 1990s, with somewhat different age groupings across blocks. This contrast is the well-documented automotive safety effect described in \citet{ooda}. The wide proliferation of cars in the mid-20th century without modern safety guidelines in this time frame led to a notable increase in automobile fatalities concentrated in younger individuals. As automotive safety improved across Europe in the 1980s and 1990s, this source of excess mortality dissipated.


\begin{figure}[h]
\includegraphics[width=\textwidth]{4way-10-10.png}
\caption{Mode of variation curve plots of the two four-way joint components. Component 1 shows stronger improvement in mortality for younger people, and component 2 shows the rise and fall of automobile fatalities.} \label{fig:mortfourwayload}
\end{figure}

Figure \ref{fig:mortthreewayload} shows the single mode of variation found as three-way shared structure between Swiss men, Spanish men, and Spanish women. Before analyzing this data, we primarily expected to see four-way shared structure and pairwise shared structure across blocks with gender or country in common; this three-way shared structure deviates from that expectation. This mode of variation indicates a contrast between the mortality rates of young adults and the rest of the population during the late 1980s and early 1990s. Considering that the automotive safety effect already appeared in the fully joint component, we suspect this mode of variation is capturing a different phenomenon. Our hypothesis given the time frame and groups affected is that this component is capturing increased mortality from HIV/AIDS in the late 20th century. The contrast focusing on young males in both countries is the main driver of this hypothesis; the corresponding effect in Spanish women seems to be concentrated in older individuals so some additional effect might be entering the mode of variation in that block. This motivates further mortality research.

\begin{figure}[h]
\includegraphics[width=\textwidth]{3way-10-10.png}
\caption{Mode of variation curve plots of three-way joint component between Swiss men, Spanish men, and Spanish women. Mode of variation contrasts mortality in young adults with others around 1990. Potentially related to the emergence of HIV/AIDS.}\label{fig:mortthreewayload}
\end{figure}





Figure~\ref{fig:mortspanishtwowayload} shows the six modes of variation found as two-way joint structure shared between the two Spanish blocks. Component 1 captures excess mortality of young people, and especially young men, during the Spanish Civil War. Component 2 is a contrast within young adults that we do not fully understand. The remaining four modes of variation seem to be harmonic components generated by the age rounding effect discussed in \citet{ooda}.

\begin{figure}[h]
\includegraphics[width=\textwidth]{2way-10-10.png}
\caption{Mode of variation curve plots of the two-way joint components between Spanish data blocks. First component driven by the Spanish Civil War in the 1930s. Last four components are primarily driven by high-frequency-harmonic age-rounding record-keeping anomalies.}\label{fig:mortspanishtwowayload}
\end{figure}

\section{Conclusions}\label{sec:conclusion}

This paper proposes DIVAS, a novel exploratory data analysis method for statistical data integration that allows for partially-shared structure between several distinct data blocks. The main contributions of DIVAS are twofold. First, we develop a rigorous, angle-based framework of statistical inference for diagnostically evaluating estimated shared structure. Second, we consider integration across both dimensions of the data blocks simultaneously and produce more thorough and higher-fidelity results as a consequence.

Future work on DIVAS could proceed in at several different directions. Methodologically, there remains room for additional refinement of noise estimation throughout the first step of DIVAS, both in the noise variance estimator and the residual matrix estimator. Structurally, DIVAS is fundamentally a linear, unsupervised statistical model. Generalizations and expansions that extend DIVAS to tackle supervised learning and nonlinear relationships between data blocks would are promising future directions. Practically, the driving force behind development of the method has always been appropriately complex data like the breast cancer omics data set that demands such methodological sophistication. Therefore we expect further improvements in DIVAS development will be found during analysis of ever more demanding data.

\section{Acknowledgement}
Jan Hannig's research was supported in part by the National Science Foundation under Grant No. DMS-1916115, 2113404, and  2210337. J.S.~Marron's research was partially supported by the National Science Foundation Grant No. DMS-2113404.

\bibliographystyle{apalike}
