1 of 56

MM Algorithms and Their Applications in Neuroimaging, Clustering and Graph Embedding

Daniel Tward, Kenneth Lange, Gary Zhou

2 of 56

Significance of Neuroimage Registration for Brain Mapping

Image registration transforms sets of image data into one coordinate system
Crucial for aligning brain images of different subjects, locations and time points
Useful for reconstructing 2D sections into 3D volumes and mapping datasets to a common coordinate reference brain atlas
Intensity-based (nonconvex) versus landmark based methods (quadratic)

Wang et al., 2020

Image registration is defined as the process of transforming a set of image data into one coordinate system. If we view the image as function I of spatial coordinates x, transforming an image I with some transformation phi means that the underlying coordinate gets transformed by the inverse of phi. This definition turns it into a group action.

Registration is a key processes in high-throughput neuroinformatics and brain imaging. It is crucial in reconstructing 2D brain histology sections into 3D volumes, aligning image series over time, and mapping datasets to a common coordinate reference brain atlas. Within this framework, registration and reconstruction elucidate how brain structures vary across population groups, time series, disease states, and imaging modalities.

Here on the left we show two adjacent coronal histology sections of an adult mouse brain sectioned at 10 micron thickness stained for Nissl and imaged with a standard bright field light microscope. Some of the structures you can see here are the CA fields and dentate Gyrus within the hippocampus, which are important for learning memory and neurodegeneration.

On the right is the widely used adult mouse brain Allen Reference Atlas, which was created as a sequence of 2D histology slices and reconstructed into a 3D volume. In practice we could either align these two sections with each other and find the optimal registration transformation, or we could find the transformation that best registers each of these slices onto the Allen reference atlas.

We divide registration models into landmark based (where a distance between corresponding feature points is minimized) and image intensity based (where a dissimilarity between two images is minimized). The first category of algorithms has the advantage of being quadratic and having closed form solutions but the disadvantage that landmarks must be identified manually or by sophisticated software. The second category possesses the opposite advantages and disadvantages. These algorithms are usually non-convex and iterative but require no manual specification of landmarks. When landmark points are correctly identified and annotated from human input based on domain knowledge, closed-form solutions for aligning landmark points result in less error and lower computational cost than intensity-based methods. However, in situations where manual annotation is insufficient or unavailable, such as high-throughput imaging, intensity-based methods are preferred. We will see in the next few slides how a minorization-maximization (MM) algorithm can realize the advantages of both approaches by minimizing a sequence of simpler problems with closed-form solutions.

3 of 56

Problem Formulation of Pairwise Rigid Registration

Minimize

R preserves norm, so it is equivalent to maximizing

Minimizers depend on the choice of interpolation method
Gaussian interpolation leads to

In the first part of my presentation I will focus on rigid pairwise registration, which is a simple but critical component as a subproblem in more complex registration schemes. For example, rigid registration can initialize non-rigid registration or facilitate joint rigid and nonrigid deformation, and pairwise registration is a frequent subproblem in sequential registration, which we will study in much more detail in our next paper. Since the serial mouse brain histology sections we study are cut and imaged using tape transfer technology, it is appropriate to consider rigid only registration here.

To formulate the problem, given two images I and J, we need to find the optimal rotation matrix R and translation vector c to minimize the objective function E, which is simply the sum of squared difference between the image J and the transformed image I. Since R preserves norm, this is equivalent to maximizing the cross integral term.

Notice that the image is only given as set of pixel intensity values on the grid points, but its necessary to evaluate the image at arbitrary non-grid points. Which means we need to choose an image interpolation method, and exact registration result (optimal transformation parameters R and c) depend on our choice of interpolation method. Here we choose Gaussian interpolation with the Gaussian kernel k of standard deviation sigma, which we will explain in detail next slide. Gaussian interpolation converts the cross term to the expression here.

4 of 56

Registration of Images under Gaussian Interpolation

5 of 56

MM Algorithm for Pairwise Rigid Registration

