Fast Volume Alignment by Frequency-Marched Newton

F ast V olume Alignmen t b y F requency-Marc hed Newton F abian Kruse 1 , V alen tin Debarnot 2 , Vinith Kishore 1 , and Iv an Dokmani ´ c 1 1 Departmen t of Mathematics and Computer Science, Universit y of Basel, Basel, 4051, Switzerland 2 INSA-Ly on, Universite Claude Bernard Ly on 1, CNRS, Inserm, CREA TIS UMR 5220, U1294, Ly on, F rance Abstract W e dev elop a fast and accurate metho d for 3D alignmen t, reco vering the rotation and transla- tion that best align a reference v olume with a noisy observ ation. Classical matc hed ﬁltering ev al- uates cross-correlation ov er a large discretized transformation space; we sho w that high-precision alignmen t can b e ac hieved far more eﬃciently b y treating p ose estimation as a con tinuous op- timization problem. Our starting point is a band-limited Wigner– D expansion of the rotational correlation, which enables rapid ev aluation and eﬃcien t closed-form gradients and Hessians. Com- bined with analytical con trol of the complexit y of trigonometric-polynomial landscap es, this mak es second-order optimization practical in a setting where it is often av oided due to noncon vexit y and noise sensitivity . W e show that Newton-type reﬁnement is stable and eﬀective when initialized at lo w angular bandwidth: a coarse lo w-resolution SO(3) searc h provides robust candidates, which are then reﬁned by iterative frequency marching and Newton steps, with translations up dated via FFT in an alternating sc heme. W e pro vide a deterministic conv ergence guarantee sho wing that, under veriﬁable sp ectral-deca y and gap conditions, the frequency-marching scheme returns a near-optimal solution whose sub optimalit y is con trolled by the Newton tolerance. On synthetic rotation-estimation b enc hmarks, the metho d attains sub-degree accuracy while substantially re- ducing runtime relative to exhaustiv e SO(3) search. In tegrated into the subtomogram-av eraging pip eline of RELION5, it matches the baseline reconstruction quality , reaching lo cal resolution at the Nyquist limit, while reducing p ose-reﬁnemen t time by more than an order of magnitude. 1 In tro duction W e consider the problem of determining the p ose—–rotation and translation—–of a particle within a three-dimensional v olume, under sev ere noise and experimental artifacts. Our primary motiv ation comes from subtomogram av eraging (ST A) in cryogenic electron tomography (cryo-ET), where highly accurate particle alignmen t is essen tial to achiev e high resolution. In practice, alignmen t is executed once p er particle (and often iterativ ely), which mak es ev en mo dest per-particle cost comp ound and motiv ates the search for metho ds that are b oth accurate and computationally eﬃcient. A classical approac h to this problem is matche d ﬁltering , which seeks the transformation that maximizes the correlation b et ween a reference template and the observ ed data. Matc hed ﬁltering is kno wn to b e optimal for detecting a deterministic signal corrupted b y additive white Gaussian noise, in the sense that it maximizes the output signal-to-noise ratio [1, 2, 3]. Let g ∈ G denote the unknown transformation to reco ver. The matched ﬁltering ob jective is the correlation C( g ) def. = ⟨ f , g ◦ h ⟩ = Z R d f ( x )( g − 1 h )( x ) dx, (1) where f , h : R d → R denote the noisy measurement and the reference template, resp ectiv ely , with d ≥ 1; in this pap er d = 3. The p ose estimation problem can then b e formulated as g ⋆ ∈ argmax g ∈G C( g ) . (2) A straigh tforward strategy is to ev aluate C o v er a suitable discretization of G . This ma y b e eﬀective for a coarse search o ver low-dimensional p ose spaces but it suﬀers from p oor scaling as accuracy 1 requiremen ts or the dimensionalit y increase. In fav orable cases—–namely when the transformation group G admits a sampling theorem on a structured grid and the correlation can b e expressed as a group conv olution–FFT-based harmonic transforms allo w fast ev aluation on that grid. This is the case for translations (via the Euclidean FFT) and for rotations on SO(3) (via Wigner/SO(3) FFTs) [4, 5]. How ev er, accuracy is fundamentally limited by grid resolution, which is tied to the frequency bandwidth of the representation. Improving accuracy th erefore requires ﬁner discretizations and higher bandwidths, leading to substan tially higher computational and memory costs despite the fa vorable asymptotic scaling of FFTs. An alternative is to treat p ose estimation as a c ontinuous optimization pr oblem o ver G . When G has a diﬀerentiable structure—–such as SO(3) or R 3 –—the correlation function admits well-deﬁned gradien ts and Hessians, enabling eﬃcien t gradien t-based and second-order optimization metho ds. This is, ho wev er, rarely done: the correlation landscap e is highly non-conv ex with many spurious lo cal maxima, and it b ecomes increasingly ill-conditioned at high resolution, in the sense that the local curv ature near prominen t maxima b ecomes highly anisotropic (i.e., the Hessian has a large condition n umber), making second-order steps sensitive to noise and n umerical error. Success thus dep ends on eﬀectiv e initialization. T o address this challenge, w e exploit the fact that smo other input volumes induce smo other corre- lation landscap es. Because correlation is a linear op eration, band-limiting the inputs to angular degree L pro duces a band-limited correlation ob jectiv e C L on SO(3). Lo wer L yields a simpler landscape with t ypically fewer critical p oin ts and low er curv ature, so the basin of attraction of the global maximizer can b e identiﬁed via a coarse search at mo derate cost. W e call the resulting pro cedure Matcha . Matc ha ﬁrst p erforms a coarse exhaustive search at lo w bandwidth L 0 to identify a set of candidate rotations. The bandwidth is then progressively in- creased through a sc hedule L 0 < L 1 < · · · < L J —a strategy called fr e quency mar ching — and at eac h stage the candidates are reﬁned via Newton steps on C L j , initialized from the previous lev el. A t lo w and intermediate bandwidths, the correlation landscap e is smo oth and well-conditioned, allowing rapid conv ergence—t ypically 1–2 Newton steps p er bandwidth suﬃce (Section 4). Because G is low- dimensional (three-dimensional for SO(3)), the Hessian is a 3 × 3 matrix and can b e inv erted in closed form. By conﬁning exhaustive search to a single lo w-frequency stage and using con tinuous optimization for reﬁnemen t, this approach achiev es sub-degree angular precision at substantially reduced compu- tational cost compared with exhaustiv e high-resolution search. By conﬁning exhaustiv e search to a single low-frequency stage and using contin uous optimization for reﬁnemen t, this approach achiev es sub-degree angular precision at substantially reduced computational cost compared with exhaustive high-resolution searc h. Crucially , the ﬁnal accuracy is determined by the contin uous reﬁnemen t, not b y the resolution of the initial grid—so Matc ha can achiev e precision far b ey ond what the coarse SO(3)- FFT grid alone would p ermit. Figure 1 illustrates the principle b ehind Matcha. At lo w frequencies, the correlation landscap e is smooth and contains few lo cal maxima, enabling reliable coarse search. As higher frequencies are introduced, the landscap e b ecomes increasingly oscillatory . By tracking domi- nan t maxima across bandwidths and reﬁning them via lo cal optimization, w e recov er the optimal p ose without exhaustive search at full resolution. 1.1 Related w ork Estimating three-dimensional rotations is a fundamental problem in imaging and signal pro cessing. In electron tomograph y , this task is t ypically addressed through exhaustiv e searc h, which b ecomes computationally prohibitive at high angular resolution, esp ecially when com bined with the estima- tion of three-dimensional translations. A standard strategy is to lev erage fast F ourier transforms on transformation spaces [6]. S2FFT [7, 5] provides accelerated harmonic transforms both on the sphere and on SO(3); more broadly , FFTs on SO(3) based on Wigner- D expansions enable fast ev al- uation of correlations ov er dense rotational grids [4]. These metho ds remain exp ensiv e, how ever, in large-scale problems such as cryo-ET subtomogram av eraging. Our implementation builds up on the ball-harmonics expansion routines of Kileel et al. [8], extending their framework with a fully GPU- in tegrated PyT orc h pipeline and accelerated NUFFT-based pro jections via cuFINUFFT [9], enabling eﬃcien t batched pro cessing on mo dern GPUs. In the context of ST A for cryo-ET, a ma jor c hallenge is the joint estimation of b oth unkno wn ori- en tations and three-dimensional translations. This problem has traditionally b een addressed through exhaustiv e sampling that com bines fast F ourier transforms o ver rotations and spatial shifts [10, 11, 12]. 2 AAAsNXiczZpLc9vIEYC5m9eKee0ml1TlAkfLKm9E06I23my8xaq1HrZlWxKtl7UiFdcAGBKzBDDwYCCRRiH/Idfkh+S35JBbKtf8hfQMQBBAj/aVHKIqm0D3N6+e7unBAHbks1hubv79nXe/9/0f/PBH7621f/yTn/7s5+9/8IvzmCfCoWcO97m4sElMfRbSM8mkTy8iQUlg+/SVPdtR+lfXVMSMh6dyEdGrgExDNmEOkSB6Nf3jOJZEvH5/fbO3qf8sfNEvLtZbxd/w9Qcbvxq73EkCGkrHJ3E86m9G8iolQjLHp1l7nMQ0Is6MTOkILkMS0Pgq1f3NrA5IXGvCBfwLpaWl1RIpCeJ4EdhABkR6cVOnhCbdKJGTT69SFkaJpKGTNzRJfEtySw3ecpmgjvQXcEEcwaCvluMRQRwJJmq3251796yDW1tsj30qxzyioc3nY0F9Mm9y0gsy40iaQjXwuC51tbBBOg5FoN1ow4OhCUEnDUphMKSqjAaRR99kepynBHykPlBJ7MQnYg6Gs7apBKNYsqBiCq4VTqWXjulcCiL4jUfZ1JNZ+nEkVYH90AHHi2lexALAyon2OKQ3Dg8CErrpWGs9Slzlsdmof5WOg8SHmeB+EoRpP0udLB3H0CXq+vwGLO/7DonidL2fZVm9rmAx8TmRxGfTEAopU1HBwqnqzim3bGppR7th0rMocTxL49rxdBnlvTUDQEMqLKCZ4gqGnUSpmpeB7kndxLrjMNCG4TmfwShjbebHbJoIGlvQX4tB7DUsPhUk8pgzhxqKSxg5GDktyt2HIVvVAnFil52sKkY2GH+mbNu14hkL464lPcoFDSDspI4wcNpGZTckXoB3grQm9jlMtKoKRlBTBGRGHaqskHsQm72t+w8IutF0Evm89Nnizmc2OM0CBsyTqKJXRZY6pqYvhgBVhu86xHe6RIB54y6En2Dzbhwt6qbO15P22KUT8CZ9l7pEzAR1s1RM7Szd7H3S3exuGpipoDRcUt3N3gMjZfsJXUEKa0Ans8V2hek9+LgLK+enffX/H7YAXloAxgr9DWB6B/3eg16/5stjBgE8VcHnMAEL6GjrKu1YgrgsgZkslUwu1FxBlToIrVRd6mWQygBWNA4rVJjeXe9/tCGg1t/eXd/6SHmQtXF386HSWZH8SNUAIaADYkwsFlppf7O79aDb6/W6Hz/YzFIgvlMrlgVesXF3TKptZUroLGBQ7XzheeRPOYzEC2DFra11S3nWFEcxTSDbcJfmNZyFTNYX/vhNQgQNu08EWVylJ/uJIrLa4jfS2RFSVAQLuMvmTS13uproRnBbw2AAt/216wvSNfHzJU2C5XL/yB3ow7GyH7h0ep49/NBSa1m7Hr8hF7DA0OAqTeB/vXLfOd27gHWUS2tgBYSFPai1XemOtZO3G68kWmsFLLwfkPn9OInus3BiFf2LYcTV3hIxBTJTvqa7JwLI39Ns3B0noasiUapFNwVXKLSKzqzxZ9b6VmbBGvtIV6Djoagsr4hH1Qp1qWxFkbmRIvMqBd02USCuUjBAEwXikoJrzChhSZwcHpcKuM7yGbdOqPIxBRRq205BaZV/a8s/a20l7FjHlPhFwZ2y4M63K3hYFjw0F1Q3y3KHRCaiLHpZFr28tWi11cta6b2y9N7SPNulaDuvMBcPS/FwSYbEBiGkmTQd29x3lR48H8Q++TAr7e1S34SBWFaxF7CJNXFK7i5BiPnVGnpMy3k8pt/Y4DDZ1gQcnqjdAPEhWELYEsOmS+dsDllUFMZZVX9L7UuBqpzOo7wUY2Ux9vXF1B6BhZANrTAJbCryOpyyCufrq/DVpsuKI0rd3JPPIXHk27jGckUgYy0rzmdzkhK0y7INlI0ox0A5iHINlIsoaqAooiYGaoKoqYGaIsozUB6imIFiiPrSQH2JqJmBmiHKN1A+3gcbqABRoYEKEcUNFEdUZKAiRL0xUG8QJQyUQNENzwQGLka1SQMlEZUYqARR1wbqGlE3BuoGUXMDNUfUwkAtEPXWQL1VVB17ZMAeocq2DdQ2onYM1A6idg3ULqL2DNQeoh4bqMeIemKgniDqqYF6iqh9A7WPqGcG6hminhuo54h6YaBeIOrAQB0g6tBAHSLqyEAdIWpooIaIemmgXiLq2EAdI+rEQJ0g6tRAnSLqzECdIercQJ0j6pWBeoWoCwN1gagvDNQXiLo0UJc4uIkfebX8bQeFECdwKhGpZDhXErhoklqIE7narTVRLcTZPIqZz8MmXIjx4mbo7VtTbw2cCYMNHAa1EOd3jkklw8kbng4RqYU4g+vtapPNpTiRJ00wwFkpRFCIIR4weFBHRi/EOJ+zJhnhjY3weJMCEcJiNsVepIV4bggaC4hwtjZ7UHKLB0UeHo2Hh+NgzDFg0AiqLcYYD+gUjVoL8VYAVThnhuxtCvBHxgDfNnj4tsnBn5gC/IkxwHdNAb5rDPA98/Ts3TI9l4beXpp6u4e5PQN2agrwU2OA7xsCfN8U4M9NAf7cGOAvjAH+whzgB8jfD7C7HyLoEENH5gA/uiXAh8jnhtiHj3GAHxsC/MQU4CfGAD/FAX5qCPAzswed3eJBQxy5Q0Pk7mBsx4ANWYxqMwT4kSnAj4wBfoHavcgDXD+Oq/dBXB241PKpuyrhomQrQBmpd2vEb6oYmRZvUepbBq3I1CFtcz0TfF478yuf35Qie63O/PSJbn70cOPBEw5zqUckrujLTPtRI01wWYxFnXceo8GwWFbGqu6aj11ErIBzuGkuTVn14aLhwsRdadXN8pDvqKzx5Gh5xHfAYvU6g4SUJ835YHGU5T8+WcRy4dPKwcvX/nWsasnGLoNm+QExkynr0V7Xyr5N1eqQCEo19kLTVZ20N/0OdapSjfw4WdXpTP4nVUrucu2tULs+HVcvanQbMFujU77LH1rr/atMtdQYYFwrp0/Vy4IHJEwqBb/6wL6RZZ2Krz5CibWq3UbrS1W7g7y0qkV+ulfV7qFn36r2McriVe0T9Kxb1T5FGbCq3UfPtlXtM5QUq9rnKA9WtS9Q4qtqD1DGq2oP0cpb1R6h5buqHaJn1Kr2JVotnK9aqU6q2hOU2qraU5TRqtoz9OxZ1Z6jZ86q9hVKLlXtBXrGrGq/QFuvqvZSZaTiSL6ee/SRMw9X593qpqK4riqul4r9VS37rq67/n6v8ha69lLZnQgSwDpQk/osipMg7z4NrxnsadT5fErnJIh8qmqB9WB0lXbaa2M2iaWgQSQXKnmBBDriTvJ393C3tqZb0G/3BrlgbaxePI+WX84M7jqJEFC/ZfMkdFk4hYt5j5JYftTlifoQAro2UJ+45KVD7tIRCR2Pi4Giuuq7EhJOfdqdMN8f2D6M487WZo7D448UibT28s5nn4F11v7/+/nwT2BN6Gsn7+yqrywMqZA8CvQbucFWJLv+8oV4WWVXDUCJYfcgPQ1BET0KaIRq9Wq8xObXNJ6xaAA5N6DzSNwbezLvz+prmzWw29im0OjKbUZXuRI8m4KPl3L9CnWZCazT/GMIq+JLZZbQTlZ8LQEZCUY1y9Lj/HdUfI9wVYMKL8zS5YyaMRKDD+cfa6SPVtdm2KfqMSx9oX/MSPGbjTR7lRaDymoQbMsiHrO8pQIcVmQ12OGCwwZILEp0p5TUQP3JQ15Buru6NvczAK/zs/RA/6yQnHHsAhuRRHLLAUcGv72q9ltFd63PKhh8fpO7NxFEUncwIX5Mu+018DobPG5w4zFZ3GsPyB2xq77dySNqbIP7CkbjVaFckYNamAuWVdHQg9ihLlwSKdX3EtAz6lqa0h+V8cjy6UQO0kXsMfi9t9nrs7A7z+/g5gELMyheLae3hIMURCKBCwjWbtnnvGUokWl7ZPkb+1tttpoVbbLKxIDFOreZrNO0WecbGa1jtFoHma1TsVvnvzRc59tbrpMps6jTlDLhVBNI1ZFXzpMvWg/uLGe+2oS+U8tqX9/obqvVTKj3rMt1Tf/aXEoe6Mv2GtefXb6lxW3sERHBzAlYOpUH2lR9g2PlSx5UAA0oN5voNbwm1Jn0q/a0ld2tTsTXEHBlKo4jEqKT0zfKQmq9jBO13qmv+mKYppkoX82rwzWPRBTyT7iwAO9lWTrIirRebBt2t8tmdrfLbcDhahtwuBQOV8JhKdxdCXcP9fmXINMpdXMzFm1wOy4puC5bWX2zsX9cCndXwt2l8BlZvT2H66V45/RxKYbrpbif9sd3xnf6xS1zqNpvvy5REJRfKNhOar8ud0J5YoYH6PfX+80vZ/HF+Vav/0nvk5e/W/98u/iq9r3Wr1u/ad1t9Vu/b33eetoats5aTmvW+nPrL62/bvxt4x8b/9z4V46++05R5pet2t/Gv/8D9m52bQ== g ω AAAsOniczZpZc9vIEYC5OVfMtZu8pCovcLSssiOaFrXRZuMtVq112Jatg9ZlrUiVawAMiVkBGHgwkEijkJ+R1+SH5I/kNW+pvOYHpGcAggB6tFfyEFXZBLq/uXq6pwcD2JHPYrm+/vf3vvf9H/zwRz9+f6X9k5/+7Oe/+ODDX57HPBEOPXO4z8WFTWLqs5CeSSZ9ehEJSgLbp6/t622lf31DRcx4eCrnEb0KyDRkE+YQCaLxOCDSE0G6nb3Zf/PB6npvXf9Z+KJfXKy2ir/hmw/Xfj12uZMENJSOT+J41F+P5FVKhGSOT7P2OIlpRJxrMqUjuAxJQOOrVHc6szogca0JF/AvlJaWVkukJIjjeWADqToZN3VKaNKNEjn59CplYZRIGjp5Q5PEtyS3lAUslwnqSH8OF8QRDPpqOR4RxJFgp3a73Xn40Dq4s8X22KdyzCMa2nw2FtQnsyYnvSAzjqQpVAOP61JXCxuk41AE2o02PBiaEHTSoBQGQ6rKaBB59G2mx3lKwFHqA5XETnwiZmA4a4tKMIolCyqm4F/hVHrpmM6kIILfepRNPZmlH0dSFdgLHfC+mOZFLACsnGiPQ3rr8CAgoZuOtdajxFVum436V+k4SHyYCe4nQZj2s9TJ0nEMXaKuz2/B8r7vkChOV/tZltXrCuYTnxNJfDYNoZAyFRUsnKrunHLLppZ2tFsmPYsSx7M0rh1Pl1HeWzMANKRiA5oprmDYSZSqeRnontRNrDsOA20YnvNrGGWszfyUTRNBYwv6azEIwIbFp4JEHnNmUENxCSMHI6dFuUcwZKtaIE7sspNVxcgG418r23at+JqFcdeSHuWCBhB2UkcYOG2jslsSz8E7QVoT+xwmWlUFI6gpAnJNHaqskHsQu35X9x8QdKPpJPJ56bPFnc9scJo5DJgnUUWviix0TE1fDAGqDN91iO90iQDzxl0IP8Fm3Tia102dryftsUsn4E36LnWJuBbUzVIxtbN0vfdJd727bmCmgtJwQXXXe5tGyvYTuoQU1oBOrudbFaa3+XEXVs5P++r/P24AvLAAjBX6G8D0Dvq9zV6/5stjBgE8VcHnMAEL6GjjKu1YgrgsgZkslUzO1VxBlToIrVRd6mWQygBWNA4rVJjeX+0/WBNQ6+/ur248UB5krd1ff6x0ViQfqBogBHRAjInFQivtr3c3Nru9Xq/78eZ6lgLxnVqxLPCKtftjUm0rU0JnDoNq5wvPE3/KYSReACtuba1byLOmOIppAtmGuzSv4Sxksr7wx28TImjYfSbI/Co92UsUkdUWv5FOkZCiIljAXTZrarnT1UQ3gtsaBgO4669dX5BuiJ8vaRIsl/tH7kAf6WwLLp2eZ48/stRa1q7Hb8gFLDA0uEoT+F+v3PdOdy9gHeXSGlgBYWEPam1XumNt5+3GS4nWWgELHwVk9ihOokcsnFhF/2IYcbW3REyBzJSvLTYDIMrG3XESuioSpVp0U3CFQqvozBp/Zq1uZBassU90BToeisryinhUrVCXypYUmRkpMqtS0G0TBeIqBQM0USAuKbjGjBKWxMnhcamA6yyfceuEKh9TQKG27RSUVvm3svizVpbCjnVMiV8U3C4Lbn+7godlwUNzQXWzKHdIZCLKopdl0cs7i1ZbvayV3i1L7y7Ms1WKtvIKc/GwFA8XZEhsEEKaSdOxzX1X6cHzQeyTj7LS3i71TRiIZRXbh52siVNydwFCzC/X0GNazuMx/cYGh8m2JuDwRO0GiA/BEsaSwKZL52wOWVQUxllWf0ftC4GqnM6ivBRjZTH29cXUHoGFkA2tMAlsKvI6nLIK5+ur8NWmy4ojSt3ck88hceTbuMZyRSBjLSrOZ3OSErTLsg2UjSjHQDmIcg2UiyhqoCiiJgZqgqipgZoiyjNQHqKYgWKI+tJAfYmoawN1jSjfQPl4H2ygAkSFBipEFDdQHFGRgYoQ9dZAvUWUMFACRTc8Exi4GNUmDZREVGKgEkTdGKgbRN0aqFtEzQzUDFFzAzVH1DsD9U5RdeyJAXuCKtsyUFuI2jZQ24jaMVA7iNo1ULuIemqgniLqmYF6hqjnBuo5ovYM1B6iXhioF4h6aaBeImrfQO0j6sBAHSDq0EAdIurIQB0hamighoh6ZaBeIerYQB0j6sRAnSDq1ECdIurMQJ0h6txAnSPqtYF6jagLA3WBqC8M1BeIujRQlzi4iR95tfxtB4UQJ3AqEalkOFcSuGiSWogTudqtNVEtxNk8ipnPwyZciPHiZujtO1NvDZwJgw0cBrUQ53eOSSXDyRueDhGphTiD6+1qk82lOJEnTTDAWSlEUIghHjB4UEdGL8Q4n7MmGeGNjfB4kwIRwmI2xV6khXhuCBoLiHC2NntQcocHRR4ejYeH42DMMWDQCKotxhgP6BSNWgvxVgBVOGOG7G0K8CfGAN8yePiWycGfmQL8mTHAd0wBvmMM8F3z9OzeMT2Xht5emnq7i7ldA3ZqCvBTY4DvGQJ8zxTgL00B/tIY4PvGAN83B/gB8vcD7O6HCDrE0JE5wI/uCPAh8rkh9uFjHODHhgA/MQX4iTHAT3GAnxoC/MzsQWd3eNAQR+7QELnbGNs2YEMWo9oMAX5kCvAjY4BfoHYv8gDXj+PqfRBXBy61fOouS7go2QpQRurdGvGbKkamxVuU+pZBKzJ1SNtczwSf1c78yuc3pcjeqDM/faKbHz3cevCEw1zqEYkr+jLTftRIE1wWY1HnncdoMCyWlbGqu+ZjFxFL4BxumktTVn24aLgwcZdadbM45Dsqazw5WhzxHbBYvc4gIeVJcz5YHGX5j0/msZz7tHLw8rV/HatasrHLoFl+QMxkynq017Wyb1O1OiSCUo290HRZJ+1Nv0OdqlQjP06WdTqT/0mVkrtceyvUrk/H1Ysa3QbM1uiU7/DH1mr/KlMtNQYY18rpU/Wy4AEJk0rBrz6wb2RZp+KrT1BirWq30PpS1W4jL61qkZ/uVrW76Nm3qn2KsnhV+ww961a1z1EGrGr30LNtVfsCJcWq9iXKg1XtPkp8Ve0BynhV7SFaeavaI7R8V7VD9Ixa1b5Cq4XzVSvVSVV7glJbVXuKMlpVe4aePavac/TMWdW+Rsmlqr1Az5hV7Rdo61XVXqqMVBzJ13OPPnLm4fK8W91UFDdVxc1CsbesZc/Vddff71XeQtdeKrsTQQJYB2pSn0VxEuTdp+ENgz2NOp9P6YwEkU9VLbAejK7STntlzCaxFDSI5FwlL5BAR9xJ/u4e7lZWdAv67d4gF6yM1Yvn0eLzmcF9JxEC6rdsnoQuC6dwMetREssHXZ6oDyGgawP1iUteOuQuHZHQ8bgYKKqrvish4dSn3Qnz/YHtwzjubaznODz+SJFIazfvfPYZWGfl/7+fj/8E1oS+dvLOLvvKwpAKyaNAv5EbbESy6y9eiJdVdtUAlBh2D9LTEBTRo4BGqFYvx0tsfkPjaxYNIOcGdBaJh2NP5v1Zfm2zAnYb2xQaXbrN6CpXgmdT8PFSrl+hLjKBdZp/DGFVfKnMEtrJiq8lICPBqK6z9Dj/HRXfI1zVoMILs3Qxo2aMxODD+cca6ZPltRn2qXoMS/f1jxkpfrORZq/SYlBZDYJtWcRjlrdUgMOKrAY7XHDYAIl5iW6XkhqoP3nIK0h3ltfmfgbgdX6WHuifJZIzjl1gI5JIbjngyOC3V9V+q+iu9VkFg89vc/cmgkjqDibEj2m3vQJeZ4PHDW49Jot77QG5I3bVtzt5RI1tcF/BaLwslCtyUAtzwaIqGnoQO9SFSyKl+l4CekZdS1P6ozIeWT6dyEE6jz0Gvw/Xe30Wdmf5HdxssjCD4tVyeks4SEEkEriAYO2Wfc5bhhKZtkeWv7G/02bLWdEmq0wMWKxzl8k6TZt1vpHROkardZDZOhW7df5Lw3W+veU6mTKLOk0pE041gVQdeek8+aK1eW8x89Um9J1aVvv6RndbrWZCvWddrGv61+ZS8kBftle4/vbyHS1uY4+ICGZOwNKpPNCm6hscK1/yoAJoQLnZRK/hNaHOpF+1p63sbnUivoGAK1NxHJEQnZy+VRZS62WcqPVOfdUXwzRdi/LVvDpc80hEIf+EcwvwXpalg6xI68W2YWerbGZnq9wGHC63AYcL4XApHJbCnaVw51CffwkynVI3N2PRBrfjkoLrspXlNxt7x6VwZyncWQhfkOXbc7heiLdPn5ZiuF6I+2l/fG98r1/cMoeq/fabEgVB+YWC7aT2m3InlCdmeID+YLXf/HIWX5xv9Pqf9D559fvVz7eKr2rfb/2m9dvW/Va/9YfW563nrWHrrOW0otafW39p/XXtb2v/WPvn2r9y9HvvFWV+1ar9rf37P3OeeI8= C L AAAsL3iczZpZc9vIEYC5OVfMtZu8pCovcLSssiOaFrXRZuMtVq112Jatg9ZlrUiVawAMiVkBGHgwkEijkD+Q1+SH5Nek8pLKa/5FegYgCKBHeyUPUZVNoPubq6d7ejCAHfksluvr/3jve9//wQ9/9OP3V9o/+enPfv6LDz785XnME+HQM4f7XFzYJKY+C+mZZNKnF5GgJLB9+tq+3lb61zdUxIyHp3Ie0auATEM2YQ6RIHq1/+aD1fXeuv6z8EW/uFhtFX/DNx+u/XrscicJaCgdn8TxqL8eyauUCMkcn2btcRLTiDjXZEpHcBmSgMZXqe5pZnVA4loTLuBfKC0trZZISRDH88AGMiDSi5s6JTTpRomcfHqVsjBKJA2dvKFJ4luSW2rYlssEdaQ/hwviCAZ9tRyPCOJIME673e48fGgd3Nlie+xTOeYRDW0+Gwvqk1mTk16QGUfSFKqBx3Wpq4UN0nEoAu1GGx4MTQg6aVAKgyFVZTSIPPo20+M8JeAd9YFKYic+ETMwnLVFJRjFkgUVU3CqcCq9dExnUhDBbz3Kpp7M0o8jqQrshQ64XEzzIhYAVk60xyG9dXgQkNBNx1rrUeIqX81G/at0HCQ+zAT3kyBM+1nqZOk4hi5R1+e3YHnfd0gUp6v9LMvqdQXzic+JJD6bhlBImYoKFk5Vd065ZVNLO9otk55FieNZGteOp8so760ZABpSAQHNFFcw7CRK1bwMdE/qJtYdh4E2DM/5NYwy1mZ+yqaJoLEF/bUYRF3D4lNBIo85M6ihuISRg5HTotwjGLJVLRAndtnJqmJkg/GvlW27VnzNwrhrSY9yQQMIO6kjDJy2UdktiefgnSCtiX0OE62qghHUFAG5pg5VVsg9iF2/q/sPCLrRdBL5vPTZ4s5nNjjNHAbMk6iiV0UWOqamL4YAVYbvOsR3ukSAeeMuhJ9gs24czeumzteT9tilE/AmfZe6RFwL6mapmNpZut77pLveXTcwU0FpuKC6671NI2X7CV1CCmtAJ9fzrQrT2/y4Cyvnp331/x83AF5YAMYK/Q1gegf93mavX/PlMYMAnqrgc5iABXS0cZV2LEFclsBMlkom52quoEodhFaqLvUySGUAKxqHFSpM76/2H6wJqPV391c3HigPstburz9WOiuSD1QNEAI6IMbEYqGV9te7G5vdXq/X/XhzPUuB+E6tWBZ4xdr9Mam2lSmhM4dBtfOF54k/5TASL4AVt7bWLeRZUxzFNIFsw12a13AWMllf+OO3CRE07D4TZH6Vnuwlishqi99I50VIUREs4C6bNbXc6WqiG8FtDYMB3PXXri9IN8TPlzQJlsv9I3egj8bKfuDS6Xn2+CNLrWXtevyGXMACQ4OrNIH/9cp973T3AtZRLq2BFRAW9qDWdqU71nbebryUaK0VsPBRQGaP4iR6xMKJVfQvhhFXe0vEFMhM+Zrungggf0+zcXechK6KRKkW3RRcodAqOrPGn1mrG5kFa+wTXYGOh6KyvCIeVSvUpbIlRWZGisyqFHTbRIG4SsEATRSISwquMaOEJXFyeFwq4DrLZ9w6ocrHFFCobTsFpVX+rSz+rJWlsGMdU+IXBbfLgtvfruBhWfDQXFDdLModEpmIsuhlWfTyzqLVVi9rpXfL0rsL82yVoq28wlw8LMXDBRkSG4SQZtJ0bHPfVXrwfBD75KOstLdLfRMGYlnF9mH7auKU3F2AEPPLNfSYlvN4TL+xwWGyrQk4PFG7AeJDsISxJLDp0jmbQxYVhXGW1d9R+0KgKqezKC/FWFmMfX0xtUdgIWRDK0wCm4q8Dqeswvn6Kny16bLiiFI39+RzSBz5Nq6xXBHIWIuK89mcpATtsmwDZSPKMVAOolwD5SKKGiiKqImBmiBqaqCmiPIMlIcoZqAYor40UF8i6tpAXSPKN1A+3gcbqABRoYEKEcUNFEdUZKAiRL01UG8RJQyUQNENzwQGLka1SQMlEZUYqARRNwbqBlG3BuoWUTMDNUPU3EDNEfXOQL1TVB17YsCeoMq2DNQWorYN1DaidgzUDqJ2DdQuop4aqKeIemagniHquYF6jqg9A7WHqBcG6gWiXhqol4jaN1D7iDowUAeIOjRQh4g6MlBHiBoaqCGiXhmoV4g6NlDHiDoxUCeIOjVQp4g6M1BniDo3UOeIem2gXiPqwkBdIOoLA/UFoi4N1CUObuJHXi1/20EhxAmcSkQqGc6VBC6apBbiRK52a01UC3E2j2Lm87AJF2K8uBl6+87UWwNnwmADh0EtxPmdY1LJcPKGp0NEaiHO4Hq72mRzKU7kSRMMcFYKERRiiAcMHtSR0QsxzuesSUZ4YyM83qRAhLCYTbEXaSGeG4LGAiKcrc0elNzhQZGHR+Ph4TgYcwwYNIJqizHGAzpFo9ZCvBVAFc6YIXubAvyJMcC3DB6+ZXLwZ6YAf2YM8B1TgO8YA3zXPD27d0zPpaG3l6be7mJu14CdmgL81Bjge4YA3zMF+EtTgL80Bvi+McD3zQF+gPz9ALv7IYIOMXRkDvCjOwJ8iHxuiH34GAf4sSHAT0wBfmIM8FMc4KeGAD8ze9DZHR40xJE7NETuNsa2DdiQxag2Q4AfmQL8yBjgF6jdizzA9eO4eh/E1YFLLZ+6yxIuSrYClJF6t0b8poqRafEWpb5l0IpMHdI21zPBZ7Uzv/L5TSmyN+rMT5/o5kcPtx484TCXekTiir7MtB810gSXxVjUeecxGgyLZWWs6q752EXEEjiHm+bSlFUfLhouTNylVt0sDvmOyhpPjhZHfAcsVq8zSEh50pwPFkdZ/uOTeSznPq0cvHztX8eqlmzsMmiWHxAzmbIe7XWt7NtUrQ6JoFRjLzRd1kl70+9QpyrVyI+TZZ3O5H9SpeQu194KtevTcfWiRrcBszU65Tv8sbXav8pUS40BxrVy+lS9LHhAwqRS8KsP7BtZ1qn46hOUWKvaLbS+VLXbyEurWuSnu1XtLnr2rWqfoixe1T5Dz7pV7XOUAavaPfRsW9W+QEmxqn2J8mBVu48SX1V7gDJeVXuIVt6q9ggt31XtED2jVrWv0GrhfNVKdVLVnqDUVtWeooxW1Z6hZ8+q9hw9c1a1r1FyqWov0DNmVfsF2npVtZcqIxVH8vXco4+cebg871Y3FcVNVXGzUOwta9lzdd3193uVt9C1l8ruRJAA1oGa1GdRnAR592l4w2BPo87nUzojQeRTVQusB6OrtNNeGbNJLAUNIjlXyQsk0BF3kr+7h7uVFd2Cfrs3yAUrY/XiebT4ZmZw30mEgPotmyehy8IpXMx6lMTyQZcn6kMI6NpAfeKSlw65S0ckdDwuBorqqu9KSDj1aXfCfH9g+zCOexvrOQ6PP1Ik0trNO599BtZZ+f/v5+M/gTWhr528s8u+sjCkQvIo0G/kBhuR7PqLF+JllV01ACWG3YP0NARF9CigEarVy/ESm9/Q+JpFA8i5AZ1F4uHYk3l/ll/brIDdxjaFRpduM7rKleDZFHy8lOtXqItMYJ3mH0NYFV8qs4R2suJrCchIMKrrLD3Of0fF9whXNajwwixdzKgZIzH4cP6xRvpkeW2Gfaoew9J9/WNGit9spNmrtBhUVoNgWxbxmOUtFeCwIqvBDhccNkBiXqLbpaQG6k8e8grSneW1uZ8BeJ2fpQf6Z4nkjGMX2IgkklsOODL47VW13yq6a31WweDz29y9iSCSuoMJ8WPaba+A19ngcYNbj8niXntA7ohd9e1OHlFjG9xXMBovC+WKHNTCXLCoioYexA514ZJIqb6XgJ5R19KU/qiMR5ZPJ3KQzmOPwe/D9V6fhd1Zfgc3myzMoHi1nN4SDlIQiQQuIFi7ZZ/zlqFEpu2R5W/s77TZcla0ySoTAxbr3GWyTtNmnW9ktI7Rah1ktk7Fbp3/0nCdb2+5TqbMok5TyoRTTSBVR146T75obd5bzHy1CX2nltW+vtHdVquZUO9ZF+ua/rW5lDzQl+0Vrj+4fEeL29gjIoKZE7B0Kg+0qfoGx8qXPKgAGlBuNtFreE2oM+lX7Wkru1udiG8g4MpUHEckRCenb5WF1HoZJ2q9U1/1xTBN16J8Na8O1zwSUcg/4dwCvJdl6SAr0nqxbdjZKpvZ2Sq3AYfLbcDhQjhcCoelcGcp3DnU51+CTKfUzc1YtMHtuKTgumxl+c3G3nEp3FkKdxbCF2T59hyuF+Lt06elGK4X4n7aH98b3+sXt8yhar/9pkRBUH6hYDup/abcCeWJGR6gP1jtN7+cxRfnG73+J71PXv1+9fOt4qva91u/af22db/Vb/2h9XnreWvYOms5Ldr6c+svrb+u/W3t72v/XPtXjn7vvaLMr1q1v7V//wf/6XOi L ▲ First lo cal maximum ▶ Second lo cal maximum ▼ Third lo cal maximum ◦ Best lo cal maxima Newton method Exhaustiv e search Bandwidth Correlation landscape Figure 1: Evolution of a 1D correlation landscap e under progressive bandwidth expansion. Each curve cor- resp onds to a diﬀerent bandwidth, with low er bandwidths pro ducing smo other ob jectives and typically fewer lo cal maxima. A coarse exhaustive search (ev aluated on a discretization depicted in green) at low bandwidth iden tiﬁes lo cal maxima. As the bandwidth increases, these lo cal maxima are trac ked across scales using fre- quency marching, while lo cal optimization reﬁnes the solution at each stage. This strategy av oids exhaustive searc h on highly oscillatory high-frequency landscapes while achieving high-precision alignment. F or ph ysical in terpretation of low-pass ﬁltering the correlation, we display corresp ondingly low-pass ﬁltered copies of a par- ticle on the left. Sev eral extensions incorp orate physical characteristics of the imaging pro cess, including the missing w edge [13] and spatially v arying densit y v ariance [14]. Chen et al. [15] introduce a fast reference-free subtomogram alignment metho d that av oids an external template by iterativ ely aligning subtomograms to one another. Con tinuous p ose reﬁnement has also b een explored in cryo-EM. Zehni et al. [16] use alternating gradien t-based optimization and ADMM for join t map–orientation estimation. This approach ad- dresses the full reconstruction problem rather than the isolated rotational correlation considered here, and do es not use second-order reﬁnement or frequency marching. Related coarse-to-ﬁne ideas do ap- p ear elsewhere in cryo-EM: Barnett et al. [17] prop ose a frequency-marching strategy for ab initio single-particle reconstruction, marching outw ard in F ourier radius from lo w to high frequency . Our use of frequency marc hing diﬀers inboth setting and scop e: we apply it to band-limited rotational cor- relation on SO(3), use it to track candidate basins via Newton reﬁnement, and provide deterministic conditions under which basin trac king succeeds. 1.2 P ap er outline Section 2 introduces the frequency-marched alignmen t algorithm, including the low-frequency ex- haustiv e search and the second-order contin uous optimization. In Section 3, we establish suﬃcien t conditions—form ulated in terms of sp ectral decay and inter-basin gaps—under which Matcha is guar- an teed to trac k the globally optimal basin through the bandwidth schedule and return a near-optimal estimate; we also quantify the adv antage of marching ov er a single bandwidth jump. In Section 4, we ev aluate accuracy and run time on syn thetic b enc hmarks, demonstrating sp eedups of several orders of magnitude ov er exhaustive searc h for isolated rotation estimation at v arying SNR, and showing that the metho d is broadly applicable to any domain requiring fast SO(3) correlation. W e then integrate Matc ha into a cry o-ET s ubtomogram av eraging pip eline. On exp erimen tal data comprising appro xi- mately 21 , 000 subtomograms, w e demonstrate that fast sub-degree rotational reﬁnement enables ST A reconstructions approaching the Nyquist limit (3 . 8 ˚ A), while reducing alignment run time by an order of magnitude compared with RELION5’s built-in reﬁnement. 3 2 Metho d W e start by in tro ducing the alignmen t problem and its sp ectral form ulation; w e then present our coarse- to-ﬁne optimization strategy . Motiv ated by subtomogram alignment in cryo-electron tomograph y , we fo cus on G = SO(3). The full alignmen t problem also includes translations; for ease of exp osition w e defer their treatment to App endix C. W e parameterize rotations g ∈ G by ZYZ Euler angles following [4], g ≡ g ( α, β , γ ) = r z ( α ) r y ( β ) r z ( γ ) , (3) where r z ( α ) =    cos α − sin α 0 sin α cos α 0 0 0 1    , r y ( β ) =    cos β 0 sin β 0 1 0 − sin β 0 cos β    , and ( α , β , γ ) ∈ [0 , 2 π ) × [0 , π ] × [0 , 2 π ). The rotation g ∈ G acts on a function f : R 3 → R as ( g ◦ f )( x ) = f  g − 1 x  . W e let B 3 ⊂ R 3 denote the unit ball and let f , h ∈ L 2 ( B 3 ) denote the observed subtomogram and the reference template, resp ectiv ely . W e deﬁne the L 2 ( B 3 ) inner pro duct by ⟨ f , h ⟩ def. = R B 3 f ( x ) h ( x ) dx . The orientation estimation problem is formulated as ˆ g ∈ argmax g ∈G C( g ) , C( g ) def. = ⟨ f , g ◦ h ⟩ . (4) 2.1 Ball-harmonic decomp osition of the rotational cross-correlation W e decomp ose f and h in a ball-harmonic basis, which separates radial and angular v ariables. Under a rotation,the radial basis functions remain unc hanged, while for each degree l , the spherical-harmonic co eﬃcien ts are mixed by the Wigner– D matrices. After summing ov er the radial index, the rotational cross-correlation can therefore b e expressed as a F ourier series on SO(3) in the Wigner– D basis. T runcating the angular conten t at degree L , we appro ximate the rotational cross-correlation (4) by C L ( g ) def. = L X l =0 X | m |≤ l X | m ′ |≤ l σ l,m,m ′ D l m,m ′ ( g ) , g ∈ G . (5) Here D l m,m ′ are the matrix co eﬃcien ts of the irreducible unitary representation of SO(3) of degree l . By the Peter–W eyl theorem, they form the analogue of the F ourier basis on SO(3). The expansion in (5) is therefore a F ourier series for the rotational cross-correlation, and truncation at degree L amounts to retaining only frequencies up to L . In particular, C L → C in L 2 (SO(3)) as L → ∞ . F or a ﬁxed rotation g , a direct ev aluation of (5) requires O ( L 3 ) op erations, since P L l =0 (2 l + 1) 2 = O ( L 3 ). This representation is also conv enient for optimization: deriv atives of C L with resp ect to the rotation parameters can b e obtained in closed form from the deriv atives of Wigner– D matrices, with the same asymptotic cost as ev aluating C L itself. In the prop osed coarse-to-ﬁne scheme, w e start from a small cutoﬀ L 0 —yielding a smo othed, low- angular-frequency ob jective—and progressiv ely increase L , thereby enric hing the angular conten t of C L ( g ). At eac h stage, the approximation at the new bandwidth is obtained simply by extending the sum ov er the precomputed co eﬃcien ts σ l,m,m ′ to include higher degrees; no recomputation of lo wer-degree co eﬃcien ts is needed. 2.2 Con tin uous optimization Because the cross-correlation C is non-conv ex, it is in general imp ossible to ﬁnd the global maximum of C directly by ﬁrst or second order contin uous optimization. It b ecomes p ossible, how ever, provided that a go od initial estimate is av ailable. Finding such an initial estimate is in turn easier for the bandlimited cross-correlation C L , since lo wer-bandwidth functions generically ha ve few er critical points and a smo other landscap e—th us a larger basin of attraction for the global maximum. 4 Once one decides to attempt a contin uous reﬁnement, a natural mo dern approac h is to use au- tomatic diﬀerentiation and gradient ascent: it only requires gradients and typically makes progress from relatively crude initializations. On landscap es such as bandlimited trigonometric p olynomials, ho wev er, second-order metho ds can b e orders of magnitude faster, since gradien t-based metho ds may require many small steps to exploit curv ature information. This b ecomes esp ecially relev an t in our setting, where the cost of ev aluating the ob jectiv e and its deriv ativ es is dominated b y the same O ( L 3 ) Wigner– D terms, making the asymptotic cost complexit y scaling in terms of L the same for b oth ﬁrst- and second-order metho ds. Moreo ver, b ecause G has a ﬁxed, low dimension (three for SO(3)), the Newton step only requires in verting a 3 × 3 matrix, which can b e done in closed form. Iteration coun t is therefore the main factor determining the diﬀerence in runtime. Gradient ascent t ypically requires many more iterations to reach high precision, whereas Newton’s metho d typically con verges in only a few iterations once initialized in the desired basin of attraction (see Figure 4c). While S O (3) admits an in trinsic Riemannian geometry (Riemannian gradien t/Hessian and a re- traction/ exp onen tial-map up date), we found that working in a lo cal coordinate chart is suﬃcient in practice. In our implementation we optimize directly in wrapp ed Euler angles θ = ( α, β , γ ) and map to rotations via g ( θ ). Since contin uous optimization is initialized with a rotation near a maximizer, iterates remain in a lo cal neighborho o d where this parameterization is well-conditioned. In this chart, the up date takes the familiar form θ t +1 = θ t − H − 1 C L ( θ t ) ∇ C L ( θ t ) , g t = g ( θ t ) , where H C L ( θ ) def. = ∇ 2 C L ( θ ) denotes the Hessian. W e remark that closed-form deriv ativ es are sub- stan tially more eﬃcien t than automatic diﬀerentiation in this setting: they av oid constructing and bac kpropagating through the Wigner- D matrices and provide direct access to the Hessian, enabling fast second-order reﬁnement. The conv ergence analysis in Section 3 is formulated in trinsically on SO(3); the Euler-angle parameterization is used only in the implementation. 2.3 Initialization via low-resolution SO(3)-FFT The success of con tinuous optimization dep ends on suitable initialization. W e obtain one b y p erforming a grid search on the low-bandwidth correlation C L 0 , which is smo oth and inexp ensiv e to ev aluate. W e use an SO(3) fast F ourier transform (SOFFT) [4, 18, 5] to ev aluate C L 0 on a grid of O (( K L 0 ) 3 ), where K ≥ 1 is the o v ersampling factor. On the natural equiangular grid, classical algorithms achiev e a complexit y of O ( L 4 0 ) [4]; with o versampling factor K and M = ( K L 0 ) 3 target rotations, nonequispaced v arian ts achiev e O ( M + L 4 0 ) [18], compared with O ( M · L 3 0 ) for direct ev aluation. W e retain the N C strongest lo cal maxima as candidate initializations; although their angular accuracy is limited by the grid spacing (a few degrees in our examples), it is suﬃcient to initialize the subsequent con tinuous optimization stages. 2.4 Coarse-to-ﬁne optimization algorithm The complete alignment pro cedure com bines a single exhaustive search at L 0 with frequency-marc hed con tinuous optimization. The ﬁnal estimate is rep orted at angular cutoﬀ L J , chosen such that C L J appro ximates the full correlation to the desired accuracy . Algorithm 1 summarizes the pro cedure. In practice, we run at most M iter Newton steps p er band and terminate early when the gradient norm, the parameter increment, or the relativ e c hange in the ob jectiv e fall b elo w a prescribed tolerance. W e use a ﬁxed bandwidth sc hedule starting from L 0 = 30, for which the exhaustive SO(3)-FFT search remains inexp ensiv e, and ending at the target bandwidth L J . Intermediate bandwidths are chosen to k eep the p er-step p erturbation within the stabilit y b ounds of our conv ergence analysis (Section 3.4). W e retain N C = 10 candidates with ov ersampling factor K = 2 for the initial SOFFT search. Newton reﬁnemen t is preferable to gradien t ascent because b oth metho ds hav e the same asymptotic p er- iteration cost, while Newton t ypically reaches high precision in substan tially fewer iterations once initialized in the correct basin (see Figure 4c). 5 Algorithm 1 Matcha, a Multi-band Angular T emplate Matching Algorithm Require: Angular cutoﬀs { L 0 , L 1 , . . . , L J } Require: Number of Newton iterations M iter Require: Number of candidates N C Require: Oversampling factor K 1: Perform exhaustiv e search with SOFFT at angular cutoﬀ L 0 with ov ersampling factor K 2: Retain the N C strongest lo cal maxima { g (0) n } N C n =1 3: for j in { 1 , . . . , J } do 4: for n in { 1 , . . . , N C } do 5: g ( j ) n ← result of M iter Newton steps on C L j initialized at g ( j − 1) n 6: end for 7: end for 8: return argmax n ∈{ 1 ,...,N C } C L J  g ( J ) n  3 Theoretical analysis The eﬃciency of Matcha (Algorithm 1) results from several factors: t ypical correlation landscap es are b enign and evolv e smo othly ov er bandwidths, the ov erall problem is lo w-dimensional with a controlled n umber of lo cal maxima to track, all ob jects inv olv ed in the second-order optimization are av ailable in eﬃcient closed form, and second order optimization is lo cally very fast. In this section we prov e that Algorithm 1 (appro ximately) recov ers the global maxim um of the correlation under suitable assumptions. Giv en these assumptions, our result is deterministic. It is p ossible that stronger guaran tees can b e given “with high probabilit y” ov er a suitable ensemble of random particles, but we lea ve this to future work. Our pro of strategy relies on the intuition that the p osition of local maxima do es not c hange muc h if the bandwidth L is increased b y a little. T urning this intuition into a frequency-marching pro of requires handling tw o main failure mo des. One is that a track ed lo cal maximum splits into t wo or more lo cal maxima. The other is that as we increase the bandwidth, a new un track ed lo cal maxim um emerges which even tually grows into a global maximum. W e will make suitable assumptions under whic h this cannot happ en. 3.1 Setting W e select a bandwidth schedule L 0 < L 1 < · · · < L J and write F j def. = C L j , j = 0 , 1 , . . . , J for the correlation function at bandwidth L j . In order to exclude the p ossibilit y that a new maxim um emerges outside the track ed union of balls S and even tually b ecomes global, we will need to assume that the global maximum at bandwidth L J is suﬃciently prominent. A technical to ol for this is the set gap. Deﬁnition 3.1. L et F : SO(3) → R and S ⊂ SO(3) b e c omp act. The set gap of F with r esp e ct to S is Γ F ( S ) def. = max g ∈ S F ( g ) − sup g / ∈ S F ( g ) . (6) Note that if Γ F ( S ) > 0 then S contains a global maximum of F . 3.2 F ailure mo de elimination The ﬁrst lemma shows that under suitable assumptions (i.e. set gap and sp ectral decay with frequency) it is imp ossible that a maxim um absent at bandwidth L 0 b ecomes global at bandwidth L J . Lemma 3.1 (Set-gap stabilit y) . L et F , e F : SO(3) → R b e c ontinuous and let S ⊂ SO(3) b e c omp act. Assume the set gap Γ F ( S ) deﬁne d in (6) satisﬁes Γ F ( S ) > 0 . If   F − e F   ∞ ≤ ε and ε < Γ F ( S ) , (7) 6 then every glob al maximizer ˜ g ⋆ ∈ argmax e F lies in S . Pr o of. Deﬁne a S := max g ∈ S F ( g ) and a S c := sup g / ∈ S F ( g ), so a S − a S c = Γ F ( S ) > 0 b y assumption. Let g S b e such that F ( g S ) = a S (it exists since S is a ﬁnite union of closed balls hence compact). F or an y g / ∈ S , e F ( g ) ≤ F ( g ) + ε ≤ a S c + ε ≤ a S − Γ F ( S ) + ε, (8) and also e F ( g S ) ≥ F ( g S ) − ε = a S − ε. Subtracting, e F ( g S ) − e F ( g ) ≥ ( a S − ε ) − ( a S − Γ F ( S ) + ε ) = Γ F ( S ) − 2 ε > 0 . Th us e F ( g S ) > e F ( g ) for g / ∈ S . In other words: any global maximizer of e F must b e in S . Since SO(3) is compact and e F is contin uous, e F attains its maxim um. In practice, this lemma is a wa y to exploit the sp ectral decay with frequency . Concretely , if the high-frequency tail ∥ C L J − C L 0 ∥ ∞ is smaller than half the set gap at bandwidth L 0 , then no newly b orn maxim um outside the trac ked union of balls S can b ecome global. Next we establish basin regularit y . Deﬁnition 3.2. Str ong c onc avity and smo othness L et U ⊂ SO(3) b e ge o desic al ly c onvex and F ∈ C 2 ( U ) . We say F is µ -str ongly c onc ave and M -smo oth on U if for every g ∈ U and every unit tangent ve ctor X ∈ T g SO(3) , − M ≤  ∇ 2 F ( g )[ X ] , X  ≤ − µ. On any Riemannian manifold suﬃcien tly small geo desic balls are geo desically conv ex. Lemma 3.2 (Second-order tw o-sided b ounds from a critical p oin t) . L et U b e ge o desic al ly c onvex, F ∈ C 2 ( U ) , and let g 0 ∈ U satisfy ∇ F ( g 0 ) = 0 . Assume F is µ -str ongly c onc ave and M -smo oth on U . Then for every g ∈ U , F ( g 0 ) − M 2 d ( g , g 0 ) 2 ≤ F ( g ) ≤ F ( g 0 ) − µ 2 d ( g , g 0 ) 2 . Pr o of. Fix g ∈ U and let γ : [0 , ℓ ] → U be the unit-sp eed minimizing geodesic from γ (0) = g 0 to γ ( ℓ ) = g , so ℓ = d ( g , g 0 ). Let h ( t ) = F ( γ ( t )). Then h ′ (0) = ⟨∇ F ( g 0 ) , ˙ γ (0) ⟩ = 0 and h ′′ ( t ) =  ∇ 2 F ( γ ( t ))[ ˙ γ ( t )] , ˙ γ ( t )  . Since ∥ ˙ γ ( t ) ∥ = 1 and F is µ -strongly concav e and M -smo oth, we hav e for all t ∈ [0 , ℓ ], − M ≤ h ′′ ( t ) ≤ − µ. In tegrate h ′′ ( t ) ≥ − M from 0 to t to get h ′ ( t ) ≥ − M t , and integrate again from 0 to ℓ : h ( ℓ ) − h (0) = Z ℓ 0 h ′ ( t ) dt ≥ Z ℓ 0 ( − M t ) dt = − M 2 ℓ 2 . Th us h ( ℓ ) ≥ h (0) − M 2 ℓ 2 , i.e. F ( g ) ≥ F ( g 0 ) − M 2 d ( g , g 0 ) 2 . The upp er b ound is identical using h ′′ ( t ) ≤ − µ . Lemma 3.3 (Basin p ersistence under a C 2 p erturbation) . L et U = B ( g ⋆ , r ) b e a close d ge o desic al ly c onvex b al l. L et Ω ⊂ SO(3) b e an op en neighb orho o d of U , and let F, E ∈ C 2 (Ω) . Assume F is µ -str ongly c onc ave on U and that g ⋆ is the ne c essarily unique maximizer of F on U . L et e F def. = F + E and set a def. = ∥∇ E ∥ L ∞ ( U ) , b def. = sup g ∈ U ∥∇ 2 E ( g ) ∥ op . Assume b < µ and a < ( µ − b ) r / 2 . Then e F is ( µ − b ) -str ongly c onc ave on U with an unique maximizer ˜ g ⋆ ∈ int( U ) such that d ( ˜ g ⋆ , g ⋆ ) ≤ a/ ( µ − b ) . 7 Pr o of. T o show strong concavit y note that for any unit tangen t X at any g ∈ U , D ∇ 2 e F ( g )[ X ] , X E =  ∇ 2 F ( g )[ X ] , X  +  ∇ 2 E ( g )[ X ] , X  ≤ − µ + b = − ( µ − b ) . By contin uity and compactness, e F attains a maxim um on U . Because e F is ( µ − b )-strongly concav e on the geo desically conv ex U , it can hav e at most one maximizer in U (if there were tw o, strict concavit y along the connecting geo desic would b e con tradicted). It remains to rule out a maximizer on ∂ U . Let γ : [0 , r ] → U b e any unit-sp eed geo desic with γ (0) = g ⋆ and γ ( r ) ∈ ∂ U . Deﬁne h ( t ) = e F ( γ ( t )). Then h ′ (0) = D ∇ e F ( g ⋆ ) , ˙ γ (0) E = ⟨∇ E ( g ⋆ ) , ˙ γ (0) ⟩ so | h ′ (0) | ≤ a . W e already sho wed h ′′ ( t ) ≤ − ( µ − b ). In tegrating h ′′ ( t ) ≤ − ( µ − b ) t wice gives h ( r ) ≤ h (0) + r h ′ (0) − µ − b 2 r 2 ≤ h (0) + ar − µ − b 2 r 2 . Since a < ( µ − b ) r/ 2, w e hav e ar − µ − b 2 r 2 < 0, hence h ( r ) < h (0). Th us no b oundary point can maximize e F , so the maximizer lies in int( U ). T o b ound the displacement, let ˜ g ⋆ b e the unique maximizer of e F in U and γ : [0 , ℓ ] → U the unit-sp eed minimizing geo desic from γ (0) = g ⋆ to γ ( ℓ ) = ˜ g ⋆ , so ℓ = d ( g ⋆ , ˜ g ⋆ ). Deﬁne ψ ( t ) = D ∇ e F ( γ ( t )) , ˙ γ ( t ) E . Then ψ ( ℓ ) = 0 because ∇ e F ( ˜ g ⋆ ) = 0, and ψ ′ ( t ) = D ∇ 2 e F ( γ ( t ))[ ˙ γ ( t )] , ˙ γ ( t ) E ≤ − ( µ − b ) . In tegrating gives 0 = ψ ( ℓ ) ≤ ψ (0) − ( µ − b ) ℓ , hence ( µ − b ) ℓ ≤ ψ (0). But ψ (0) = D ∇ e F ( g ⋆ ) , ˙ γ (0) E = ⟨∇ E ( g ⋆ ) , ˙ γ (0) ⟩ ≤ ∥∇ E ( g ⋆ ) ∥ ≤ a . Therefore ℓ ≤ a/ ( µ − b ). 3.3 Main result W e now analyze the multi-band marching scheme in Algorithm 1: w e start with a lo w-bandwidth exhaustiv e searc h at L 0 , keep P = N C candidate basins, and for j = 1 , . . . , J reﬁne each candidate by M iter Newton steps for F j = C L j , then select the b est candidate at L J . The key pro of obligations are: (i) basins do not split internally as j increases; (ii) no untr acke d maxim um can b ecome a global maxi- m um at L J ; (iii) b ecause Newton is run for a ﬁnite n umber of steps, the output is only approximately optimal. Assumption 3.1. W e make the following assumptions (A) T rack ed basins at L 0 and con vexit y radius. Fix r > 0 such that ev ery closed geo desic ball B ( g, ρ ) with ρ ≤ r is geo desically conv ex. There exist p oin ts g p, 0 ∈ SO(3), p = 1 , . . . , P , suc h that the closed balls U p def. = B ( g p, 0 , r ) ⊂ SO(3) , p = 1 , . . . , P , are pairwise disjoint, and g p, 0 is the unique maximizer of F 0 on U p . W e let S def. = S P p =1 U p . (B) Small p er-step C 2 c hanges on S with non-degeneracy budget. F or j ≥ 1, let E j def. = F j − F j − 1 and a j def. = ∥∇ E j ∥ L ∞ ( S ) , b j : = sup g ∈ S ∥∇ 2 E j ( g ) ∥ op , δ j def. = ∥ E j ∥ ∞ . W e hav e J X j =1 b j ≤ µ 2 , a j < µr 8 , J X j =1 2 a j µ ≤ r 2 , (9) and F 0 is µ -strongly concav e on each U p . (C) Set-gap prev ents late-b orn global maximizers outside S . W e hav e Γ F 0 ( S ) > 0 and 2 J X j =1 δ j < Γ F 0 ( S ) . (10) It thus follows ∥ F J − F 0 ∥ ∞ ≤ P j δ j < Γ F 0 ( S ) / 2. 8 (D) Finite Newton accuracy . F or each j = 1 , . . . , J and each p , after M iter Newton steps on F j initialized from the previous lev el, the returned iterate g ( j ) p satisﬁes d  g ( j ) p , m p,j  ≤ τ , τ ≤ r 2 (11) where m p,j denotes the unique maximizer of F j on U p (whose existence/uniqueness is prov ed b elo w). (E) Final in ter-basin winner is unique and separated. Deﬁne p ⋆ so that F J ( m p ⋆ ,J ) = max p =1 ,...,P F J ( m p,J ) , and assume this winner is unique. Let γ def. = F J ( m p ⋆ ,J ) − max p  = p ⋆ F J ( m p,J ) > 0 . (F) Curv ature b ound conv erts distance-to-optimum in to v alue-to-optim um. Let M J def. = sup g ∈ S   ∇ 2 F J ( g )   op < ∞ . (ﬁnite b ecause F J ∈ C 2 and S is compact). The Newton tolerance satisﬁes M J 2 τ 2 ≤ γ 4 . (12) Theorem 3.4 (Main result) . The output of Matcha (A lgorithm 1), ˆ g ∈ argmax p =1 ,...,P F J  g ( J ) p  is ne ar-optimal; namely, d ( ˆ g , m p ⋆ ,J ) ≤ τ and 0 ≤ max g ∈ SO(3) F J ( g ) − F J ( ˆ g ) ≤ M J 2 τ 2 . Pr o of. P art (i): W e start by proving existence/uniqueness of basin maximizers for all j , without splitting. Fix a basin index p . W e prov e by induction on j that F j has a unique maximizer m p,j ∈ in t( U p ) and that d ( m p,j , m p, 0 ) ≤ j X k =1 2 a k µ . (13) Base j = 0 holds by Assumption (A). Induction step: assume the claim for j − 1. Let µ j − 1 def. = µ − P j − 1 k =1 b k . By Assumption (9), µ j − 1 ≥ µ/ 2. Because F j − 1 = F 0 + P j − 1 k =1 E k and   ∇ 2 E k   op ≤ b k on S , we hav e on U p ⊂ S the concavit y b ound  ∇ 2 F j − 1 ( g )[ X ] , X  ≤ − µ + j − 1 X k =1 b k = − µ j − 1 , so F j − 1 is µ j − 1 -strongly concav e on U p . In particular m p,j − 1 is the unique maximizer of F j − 1 on U p . W e will apply Lemma 3.3 on the smaller conv ex ball U ′ def. = B ( m p,j − 1 , r / 2) . First, U ′ is geo desically con vex b ecause r / 2 is b elo w the conv exit y radius assumed in (A). Second, U ′ ⊂ U p : indeed by (13) at level j − 1 and (9), d ( m p,j − 1 , m p, 0 ) ≤ j − 1 X k =1 2 a k µ ≤ r 2 , 9 hence for any g ∈ U ′ , d ( g , m p, 0 ) ≤ d ( g , m p,j − 1 ) + d ( m p,j − 1 , m p, 0 ) ≤ r 2 + r 2 = r , so g ∈ U p . On U ′ ⊂ S we ha ve ∥∇ E j ∥ ∞ ≤ a j and ∥∇ 2 E j ∥ op ≤ b j . T o apply Lemma 3.3 with strong concavit y constan t µ j − 1 and radius r / 2, we need b j < µ j − 1 and a j < ( µ j − 1 − b j ) r 4 . But µ j − 1 − b j = µ − j X k =1 b k ≥ µ − J X k =1 b k ≥ µ 2 , so in particular b j < µ j − 1 , and ( µ j − 1 − b j ) r 4 ≥ µr 8 . By (9), a j < µr / 8, hence the conditions of Lemma 3.3 are satisﬁed. Lemma 3.3 therefore yields a maximizer m p,j ∈ int( U ′ ) ⊂ int( U p ) of F j on U ′ , and in particular on U p . T o see that m p,j is in fact the unique maximizer of F j on all of U p , note that for an y g ∈ U p and an y unit tangent X ∈ T g SO(3),  ∇ 2 F j ( g )[ X ] , X  =  ∇ 2 F 0 ( g )[ X ] , X  + j X k =1  ∇ 2 E k ( g )[ X ] , X  ≤ − µ + j X k =1 b k . By (9), P j k =1 b k ≤ µ/ 2, so F j is ( µ/ 2)-strongly conca ve on U p and therefore has at most one maximizer on the geo desically conv ex U p . Since m p,j ∈ U ′ ⊂ U p is already a maximizer, it is the unique maximizer of F j on U p . Lemma 3.3 gives d ( m p,j , m p,j − 1 ) ≤ a j µ j − 1 − b j ≤ a j µ/ 2 = 2 a j µ , so (13) follows by summing in j . P art (ii): W e now show that global maximizers of F J all lie in S (and are th us trac ked combined with part (i)). W e hav e ∥ F J − F 0 ∥ ∞ =       J X j =1 E j       ∞ ≤ J X j =1 ∥ E j ∥ ∞ = J X j =1 δ j . Assumption (10) gives 2 ∥ F J − F 0 ∥ ∞ < Γ F 0 ( S ). Apply Lemma 3.1 with F = F 0 , e F = F J , and this S . Then every global maximizer of F J lies in S . This pro ves (ii). P art (iii): By (ii), max SO(3) F J = max g ∈ S F J ( g ) = max p max g ∈ U p F J ( g ) = max p F J ( m p,J ). Let p ⋆ b e the unique maximizer index and γ > 0 the inter-basin gap. W e ﬁrst verify that the returned Newton p oin t lies in the track ed basin. By (11) and the drift b ound (13) (with j = J ), d ( g ( J ) p , m p, 0 ) ≤ d ( g ( J ) p , m p,J ) + d ( m p,J , m p, 0 ) ≤ τ + J X k =1 2 a k µ ≤ r 2 + r 2 = r . Hence g ( J ) p ∈ U p . Now ﬁx any p . Because m p,J maximizes F J on U p and g ( J ) p ∈ U p , F J ( g ( J ) p ) ≤ F J ( m p,J ) . (14) On the other hand, for every g ∈ U p and every unit tangent v ector X ∈ T g SO(3), − M J ≤ ⟨∇ 2 F J ( g )[ X ] , X ⟩ ≤ − µ/ 2 , 10 where the upp er b ound is the ( µ/ 2)-strong concavit y established in P art (i) and the low er b ound holds b y deﬁnition of M J . Hence F J is ( µ/ 2)-strongly concav e and M J -smo oth on U p . Applying Lemma 3.2 on U p with g 0 = m p,J (noting ∇ F J ( m p,J ) = 0 since m p,J ∈ int( U p )) and using (11) giv es the lo wer b ound F J ( g ( J ) p ) ≥ F J ( m p,J ) − M J 2 τ 2 . (15) F or the winning basin p ⋆ , (15) gives F J ( g ( J ) p ⋆ ) ≥ F J ( m p ⋆ ,J ) − M J 2 τ 2 . F or any other basin p  = p ⋆ , (14) and the deﬁnition of γ give F J ( g ( J ) p ) ≤ F J ( m p,J ) ≤ F J ( m p ⋆ ,J ) − γ . Therefore F J ( g ( J ) p ⋆ ) − F J ( g ( J ) p ) ≥ γ − M J 2 τ 2 . By (12), ( M J / 2) τ 2 ≤ γ / 4, hence the diﬀerence is at least 3 γ / 4 > 0. So the algorithm’s ﬁnal argmax o ver p must select p ⋆ . Thus ˆ g = g ( J ) p ⋆ ∈ U p ⋆ and d ( ˆ g , m p ⋆ ,J ) ≤ τ by (11). P art (iv): Finally , since m p ⋆ ,J is a global maximizer of F J (from (ii) and the deﬁnition of p ⋆ ), max SO(3) F J − F J ( ˆ g ) = F J ( m p ⋆ ,J ) − F J ( ˆ g ) . Applying (15) with p = p ⋆ and ˆ g = g ( J ) p ⋆ w e get F J ( ˆ g ) ≥ F J ( m p ⋆ ,J ) − M J 2 τ 2 implying 0 ≤ F J ( m p ⋆ ,J ) − F J ( ˆ g ) ≤ M J 2 τ 2 . 3.4 Beneﬁts of marching The assumptions of Theorem 3.4 are stated in terms of the abstract quantities a j , b j , δ j . T o con- nect these with the bandwidth schedule, we use the fact that bandlimited functions on SO(3) satisfy Bernstein-t yp e deriv ative b ounds: higher deriv atives are controlled b y the bandwidth and the sup-norm of the function itself. Lemma 3.5 (Bernstein inequality on SO(3)) . L et m b e a ﬁxe d bi-invariant Riemannian metric on SO(3) , and let Π L : = span { D ℓ m,m ′ : 0 ≤ ℓ ≤ L, | m | ≤ ℓ, | m ′ | ≤ ℓ } . Then ther e exist c onstants B 1 , B 2 > 0 , dep ending only on the metric normalization, such that every F ∈ Π L satisﬁes ∥∇ F ∥ L ∞ (SO(3)) ≤ B 1 (1 + L ) ∥ F ∥ L ∞ (SO(3)) , and ∥∇ 2 F ∥ L ∞ op (SO(3)) ≤ B 2 (1 + L ) 2 ∥ F ∥ L ∞ (SO(3)) . Pr o of. Cho ose a ﬁxed orthonormal basis of left-inv arian t vector ﬁelds X 1 , X 2 , X 3 for the metric m , and set L : = − ( X 2 1 + X 2 2 + X 2 3 ) . F or a bi-inv arian t metric on SO(3), the op erator L agrees with the Laplace–Beltrami op erator up to a constant dep ending only on the metric normalization. Hence there exists κ n > 0 such that each Wigner function D ℓ m,m ′ is an eigenfunction of L with eigenv alue κ g ℓ ( ℓ + 1). Therefore Π L ⊂ E κ g L ( L +1) ( L ) , 11 where E λ ( L ) denotes the span of the eigenspaces of L with eigenv alue at most λ . By Pesenson’s Bernstein inequality on compact homogeneous manifolds, this implies max i ∥ X i F ∥ ∞ ≤ C p L ( L + 1) ∥ F ∥ ∞ , max i,j ∥ X i X j F ∥ ∞ ≤ C L ( L + 1) ∥ F ∥ ∞ , with a constant C dep ending only on the ﬁxed metric normalization. Since X 1 , X 2 , X 3 is an orthonormal frame, w e hav e p oin twise |∇ F | 2 = 3 X i =1 | X i F | 2 , so ∥∇ F ∥ ∞ ≤ √ 3 max i ∥ X i F ∥ ∞ ≲ (1 + L ) ∥ F ∥ ∞ . F or the Hessian, in the frame ( X i ) one has Hess F ( X i , X j ) = X i X j F − ( ∇ X i X j ) F . Because the metric is bi-in v ariant, ∇ X i X j = 1 2 [ X i , X j ] , and the brac ket [ X i , X j ] is a ﬁxed linear combination of X 1 , X 2 , X 3 with co eﬃcien ts dep ending only on the metric normalization. Hence | Hess F ( X i , X j ) | ≲ max p,q | X p X q F | + max p | X p F | . Since the op erator norm of a bilinear form in a ﬁxed orthonormal basis is controlled by the maxim um of its matrix entries, it follows that ∥∇ 2 F ∥ L ∞ op ≲ max p,q ∥ X p X q F ∥ ∞ + max p ∥ X p F ∥ ∞ ≲ (1 + L ) 2 ∥ F ∥ ∞ . This prov es the claim. A natural question is whether intermediate bandwidths are b eneﬁcial and if so, by how muc h, compared with making a single big jump from L 0 to L J . Indeed, frequency marching controls the p erturbation induced b y higher angular frequencies at the scale of each incremen t rather than the ﬁnal bandwidth. A single jump from L 0 to L J requires to b ound the full sp ectral tail b y the curv ature of the ﬁnal high-bandwidth ob jectiv e; as we show, marching distributes this control ov er smaller increments. Since each increment E j = F j − F j − 1 has bandwidth at most L j , Lemma 3.5 yields a j ≤ B 1 L j δ j , b j ≤ B 2 L 2 j δ j , (16) where we recall that δ j = ∥ E j ∥ ∞ . This allows to write the assumptions in terms of the quantities ( L j , δ j ) J j =1 . W e let ≲ absorb the Bernstein constan ts B 1 , B 2 from Lemma 3.5. Supp ose w e skip all in termediate lev els and apply the basin p ersistence argument directly to the single p erturbation E tot = F J − F 0 . Since E tot has bandwidth L J w e hav e by Lemma 3.5 ∥∇ 2 E tot ∥ ∞ ≲ L 2 J ∥ E tot ∥ ∞ . Preserving strong concavit y of the track ed basins therefore is guaran teed by ∥ E tot ∥ ∞ ≲ µ L 2 J . (17) By con trast, when marching we need to apply Lemma 3.5 to eac h increment E j separately . Us- ing (16), the condition P j b j ≤ µ/ 2 from Assumption (B) b ecomes B 2 J X j =1 L 2 j δ j ≤ µ 2 . 12 One-shot: L 0 → L J L J true wrong max jump L 0 g ⋆ 0 Marc hing: L 0 → L 1 → · · · → L J L 0 L 1 L 2 L 3 L J ✓ land in basin climb to basin p eak Figure 2: One-shot v ersus frequency marching on a sc hematic 1D correlation landscape. Left: jumping directly from L 0 to L J , Newton’s metho d c an conv erge to a spurious lo cal maximum that emerged at high bandwidth. Righ t: under the conditions of Theorem 3.4, marching through in termediate bandwidths keeps the track ed maximizer in the correct basin at each level; Newton reﬁnes within the basin b efore pro ceeding to the next bandwidth. W rite δ tot = P j δ j and deﬁne the energy-weigh ted eﬀective bandwidth L 2 2 def. = P J j =1 L 2 j δ j δ tot , (18) so that the marc hing condition reads δ tot ≲ µ/L 2 2 . Note that δ tot ≥ ∥ E tot ∥ ∞ b y the triangle inequalit y , so the marching analysis works with a potentially larger n umerator. The gain comes from the denom- inator: since L 2 ≤ L J , with equalit y only when all energy is concen trated at the ﬁnal level, marching comp ensates by replacing the worst-case L 2 J w eight with the eﬀective weigh t L 2 2 . The denominator impro ves from L 2 J to L 2 2 , but this can b e oﬀset by the larger numerator δ tot = P j δ j . Hence march- ing is most adv antageous when the sp ectral increments are concentrated at low bandwidths (which is the typical regime for smo oth signals or noise-dominated high-frequency conten t) and the cancel- lation factor δ tot / || E tot || ∞ is not too large. An analogous comparison holds for the gradien t-based conditions. The one-shot approac h is ensured by ∥ E tot ∥ ∞ ≲ µr /L J , whereas marching imp oses tw o p er-step conditions: L j δ j ≲ µr (to apply basin p ersistence at each level) and P j L j δ j ≲ µr (to con trol the cumulativ e drift of the track ed maxima). Both are easier to satisfy when the increments δ j deca y with j . A to y mo del. T o quantify the gain, consider an example of uniform bandwidth steps L j = L 0 + ( j − 1)∆ with geometrically decaying increments δ j = δ q j − 1 , q ∈ (0 , 1). A direct computation in the J → ∞ limit gives L 2 2 = L 2 0 + 2 L 0 ∆ q 1 − q + ∆ 2 q (1 + q ) (1 − q ) 2 . (19) The key observ ation is that this expression is indep endent of L J : no matter how large the ﬁnal bandwidth, the eﬀectiv e penalty dep ends only on L 0 and the deca y rate q . When q ≪ 1 (rapid sp ectral decay), L 2 ≈ L 0 , and the (a priori) marc hing requiremen t δ tot ≲ µ/L 2 0 is generically muc h milder than the one-shot requirement (17). Note that the ﬁnal sub optimalit y b ound in Theorem 3.4 still scales with L 2 J , since it reﬂects the curv ature of the ﬁnal ob jective rather than the cumulativ e p erturbation. 13 (a) SNR = 20dB. (b) SNR = 0dB. (c) SNR = − 20dB. Figure 3: Slice of a rib osome volume corrupted by Gaussian noise at the indicated SNR levels. 4 Exp erimen ts W e ﬁrst v alidate the proposed method on syn thetic data, demonstrating substantial gains in both sp eed and accuracy ov er exhaustive rotational search. W e then sho w that our approach can b e used as a drop-in replacemen t in the subtomogram a veraging (ST A) pip eline of RELION, accelerating rotational alignmen t by an order of magnitude while reﬁning a rib osome to 4 ˚ A, reaching the lo cal Nyquist limit of 3 . 8 ˚ A. The exp erimen ts were conducted either on one A100-80GB GPU (Section 4.1) or on four R TX4090 GPUs (Section 4.2). They can be repro duced using co de av ailable at https://gith ub.com/swing- researc h/Matcha 4.1 Rotational searc h: syn thetic examples Let h ∈ R N × N × N denote a rib osome density at 13 ˚ A resolution (EMDB:EMD-3228), obtained from the EMPIAR-10045 dataset 1 . W e synthetically generate measurements f = g ⋆ ◦ h + η , where g ⋆ ∈ SO(3) is the ground truth rotation (assumed unknown in the search) and η is additive white Gaussian noise. W e consider three signal-to-noise ratios: 20 dB, 0 dB, and − 20 dB. Figure 3 sho ws a represen tative volume slice at each level; at 0 dB the signal and noise hav e equal energy , and in cryogenic electron microscopy it is common to work at or b elo w − 20 dB. W e rep ort the results at 0 dB in the main text and defer the remaining SNR levels to App endix B.2. A t eac h SNR, we generate 1000 randomly rotated volumes of size N = 200 and compare the p erformance of Algorithm 1 with SOFFT [4, 18, 7]. SOFFT ev aluates the cross-correlation on a discrete grid of rotations whose resolution is con trolled by an ov ersampling factor K : increasing K reﬁnes the grid at the cost of ev aluating ( K − 1) 3 additional p oin ts. All SOFFT exp erimen ts w ere run with a batch size of one, as larger batches exceed GPU memory at higher bandwidths; even so, SOFFT frequently pro duces out-of-memory errors at large ov ersampling factors (see Figure 4b). F or Matc ha, we select 10 candidate rotations as the highest lo cal maxima of the cross-correlation computed by SOFFT at a coarse bandwidth L 0 = 30 with ov ersampling factor 2. Candidates are then reﬁned across bandwidth L j ∈ { 30 , 40 , 60 , L max } , skipping an y L j that exceed L max : at each L j , a single Newton step is applied to ev ery candidate, and reﬁned estimates are propagated to the next band. After optimization at L max , the candidate with the highest correlation score is retained. Figure 4a rep orts the w all-clo c k time and 90th-p ercen tile alignment error for a ﬁxed bandwidth L max = 40, whic h yields the fastest inference but not the best resolution—especially for SOFFT, whose accuracy is fundamen tally limited by its grid spacing. A t the highest o versampling factor that do es not cause out-of-memory errors ( K = 18), SOFFT achiev es 0 . 13 ◦ accuracy but requires ov er a hundred seconds. F or K = 8 and ab ov e, wall-clock times range from ten seconds to several minutes. By con trast, Matcha ac hieves 0 . 03 ◦ accuracy–o ver four times more precise–within a few seconds, despite starting from a coarse initial grid. 1 https://www.ebi.ac.uk/emdb/EMD-3228 14 This gain stems from a structural adv antage: unlik e SOFFT, whose accuracy is fundamentally limited by the grid spacing at a giv en bandwidth, Matcha uses the grid only to identify candidate basins and then reﬁnes con tinuously to a precision that is indep enden t of the initial grid resolution. Matc ha also uses memory more eﬃciently , whic h enables larger batc h sizes when processing m ultiple v olumes. Under identical hardware constraints, SOFFT is limited to a batch size of 10 for K ≤ 8, a batc h size of 4 for K = 10, and a batc h size of 1 for K = 18, whereas Matcha supp orts batch sizes of up to 100. The maximum bandwidth L max is a key hyperparameter of b oth metho ds: small v alues yield fast computation but introduce bias that limits accuracy , while larger v alues are slow er but more accurate. Figure 4b conﬁrms that higher L max consisten tly improv es estimation accuracy . The eﬀect is mo derate for Matcha and more pronounced for SOFFT at all o versampling factors. Bey ond L max = 100, im- pro vemen ts for Matcha are marginal, whereas SOFFT b ecomes infeasible due to memory constraints. W e also compare the Newton reﬁnement used in Matcha with a ﬁrst-order gradient metho d. Since b oth methods require ev aluating the same band-limited ob jective, the practical diﬀerence is the num b er of reﬁnemen t steps needed to reach high precision. Figure 4c shows that Newton reﬁnement reaches close to optimal accuracy after v ery few steps, whereas the ﬁrst-order metho d requires muc h more steps for the desired accuracy . This supp orts our use of one Newton step p er bandwidth level in the syn thetic exp erimen ts. 4.2 Nyquist resolution in cryo-ET using ST A W e in tegrate Matcha into the subtomogram av eraging pip eline of RELION5 [19] as a drop-in re- placemen t for the rotational alignment p erformed during 3D auto-r eﬁne . All other pro cessing steps (classiﬁcation, CTF-reﬁnemen t, Ba yesian p olishing) remain unchanged and are executed using RE- LION’s own implemen tation throughout; since these steps are shared b et ween the baseline and our metho d, w e exclude them from the run time comparison and report only the cumulativ e wall-clock time sp en t on alignment. W e follow the ST A tutorial 2 for the Chlamy dataset (EMPIAR-11830), originally published b y [20]. The dataset contains 33 bin4 tomograms, from which 20 , 990 particles were extracted using p ytom-match-pic k [21]. The RELION baseline reﬁnes progressively through three binning levels (bin4, bin2, bin1), reac hing 3 . 8 ˚ A lo cal resolution at bin1—the Nyquist limit—under gold-standard FSC at a threshold of 0 . 143. The cumulativ e alignment time across all three stages is approximately 3 . 5 hours on four R TX 4090 GPUs. W e replace 3D auto-r eﬁne at the bin2 and bin1 stages, retaining RELION’s bin4 alignmen t as initialization. Matc ha estimates b oth rotations and translations using the same half-set split as RE- LION, with translations reco vered via FFT-based correlation (details in Appendix C). At eac h binning lev el, we p erform T = 3 outer rotation–translation alternations. Within each alternation, we retain N C = 10 rotation candidates from a lo cal searc h at angular cutoﬀ L 0 = 30 with o versampling factor 2, and reﬁne with M = 5 Newton steps—more conserv ative than the 1–2 steps suﬃcient in the synthetic setting (Section 4.1), to ensure conv ergence under the noisier conditions of experimental data—through bandwidths { 30 , 40 , L max } , where L max = 45 (bin2) and L max = 60 (bin1). Shifts are estimated via FFT-based correlation with subpixel reﬁnemen t [22]. Our bin1 reconstruction attains 3 . 8 ˚ A lo cal resolution (Nyquist), identical to the RELION baseline, and 4 . 0 ˚ A global resolution under the same FSC criterion. The cumulativ e alignment time drops from appro ximately 3 . 5 hours (RELION) to roughly 20 minutes (Matcha), a speedup of o ver 10 × (Figure 5). Since the non-alignmen t pro cessing steps add a ﬁxed ov erhead of approximately 90 minutes to b oth pip elines, the end-to-end sp eedup is smaller but alignment is the dominan t cost in standard ST A w orkﬂows. 2 https://tomoguide.gith ub.io 15 1 10 100 0 1 2 K=1 2 . 34 ◦ K=2 1 . 16 ◦ K=3 0 . 78 ◦ K=4 0 . 59 ◦ K=8 0 . 3 ◦ K=10 0 . 24 ◦ K=18 0 . 13 ◦ Matcha B10 0 . 03 ◦ Matcha B100 0 . 03 ◦ Execution time (sec) 90th p ercen tile error (degrees) (a) Accuracy vs. execution time ( L max = 40). 40 60 80 100 0 1 2 3 K=1 3 . 12 ◦ 0 . 89 ◦ K=3 1 . 05 ◦ 0 . 31 ◦ K=10 0 . 31 ◦ 0 . 14 ◦ K=18 Matcha 0 . 08 ◦ 0 . 03 ◦ Bandwidth L max 90th p ercen tile error (degrees) (b) Eﬀect of L max for v arying oversampling K . 1 2 3 4 5 10 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 Number of steps Newton’s metho d Gradient ascent (c) Newton’s metho d vs. gradient ascent ( L max = 40). Figure 4: Hyp erparameter study for Matcha and SOFFT at SNR = 0 dB. In (a), dot size enco des alignment error; Matc ha ac hieves 0 . 03 ◦ accuracy in seconds, whereas SOFFT at its highest feasible o versampling ( K = 18) requires ov er a hundred seconds for 0 . 13 ◦ . In (b), increasing L max impro ves b oth metho ds, but the eﬀect is more pronounced for SOFFT; b ey ond L max = 100, SOFFT b ecomes memory-prohibitive at large K . In (c), a single Newton step matches the accuracy that gradient ascent reaches only after several iterations, supp orting the use of one Newton step p er band in (a) and (b). 16 8.1Å 4.5Å 3.8Å Ours 5min ⏲ : 60m in 10min 22min ⏲ ⏲ ⏲ B as eline 5m in ⏲ 75m in 210m in ⏲ ⏲ ⏲ ⏲ ⏲ ⏲ Figure 5: Subtomogram av eraging on the Chlamy dataset (EMPIAR-11830) using RELION5. Matcha replaces 3D auto-r eﬁne at the bin2 and bin1 stages, matching the baseline resolution (3 . 8 ˚ A at Nyquist) while reducing cum ulative alignmen t time by ov er 10 × . Intermediate pro cessing steps (CTF-reﬁnement, classiﬁcation, tem- plate updates) are shared b et ween b oth pip elines and excluded from the reported timings. 5 Conclusion W e presented Matc ha, a frequency-marc hing algorithm for fast rotational alignmen t that replaces exhaustiv e high-resolution search with a combination of coarse SO(3)-FFT initialization and Newton reﬁnemen t across progressively ﬁner angular bandwidths. The key algorithmic insight is that band- limited correlation landscap es are smo oth enough for second-order optimization to conv erge in very few steps, while the low-dimensional structure of SO(3) keeps the p er-step cost negligible. The k ey theoretical insight is that frequency marching—reﬁning through intermediate bandwidths rather than jumping directly to the target resolution—relaxes the sp ectral-deca y conditions required for prov able basin tracking from a worst-case L 2 J scaling to an eﬀective L 2 2 scaling that is indep enden t of the ﬁnal bandwidth. On synthetic b enc hmarks, Matcha achiev es sub-degree angular accuracy several orders of magni- tude faster than exhaustive searc h. In tegrated in to the RELION5 subtomogram av eraging pip eline, it reduces alignmen t run time by ov er 10 × while matching Nyquist-limited reconstruction qualit y , demon- strating that the approach is practical for large-scale cry o-ET datasets. Sev eral directions remain op en. First, our con vergence analysis assumes a deterministic setting with known signal structure; it would b e v aluable to establish probabilistic guarantees—for instance, sho wing that the spectral-decay and gap conditions hold with high probabilit y o ver a suitable ensem ble of noisy particles. Second, the curren t alternating rotation–translation scheme inherits the limitations of blo c k co ordinate ascent. A joint con tinuous optimization ov er the full six-dimensional p ose space, exploiting the semi-direct pro duct structure of SE(3), could improv e b oth accuracy and con vergence sp eed. Third, the image formation mo del in cryo-ET includes degradations b ey ond additiv e noise— notably the contrast transfer function (CTF) and the missing w edge—which are curren tly handled b et ween alignment iterations rather than within them. Incorporating these eﬀects directly into the band-limited correlation ob jective could p oten tially reduce the num b er of outer iterations needed and further accelerate the pip eline. Finally , in practice the reference template is itself unknown and m ust b e iteratively estimated from the data, and the particle p opulation ma y contain multiple conforma- tional states. The sp eed of Matc ha mak es it feasible to ev aluate alignment against multiple candidate templates p er iteration, opening a path to ward faster heterogeneity analysis in subtomogram av eraging. 17 A Ball harmonics, rotational correlation, and radial sparsit y The ball harmonics provide an orthonormal, rotation-steerable basis for the function space L 2 ( B 3 ) that is naturally ordered by Laplacian eigen v alues. This structure mak es the representation particularly con venien t for rotational alignment via cross-correlation, and it induces a characteristic r adial sp arsity pattern under the eigenv alue-based truncation used as describ ed in [8]. In this part, we give more precise deﬁnition of the ball harmonics, their prop erties and we pro vide references for the main results. A.1 Ball harmonics expansion and eigen v alue-based truncation W e let B 3 denote the unit ball in R 3 and let u : B 3 → R b e supp orted on B 3 . It admits an expansion of the form u ( x ) = X ( k,ℓ,m ) ∈I Λ ˆ u k,ℓ,m ψ k,ℓ,m ( x ) , (20) where the basis { ψ k,ℓ,m } k,l,m is known as the ball harmonics basis. These basis functions are separable in to radial (indexed by k ) and angular (indexed by ( l, m )) factors and form an orthonormal basis of L 2 ( B 3 ), see [8] for con ven tions. The basis functions { ψ k,ℓ,m } k,ml are Diric hlet Laplacian eigenfunctions on B 3 with eigenv alues/frequencies λ ℓk = j ℓ,k , the k th p ositiv e ro ot of the spherical Bessel function j ℓ , and satisfy the eigenv alue equation − ∆ ψ k,ℓ,m = λ 2 ℓk ψ k,ℓ,m in B 3 , with ψ k,ℓ,m = 0 on ∂ B 3 . In practice, we k eep only frequencies λ lk that are b elo w a certain threshold (frequency budget) Λ > 0, and retain the index set I Λ := { ( k , ℓ, m ) : λ ℓk ≤ Λ , ℓ ∈ Z ≥ 0 , m ∈ {− ℓ, . . . , ℓ } , k ∈ Z > 0 } . F or a giv en angular degree ℓ , deﬁne the retained radial set K ℓ := { k ∈ Z > 0 : λ ℓk ≤ Λ } . The function ( ℓ, k ) 7→ λ ℓk is increasing in b oth ℓ and k , so the cardinality | K ℓ | decrease monotonically with ℓ . Lo w-rank approximation can b e emplo yed to reduce the num b er of co eﬃcien ts in the set I Λ , see App endix B.1. A.2 Computational complexit y A key practical prop ert y is that the truncated expansion in Equation (20) can b e computed in nearly- linear time (up to p olylogarithmic factors) in the num b er of vo xels. Lemma A.1 (Ball harmonics transform complexity) . Given a volume sample d on an N × N × N Cartesian grid, the b al l harmonics de c omp osition up to a b andwidth Λ = O ( N ) c an b e c ompute d with c omplexity O ( N 3 log 2 N ) while achieving numeric al ac cur acy ε = O (1 / N ) . Pr o of. This follows from [8, Theorem 3.1]. Cho osing ε = 1 / N giv es | log ε | = Θ(log N ), so the stated complexit y b ound reduces to O ( N 3 (log N ) 2 ) for Λ = O ( N ). A.3 Wigner– D matrices An orthogonal basis for L 2 (SO(3)) is given by the Wigner D -matrices. As a result, rotation of ball harmonic functions can b e expressed as a linear com bination of other ball harmonic functions and Wigner D -matrices with angular indices up to the same order. This prop ert y will b ecome handy to eﬃcien tly ev aluate the cross-correlation for a large n umber of rotations, as well as for calculating their ﬁrst and second order deriv atives. More precisely , for any g ∈ SO(3), we hav e ψ k,ℓ,m ( g − 1 x ) = ℓ X m ′ = − ℓ D ℓ mm ′ ( g ) ψ k,ℓ,m ′ ( x ) , ∀ x ∈ B 3 . 18 As deﬁned in Equation 3, an y rotation g ∈ SO(3) can b e parameterized by ZYZ Euler angles ( α, β , γ ) ∈ [0 , 2 π ) × [0 , π ] × [0 , 2 π ). F or eac h angular degree ℓ ∈ Z ≥ 0 , the Wigner D -matrix D ℓ ( g ) ∈ C (2 ℓ +1) × (2 ℓ +1) , has entries indexed by order m, m ′ ∈ {− ℓ, . . . , ℓ } and is deﬁned by D ℓ mm ′ ( α, β , γ ) = e − imα d ℓ mm ′ ( β ) e − im ′ γ , (21) where d ℓ mm ′ are the Wigner d -functions . The Wigner d -functions dep end only on the p olar angle β and admit explicit expressions in terms of ﬁnite sums [23] and can be ev aluated eﬃcien tly . The collection { D ℓ mm ′ } ℓ,m,m ′ forms an orthogonal basis of L 2 (SO(3)) with resp ect to the Haar measure dg , normalized so that Z SO(3) D ℓ mm ′ ( g ) D ℓ ′ nn ′ ( g ) dg = 1 8 π 2 δ ℓℓ ′ δ mn δ m ′ n ′ . (22) The Wigner D -matrices admit clos ed-form expressions not only for their v alues but also for their deriv atives with respect to the rotation parameters. In particular, deriv atives with resp ect to the Euler angles ( α, β , γ ) can b e written analytically in terms of the same Wigner D -matrices, see [23]. Diﬀeren tiation with resp ect to α and γ acts diagonally , ∂ α D ℓ mm ′ = − im D ℓ mm ′ , ∂ γ D ℓ mm ′ = − im ′ D ℓ mm ′ , while deriv atives with resp ect to β can be expressed as ﬁnite linear combinations of d ℓ mm ′ (cos β ) with shifted indices and known co eﬃcien ts. As a result, gradients and higher-order deriv atives of Wigner D -matrices are av ailable analytically and can b e ev aluated eﬃciently . A.4 Rotational cross-correlation in the ball harmonics basis Rotations act only on the angular index via Wigner– D matrices. Consequently , the rotational cross- correlation of f , h ∈ L 2 ( B 3 ), C( g ) = Z B 3 f ( x ) h ( R − 1 g x ) dx, g ∈ SO(3) , admits a natural Wigner– D expansion. In a band-limited setting with angular cutoﬀ L , we write C L ( g ) = L X ℓ =0 ℓ X m = − ℓ ℓ X m ′ = − ℓ σ ℓmm ′ D ℓ mm ′ ( g ) , where D ℓ mm ′ ( g ) are the Wigner– D matrix entries deﬁned in Equation (21). The co eﬃcien ts are deter- mined directly by the ball harmonic co eﬃcien ts of f and h : σ ℓmm ′ = X k ∈ K ℓ ˆ f k,ℓ,m ˆ h k,ℓ,m ′ . W e will see that these co eﬃcien ts admit a low-dimension decomp osition. A.5 Lo w-rank structure induced by radial sparsit y F or each ﬁxed ℓ , deﬁne the (2 ℓ + 1) × (2 ℓ + 1) matrix A ℓ =  σ ℓmm ′  ℓ m,m ′ = − ℓ . Using (A.4), A ℓ admits the factorization A ℓ = X k ∈ K ℓ a ℓk b ∗ ℓk , ( a ℓk ) m = ˆ f k,ℓ,m , ( b ℓk ) m = ˆ h k,ℓ,m , (23) 19 and therefore rank( A ℓ ) ≤ | K ℓ | . (24) Since eigenv alue-based truncation enforces | K ℓ | to decrease with ℓ , the co eﬃcien ts blo c ks A l of the correlation b ecome progressively low-rank at higher angular degrees. In practice, the eﬀectiv e rank is often ev en smaller than | K l | due to additional dep endencies among the vectors { a ℓk } k ∈ K ℓ and { b ℓk } k ∈ K ℓ . In particular, rank( A l ) do es not dep end on the rotation of h : Let h g ∗ = h ( R − 1 g ∗ x ) b e the volume h rotated by g ∗ ∈ SO(3) with ( b ℓk ) ( g ∗ ) m and A ( g ∗ ) l the corresp onding factorization of the ball harmonics co eﬃcien ts of h g ∗ . Since A ( g ∗ ) l = X k a lk b ∗ lk ( g ∗ ) = X k a lk b ∗ lk D l ( g ∗ ) † = A l D l ( g ∗ ) † . and b ecause the matrices D l ( g ) are unitary , and hence inv ertible, for an y g ∈ SO(3), w e hav e rank( A ( g ∗ ) l ) = rank( A l D l ( g ∗ ) † ) = rank( A l ) ≤ | K l | . B Additional Exp erimen ts B.1 Empirical eﬀectiv e rank of Wigner blo c ks v ersus the theoretical b ound In App endix A.5 we remark that, for eac h angular degree l , the Wigner blo c k A l ∈ C (2 l +1) × (2 l +1) of the rotational correlation admits a factorization such that rank( A l ) ≤ | K l | , where K l is the set of radial indices retained at degree l b y the eigenv alue cutoﬀ Λ. This provides a the or etic al (worst-case) rank upper b ound determined purely by the truncation scheme. In practice, ho wev er, w e observe substan tially low er rank b eha vior. F or selected degrees l ∈ { 20 , 30 , 50 , 80 } , see Figure 6a, the singular v alues exhibit a pronounced sp ectral gap: only a small num b er of mo des carry signiﬁcan t mass, while the remaining singular v alues rapidly decay to a near-zero ﬂo or. Motiv ated b y this b eha vior, we deﬁne the eﬀective rank as r eﬀ ( l ) def. = min { k : s k +1 ( A l ) s k ( A l ) < ρ } , where s 1 ( A l ) ≥ s 2 ( A l ) ≥ · · · are the singular v alues of A l and ρ is a ﬁxed numerical tolerance. Figure 6b shows r eﬀ ( l ) alongside the upp er-bound | K l | across degrees for ρ = 0 . 2. Across degrees, r eﬀ ( l ) sta ys w ell below | K l | with high- l blo c ks approac hing rank one. This indicates that the correlation is eﬀectiv ely supp orted on a low-dimensional subspace relative to the ambien t (2 l + 1) × (2 l + 1) blo c k size. B.2 Ablation study for other noise regimes W e follow the exp erimen tal setting presented in Section 4.1 and display the results for the extremely noisy setting, SNR=-20dB, and the easier setting, SNR=20dB. A key h yp erparameter of Matc ha and SOFFT is the maximum bandwidth L max used to ev aluate the truncated correlation in Equation (5). This parameter controls the highest angular frequencies included in the alignment score and trades high-resolution information for computational time. With Matcha, the exp erimen ts indicate that accurate alignment can already b e achiev ed at rela- tiv ely low bandwidth when the signal-to-noise ratio is mo derate or high. As shown in Figures 4b, the alignmen t accuracy remains high across the entire range of tested L max , suggesting that low-frequency comp onen ts alone are suﬃcien t to reliably localize the correct rotation in these regimes. A t very low SNR, e.g. SNR = − 20 dB, Figure 7b, the ov erall accuracy of b oth metho ds degrades, reﬂecting the increased diﬃcult y of the alignmen t problem. Nevertheless, b oth metho ds contin ue to b eneﬁt from increasing the bandwidth. F or Matcha, higher v alues of L max lead to a gradual impro vemen t in accuracy even in this strongly noisy setting, although the gains b ecome progressively smaller as additional high-frequency comp onen ts are increasingly dominated by noise. A similar trend is observed for SOFFT. 20 0 50 100 150 10 − 12 10 − 7 10 − 2 10 3 Singular v alue index Singular v alue amplitude (log) l = 20 l = 30 l = 50 l = 80 (a) Singular values of A l for selected degrees. 0 50 0 10 20 30 Degree l Rank Eﬀective rank Upper-b ound | K l | (b) Eﬀective rank vs. theoretical upp er b ound | K l | . Figure 6: Spectral prop erties of the Wigner blo c ks A l . In (a), the singular v alues exhibit a sharp sp ectral gap at each degree, with only a few mo des carrying signiﬁcant energy . In (b), the eﬀective rank (deﬁned by the ratio criterion s k +1 /s k < ρ with ρ = 0 . 2) stays well b elo w the theoretical b ound | K l | across all degrees, with high- l blo c ks approaching rank one. This low-rank structure can b e exploited to accelerate the ev aluation of C L . Ov erall, these results sho w that, as exp ected, increasing the bandwidth is b eneﬁcial for b oth meth- o ds across all noise levels. At the same time, Matcha exhibits strong robustness to bandwidth trun- cation at mo derate and high SNR, ac hieving high accuracy even at lo w L max , while still being able to exploit additional frequency information when av ailable. Figure 7a and Figure 7b present the execution times for Matc ha and SOFFT with v arying o ver- sampling factor K at SNR 20dB and -20dB resp ectiv ely . In general, execution time for b oth Matcha and SOFFT do es not dep end on the noise level. 40 60 80 100 0 1 2 3 K=1 3 . 12 ◦ 0 . 89 ◦ K=3 1 . 05 ◦ 0 . 3 ◦ K=10 0 . 3 ◦ 0 . 13 ◦ K=18 Matcha 0 . 07 ◦ 0 ◦ Bandwidth L max 90th p ercen tile error (degrees) (a) SNR= 20dB. 40 60 80 100 0 1 2 3 K=1 3 . 18 ◦ 0 . 93 ◦ K=3 1 . 11 ◦ 0 . 43 ◦ K=10 0 . 47 ◦ 0 . 36 ◦ K=18 Matcha 0 . 34 ◦ 0 . 32 ◦ Bandwidth L max 90th p ercen tile error (degrees) (b) SNR= − 20dB. Figure 7: Alignment accuracy of Matcha and SOFFT with v arying ov ersample factor K for diﬀerent L max for diﬀeren t SNR. C Alternating searc h o ver rotations and translations In subtomogram alignment, the relativ e p ose b et ween t wo volumes is t ypically describ ed by a rotation and a translation. Given t wo v olumetric densities f , h ∈ L 2 ( B 3 ), we consider the 6-dimensional alignment problem ( ˆ τ , ˆ g ) ∈ argmax τ ∈ R 3 , g ∈ SO(3) C( τ , g ) , C( τ , g ) def. =  f , g ◦ S τ h  , (25) 21 1 10 100 0 1 2 K=1 2 . 34 ◦ K=2 1 . 16 ◦ K=3 0 . 77 ◦ K=4 0 . 58 ◦ K=8 0 . 29 ◦ K=10 0 . 24 ◦ K=18 0 . 12 ◦ Matcha B10 0 . 01 ◦ Matcha B100 0 . 01 ◦ Execution time (sec) 90th p ercen tile error (degrees) (a) SNR: 20dB 1 10 100 0 1 2 K=1 2 . 37 ◦ K=2 1 . 22 ◦ K=3 0 . 87 ◦ K=4 0 . 7 ◦ K=8 0 . 46 ◦ K=10 0 . 41 ◦ K=18 0 . 38 ◦ Matcha B10 0 . 33 ◦ Matcha B100 0 . 33 ◦ Execution time (sec) 90th p ercen tile error (degrees) (b) SNR: − 20dB Figure 8: Alignment accuracy of Matcha and SOFFT versus execution time for L max = 40 and v arying SNR. Large dot sizes corresp ond to larger alignment error. Matc ha is one order of magnitude more precise and sev eral orders of magnitude faster than exhaustiv e search SOFFT. where ( S τ h )( x ) = h ( x − τ ) denotes the translation op erator. It is unrealistic to exhaustively discretize the parametrization of the cross-correlation as deﬁned in Problem (25) as it lives in six dimensions. A common approach is to alternate b et ween optimization o ver shifts and ov er the rotations, e.g. [15]. F or ﬁxed shift τ , optimization ov er g can b e p erformed eﬃcien tly using Algorithm 1. F or ﬁxed rotation g , the optimal translation can b e estimated eﬃciently on a ﬁne grid b y ev aluating the cross-correlation with 3D FFT. More precisely , w e start from zero shift and estimate the rotation using Algorithm 1. Then, giv en the estimated rotation, we estimate the translation of the rotated volume via a 3D FFT-based correlation. This alternating pro cess is rep eated for a small num b er of iterations, e.g. T = 3 in the exp erimen t of Section 4.2. Remarks. (i) Eac h translation up date costs O ( N 3 log N ) on an N × N × N grid (or less if restricted to a b ounded search window), and is typically inexp ensiv e compared to the rotation reﬁnement at higher angular cutoﬀs. (ii) The outer lo op is a blo c k co ordinate ascent metho d and is not guaran teed to ﬁnd the global maximizer of the 6-dimensional ob jective; nevertheless, empirically , it has b een emplo yed successfully and is known to b e robust if relatively accurate initial translation is provided. (iii) If desired, multi-resolution translation reﬁnement can b e employ ed to achiev e sub-pixel accuracy , for example using the eﬃcient upsampled phase-correlation scheme of [22]. (iv) T ranslational search can b e restricted to a lo cal neighborho o d if an accurate initial translation is av ailable; as is the case for the hierarchical binning reﬁnemen t data typically emplo yed in ST A, see Exp erimen t 4.2. References [1] George L. T urin. An introduction to matched ﬁlters. IRE T r ansactions on Information The ory , 6(3):311–329, 1960. [2] Steven M. Ka y . F undamentals of Statistic al Signal Pr o c essing, V olume II: Dete ction The ory . Pren tice Hall, 1998. [3] V alentin Debarnot and Pierre W eiss. Blind in verse problems with isolated spikes. Information and Infer enc e: A Journal of the IMA , 12(1):26–71, 2023. [4] Peter J. Kostelec and Daniel N. Ro c kmore. FFTs on the Rotation Group. Journal of F ourier A nalysis and Applic ations , 14(2):145–179, April 2008. [5] Matthew A. Price and Jason D. McEwen. Diﬀerentiable and accelerated spherical harmonic and Wigner transforms. Journal of Computational Physics , 510:113109, August 2024. 22 [6] Julio A Ko v acs and Willy W riggers. F ast rotational matc hing. Biolo gic al Crystal lo gr aphy , 58(8):1282–1286, 2002. [7] Jason D McEwen, Martin B ¨ uttner, Boris Leistedt, Hirany a V Peiris, and Yv es Wiaux. A nov el sampling theorem on the rotation group. IEEE Signal Pr o c essing L etters , 22(12):2425–2429, 2015. [8] Jo e Kileel, Nicholas F Marshall, Oscar Mic kelin, and Amit Singer. F ast expansion into harmonics on the ball. SIAM Journal on Scientiﬁc Computing , 47(2):A1117–A1144, 2025. [9] Y u-hsuan Shih, Garrett W right, Joakim And ´ en, Johannes Blaschk e, and Alex H Barnett. cuﬁnuﬀt: a load-balanced gpu library for general-purp ose non uniform ﬀts. In 2021 IEEE International Par al lel and Distribute d Pr o c essing Symp osium Workshops (IPDPSW) , pages 688–697. IEEE, 2021. [10] Julio A Kov acs, P ablo Chac´ on, Y ao Cong, Essam Metw ally , and Willy W riggers. F ast rotational matc hing of rigid bo dies b y fast fourier transform acceleration of ﬁv e degrees of freedom. Biolo gic al Crystal lo gr aphy , 59(8):1371–1376, 2003. [11] Leandro F arias Estrozi and Jorge Nav aza. F ast pro jection matching for cryo-electron microscop y image reconstruction. Journal of Structur al Biolo gy , 162(2):324–334, 2008. [12] A. Bartesaghi, P . Sprechmann, J. Liu, G. Randall, G. Sapiro, and S. Subramaniam. Classiﬁcation and 3D av eraging with missing wedge correction in biological electron tomography . Journal of Structur al Biolo gy , 162(3):436–450, June 2008. [13] Maxim Shatsky , P ablo Arb elaez, Rob ert M Glaeser, and Steven E Brenner. Optimal and fast rotational alignment of v olumes with missing data in fourier space. Journal of Structur al Biolo gy , 184(2):345–347, 2013. [14] Min Xu, Martin Beck, and F rank Alb er. High-throughput subtomogram alignment and clas- siﬁcation by fourier space constrained fast v olumetric matc hing. Journal of structur al biolo gy , 178(2):152–164, 2012. [15] Y uxiang Chen, Stefan Pfeﬀer, Thomas Hrab e, Jan Michael Sch uller, and F riedrich F¨ orster. F ast and accurate reference-free alignment of subtomograms. Journal of Structur al Biolo gy , 182(3):235– 245, June 2013. [16] Mona Zehni, Laur` ene Donati, Emman uel Soubies, Zhizhen Zhao, and Michael Unser. Join t an- gular reﬁnemen t and reconstruction for single-particle cryo-em. IEEE T r ansactions on Image Pr o c essing , 29:6151–6163, 2020. [17] Alex Barnett, Leslie Greengard, Andras P ataki, and Marina Spiv ak. Rapid solution of the cryo-em reconstruction problem by frequency marching. SIAM Journal on Imaging Scienc es , 10(3):1170– 1195, 2017. [18] Daniel Potts, J ¨ urgen Prestin, and An tje V ollrath. A fast algorithm for nonequispaced F ourier transforms on the rotation group. Numeric al Algorithms , 52(3):355–384, No vem b er 2009. [19] Alister Burt, Bogdan T oader, Rangana W arshamanage, Andriko von K ¨ ugelgen, Euan Pyle, Jasenk o Ziv anov, Dari Kimanius, T anmay AM Bharat, and Sjors HW Scheres. An image pro cess- ing pip eline for electron cryo-tomograph y in relion-5. FEBS op en bio , 14(11):1788–1804, 2024. [20] Ron Kelley , Sagar Kha vnek ar, Ricardo D Righetto, Jessica Heebner, Martin Obr, Xianjun Zhang, Saik at Chakrab ort y , Grigory T agiltsev, Alicia K Michael, Soﬁe V an Dorst, et al. T ow ard comm unity-driv en visual proteomics with large-scale cry o-elec tron tomography of c hlamydomonas reinhardtii. Mole cular Cel l , 86(1):213–230, 2026. [21] Marten L. Chaillet, Sander Ro et, Remco C. V eltk amp, and F riedrich F¨ orster. pytom-matc h-pick: A tophat-transform constraint for automated classiﬁcation in template matching. Journal of Structur al Biolo gy: X , 11:100125, June 2025. [22] Manuel Guizar-Sicairos, Sam uel T. Th urman, and James R. Fienup. Eﬃcient subpixel image registration algorithms. Optics L etters , 33(2):156–158, January 2008. 23 [23] V. K. Khersonskii, A. N. Mosk alev, and D. A. V arshalovic h. Quantum The ory Of Angular Mo- mentum . W orld Scientiﬁc Publishing Company , 1988. Accepted: 2021-08-20T11:39:57Z. 24

Fast Volume Alignment by Frequency-Marched Newton

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment