SHIFT- AND STRETCH-INVARIANT NON-NEGATIVE MATRIX FACTORIZATION WITH AN APPLICATION TO BRAIN TISSUE DELINEATION IN EMISSION TOMOGRAPHY DATA Anders S. Olsen1∗ , Miriam L. Navarro1,2 , Claus Svarer1 , Jesper L. Hinrich3 , Morten Mørup3 , Gitte M. Knudsen1,4 1
Neurobiology Research Unit, Copenhagen University Hospital Rigshospitalet, Copenhagen, Denmark Department of Neuroscience, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark 3 Department of Applied Mathematics and Computer Science, Technical University of Denmark, Kgs. Lyngby, Denmark 4 Department of Clinical Medicine, Faculty of Health and Medical Sciences, University of Copenhagen, Copenhagen, Denmark ∗ Corresponding author: [email protected]
arXiv:2604.08161v1 [cs.LG] 9 Apr 2026
2
ABSTRACT
[4], [5], [6], [7], see [8] for a recent review. While effective in certain contexts, these approaches generally assume that components Dynamic neuroimaging data, such as emission tomography are fixed in time and therefore fail to account for delays or dispermeasurements of radiotracer transport in blood or cerebrospinal sion of TACs across regions, evident as relative shifts and stretching fluid, often exhibit diffusion-like properties. These introduce of TACs, which can result in misclassification or blurred regional distance-dependent temporal delays, scale-differences, and stretchboundaries. Fig. 1 shows example data, where a radiotracer is ining effects that limit the effectiveness of conventional linear modfused in the cisterna magna (CM) and subsequently disperses in the eling and decomposition methods. To address this, we present CSF and into the gray matter. This causes TACs to differ in onset the shift- and stretch-invariant non-negative matrix factorization and tracer movement speed, leading to delays and stretching even framework. Our approach estimates both integer and non-integer between areas of the same tissue type. temporal shifts as well as temporal stretching, all implemented in the To overcome this limitation, we propose a shift- and stretchfrequency domain, where shifts correspond to phase modifications, invariant extension of NMF. Building on shift-invariant NMF and where stretching is handled via zero-padding or truncation. [9], which estimates integer (i.e., whole samples) and non-integer The model is implemented in PyTorch (https://github.com/anders(sub-sample) temporal delays, our model additionally incorporates s-olsen/shiftstretchNMF). We demonstrate on synthetic data and stretch-invariance to capture differences in the speed of tracer kibrain emission tomography data that the model is able to account for netics. This flexibility enables more accurate alignment of TACs stretching to provide more detailed characterization of brain tissue across tissues exhibiting heterogeneous dynamics. We validate the structure. method on synthetic data and animal brain SPECT experiments, and Index Terms— Non-negative matrix factorization, shift-invariance, provide an open-source PyTorch toolbox (https://github.com/anderss-olsen/shiftstretchNMF). stretch-invariance, SPECT data. 1. INTRODUCTION
2. METHODS
Dynamic medical imaging, e.g., positron emission tomography (PET) or single-photon emission computed tomography (SPECT) enables tracking of radiotracer kinetics and physiological processes across space and time. These measurements can provide valuable insights into brain function including movement of cerebrospinal fluid (CSF) or binding patterns of neuroreceptors. However, accurate analysis of tracer-based imaging data is hindered by several challenges. Extracting regional signals typically requires aligning emission data with high-resolution anatomical images (MRI or CT). This registration step is error-prone due to subject motion and physiological variability between scans. Furthermore, partial volume effects blur tissue boundaries, degrading the precision of regional estimates [1]. An alternative strategy is to model the latent composition of the emission data, bypassing anatomical information. By clustering or decomposing raw time–activity curves (TACs), one can identify distinct regions with unique kinetic signatures and thus mitigate alignment and resolution issues. Various clustering and decomposition techniques have been applied to this task, including K-means, fuzzy C-means, Gaussian mixture models, and non-negative matrix factorization (NMF) [2], [3],
2.1. Shift- and stretch-invariant non-negative matrix factorization Given a data set X ∈ RP ×N , where j = 1, . . . , P indexes channels (e.g., voxels in imaging data) and t = 1, . . . , N indexes temporal samples, the standard bilinear factor analysis model is xj (t) ≈
X
aj,k sk (t),
(1)
k
yielding a set of K channel maps A ∈ RP ×K and temporal profiles S ∈ RK×N . Additional constraints may be imposed on the decomposition, such as statistical independence in independent component analysis [10] or non-negativity in NMF [11], the latter being particularly suitable when the data itself is non-negative. The shift-invariant factor analysis extends this framework by allowing temporal profile shifts (delays) for each channel [12]: xj (t) ≈
X k
aj,k sk (t − τj,k ) ,
(2)
time-scaling factor. Consequently, working in the frequency domain makes model estimation more computationally efficient: x̃j (f ) ≈
X
f
aj,k s̃k (rj,k f ) e−i2π N τj,k ,
(5)
k
where s̃√ k and x̃j are one-sided discrete Fourier transforms (DFT), i = −1, and f = 0, . . . , NF F T . Of note, in the continuous Fourier transform, the time-frequency scaling property requires amplitude modulation to satisfy Parseval’s theorem. In the DFT, no closed-form analogue exists for non-integer r [15]. Instead, we rescale the coefficients of truncated spectra/signals empirically (see Sec. 2.3). Phase shifts are frequency-independent, hence frequency scaling only applies to s̃k (f ). In condensed form, the loss function to be optimized is P
L=
1 X w⊙ N j=1
x̃j −
X (τj,k ) (rj,k ) ⊙ s̃k ãj,k
!
2
,
(6)
k
f (τj,k ) where ⊙ is element-wise multiplication, ãj,k = aj,k e−i2π N τj,k (rj,k ) contains frequency-scaled comis a vector of size NF F T and s̃k ponent profiles. The vector w = [ √12 , 1, . . . , 1, √12 ] ensures correct weighting of DC and Nyquist terms for comparison between timedomain and one-sided frequency domain loss functions according to Parseval’s identity, which establishes a direct relation between the time and frequency domain signal energy.
2.2. Estimation of shifts
Fig. 1. SPECT data from an example pig after 99m Tc-DTPA radiotracer injection in the cisterna magna (CM). Time-activity curves in the time-series plot - from selected voxels and timepoints indicated by squares - were normalized to unit norm for visualization.
Non-linear optimization of the shift-parameter τj,k in previous studies using the Newton-Raphson method found that estimated shifts seldom exceeded one temporal index due to a large distance between minima in the parameter landscape [9], [16]. Instead, shifts were proposed to be estimated first using the cross-correlation function for integer shifts and subsequently non-linear finetuning for obtaining non-integer shifts. For each component k, we compute the residual spectrum with all other sources projected out: g̃j,k (f ) = x̃j (f ) −
with τj,k ∈ R denoting the per-channel, per-component shift (temporal delay in TACs). Shift-invariant NMF estimation has been proposed using multiplicative updates [9] and a time-frequency gradient method [13]. A different generalization of NMF is to account for stretch (dispersion) by expanding or compressing the temporal profiles (time-scaling): X xj (t) ≈ aj,k sk (t/rj,k ) , (3) k
with rj,k ∈ R+ denoting stretch (r > 1) or compression (r < 1). This idea was recently explored using spline interpolations in the temporal domain [14]. Here, we propose combining both shift- and stretch-invariance in the same model such that for each channel j and component k, both a shift τj,k and a stretch rj,k is learned: X t − τj,k xj (t) ≈ aj,k sk . (4) rj,k k
In the frequency domain, temporal shifts correspond to constant phase modulations and stretching to frequency scaling by the inverse
X (τj,k′ ) ãj,k′ (f ) s̃k′ (f ).
(7)
k′ ̸=k
The cross-spectrum is h̃j,k (f ) = g̃j,k (f )∗ s̃k (f ), where (·)∗ is the complex conjugate. The inverse DFT gives the cross-correlation function where the optimal delay is lj,k = arg max hj,k (t), t
τj,k = lj,k − N,
(8)
assuming zero-based indexing. We note that computing the crosscorrelation function can be avoided by estimating the slope of the phase of the cross-spectrum ∠h̃j,k (f ). However, this requires “unwrapping” the phase to avoid discontinuities at ±π and combined with the slope estimation, this turns out to actually be slower than the inverse DFT of the cross-spectrum. 2.3. Estimation of stretches Here we extend the above cross-correlation method to also simultaneously estimate stretches. We start by constructing a library of (b) spectral profiles s̃k (f ), b = [−NF F T /2+1, . . . , 0, . . . , NF F T /2− 1] by truncating the spectrum (b < 0) or zero-padding (b > 0). To obtain the same NF F T for all profiles we computed the inverse
DFT and zero-padded in the time-domain (b < 0) or truncated the temporal representation (b > 0) followed by the forward DFT. In the case of truncating in either domain, the retained coefficients were rescaled by the square root of the ratio of the original energy to the truncated energy to ensure the total energy was conserved. During optimization, the library is constructed for each update of S. Of note, the integer index variable b is related to the continuous frequency axis scaling parameter through r = 1 + NFbF T . Residuals are then computed as X (τj,k′ ) (b ′ ) ãj,k′ (f ) s̃k′ j,k (f ),
g̃j,k (f ) = x̃j (f ) −
(9)
k′ ̸=k
where bj,k indexes into the library of spectral profiles. The (b) (b) cross-spectrum is computed over h̃j,k (f ) = g̃j,k (f )∗ s̃k (f ) such that for each channel j and component k, it is now a 2D maximization problem over dimensions b and t. (b)
[lj,k , bj,k ] = arg max hj,k (t), t,b
τj,k = lj,k − N.
(10)
Fig. 2. Two-component NMF models applied to two-component synthetic data with random shifts and stretches. Model performance is measured using matched correlation of the A-matrix to the ground truth design matrix (variance is over 25 repeated model estimates), while the learned profiles (bottom) correspond to the S-matrix of the estimate with the lowest loss.
2.4. Estimation of A and S Given the estimated shifts and stretches, an analytical estimate of the matrix A is given in closed form (extended from [9] to also include stretching): (b
aj,k =
)
hj,kj,k (lj,k )
, (11) s̃H k s̃k i.e., the value of the maximal cross correlation scaled by the component energy, where (·)H is the Hermitian transpose operator. Any negative elements in A are clipped to zero. We propose non-linear optimization of S in PyTorch, which has repeatedly shown to produce comparable results to models estimated via analytically derived update rules for various non-convex unsupervised problems with appropriate initialization while being considerably faster [17], [18]. Using this approach, we optimize S̄ (S-bar) unconstrained while the softplus function (which is differentiable everywhere) is used to impose element-wise non-negativity prior to calculating shifts, stretches, and loss, i.e., sk,t = ln 1 + es̄k,t . Non-negativity in the time-domain does not have a direct equivalent in the frequency domain, so S̄ is in the time-domain and the spectral library of stretches is established for each iteration. We used ADAM (learning rate 0.1) in all experiments. The stopping criterion searched the latest 50 iterations, selected the lowest two losses, and stopped optimization if the relative change between the two losses was below 10−10 or if the loss only increased. 2.5. Data We used experimental data from five pigs infused with 99m Tc-DTPA in the CM immediately followed by N ≈ 100 SPECT-acquisitions, each accumulating radiation counts over six minutes (see an example subject in Fig. 1). These data were acquired to investigate radiotracer flow from CSF to the brain parenchyma to describe the glymphatic system in gyrated brains. The poor spatial resolution of SPECT data renders it particularly susceptible to partial volume effects [1] and radiotracer movement is slow [19], making a shiftand stretch-invariant model particularly viable to extract whole-brain tissue compartments. Due to vast scale differences between tissue compartments, the data in each voxel was normalized to unit
norm. To avoid inflating noise channels, only voxels with summed counts higher than 5·106 Bq/cc were included. Moreover, the spinal cord was excluded, totaling P ≈ 20.000 brain voxels. Representing shifts in the frequency domain assumes circular shifting, and to avoid late data entering as early data when shifted, we zero-padded the TACs by 20% of N in the time-domain. 2.6. Experiments Combined, we tested four models: (1) NMF, (2) Integer-shift NMF, (3) Non-integer-shift NMF (initialized from (2)), and (4) Shiftstretch NMF (initialized from (2)). For all models, we used K-shape (K-means with the cross-correlation distance [20]) for initialization, which we found to perform marginally better than other typical NMF-initializers (not shown). This procedure only initialized S while an initial guess for A was entered using a least-squares solution. Models (1) and (3) also optimized A non-linearly in PyTorch, since these models had no access to the cross-correlation function. 3. RESULTS 3.1. Synthetic data We generated synthetic data by applying random integer shifts and stretches between [−N/4, N/4] to two true components: One was half of a cosine period (soft peak), and the other was a Laplace probability distribution (sharp hump) (Fig. 2). Both ground truth components were pre- and appended by zeros to avoid the applied shifts and stretches to cause data to extend outside the data domain. 100 synthetic channels were generated for each component and model performance measured using matched correlation between the groundtruth and the learned A. The Shift-stretch NMF outperformed the other models in identifying the true labels (Fig. 2). The matched correlation coefficients were not precisely equal to one, which we attribute to the frequency domain loss function causing a loss of precision in the learned temporal profiles. The learned profiles were closer to the ground truth for the Shift-stretch model, while the shift models were not able to reconstruct particularly the half-cosine component, which had two
Fig. 3. Effect of model order K on variance explained (shaded area is variance over pigs).
humps instead of one. In comparison, the conventional NMF did not yield satisfactory component profiles. 3.2. Pig brain SPECT data Using brain SPECT data from five pigs, we observed that the Shiftstretch NMF model (applied to each pig separately) outperformed the other models across model orders, particularly for K = 1 and K = 2, indicating that the added stretching parameters are mostly beneficial for low numbers of components (Fig. 3). With more components, Shift-models and even the conventional NMF are able to fit the data comparably to Shift-stretch NMF. Example model fits for K = 3 components for one pig (pig5, the same as in Fig. 1) are shown in Fig. 4. For all models, one component corresponds roughly to CM - the infusion site, the second one to a broad CSF component, and the third one to gray matter. While the segmentation maps appear similar, the boundaries between clusters appear slightly sharper in NMF and Non-integershift NMF, i.e., the models where both A and S were optimized non-linearly. The temporal profiles were more different between models: The shift and stretch models had a significantly narrower CM component profile, indicating that the conventional NMF averages over many delayed channels. The shift-stretch NMF profiles are more noisy while still contributing to higher explained variance (Fig. 3), which we attribute to the manipulation of high-frequency components via zero-padding of spectra (stretching in time domain), such that high-frequency components are non-zero for fewer channels. To avoid this phenomenon, profiles could be reconstructed from only the first NF F T /2 frequencies. Fig. 4(right) shows the reconstructed example time-activity curves from Fig. 1 computed as the right-hand-side of Eqs. (1)-(4). The reconstructions are clearly better for the shift and shift-stretch models. 4. DISCUSSION We proposed a novel shift- and stretch-invariant version of NMF and applied it to synthetic data as well as SPECT data displaying slow dispersion over multiple tissue compartments (CSF, gray matter). Previous methods for such data did not account for shift and stretch-differences across the brain and would therefore erroneously assume distance-independent data. The method is applicable to data with similar shapes across channels and may therefore also be applied to, e.g., time-locked EEG or physiological measurements. The proposed stretching mechanism via truncating/zero-padding in the frequency domain is novel for matrix decompositions, to the best of our knowledge. A previous stretched NMF implementation [14] stretched via spline interpolations in the time-domain. We argue
Fig. 4. Extracted channel maps (left), temporal profiles (mid), and selected channel reconstructions (right) for pig-5 using K = 3 components for the four tested models. CSF: cerebrospinal fluid.
that our approach is simpler and has an easier optimization framework. The combination between shifting and stretching results in a versatile matrix decomposition framework, and the implementation in PyTorch means that the framework is very easily altered to related linear decomposition methods such as principal/independent component analysis or sparse coding by replacing non-negativity with orthonormality/independence/sparsity constraints on S. Some complications with the model should be mentioned. Shifting is circular, but stretching is not, and the circularity property of phase shifting in the frequency domain is likely seldom desired in real-world scenarios. We propose future studies to approach shifting by constructing a library of candidate temporally shifted profiles, similarly to what we did for stretching. Furthermore, the data to which the model is applied should be smooth to avoid large transients, which cause Gibbs ringing when manipulating the frequency components of component profiles - the presented SPECT data did not conform to this. In Fig. 4, component number three had a small hump in the beginning for most models, which we attribute to ghosting, i.e., the lack of complete uniqueness in the model. Future studies could enforce unimodality by penalizing sign-switching of gradients in S, which could also remedy the observed noise in component profiles. While the components were estimated by shifting and stretching each channel, the estimated temporal profiles for different components cannot be assumed to be aligned to each other (see Fig. 2), complicating, e.g., kinetic modeling. In future research, this problem may be approached by selecting a reference region/channel to which other channels are shifted/stretched before computing a weighted average via the estimated A-matrix. Future modeling could also focus on non-linear optimization of A in stretch-invariant models to further refine channel maps.
Compliance with ethical standards All pig experiments were performed in accordance with the European Communities Council Resolves of 22 September 2010 (2010/63/EU) and approved by the Danish Veterinary and Food Administration’s Council for Animal Experimentation (Journal No. 2022–15–2934–00156), and followed the ARRIVE guidelines.
[12]
R. A. Harshman, S. Hong, and M. E. Lundy, “Shifted factor analysis—Part I: Models and properties,” Journal of Chemometrics, 2003. DOI: 10.1002/cem.808
[13]
K. H. Madsen, L. K. Hansen, and M. Mørup, “Time and Frequency Domain Optimization with Shift, Convolution and Smoothness in Factor Analysis Type Decompositions,” Report, 2009.
[14]
R. Gu, Y. Rakita, L. Lan, Z. Thatcher, G. E. Kamm, D. O’Nolan, B. Mcbride, A. Wustrow, J. R. Neilson, K. W. Chapman, Q. Du, and S. J. L. Billinge, “Stretched non-negative matrix factorization,” npj Computational Materials, 2024. DOI: 10.1038/s41524- 02401377-5
[15]
S. A. Talwalkar and S. L. Marple, “Time-frequency scaling property of Discrete Fourier Transform (DFT),” in 2010 IEEE International Conference on Acoustics, Speech and Signal Processing, 2010. DOI: 10.1109/ICASSP.2010.5495902
[16]
M. Mørup, L. K. Hansen, S. M. Arnfred, L.-H. Lim, and K. H. Madsen, “Shift-invariant multilinear decomposition of neuroimaging data,” NeuroImage, 2008. DOI: 10.1016/j.neuroimage. 2008.05.062
[17]
A. S. Olsen, A. Brammer, P. M. Fisher, and M. Moerup, Uncovering dynamic human brain phase coherence networks, 2024. DOI: 10 . 1101/2024.11.15.623830
[18]
A. S. Olsen, J. D. Nielsen, and M. Mørup, “Coupled Generator Decomposition for Fusion of Electro- and Magnetoencephalography Data,” in 2024 32nd European Signal Processing Conference (EUSIPCO), 2024. DOI: 10 . 23919 / EUSIPCO63174 . 2024 . 10715032
[19]
C. Nilsson, F. Ståhlberg, C. Thomsen, O. Henriksen, M. Herning, and C. Owman, “Circadian variation in human cerebrospinal fluid production measured by magnetic resonance imaging,” The American Journal of Physiology, 1992. DOI: 10 . 1152 / ajpregu . 1992.262.1.R20
[20]
J. Paparrizos and L. Gravano, “K-Shape: Efficient and Accurate Clustering of Time Series,” SIGMOD Rec., 2016. DOI: 10.1145/ 2949741.2949758
Funding ASO and MLN were supported by JPND research, a Horizon 2020 supported EU joint programme. JLH and MM were supported by the Independent Research Fund Denmark, grant no. 10.46540/203500294B.
References [1]
H. Marquis, K. Willowson, and D. Bailey, “Partial volume effect in SPECT & PET imaging and impact on radionuclide dosimetry estimates,” Asia Oceania Journal of Nuclear Medicine and Biology, 2023. DOI: 10.22038/AOJNMB.2022.63827.1448
[2]
P. D. Acton, L. S. Pilowsky, H. F. Kung, and P. J. Ell, “Automatic segmentation of dynamic neuroreceptor single-photon emission tomography images using fuzzy clustering,” European Journal of Nuclear Medicine, 1999. DOI: 10.1007/s002590050425
[3]
J. Ashburner, J. Haslam, C. Taylor, V. J. Cunningham, and T. Jones, “CHAPTER 59 - A Cluster Analysis Approach for the Characterization of Dynamic PET Data,” in Quantification of Brain Function Using PET, 1996. DOI: 10.1016/B978-012389760-2/50061X
[4]
A.-E.-O. Boudraa, J. Champier, L. Cinotti, J.-C. Bordet, F. Lavenne, and J.-J. Mallet, “Delineation and quantitation of brain lesions by fuzzy clustering in Positron Emission Tomography,” Computerized Medical Imaging and Graphics, 1996. DOI: 10 . 1016 / 0895 6111(96)00025-0
[5]
M. K. Jaakkola, M. Rantala, A. Jalo, T. Saari, J. Hentilä, J. S. Helin, T. A. Nissinen, O. Eskola, J. Rajander, K. A. Virtanen, J. C. Hannukainen, F. López-Picón, and R. Klén, “Segmentation of Dynamic Total-Body [18F]-FDG PET Images Using Unsupervised Clustering,” International Journal of Biomedical Imaging, 2023. DOI: 10. 1155/2023/3819587
[6]
K. M. Kim, H. Watabe, M. Shidahara, J. Y. Ahn, S. Choi, N. Kudomi, K. Hayashida, Y. Miyake, and H. Iida, “Noninvasive estimation of cerebral blood flow using image-derived carotid input function in H/sub 2//sup 15/O dynamic PET,” in 2001 IEEE Nuclear Science Symposium Conference Record (Cat. No.01CH37310), 2001. DOI: 10.1109/NSSMIC.2001.1008569
[7]
J. S. Lee, D. Lee, S. Choi, K. S. Park, and D. S. Lee, “Non-negative matrix factorization of dynamic images in nuclear medicine,” in 2001 IEEE Nuclear Science Symposium Conference Record, 2001. DOI : 10.1109/NSSMIC.2001.1009222
[8]
O. Rainio, M. K. Jaakkola, and R. Klén, “Quantitative evaluation of unsupervised clustering algorithms for dynamic total-body PET image analysis,” Journal of Medical Engineering & Technology, 2025. DOI : 10.1080/03091902.2025.2466834
[9]
M. Mørup, K. H. Madsen, and L. K. Hansen, “Shifted Non-Negative Matrix Factorization,” in 2007 IEEE Workshop on Machine Learning for Signal Processing, 2007. DOI: 10 .1109 / MLSP .2007 . 4414296
[10]
A. J. Bell and T. J. Sejnowski, “An information-maximization approach to blind separation and blind deconvolution.,” Neural computation, 1995. DOI: 10.1162/neco.1995.7.6.1129
[11]
D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, 1999. DOI: 10 . 1038 / 44565