Gradient based method common in this setting, but its descent is not guaranteed, and speed depends on the careful selection of step size.
MM is a method for constructing optimization algorithms, by exploiting problem structure and replacing a hard problem by easier surrogate problems.
MM stands for either Majorization-Minimization or Minorization-Maximization
Some advantages include:
Avoids need for line-search and for adjusting descent step sizes
Separates parameters (for parallelization)
Linearizes problem (for analytical solutions)
Enforces constraints without parametrization

6 of 56

MM Algorithm

To minimize f(x), we can majorize it with a surrogate function, g(x| . ), anchored at its current iterate
g satisfies both tangency and domination conditions: and for every x
Minimizing the surrogate function to get the next iterate makes sure that f(x) is going downhill:

7 of 56

An Example for Minimizing f(x)=cos x

8 of 56

How to Find Surrogate Functions

They need to be easier to optimize
We usually use inequalities and convexity to find them:
Jensen’s inequality
Cauchy Schwarz inequality
von Neumann trace inequality
This step requires mathematical ingenuity (trial and error)

9 of 56

Optimizing Gaussian-Interpolated Cross Term Integral

(We use the fact that integral of product of Gaussians is a Gaussian again)

Pixel values

distance squared

10 of 56

Simplified Objective and its MM Surrogate

We derived the surrogate function as a weighted sum of square distance

11 of 56

Now these Functions are Easy to Optimize

We can now use the point-set rigid optimization method, which amounts to calculating the SVD of a 2x2 or 3x3 matrix.
This method is efficient because we turned an image problem into a point-set problem, which can be solved via linear algebra

12 of 56

Reference: Least-Squares Rigid Motion Using SVD by Olga Sorkine-Hornung and Michael Rabinovich

13 of 56

Issue of Truncation

