The Procrustes problem is aimed at matching matrices using similarity transformations by minimizing their Frobenius distance. It allows comparison of matrices with dimensions defined in an arbitrary coordinate system. This method raised the interest of applied researchers hence highlighting its potentiality through a plethora of applications in several fields, such as ecology (Saito et al., Reference Saito, Fonseca-Gessner and Siqueira2015), biology (Rohlf & Slice, Reference Rohlf and Slice1990), analytical chemometrics (Andrade et al., Reference Andrade, Gómez-Carracedo, Krzanowski and Kubista2004), and psychometrics (Green, Reference Green1952; McCrae et al., Reference McCrae, Zonderman, Costa, Bond and Paunonen1996).
The interest of a large audience from applied fields stimulates, in parallel, the growth of a vast body of literature. Despite this, essentially all applications comprise spatial coordinates (i.e., two- or three-dimensional). Haxby et al. (Reference Haxby, Guntupalli, Connolly, Halchenko, Conroy, Gobbini, Hanke and Ramadge2011) first introduced the use of this approach into a different context: align functional Magnetic Resonance Images (fMRI). The coordinates are hence substituted by voxels (i.e., three-dimensional pixels), and the problem becomes inherently high-dimensional. The approach rapidly grew in popularity in the neuroimaging community because of its effectiveness. However, the proposed solution is naive; the extension from the spatial context to a more general and high-dimensional one is a theoretical challenge that needs adequate attention.
The most serious concern is the results’ interpretability. In most cases, Procrustes methods turn into an ill-posed problem. It is a barely noticeable problem with spatial coordinates because the solution is unique up to rotations; hence, the user has the freedom to choose the point of view that provides the nicest picture. When the dimensions do not have a spatial meaning, any rotation completely changes the interpretation of the results.
To tackle this problem, we revise the perturbation model (Goodal, Reference Goodall1991), which rephrases the Procrustes problem as a statistical model. The matrices are defined as a random perturbation of a reference matrix plus an error term. The perturbation is expressed by rotation, scaling, and translation, and the matrix normal distribution (Gupta & Nagar, Reference Gupta and Nagar2018) is assumed for the error terms. Like Green and Mardia (Reference Mardia, Fallaize, Barber, Jackson and Theobald2013) and Mardia et al. (Reference Green and Mardia2006), we assume that the orthogonal matrix parameter follows the von Mises–Fisher distribution. We prove that the proposed prior distribution is conjugate, making the estimation process quite fast. Indeed, the maximum a posterior estimate of the orthogonal matrix parameters results in a minor modification of the original solution because the prior information enters into the pairwise cross-product of the matrices to be aligned. The prior distribution plays the role of regularizing term, resolving the non-identifiability of the orthogonal matrix parameter. In the application to fMRI data in Sect. 4, we further show that specification of a prior distribution permits the integration of functional and topological aspects, which largely improves the results’ interpretability. We then propose a comprehensive approach to the Procrustes problem: the ProMises (Procrustes von Mises–Fisher) model.
The second problem raised by the extension to the high-dimensional framework is computational. The estimation algorithm of a Procrustes-based method involves a series of singular value decompositions of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m\times m$$\end{document}
matrices (m dimensions). In a typical fMRI data set, the subjects (i.e., the matrices to be aligned) have a few hundred (observations/rows) n and hundreds of thousands of voxels (dimensions/columns) m. We prove that the minimization problem can be solved by a series of singular value decompositions of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n\times n$$\end{document}
matrices, reducing the computation burden and making the ProMises model applicable to matrices with virtually any number of columns. We denote this approach as the Efficient ProMises model.
We emphasize here that the problem of aligning fMRI data is not three-dimensional as it could appear at first glance, although it is high-dimensional: Each voxel is one dimension of the problem. Nevertheless, in such conditions, three critical issues arise: The first is that the Procrustes method combines any voxel inside the brain without distinguishing between adjacent and distant anatomical locations. This can be questionable because the voxels have a spatial organization, and, despite the inter-subject variability, we expect some degree of spatial similarity between subjects in their functional organization. Therefore, we want the subjects’ voxels of a given location to be more likely to contribute to the construction of voxels with the same location in the common space. The second issue revolves around the non-identifiability of orthogonal transformations
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
, where i indexes the subjects. The Procrustes method does not return a unique solution of the maximum likelihood estimate for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
: Given a solution (i.e., a three-dimensional image), any linear combination that mixes the voxels’ values is an equivalent solution to the problem. The solutions are equivalent only from a mathematical point of view because the practical consequence is the loss of results’ topological interpretability. The third issue is the computational load: applying the Procrustes-based alignment to the whole brain implies the decomposition of many square matrices of dimensions roughly equal to 200.000 (i.e., the number of voxels.)
The Efficient ProMises model resolves all these three issues. The use of a properly chosen prior shrinks the estimate to the anatomical solution (i.e., no rotation), hence making the solution unique and interpretable from an anatomical point of view. Finally, as mentioned before, the Efficient implementation permits performing functional alignment on high-dimensional data such as fMRI data.
The paper is organized as follows. Section 1 introduces the perturbation model (Goodall, Reference Goodall1991), stressing its critical issues. Section 2 illustrates the ProMises model and its challenges: identifiability, interpretability, and flexibility. Section 3 defines the Efficient version of the ProMises model, which permits application of the functional alignment to high-dimensional data. Finally, the presented model is evaluated by analyzing task-related fMRI data in Sect. 4. The entire code used is available in https://github.com/angeella/ProMisesModel using the programming language Python (Van Rossum and Drake Jr, Reference Van Rossum and Drake1995) and in https://github.com/angeella/alignProMises using the R (R Core Team, 2018) package alignProMises. We report the proofs of the main lemmas here, whereas the remaining proofs are included in the supplementary material.
1. Perturbation Model
1.1. Background
Let
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_{i} \in \mathrm{I\!R}^{n \times m} \}_{i = 1,\dots ,N}$$\end{document}
be a set of matrices to be aligned. The Procrustes-based method uses similarity transformations to match each matrix to the target one as closely as possible, according to the Frobenius distance.
Each matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
could be then assumed to be a similarity transformation of a shared matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}} \in \mathrm{I\!R}^{n \times m}$$\end{document}
, which contains the common reference space’s coordinates, plus a random error matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{E}}_i \in \mathrm{I\!R}^{n \times m}$$\end{document}
. The perturbation model proposed by Goodall (Reference Goodall1991) is then reported.
Definition 1
(Goodall, Reference Goodall1991) Let
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_{i} \in \mathrm{I\!R}^{n \times m} \}_{i = 1,\dots ,N}$$\end{document}
be a set of matrices to be aligned and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {O}}(m)$$\end{document}
the orthogonal group in dimension m. The perturbation model is defined as
where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{E}}_i \sim \mathcal{MN}\mathcal{}_{n,m}(0,\varvec{\Sigma }_n,\varvec{\Sigma }_m)$$\end{document}
—i.e., the matrix normal distribution with
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n \in \mathrm{I\!R}^{n \times n}$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m \in \mathrm{I\!R}^{m \times m}$$\end{document}
scale parameters—
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}} \in \mathrm{I\!R}^{n \times m}$$\end{document}
is the shared matrix,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _i \in \mathrm{I\!R}^{+}$$\end{document}
is the isotropic scaling,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{t}}_i \in \mathrm{I\!R}^{1 \times m}$$\end{document}
defines the translation vector, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbf {1}}_n \in \mathrm{I\!R}^{1 \times n}$$\end{document}
is a vector of ones.
To simplify the problem, the column-centered
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
are considered. The distribution is
where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{C}}_n = {\varvec{I}}_n - \frac{1}{n} {\varvec{J}}_n$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{I}}_n \in \mathrm{I\!R}^{n \times n}$$\end{document}
is the identity matrix,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{J}}_n$$\end{document}
is a
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n \times n$$\end{document}
matrix of ones, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\text {vec}({\varvec{A}})$$\end{document}
the vectorization of the matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{A}}$$\end{document}
. Let the singular value decomposition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{C}}_n = \varvec{\Gamma } \varvec{\Delta } \varvec{\Gamma }^\top $$\end{document}
, where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Gamma } \in \mathrm{I\!R}^{n \times (n-1)}$$\end{document}
, then:
The
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Gamma }^\top $$\end{document}
transformation leads to independence between the rows of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{E}}_i$$\end{document}
without affecting the estimation of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
.
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Gamma }^\top {\varvec{C}}_n {\varvec{X}}_i$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Gamma }^\top {\varvec{C}}_n {\varvec{M}}$$\end{document}
have now
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n-1$$\end{document}
rows; however, we can simply re-project them on
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mathrm{I\!R}^{n \times m}$$\end{document}
using
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Gamma }$$\end{document}
. Because we can always write
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tilde{{\varvec{X}}}_i = \varvec{\Gamma } \varvec{\Gamma }^\top {\varvec{C}}_n {\varvec{X}}_i$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tilde{{\varvec{M}}} = \varvec{\Gamma } \varvec{\Gamma }^\top {\varvec{C}}_n {\varvec{M}}$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\tilde{\varvec{\Sigma }}_n = \varvec{\Gamma }\varvec{\Gamma }^\top {\varvec{C}}_n \varvec{\Sigma }_n {\varvec{C}}_n^\top \varvec{\Gamma }\varvec{\Gamma }^\top $$\end{document}
without loss of generality, we can re-write Definition 1 as follows:
Definition 2
where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{E}}_i \sim \mathcal{MN}\mathcal{}_{n m}(0,\varvec{\Sigma }_n,\varvec{\Sigma }_m)$$\end{document}
. This way, we obtain the following:
The following notation is also adopted:
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$||\cdot ||$$\end{document}
to indicate the Frobenius norm and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$<\cdot ,\cdot>$$\end{document}
for the Frobenius inner product (Golub & van Loan, Reference Golub and Van Loan2013).
The main objective of this work is comparing the shapes
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
instead of the form’s analysis of the matrices. For that, the parameters of interest are
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _i$$\end{document}
, whereas
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
are considered as nuisance parameters for each
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1, \dots , N$$\end{document}
. The estimation of the unknown parameters changes if these nuisance parameters are known. Section 1.2 initially presents the estimates under this assumption and then provides the more realistic case of unknown nuisance parameters.
1.2. Estimation of the Perturbation Model
We formalize some results from Theobald and Wuttke (Reference Theobald and Wuttke2006) in the case of known nuisance parameters
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
, with
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
positive definite matrices by the following theorem:
Theorem 1
(Theobald & Wuttke Reference Theobald and Wuttke2006) Consider the perturbation model described in Definition 2, and the singular value decomposition
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1}= {\varvec{U}}_i {\varvec{D}}_i {\varvec{V}}_i^\top $$\end{document}
. The maximum likelihood estimators equal
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{{\varvec{R}}}_i = {\varvec{U}}_i {\varvec{V}}_i^\top $$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{\alpha _i}_{\varvec{{\hat{R}}}_i} = ||\varvec{\Sigma }_m^{-1/2} \hat{{\varvec{R}}}_i^\top {\varvec{X}}_i^\top \varvec{\Sigma _n}^{-1/2}||^2/tr{({\varvec{D}}_i)}$$\end{document}
.
Now consider N independent observations
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_1, \dots , {\varvec{X}}_N$$\end{document}
. The joint log-likelihood is simply the sum of N log-likelihoods.
In the case of unknown nuisance parameters
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
, the joint likelihood cannot be written as product of separated likelihoods, one for each
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
, because each of the unknown parameters is a function of the others. The solution must be found by an iterative algorithm. In particular, the two covariance matrices
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
can be estimated by a two-stage algorithm defined in Dutilleul (Reference Dutilleul1999), where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{\Sigma }}}_n = \{\sum _{i=1}^{N}({\varvec{X}}_i - \varvec{{\hat{M}}}) \varvec{{\hat{\Sigma }}}_m^{-1} ({\varvec{X}}_i - \varvec{{\hat{M}}})^\top \}/Nm$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{\Sigma }}}_m = \{\sum _{i=1}^{N}({\varvec{X}}_i - \varvec{{\hat{M}}})^\top \varvec{{\hat{\Sigma }}}_n^{-1} ({\varvec{X}}_i - \varvec{{\hat{M}}})\}/Nn$$\end{document}
are maximum likelihood estimators.
The necessary and sufficient condition for the existence of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{\Sigma }}}_n$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{\Sigma }}}_m$$\end{document}
is
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N \ge \frac{m}{n} + 1$$\end{document}
, assuming
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
are positive definite matrices. In real applications, this assumption could be problematic. For example, in fMRI data analysis, m roughly equals 200, 000, and n approximately equals 200; therefore, the researcher would have to analyzing at least 1, 001 subjects, which is virtually impossible because fMRI is costly.
Various solutions can be found in the literature: Theobald and Wuttke (Reference Theobald and Wuttke2006) proposed a regularization for the covariance matrix, whereas Lele (Reference Lele1993) estimated
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
using the distribution of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i {\varvec{X}}_i^\top $$\end{document}
. In this work, we maintain a general formulation of the estimator for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{{\hat{\Sigma }}}}_m = g(\varvec{{\hat{\Sigma }}}_n, \varvec{{\hat{M}}}, {\varvec{X}}_i)$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{{\hat{\Sigma }}}}_n = g(\varvec{{\hat{\Sigma }}}_m, \varvec{{\hat{M}}}, {\varvec{X}}_i)$$\end{document}
, because we aim to find a proper estimate of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
rather than
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
. The shared matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
is estimated by the element-wise arithmetic mean of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\hat{{\varvec{X}}}_i\}_{i = 1, \dots , N}$$\end{document}
, where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{{\varvec{X}}}_i = {\hat{\alpha }}_{\hat{{\varvec{R}}}_i}^{-1} {\varvec{X}}_i \hat{{\varvec{R}}}_i$$\end{document}
. We then modified the iterative algorithm of Gower (Reference Gower1975) (i.e., generalized Procrustes analysis), to estimate
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
.
Groisser (Reference Groisser2005) proved the convergence of the generalized Procrustes analysis algorithm. However, it leads to non-identifiable estimators of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1, \dots , N$$\end{document}
, proved by the following lemma:
Lemma 1
Let
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\varvec{{\hat{R}}}_i\}_{i=1,\dots ,N}$$\end{document}
be the maximum likelihood solutions for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{R}}_i\}_{i=1,\dots ,N}$$\end{document}
with
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
as unknown parameters. If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Z}} \in {\mathcal {O}}(m)$$\end{document}
, then
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\varvec{{\hat{R}}}_i {\varvec{Z}}\}_{i = 1, \dots , N}$$\end{document}
are still valid maximum likelihood solutions for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\varvec{{\hat{R}}}_i\}_{i = 1, \dots , N}$$\end{document}
.
Proof.
Consider
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Z}} \in {\mathcal {O}}(m)$$\end{document}
, and
Consider the proof of Theorem 1 placed in the supplementary material, we have
Since
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Z}}^\top {\varvec{Z}} = {\varvec{Z}} {\varvec{Z}}^\top = {\varvec{I}}_m$$\end{document}
, the solutions
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{\hat{{\varvec{R}}}_i {\varvec{Z}}\}_{i = 1, \dots , N}$$\end{document}
are still valid solutions for the maximization (3).
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\square $$\end{document}
To sum up, in the more realistic case, when we must estimate the nuisance parameters, the Procrustes solutions are infinite in general, in both in the high-dimensional case and the low-dimensional by Lemma 1. We emphasize here that, to resolve the non-identifiability of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
, the proposed ProMises model imposes a prior distribution for the parameter
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
.
2. ProMises Model
2.1. Background
In the previous section, we justified how the perturbation model could be problematic because Lemma 1 proves the non-identifiability of the parameter
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
in the realistic case of unknown nuisance parameters. This scenario returns to be critical in several applications, where the m columns of the matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
do not express the three-dimensional spatial coordinates, such as in the fMRI data framework illustrated in Sect. 4. In the high-dimensional case, the final orientations of the aligned data can be relevant for interpretation purposes.
For that, we propose its Bayesian approach: the ProMises model. We stress here that a proper prior distribution leads to a closed-form and interpretable, unique point estimate of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
. The specification of the prior parameters is essential, especially in our high-dimensional context.
2.2. Interpretation of the Prior Parameters
Because
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i \in {\mathcal {O}}(m)$$\end{document}
, a proper prior distribution must take values in the Stiefel manifold
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_{m}(\mathrm{I\!R}^m)$$\end{document}
. The matrix von Mises–Fisher distribution is a non-uniform distribution on
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$V_{m}(\mathrm{I\!R}^m)$$\end{document}
, which describes a rigid configuration of m distinct directions with fixed angles. It was proposed by Downs (Reference Downs1972) and investigated by many authors (e.g., Chikuse, Reference Jupp and Mardia1979; Jupp & Mardia Reference Chikuse2003).
We report below the formal definition of the von Mises–Fisher distribution.
Definition 3
(Downs, Reference Downs1972) The von Mises–Fisher distribution for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i \in {\mathcal {O}}(m)$$\end{document}
is
where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$C({\varvec{F}}, k)$$\end{document}
is a normalizing constant,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}} \in \mathrm{I\!R}^{m \times m}$$\end{document}
is the location matrix parameter, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k \in \mathrm{I\!R}^{+}$$\end{document}
is the concentration parameter.
The parameter k defined in (4) balances the amount of concentration of the distribution around
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
. If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k \rightarrow 0$$\end{document}
, the prior distribution is near a uniform distribution (i.e., unconstrained). If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k \rightarrow +\infty $$\end{document}
, the prior distribution tends toward a Dirac distribution (i.e., maximum constraint).
A proper specification of the prior distribution leads to improved estimation of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
. Therefore, the core of the ProMises model is the specification of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
defined in (4). Consider the polar decomposition and singular value decomposition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}= {\varvec{P}} {\varvec{K}} = {\varvec{L}} \varvec{\Sigma }{\varvec{B}}^\top = {\varvec{L}} {\varvec{B}}^\top {\varvec{B}} \varvec{\Sigma } {\varvec{B}}^\top $$\end{document}
, where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{P}}, {\varvec{L}}, {\varvec{B}} \in {\mathcal {O}}(m)$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{K}} \in \mathrm{I\!R}^{m \times m}$$\end{document}
are symmetric positive semi-definite matrices, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma } \in \mathrm{I\!R}^{m \times m}$$\end{document}
diagonal matrix with non-negative real numbers on the diagonal. The mode of the density defined in (4) equals
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{P}}$$\end{document}
(Jupp & Mardia, Reference Jupp and Mardia1979), so the most plausible rotation matrix depends on the orientation characteristic of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
. Merging the two decompositions,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{P}} = {\varvec{L}} {\varvec{B}}^\top $$\end{document}
describes the orientation part of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{K}} = {\varvec{B}} \varvec{\Sigma } {\varvec{B}}^\top $$\end{document}
defines the concentration part. The mode is specified by the product of the left and right singular vectors of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
. These decompositions are useful to understand when the density (4) is uni-modal. If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
has full rank,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }$$\end{document}
does too, the polar decomposition is unique, and thus the mode of the density (i.e.,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{P}}$$\end{document}
is the global maximum). Let
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
be a full rank matrix, then the maximum equals
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\max _{{\varvec{R}}_i \in {\mathcal {O}}(m)} \,\, tr({\varvec{F}} {\varvec{R}}_i^\top ) = tr\{{\varvec{L}} {\varvec{B}}^\top {\varvec{B}} \varvec{\Sigma } {\varvec{B}}^\top ({\varvec{L}} {\varvec{B}}^\top )^\top \} = tr(\varvec{\Sigma })$$\end{document}
.
To sum up, the prior specification allows us to include a priori information about the optimal orientation in the perturbation model. We anticipate here the result of Lemma 3: If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
is defined as a full rank matrix, the maximum a posteriori solution
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{R}}^{'}}_i$$\end{document}
will be unique with the orientation structure of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
.
Two simple examples of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
are delineated below.
Example 1
The most simple definition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
is
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{I}}_m$$\end{document}
(Lee, Reference Lee2018). The eigenvalues are all 1, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{L}}$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{B}}$$\end{document}
are equal to
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{e}}_1,\dots , {\varvec{e}}_m$$\end{document}
, where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{e}}_i$$\end{document}
is the standard basis forming an orthonormal basis of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\mathrm{I\!R}^{m}$$\end{document}
. The prior distribution shrinks the possible solutions for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
toward orthogonal matrices that consider only the combination of variables with the same location.
Alternatively, considering the fMRI scenario, the hyperparameter
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
can be defined as an Euclidean similarity matrix using the 3D anatomical coordinates x, y, and z of each voxel:
where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i,j = 1,\dots m$$\end{document}
. In this way,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
is a symmetric matrix with ones in the diagonal, which means that voxels with the same spatial location are combined with weights equalling 1, and the weights decrease as the voxels to be combined are more spatially distant.
Example 2
Consider N matrices, one for each plant, describing the three-dimensional spatial trajectories of a climbing plant, Pisum sativum, having wooden support as a stimulus (Guerra et al., Reference Guerra, Peressotti, Peressotti, Bulgheroni, Baccinelli and D’Amico2019) across time. The spatiotemporal trajectories of the plants are analyzed until they come to grasp the stick. The aim is to functionally align the three time series, one for each coordinate (x, y, z); then, we have
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m=3$$\end{document}
. In this case, we could suppose that the rotation along the z axis is not of interest because the functional misalignment between plants can be along the x-axis and y-axis with the z-axis being the one that reflects the growth of the plants and the x-axis and y-axis describing the elliptical movement (circumnutation) of the plants (Guerra et al., Reference Guerra, Peressotti, Peressotti, Bulgheroni, Baccinelli and D’Amico2019). Therefore, the
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
location matrix parameter can be described as follows:
The axes x and y have the same probability of entering in the calculation of the first two dimensions of the final common space, whereas the z axis is not considered.
We use the kinematic plant data from Guerra et al. (Reference Guerra, Peressotti, Peressotti, Bulgheroni, Baccinelli and D’Amico2019), consisting of five matrices/plants. Figure 1 shows the elliptical movement expressed by the axes x and y, in the case of unaligned and aligned plant trajectories. The rotation transformations are estimated by the ProMises model with
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
expressed as (5). We do not go into detail about the meaning of the results because we have introduced this example to explain the usefulness of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
. However, we can note how the ProMises model aligns the final coordinates of the tendrils (i.e., when the plant touches the wooden support).
Left panel: Unaligned spatial trajectories of the tendrils of two plants. Right panel: Aligned spatial trajectories of the tendrils of two plants.
2.3. Von Mises–Fisher Conjugate Prior
The von Mises–Fisher distribution (4) was proved by Khatri and Mardia (Reference Khatri and Mardia1977) to be a member of the standard exponential family (Barndorff–Nielsen, Reference Barndorff-Nielsen2014). Green and Mardia (Reference Green and Mardia2006) mentioned that the von Mises–Fisher distribution is a conjugate prior for the matrix normal distribution, which we formally prove in the following lemma under the perturbation model’s assumptions.
Lemma 2
Consider the perturbation model of Definition 2, with
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
distributed according to (4), then the posterior distribution
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f({\varvec{R}}_i| k, {\varvec{F}}, {\varvec{X}}_i)$$\end{document}
is a conjugate distribution to the von Mises–Fisher prior distribution with location posterior parameter equalling the following:
The posterior location parameter is the sum of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1}$$\end{document}
and the prior location parameter
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
multiplied by k. Consider the singular value decomposition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1}$$\end{document}
:
The right part of (7)
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{V}}_i {\varvec{D}}_i {\varvec{V}}_i^\top $$\end{document}
is the elliptical part of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1}$$\end{document}
, which is a measure of variation relative to the decomposition
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{U}}_i {\varvec{V}}_i^\top $$\end{document}
(i.e., the maximum likelihood estimator of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
). Focus on the right part of (6),
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}={\varvec{P}} {\varvec{K}}$$\end{document}
, which is the polar decomposition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
, where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{P}}$$\end{document}
is the mode of the von Mises–Fisher distribution, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{K}}$$\end{document}
is its measure of variation. Therefore,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}^\star $$\end{document}
is expressed as a combination of the maximum likelihood estimate
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{R}}}_i = {\varvec{U}}_i {\varvec{V}}_i^\top $$\end{document}
and the prior mode
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{P}}$$\end{document}
, multiplied by corresponding measures of variation.
Thanks to the conjugacy, the estimation process remains simple; with a small modification we keep the previous algorithm without increasing the computational burden.
2.4. Estimation of the ProMises Model
This section delineates the estimation process for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
using the ProMises method. First of all,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f(\varvec{X_i} | \alpha _i, \varvec{R_i})$$\end{document}
depends only on the product
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _i \varvec{R_i}$$\end{document}
, we thus refer to the distribution
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f(\varvec{X_i} | \alpha _i \varvec{R_i})$$\end{document}
instead of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$f(\varvec{X_i} | \alpha _i, \varvec{R_i})$$\end{document}
defined in (2). The following density is then considered as prior distribution for the product
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _i \varvec{R_i}$$\end{document}
:
The following theorem delineates the estimation of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
with known nuisance parameters:
Theorem 2
The ProMises model is defined as the perturbation model specified in Definition 2 imposing the prior distribution (8) for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\alpha _i {\varvec{R}}_i$$\end{document}
. Let the singular value decomposition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1} + k {\varvec{F}}$$\end{document}
be
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{U}}_i {\varvec{D}}_i {\varvec{V}}_i^\top $$\end{document}
. Then, the maximum a posteriori estimators equal
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{R}}^{'}}_i = {\varvec{U}}_i {\varvec{V}}_i^\top $$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{\alpha _i}_{\varvec{{\hat{R}}^{'}}_i}^{'}=||\varvec{\Sigma }_m^{-1/2} \varvec{{\hat{R}}^{'\top }}_i {\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1/2}||^2/tr({\varvec{D}}_i)$$\end{document}
.
The prior information about
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
’s structure is directly entered in the singular value decomposition step; the maximum a posteriori estimator turns out to be a slight modification of the solution given in Theorem 1. We decompose
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1} + k {\varvec{F}}$$\end{document}
instead of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1}$$\end{document}
.
Let
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_{i} \in \mathrm{I\!R}^{n \times m} \}_{i = 1,\dots ,N}$$\end{document}
be a set of independent matrices. Then, the joint posterior distribution is simply the product of the single posterior distribution.
If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
are unknown, the maximization problem has no closed-form solution, like in Sect. 1.2. Because we proved that the prior specification modifies only the singular value decomposition step of Theorem 1, we then modify the Line 5 of Algorithm 1, as follows:
2.5. On the Choice of the Parameter of the von Mises–Fisher Distribution
We choose the von Mises–Fisher distribution as prior distribution for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
for its useful and practical properties. First, as shown, it is a conjugate prior distribution, leading to a direct calculation and interpretation of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\hat{R}}^{'}}_i$$\end{document}
. Second, it expresses the orthogonality constraint imposed by the Procrustes problem. Finally, the definition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
does not require strong assumptions (Downs, Reference Downs1972); nevertheless, if we specify it as a full-rank matrix, we guarantee solution’s uniqueness. This permits formulation of the below lemma.
Lemma 3
If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
has full rank, the maximum a posteriori estimates for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
given by Theorem 2 are unique.
Proof.
Consider the proof of Lemma 1. Multiplying by
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Z}}$$\end{document}
leads to the following maximization:
since the cyclic permutation invariance property of the trace does not work as Lemma 1 having the additional term
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$k tr( {\varvec{F}} {\varvec{R}}_i)$$\end{document}
.
In addition, recalling Lemma 4, the solution for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
is unique if and only if
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i^\top \varvec{\Sigma }_n^{-1} {\varvec{M}} \varvec{\Sigma }_m^{-1}+ k {\varvec{F}}$$\end{document}
has full rank. If
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
is defined with full rank, so
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{{\tilde{X}}}_i^\top \varvec{{\tilde{M}}} + k {\varvec{F}}$$\end{document}
, and the solution for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$R_i$$\end{document}
is unique. Furthermore, recalling Jupp and Mardia (Reference Jupp and Mardia1979), the mode of the von Mises–Fisher is the orientation part of location matrix parameter. Because the polar decomposition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}^\star $$\end{document}
is unique, the maximum a posteriori estimate is unique.
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\square $$\end{document}
To sum up, the ProMises model enables resolving the non-identifiability of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{R_i}$$\end{document}
that characterizes the perturbation model. The prior information inserted in the model permits guidance of the estimation process, computing a unique and interpretable data orthogonal transformation. Finally, all these properties are reached without complicating the estimation process of the perturbation model; we only modify the singular value decomposition step of Algorithm 1.
3. Efficient ProMises Model
The framework depicted above can be applied both in low- and high-dimensional settings. However, the extension to the high-dimensional case does not come for free if the perturbation model is used. When
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n<m$$\end{document}
, rank equals n, the identifiability of the solution is lost even when the nuisance parameters are known. The lemma below formally states:
Lemma 4
Consider
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i \in \mathrm{I\!R}^{n \times m}$$\end{document}
, if
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n < m$$\end{document}
, then the maximum likelihood estimate for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
defined in Theorem 1 is not unique.
Although the ProMises model provides unique solutions even in high-dimensional frameworks, a second issue remains prominent: the computational load. At each step, the presented algorithms perform N singular values decompositions of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m\times m$$\end{document}
matrices, which have a polynomial-time complexity
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$O(m^3)$$\end{document}
. When m becomes large, as in fMRI data where m is a few hundred thousands, the computation runtime, and the required storing memory, becomes inadmissible.
This section proposes the Efficient ProMises model, which resolves the two above points. The method is efficient in terms of space and time complexity and fixes the non-identifiability of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
. The algorithm allows a faster and more accessible shape analysis without loss of information in the case of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n \ll m$$\end{document}
. It essentially merges the thin singular value decomposition (Bai et al., Reference Bai, Demmel, Dongarra, Ruhe and van der Vorst2000) with the Procrustes problem.
In practice, the Efficient ProMises approach projects the matrices
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
into an n-lower-dimensional space using a specific semi-orthogonal transformation (Abadir & Magnus, Reference Abadir and Magnus2005; Groß et al., Reference Groß, Trenkler and Troschke1999)
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Q}}_i$$\end{document}
, with dimensions
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m \times n$$\end{document}
, which preserve all the data’s information. It aligns, then, the reduced
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n\times n$$\end{document}
matrices
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_i {\varvec{Q}}_i \in \mathrm{I\!R}^{n \times n }\}_{i = 1, \dots , N}$$\end{document}
by the perturbation or ProMises model. Finally, it projects the aligned matrices back to the original
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n\times m$$\end{document}
-size matrices
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_i \in \mathrm{I\!R}^{n \times m }\}_{i = 1, \dots , N}$$\end{document}
using the transpose of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{Q}}_i\}_{i = 1, \dots , N}$$\end{document}
.
The following theorem proves that the maximum defined in Eq. (3) using
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_i {\varvec{Q}}_i \in \mathrm{I\!R}^{n \times n }\}_{i = 1, \dots , N}$$\end{document}
equals the original maximum because the Procrustes problem analyzes the first
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n \times n$$\end{document}
dimensions of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
. The maximum remains the same if we multiply
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\{{\varvec{X}}_i \in \mathrm{I\!R}^{n \times m }\}_{i = 1, \dots , N} $$\end{document}
by
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Q}}_i$$\end{document}
.
Theorem 3
Consider the perturbation model in Definition 2 with
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m = \sigma ^2 {\varvec{I}}_m$$\end{document}
and the thin singular value decompositions of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i = {\varvec{L}}_i {\varvec{S}}_i {\varvec{Q}}_i^\top $$\end{document}
for each
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$i = 1, \dots , N$$\end{document}
, where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Q}}_i$$\end{document}
has dimensions
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m \times n$$\end{document}
. The following holds
Additionally, the condition for the existence of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\hat{\varvec{\Sigma }}_n$$\end{document}
by Dutilleul (Reference Dutilleul1999) is satisfied because
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
now has dimensions
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n \times n$$\end{document}
, and the functional alignment needs at least
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$N = 2$$\end{document}
observations.
So, Theorem 3 is used to define an Efficient version of the ProMises model.
Lemma 5
Consider the assumptions of Theorem 3, then
where
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{F}} \in \mathrm {I\!R}^{m \times m}$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${{F}}^\star = {{Q}}_i^\top {{F}} {{Q}}_j \in \mathrm {I\!R}^{n \times n}$$\end{document}
.
Proofs of Theorem 3 and Lemma 5 are shown in the supplementary material, whereas here, we make some further considerations about the proposed method.
At first glance, the assumption
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m = \sigma ^2 {\varvec{I}}_m$$\end{document}
may dilute the resul’s impact. However, this assumption does not imply that the data are column-wise independent because this dependence is modeled by
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
. Additionally, the joint and accurate estimate of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i$$\end{document}
requires a large number of observations. So, it is common in real applications to set
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_m$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$\varvec{\Sigma }_n$$\end{document}
to be proportional to the identities in Procrustes-like problems (Haxby et al., Reference Haxby, Guntupalli, Nastase and Feilong2020). When the model is high-dimensional, this problem becomes even more pronounced because of the huge number of parameters to be estimated.
The Efficient ProMises approach reaches the same maximum while working in the reduced space of the first n eigenvectors, which contains all the information, instead of the full data set. Therefore, the original problem estimates orthogonal matrices of size
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m\times m$$\end{document}
:
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i \in {\mathcal {O}}(m)$$\end{document}
, whereas the Efficient solution provides a set of orthogonal matrices of size
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n\times n$$\end{document}
:
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{R}}_i^{\star } \in {\mathcal {O}}(n)$$\end{document}
. Even when the solution is projected back into the
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m\times m$$\end{document}
space through
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{Q}}_i{\varvec{R}}_i^{\star }{\varvec{Q}}_i^\top $$\end{document}
, the rank remains n, whereas the matrices of the original solutions have rank m. This should clarify that the Efficient approach reaches the same fit to the data under a different set of constraints that is
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n\times n$$\end{document}
orthogonal matrices instead of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$m\times m$$\end{document}
matrices; hence, the solutions of the two algorithms will not be identical.
Then, we add on Algorithm 1 the lines used to reduce the dimensions of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{X}}_i$$\end{document}
, and we modify Line 5 to insert the prior information:
The Efficient approach reduces the time complexity from
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {O}}(m^3)$$\end{document}
to
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {O}}(m n^2)$$\end{document}
and the space complexity from
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {O}}(m^2)$$\end{document}
to
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathcal {O}}(m n)$$\end{document}
.
4. Functional Magnetic Resonance Imaging Data Application
4.1. Motivation
The alignment problem is recognized in fMRI multi-subject studies because the brain’s anatomical and functional structures vary across subjects. The most used anatomical alignments are the Talairach normalization (Talairach & Tournoux, Reference Talairach and Tournoux1988) and the Montréal Neurological Institute (MNI) space normalization (Jenkinson et al., Reference Jenkinson, Bannister, Brady and Smith2002), where the brain images are aligned to an anatomical template by affine transformations using a set of major anatomical landmarks. However, this alignment does not explore the between-subjects variability in anatomical positions of the functional loci. The functional brain regions are not consistently placed on the anatomical landmarks defined by the Talairach and MNI templates. The anatomical alignment is then an approximate inter-subject registration of the functional cortical areas. Haxby et al. (Reference Haxby, Guntupalli, Connolly, Halchenko, Conroy, Gobbini, Hanke and Ramadge2011) proved that functional brain anatomy exhibits a regular organization at a fine spatial scale shared across subjects.
Therefore, we can assume that anatomical and functional structures are subject-specific (Conroy et al., Reference Conroy, Singer, Haxby and Ramadge2009; Sabuncu et al., Reference Sabuncu, Singer, Conroy, Bryan, Ramadge and Haxby2010) and that the neural activities in different brains are noisy rotations of a common space (Haxby et al., Reference Haxby, Guntupalli, Connolly, Halchenko, Conroy, Gobbini, Hanke and Ramadge2011). Functional alignments (e.g., Procrustes methods) attempt to rotate the neural activities to maximize similarity across subjects.
Specifically, each subject’s brain activation can be represented by a matrix, where the rows represent the stimuli/time points, and the columns represent the voxels. The stimuli are time-synchronized among subjects, so we have correspondence among the matrices’ rows. However, the columns are not assumed to be in correspondence among subjects, as explained before. Each time series of brain activation (i.e., each of the matrices’ columns) represents the voxels’ functional characteristics that the anatomical normalization fails to align. We aim to represent the neural responses to stimuli into a common high-dimensional space, rather than in a canonical anatomical space that does not consider the variability of functional topographies loci.
Figure 2 shows three voxels’ neural activities (i.e.,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_1$$\end{document}
,
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_2$$\end{document}
, and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_3$$\end{document}
) three columns of the data matrix in two subjects recorded across time. The functional pattern of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_2$$\end{document}
is equal across subjects, whereas
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_1$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_3$$\end{document}
are swapped. A rotation matrix can resolve this misalignment, with the swap being a particular case of the rotation matrix. For further details about the motivation in using functional Procrustes-based alignment in fMRI studies, see Haxby et al. (Reference Haxby, Guntupalli, Nastase and Feilong2020).
Illustration of functional misalignment between fMRI images, where three voxels’ time series are plotted considering two subjects. The time series of voxels
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_1$$\end{document}
and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$v_3$$\end{document}
of the second subject are swapped with respect to the first subject.
4.2. Data Description
We apply the proposed method to data from Pernet et al. (Reference Pernet, McAleer, Latinus, Gorgolewski, Charest, Bestelmeyer, Watson, Fleming, Crabbe, Valdes-Sosa and Belin2015), available at https://openneuro.org/datasets/ds000158/versions/1.0.0. The study consists of neural activations of 218 subjects passively listening to vocal (i.e., speech) and nonvocal sounds. Because the application has had a mere illustrative purpose, we choose to use a small number of subjects (18) to facilitate the example’s reproducibility by the readers. We preprocessed the data using the Functional MRI of the Brain Software Library (FSL) (Jenkinson et al., Reference Jenkinson, Beckmann, Behrens, Woolrich and Smith2012) using a standard processing procedure (i.e., high-pass filtering, brain extraction, spatially smoothing, registration to standard MNI space, dealing with motion and differences in slice acquisition time). Anatomical and functional alignment (based on the ProMises model) is compared, having images preprocessed in the same way, but in one case, the functional alignment is applied, while in the other case not. For details about the experimental design and data acquisition, please see Pernet et al. (Reference Pernet, McAleer, Latinus, Gorgolewski, Charest, Bestelmeyer, Watson, Fleming, Crabbe, Valdes-Sosa and Belin2015).
4.3. Functional Connectivity
We performed region of interest and seed-based correlation analysis (Cordes et al., Reference Cordes, Haughton, Arfanakis, Wendt, Turski, Moritz, Quigley and Meyerand2000). The seed-based correlation map shows the level of functional connectivity between a seed and every voxel in the brain, whereas the region of interest analysis expresses the functional correlation between predefined regions of interest coming from a standard atlas. The analysis process is defined as follows: First, the subject images are aligned using Algorithm 3, then the element-wise arithmetic mean across subjects is calculated, and finally, the functional connectivity analysis is developed on this average matrix.
We take the frontal pole as seed, being a region with functional diversity (Liu et al., Reference Liu, Qin, Li, Fan, Wang, Jiang and Yu2013). The anatomical alignment considered here refers to the MNI space normalization (Jenkinson et al., Reference Jenkinson, Bannister, Brady and Smith2002). Figure 3 shows the correlation values between the seed and each voxel in the brain using data without functional alignment (top of Fig. 3) and with functional alignment using the Efficient ProMises model (bottom of Fig. 3). The first evidence is that the functional alignment produces more interpretable maps, where the various regions, such as the superior temporal gyrus, are delineated by marked spatial edges, while the non-aligned map produces more spread regions, hence being less interpretable. It is interesting to evaluate the regions more correlated with the frontal pole, for example, the superior temporal gyrus. This region is associated with the processing of auditory stimuli. The correlation of the superior temporal gyrus with the seed is clear in the bottom part of Fig. 3, where functionally aligned images are used.
Seed-based correlation map for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
, using data only aligned anatomically (top figure), and data also functionally aligned by the Efficient ProMises model (bottom figure). The black point refers to the seed used (i.e., frontal pole with MNI coordinates (0, 64, 18)). So, the brain map indicates the level of correlation between each voxel and the frontal pole.
In contrast, the region of interest correlations analysis shows the integration mechanisms between specialized brain areas. Figure 4 indicates the correlation matrices of time-series extracted from the 39 main regions of the atlas of Varoquaux et al. (Reference Varoquaux, Gramfort, Pedregosa, Michel and Thirion2011). Using functionally aligned data (right side of Fig. 4), we can see delineated blocks of synchronized regions that can be interpreted as large-scale functional networks. Instead, using data without functional alignment (left side of Fig. 4) the distinctions between blocks are clearly worse. Using functionally aligned data, the left and right visual systems, composed of the dorsolateral prefrontal cortex (DLPFC), frontal pole (Front pol), and parietal (Par), are clearly visible, whereas in the analysis using functionally nonaligned data, this distinction is hidden by noise.
Correlation matrix for
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{M}}$$\end{document}
, using data only aligned anatomically (left figure) and data also functionally aligned by the Efficient ProMises model (right figure). The cells of the matrix represent the correlation between the regions (represented by the row/column labels) of the Varoquaux et al. (Reference Varoquaux, Gramfort, Pedregosa, Michel and Thirion2011)’s atlas.
The preprocessed data are available on the GitHub repository: http://github.com/angeella/fMRIdata, as well as the code used to perform functional connectivity: http://github.com/angeella/ProMisesModel/Code/Auditory.
5. Discussion
The ProMises model provides a methodologically grounded approach to the Procrustes problem allowing functional alignment on high-dimensional data in a computationally efficient way. The issues of the perturbation model (Goodall, Reference Goodall1991)—non-uniqueness, critical interpretation, and inapplicability when
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$n \ll m$$\end{document}
—are completely surpassed thanks to our Bayesian extension. Indeed, the ProMises method returns unique and interpretable orthogonal transformations, and its efficient approach extends the applicability to high-dimensional data. The presented method is particularly useful in fMRI data analysis because it allows the functional alignment of images having roughly
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$$200 \times $$\end{document}
200,000 dimensions, obtaining a unique representation of the aligned images in the brain space and a unique interpretation of the related results.
In the application example presented in Sect. 4, a subsample was analyzed. However, the algorithm has a linear growth in N, and therefore, it is not a problem to work with larger samples. Also, the algorithm permits a parallel computation for the subjects.
The Bayesian framework gives the user the advantage and the duty to insert prior information into the model through k and
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
. The parameter k plays the role of regularization parameter, which is rarely known a priori. We estimated it by cross-validation, although it may be interesting to adapt approximations-based methods (e.g., generalized cross-validation) to reduce the computational burden. Alternatively, we could assume a prior distribution taking values in
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\mathbb {R}}^{+}$$\end{document}
for the regularization parameter k and proceed to jointly estimate this parameter as well. More interestingly, the matrix
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
addresses the estimate of the optimal rotations, which is favorable in the analysis of fMRI data because, in this context, the variables have a spatial anatomical location. In the example in Sect. 4, our definition of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
favors the combination of voxels with equal location. However, a more thoughtful specification can entirely exploit the voxels’ specific spatial position in the anatomical template. This opens up the possibility to explore various specifications of
\documentclass[12pt]{minimal}\usepackage{amsmath}\usepackage{wasysym}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{amsbsy}\usepackage{mathrsfs}\usepackage{upgreek}\setlength{\oddsidemargin}{-69pt}\begin{document}$${\varvec{F}}$$\end{document}
and will be the subject of further research.
Acknowledgements
The authors thank the Editor, the Associate Editor as well as the two anonymous referees for helpful comments that greatly improved the paper. Angela Andreella gratefully acknowledges funding from the grant BIRD2020/SCAR_ASEGNIBIRD2020_01 of the University of Padova, Italy, and PON 2014-2020/DM 1062 of the Ca’ Foscari University of Venice, Italy. Some of the computational analyses done in this manuscript were carried out using the University of Padova Strategic Research Infrastructure Grant 2017: “CAPRI: Calcolo ad Alte Prestazioni per la Ricerca e l’Innovazione”, http://capri.dei.unipd.it. The authors thank Prof. Umberto Castiello and Dr. Silvia Guerra for sharing the cinematic plant data.
Funding
Open access funding provided by Università degli Studi di Padova within the CRUI-CARE Agreement.
Open access