A lot of terms in the sum are small due to negative exponentials and do not contribute significantly.
The number of summands is (# of I pixels) * (# of J pixels), which can be very large and computationally expensive (quadratic complexity).
We truncate these terms and keep only coefficients that are large enough.

14 of 56

Details of Truncation

15 of 56

Error Bound for Truncated SVD

We bounded the error of the optimal rotation and translation of our truncation
Error is determined by the smallest singular value of the SVD matrix
Given two matrices S and S’ (the covariance matrix) such that their difference |S-S’| is small, we showed that the error |UV-U’V’| is bounded by twice the reciprocal of the smallest singular value times their difference.

Experiments show that the smallest singular value is usually large

16 of 56

Smallest Singular Value of SVD for Registering A Pair of Real Images with MM algorithm is large

17 of 56

Performance Comparison with GD methods

MM with Nesterov acceleration compared with Nesterov-accelerated Block Gradient Descent with fine-tuned step-size parameters.
MM holds the advantage of being faster, more robust and parameter-free.
Advantages have been demonstrated for both simulated and real images under noisy conditions.

18 of 56

Real Images of Mouse Brain Microscopy Slices

19 of 56

Performance Comparison on Simulated Images (Blue-MM, Beige-GD)

20 of 56

Statistical Analysis of Registering 102 Pairs of Real Images

90% of the times the MM algorithm outperforms block descent in speed at 200th iteration.
Block descent method occasionally (8% of the times) converges to a local maximum lower than the maximum achieved by MM, while the reverse never happens.

Difference in Cost= MM Cost-GD Cost (blue-MM, Red-GD)

21 of 56

Possible Extensions

We added a scaling factor (replace R with sR)
Making the algorithm more robust to imaging artifacts.
Robust sequential registration (focus of our next paper)
Consider a general affine transformation (replace R with some invertible matrix A)

22 of 56

Robust Alignment of Sequential Images

Important for 3D reconstruction and common coordinate mapping
Due to many imaging artifacts, algorithm has to be robust
Sequential pairwise registration leads to accumulation of errors

One important application of sequential neuroimage registration is in 3D brain histology reconstruction. Here the goal to is visualize the spatial relationships of neurons and brain regions. High-quality 3D histology reconstruction enables: (a) examination of neural structures with respect to their environment in 3D with or without the help of 3D imaging, (b) correlation of post mortem microscopy data with non-invasive brain imaging data to improve signal interpretation, and (c) the creation of digital brain reconstructed volumes called atlases that assist in 2D and 3D visualization and provide quantitative measures of brain anatomy.

Nissl stained, 10 micron thickness coronal histology sections of an adult mouse brain under a standard bright field light microscope. We see the x-z and y-z cross section views the original and aligned slices. Note the hippocampus formation (elongated arcs) on the Y-Z and X-Z views wouldn’t be visible unless we sequentially register the sections correctly.

23 of 56

Starting Point of our Registration Problem

Goal: rigidly align a sequence of brain microscopy slices to minimize the sum squared differences with respect to L2 norm.
This non-linear intensity-based problem requires an iterative optimization approach. Rigid usually an initialization step before non-rigid alignment.
Efficiency and robustness become important as datasets scale up.

24 of 56

Rigid Registration of a Sequence of Images

Reminder we need to find a set of rigid transformations R_i such that we minimize the objective function

We proposed (avg. of aligned neighbors)

then derived MM surrogate function below via Jensen’s inequality

25 of 56

Another Model of Sequential Registration

Equivalent when registration is rigid
Non-equivalent if transformations are deformable
We need both parameters and a reconstructed volume (atlas)
Decouples all pairs of transformations and easy to add robustness

26 of 56

Robust Objective Function

We control a special robust loss function with the parameter c

27 of 56

Transformations as Diffeomorphisms

Transformations and their inverses smooth/differentiable
Preserves topology and physical realism
We remove noise in the reconstructed volume (atlas) in 3D (possibly anisotropically) and control it with the parameter a

28 of 56

LDDMM as a Deformable Pairwise Registration Method:

Large Deformation Diffeomorphic Metric Matching

LDDMM as an ODE-constrained optimization
Guarantees diffeomorphism
Hilbert gradient descent or geodesic shooting
Computationally expensive
Stationary Velocity Field (SVF)
Scaling and Squaring method solves SVF

Now that we have robustness and smoothness of the atlas, how do we ensure diffeomorphic transformations? The answer is the so called LDDMM. LDDMM can be viewed as an optimal control problem where the control variable is infinite dimensional time-dependent vector field with a special norm V. Specifically we need to find the minimizer vector field V, and solve the ODE problem of V at every possible initial point to get the flow of the vector field at time t=1 phi1 so that the minimal vector field minimizes the sum of the special norm of v and the L2 difference between the image I1 and the image I0 transformed by the flow of that vector field at t=1. The resulting flow has the nice property of being diffeomorphic, and this problem is usually solved with either Hilbert gradient descent or geodesic shooting. As you can see, this can be computationally expensive, so in our problem we replace the time varying vector field with the time independent vector field, so that the flow at time 1 of this vector field can now be efficiently estimated with the so called scaling and squaring method, which we will explain in the next slide. Now we get the final form of the objective function F(I,phi)

29 of 56

Scaling and Squaring to Solve SVF Registration

We use repeated interpolations and function compositions to estimate the flow of a time-independent vector field at t=1

30 of 56

Optimizing the Objective Function

First use MM to reduce robust function into weighted least squares

31 of 56

Another MM for Updating Reconstructed Volumes (Atlases)

Inverting (LL+K) is difficult, but inverting (LL+id) with FFT is easy
Design a surrogate function that can be easily solved with FFT

32 of 56

Applications

original Nissl stained 10 micron adult mouse brain histology sections unaligned and noisy with imaging artifacts (discoloration and missing tissue due to slice tearing), we then have the aligned images, where the hippocampal formations become visible in the x-z and y-z views. dentate gyrus (dark stained semicircular curves) We also have the robustly smoothed atlas that preserves key structures of the mouse brain and significantly decreases outliers and imaging defects. Note that the x-y plane of the atlas appears less smooth than the other two planes because the anisotropic Laplacian coefficient along this axis is smaller than the Laplacian coefficient along the other two axes. At last we have the pixelwise robustness weight ws that measure possible pixel outliers. The fourth row shows the robustness weights computed for each plane view, where a higher weight represents better reliability of a pixel's intensity and a lower chance of it being an outlier.

33 of 56

Koay et al. 2016

The third dataset consists of video frames collected with a two-photon microscope on a mouse expressing a fluorescent protein in the parietal cortex. In this calcium two-photon video, the sudden appearance of bright spots represents neurons firing. Therefore, an area with low weights, which depend on the differences between the original image and the reconstructed volume image, should represent neuron firing spikes rather than imaging defects or noise. The first row shows the original unaligned frames with random motion between each frame due to \textit{in vivo} imaging with several bright neurons and a particularly bright (firing) neuron at the lower left corner of the x-y plane. The second row shows the same frames after deformation and rigid alignment so that the pixels of the same firing neuron are kept in the same location represented by the vertical bright bands. The third row represents a smoothed volume reconstruction, where the t-axis appears more smoothed than the x-y axes due to the aforementioned anisotropy factor. Note that the bright spot (firing neuron) in the lower left corner is missing in the reconstructed volume because it is a transient signal with a low robustness weight, while the more durable signals of several less bright and longer-lasting neurons are preserved. The last row shows images of robustness weights. For example, we see a dark spot in the lower left corner that represents a rapid neuron firing.

34 of 56

Based on the phantom dataset with known ground truth, we now study the effects of the Laplacian scaling constant $a$ of equation (\ref{eq:Smooth operator}) and the robustness constant $c$ of equation (\ref{objective}) on algorithm speed and registration accuracy. The first two rows of Figure \ref{quantitative} show that faster convergence occurs as $c$ increases or $a$ decreases. Overall, the different combinations $(a,c)$ lead to little discernible differences in the accuracy of image denoising and registration. A slight increase in error of the converged translation vector can be observed in the top right panel of Figure \ref{quantitative}. This error is compensated by a larger and slowly changing reduction in error of the rotation angle in the middle right panel. Row 3 of Figure \ref{quantitative} displays the noise power spectrum of the realigned calcium two-photon imaging frames under different pairs $(a,c)$. We define the noise power spectrum as the average of the x-axis and the y-axis norm squared Fourier transform of the noise of the frames. The realigned frames exhibit higher power at high frequencies than the original frames, indicating improved resolution. Conversely, averaging unaligned sections leads to blur, decreasing power at high frequencies for the original frames.

35 of 56

he first column depicts a failure under an extremely low value of $c$. As seen in Figure \ref{rob 02}, an increase in $c$ leads to a rapid decrease in the weights as a function of the residual error $s$. In other words, even a small error between the reconstructed volume and the transformed image can result in a near-zero weight under a small $c$. Conversely, a large $c$ implies that the weights are all close to $1$ regardless of the size of the residual error $s$. Thus, matching terms become essentially unweighted. The changing reconstructed volume images from left to right are formed by gradually equalizing weights in matching. This relationship is confirmed in the row of the robustness weight images, where the leftmost image depicts an extreme case where all weights are nearly $0$. In short, lower values $c$ correspond to higher sensitivity to errors. Thus, the outliers detected by low robustness weights become less pronounced as $c$ increases. The failure when $c$ is nearly $0$ also leads to an imperfect alignment of the original images because all the matching terms in the surrogate function have a coefficient near $0$.

36 of 56

37 of 56

38 of 56

Testing Robustness with Imaging Artifacts

39 of 56

Multimodality Robustness with MIND Preprocessing

40 of 56

Analysis of Spatial Transcriptomics and Brain Regional Network:

Sylvester Equation Connection

Spatial transcriptomics combines both spatial and RNA gene expression data.
One central task in ST is to learn cell types
Reasonable assumptions lead to regularized clustering
Brain regional networks for neuroscience and clinical applications
fMRI connectivity data is noisy, high dimensional and variable
Link to Sylvester equation via graph theory

Yao et al. (2023)

We wanted to contribute to two important problems in computational biology today, analyzing spatial transcriptomics and functional MRI data.Spatial transcriptomics is a cutting-edge molecular profiling technique that maps gene expression within the structural context of intact tissues. Unlike standard bulk or single-cell RNA sequencing, which destroy a tissue's physical structure, this method reveals exactly where genes are active, allowing researchers to study cell-to-cell interactions and complex tissue architecture. One central task in ST is to learn cell types by organically combining both types of information, which we will next model in our paper as a constrained regularized graph-based clustering problem if we make some reasonable assumptions. Another problem we want to explore is the analysis/detection of brain regional networks from fMRI brain scan data, which measures the blood oxygen concentration in real time.
Why study brain regional networks? (importance for neuroscience & clinical applications) Many neurological conditions are understood as disruptions of functional neural networks. Default Mode Network (DMN), Dorsal Attention Network (DAN), Salience Network, Executive Control Network
How do we study it? (fMRI functional connectivity, or DTI structural connectivity)
FC measures similarity between brain signals from two regions with Pearson’s correlation.
Challenges: high-dimensional data, multiple subjects, variability across individuals and conditions.
A potential solution to this challenge is to use a Graph-based approach to represent brain regions as nodes and connectivity as edges and jointly embed multiple graphs into a continuous lower dimensional vector space.
In both cases to address these problems in computational problems we wanted to develop unsupervised learning models to provide low dimensional embeddings which can be used for applications like clustering, classification and visualization.. Interestingly we were able to show that despite their very different fields solving these models reduces via a graph theory perspective to solving a specific type of matrix equation called the sylvester equation or its modified and transposed versions.
Left figure depicts the fMRI scan of the Default mode network in the brain the dorsal medial prefrontal cortex, the posterior cingulate cortex, the precuneus, and the angular gyrus. The figure on the right shows a slice of the Allen mouse brain spatial transcriptomics project, where each dot represents a cell (most of them are neurons ) and color of the dot represents a cell type.

41 of 56

Sylvester and Transposed Sylvester Equation

AX+XB=C and AX+XᵀB=C
Used in control theory and factor analysis
Vectorization inefficient
Established algorithms: Bartels-Stewart and conjugate gradient
Focus: A large, possibly sparse and B small and dense for Sylvester (A tall B wide for the transposed Sylvester)
Example: previous spatial transcriptomics problem
Brain regional functional network detection is another
We propose MM algorithms for these equations
Then we focus on these two applications

Here are the standard Sylvester equation and the transposed Sylvester equation, where A,B,C are matrices of compatible dimensions (A, B square matrices in Sylvester equations, while dimensions of A are equal to the dimensions of B transposed in the transposed Sylvester equation). They are commonly encountered in control theory factor analysis and many other fields. Because of linearity, they can be vectorized to form a linear equation whose dimension is the product of the dimension of A and B, so solving with vectorization is inefficient and computational expensive. Common established methods for solving it include B-S and CG. BS uses Schur decomposition of A and B to decompose the equation into many triangular matrix equations, however, it is a direct solver with computational complexity of the dim of the larger matrix cubed, also note the complexity doesn’t depend on the sparsity of either matrix, so does become more efficient with sparse matrices. CG is an iterative method, which is an modification of the original conjugate gradient that requires positive semidefinite A and its convergence sparsity-dependent, with per iteration complexity of big O of the density of the larger matrix. Here we focus on the cases where A is large, possibly sparse and B is small and dense for the Sylvester equation, and A is tall, B wide for the transposed equation. This assumption is not necessary for the MM algorithms we will derive, but when the assumption holds our algorithms achieve high computational efficiency by exploiting the dimensionality asymmetry. The previous spatial transcriptomics problem is one such example, where the dim of A iis the number of cells, dimension of B is the number of cell types. For the brain regional functional network detection problem is another such example, where the joint graph embedding model gives a modified transposed Sylvester equation where the row dimension of A is number of regions or nodes in the brain, column dim of A is the dimension of the latent vector space, so A is tall and B is wide.

42 of 56

Problem Formulation for Spatial Transcriptomics

Row of Y: features attached to each node (cell)
Row of X: cluster centers
Row of A: admixture coefficients for each cell.
Penalty: smoothness on the admix coefficients of neighboring cells
How is this related to the Sylvester equation? Graph Laplacian L
# of cell types << # of cells, so dim L >> dim XXᵀ
Efficiently solve this type of Sylvester equations

Here we explain the regularized cell type clustering formulation of the spatial transcriptomics problem in detail. F is the objective function which we optimize with respect to the two matrices A and X. Each row of Y represents the RNA gene expression profile of each cell, each row of X represents the gene expression centroids of each cell type, each row of A represents the proportions or probabilities of that cell to belong to each of the cell types. The penalty term simply means that we discourage cells in the same neighborhood to be too dissimilar or have very different cell type compositions or probability distributions. Here the neighborhood Ni of the cell i can be defined according to the context of the problem, but here we define it to consist of all cells j that are within a certain distance from cell i. How is this related to the Sylvester equation then? If we write the neighbor structure as a graph, where the nodes are the cells and edges are drawn according to the neighborhood defined, it turns out that we can write the second regularization/penalty term in terms of the graph Laplacian of that graph, which is defined as the degree matrix minus the adjacency matrix of the graph, so the second terms becomes the trace of the product ALA. Estimating X with A fixed is simply ordinary least square, whereas estimating A with X fixed requires us to take the gradient of f with respect to the matrix A, which gives us this expression. So we see that this becomes a Sylvester equation in A, and two given matrices to the left and right of A are the graph Laplacian L and XX transposed. Note the number of cell is a lot larger than the number of cell types, so dim L>>dim XX. So we focus the efficiently solving this type of Sylvester equations with dimension asymmetry and possible sparsity.

43 of 56

MM Algorithms for General Sylvester Equation

Matrix A large and possibly sparse, B small and dense (m>>n)
Eigendecomposition feasible for B, but not for A
Power iteration feasible to obtain largest eigenvalue of A

44 of 56

MM for Regularized Clustering

45 of 56

MM for Transposed Sylvester Equation

Bartels-Stewart method does not work here
Assume A tall, B wide

46 of 56

Performance Comparisons with Numerical Experiments

47 of 56

Regularized Clustering in Spatial Transcriptomics

Mouse brain data set with 5210 cells and 451 genes into 14 cell types
Ill-defined unless A is constrained
Iterative projection onto the unit simplex
Same # of cell types as Allen’s

Yao et al. (2023)

48 of 56

Clustering Results

Clustering similarity measured with adjusted Rand index
Easy to select the top r genes driving the clustering
Number r can be selected based on BIC curve

Yao et al. (2023)

49 of 56

Joint Graph Embedding

Joint Graph Embedding generalizes Adjacency Spectral Embedding
Node i has a latent representation xi of lower dim n
P(i,j)= xiᐧxj
ASE is a low-rank matrix decomposition
JGE is ASE for m graphs with s common vertices
Each graph has its own weights (loading vector)
Equivalent to a certain tensor decomposition

50 of 56

51 of 56

JGE for Brain Regional Network Detection from fMRI data

Correlations on all pairs of 112 brain regions from ABIDEII (Autism Brain Imaging Data Exchange)
Thresholded to adjacency matrices for 21 subjects
High dimensional, noisy and variable
Embedding: dimension reduction, interpretability and visualization
JGE proposed by Wang et al. (2021) finds a shared embedding space for multiple graphs across subjects or time points

52 of 56

Tward and He, (2023)

53 of 56

MM Performance for Simulated Graphs JGE

54 of 56

MM Performance for brain fMRI JGE

(GD used in Wang et al. (2021))

55 of 56

Future Directions

Directed and weighted graphs
Smoothing of loading vectors in JGE
Larger data sets
Generalized Sylvester

56 of 56

Thank You for Your Time!

References

Zhou G, Tward D, Lange K. A Majorization-Minimization Algorithm for Neuroimage Registration. SIAM J Imaging Sci. 2024;17(1):273-300. doi: 10.1137/22m1516907. PMID: 38550750; PMCID: PMC10977051.

Zhou G, Tward D, Lange K. A Robust MM Algorithm for Sequential Neuroimage Registration. Submitted to SIAM J Imaging Sci

Zhou G, Tward D, Lange K. MM Algorithms for Sylvester Equations with Applications in Clustering and Graph Embedding. Preprint