Alberto Pepe Dissertation - Full

• Original PDF: https://bmva-archive.org.uk/theses/2025/2025-pepe.pdf

Building the Dialog

I want to make a solveit dialogue for Alberto Pepe's dissertation, which I originally had as PDF.

Then after creating a free account on DataLab.io, I used that to generate Markdown and HTML versions of it, and uploaded them to solveit as pepe-diss.md and pepe-diss.html. The HTML version includes figures!

Hi, can you help me build a solveit dialogue so I can do a deep reading of this dissertation? I'd like any section headings or headers to signal subdivision into new separate blocks. I'd like the figures to be included, along with all charts and equations.

Data Lab already did a good job of converting the math. So there's one set of it in the Markdown, another set in the HTML. I suspect since you're mostly going to be building Markdown,... You could just import using Markdown for the most part and then only turn to the HTML to find the figures to fill in.

That makes sense? Do you think he could get started? So you're going to create many, many, many new note cells. And I suggest that you always append to the end of the dialogue so that they don't end up in some weird other order.

May take so long to do that you might hit time out cancellation of your operations in which the task will be only partially completed

So you'll need some way to resume from your previous progress based on which notes have already been added to this dialogue.-- We don't need to re-add the same note twice, And we can skip over any notes that have already been added and get to the main good stuff.

I suggest you first output a description of your strategy for how you're going to tackle this. It is a long dissertation. With many figures and tons of mathematics. PDF version is 228 pages. So, write out your strategy first and then pause, and then I will read your strategy and tell you when to execute it.

As I say, Data Lab has already done the heavy lifting for you. You don't need to do any OCR. You can go directly to the Markdown and or HTML. Again, I think favoring the markdown for most things and only using the HTML for the figures may be a good plan, but it's up to you.

Great, given the kinds of timeouts that you're subject to, we might want to proceed by having you do each phase one at a time. Do you want to try to parse the HTML yourself, or do you want to write some Python code to parse it and save the images, and then I'll execute that Python? I think I get longer timeout limits than you do when executing code.

Wow, that executed instantaneously. Makes me a little suspicious.

Okay, I guess we're good to go. Ready for phase two?

Uh-oh. What do you think went wrong? I think using dialogue helper is the right way to go, but perhaps your uppercase n isn't right or something else. Do you wanna maybe review the list of available dialog helper functions?

looks like that's not the right way to do that. Can you maybe double check that?

Okay, that looks like it actually worked. it looks pretty good. I see the math, I see the figures, lots of little subdivided individual notes. Are we done or were there additional phases still to be done?

Machine Learning with Geometric Algebra: Multivectors for Modelling, Understanding and Computing

The crest of the University of Cambridge, featuring a shield divided into four quadrants by a cross. The quadrants contain four lions passant guardant in pale, and the center of the cross contains a closed book.

University of Cambridge crest

Alberto Pepe

Supervisor: Prof. Joan Lasenby

Department of Engineering
University of Cambridge

This dissertation is submitted for the degree of
Doctor of Philosophy

Darwin College

August 2025

To my nephew Gian Maria, for the light and colour he brought into our lives

Declaration

This thesis is the result of my own work and includes nothing which is the outcome of work done in collaboration except as declared in the preface and specified in the text. It is not substantially the same as any work that has already been submitted, or, is being concurrently submitted, for any degree, diploma or other qualification at the University of Cambridge or any other University or similar institution except as declared in the preface and specified in the text. It does not exceed the prescribed word limit for the relevant Degree Committee.

Alberto Pepe
August 2025

Acknowledgements

When I began my PhD, I saw it as a marathon: every day, I aimed to move forward, however small the step. Some days flew by; others dragged. Now, approaching the finish line, I realise I could never have done this alone. I may not get another chance like this thesis to properly thank everyone, so I'll do my best not to leave anyone out.

First and foremost, I thank my supervisor, Joan Lasenby. It's often said a PhD depends greatly on who guides you, and I couldn't agree more. Thank you for your contagious optimism, the many hours spent discussing our research, the freedom to explore whatever sparked my interest, your support for my work in industry, and your patient advice as you read through countless and often messy drafts. This opportunity was truly once in a lifetime, and all my growth as a researcher I owe it to your guidance.

To my mum, Silvana: thank you for always supporting my choices and encouraging me to explore the world. Your strength and courage will stay with me forever.

To my sister Alessia and my brother Carmine: I've always admired your brilliance and drive. No matter how far apart we are, I carry a piece of you with me.

To my entire family, my roots and constant center of gravity: my dad Salvatore, Nonna Anna, aunts Gerardina, Alda, and Maria Pia; Andrea and Lorenza; Ilaria and Gian Maria; Francesca, Filippo, Francesco, Federico, Angelo; zio Giuseppe and zia Mariaconsiglia; zio Vito, zio Raffaele, and zio Mino.

Thank you to B. R. for your guidance, the long research discussions, and all the in-between moments I deeply treasure.

To Antonio V. and Alda, my chosen siblings since birth. Few are lucky to have such enduring friendship, and I'm proud to be among them.

To my high school friends: Adua, Camilla, Manuela, Margherita, Mariachiara, Michela, Luca, and Valentina; thank you for the memories, the support, and your constant belief in me. You are my safe haven.

To Antonio A. and Eleonora, thank you for the laughter and adventures since Shanghai, and for all those still to come.

To Piera and Gianluca, from my Turin days, thank you for the lasting bond, built on mutual admiration and affection.

To Deesha, Cindy, Fariza, and Maria, some of the best friends I've made here, thank you for the unforgettable moments and for letting me be fully myself.

Thank you to Cambridge, and to the many brilliant, inspiring students who've shaped my perspective on the world. It's been an honour to be part of such a unique place.

To Darwin College, my second home, and especially Sabrina: thank you for always brightening the day.

To Qian, Haris, Joana, Yuxin, Zhening, Gabriel, Pooya, Ralph, and the entire SigProc Lab. Thank you for sharing your ideas and for making the Engineering Department feel less gray.

To Jasivan: our summer in the Evergreen State made my Microsoft experience unforgettable. Thank you for the baseball, Vegas, and everything in between.

To every book and every piece of music I turned to over these years that has kept me grounded and focused.

And finally, to Namaste Cambridge and their legendary £9.49 chicken biryani: you got me through more than a few all-nighters, and for that, you have my gratitude.

Abstract

Geometric Algebra (GA) has been successfully applied in several fields, including physics, graphics, and robotics, but its potential in Machine Learning (ML) and Deep Learning (DL) remains largely unexplored. This thesis addresses that gap by investigating the application of GA to a variety of ML tasks that are inherently geometric in nature. Our premise is simple yet powerful: if, as it is often said, ML is a clever rebranding of linear algebra, then geometric problems in ML deserve to be tackled with GA, an extension of linear algebra designed to represent geometric objects and perform transformations on them naturally and compactly.

The GA framework provides a unified language to represent, understand, and manipulate geometric entities using multivectors, rotors, and sandwich products. The result is a set of tools that are not only mathematically elegant but also highly effective in practice, spanning a broad range of domains: Chapter 2 tackles regression tasks on rotation groups and molecular geometry optimization; Chapter 3 explores protein modelling and structure prediction; Chapter 4 addresses 3D camera pose estimation and 3D line alignment; Chapter 5 focuses on the solution of partial differential equations for computational fluid dynamics and electromagnetism.

We demonstrate that GA can be employed as a versatile, practical and principled framework for building geometry-aware ML systems independently of the employed architecture. By embedding geometric priors directly into the model architecture via GA, we unlock several advantages, including lower regression errors, robustness to noise and transformations and interpretability of intermediate computations. More broadly, this thesis is an invitation to rethink how neural networks should represent, model, and transform geometric data, grounding these operations in an algebra that reflects the structure of the problem itself. As ML continues to engage with increasingly complex and structured data, the need for such expressive representations will only grow, and GA offers a compelling framework to meet that need.

Table of contents

• 4.2.3 Method . . . . . 87
  • 4.2.4 Poses as motors in 1d-Up CGA . . . . . 87
  • 4.2.5 Architecture: InceptionV3 + sandwich product layers . . . . . 90
  • 4.2.6 Experiments . . . . . 91
  • 4.2.7 Datasets: Cambridge Landmarks and 7 Scenes . . . . . 91
  • 4.2.8 Error metrics . . . . . 92
  • 4.2.9 Implementation details . . . . . 92
  • 4.2.10 Results . . . . . 94
  • 4.2.11 Geometric interpretation of intermediate outputs . . . . . 98
• 4.3 3D Line Alignment with Define, Refine, Align . . . . . 100
  • 4.3.1 Related work on partial-to-partial registration . . . . . 102
  • 4.3.2 Method . . . . . 103
  • 4.3.3 3D Lines as bivectors in 1d-Up CGA . . . . . 103
  • 4.3.4 Architecture: attentional, rotational and equivariant modules . . . . . 106
  • 4.3.5 Datasets: Structured3D and Semantic3D . . . . . 109
  • 4.3.6 Error metrics . . . . . 110
  • 4.3.7 Implementation details . . . . . 114
  • 4.3.8 Results . . . . . 114
• 4.4 Conclusion . . . . . 116

5 Partial Differential Equations (PDEs) 119

• 5.1 Introduction . . . . . 119
• 5.2 Background . . . . . 120

Part I: Computational Fluid Dynamics 122

• 5.3 GA-ReLU: an activation function for Navier-Stokes PDEs . . . . . 123
  • 5.3.1 Problem definition . . . . . 123
  • 5.3.2 Data generation . . . . . 127
  • 5.3.3 Training details . . . . . 127
  • 5.3.4 Experiments . . . . . 127
  • 5.3.5 Results . . . . . 128
• 5.4 Fengbo: a neural operator for 3D PDEs over irregular geometries . . . . . 129
  • 5.4.1 Method . . . . . 130
  • 5.4.2 Architecture . . . . . 130
  • 5.4.3 Experiments . . . . . 135
  • 5.4.4 Implementation details . . . . . 139
  • 5.4.5 Training details . . . . . 143

5.4.6 Ablation study . . . . . 144

5.4.7 Notes on computational complexity . . . . . 149

Part II: Electromagnetism 153

5.5 STAResNet: working in STA to solve Maxwell’s PDEs . . . . . 154

5.5.1 Problem definition . . . . . 154

5.5.2 Architecture: a ResNet in STA . . . . . 158

5.5.3 Experiments in 2D . . . . . 161

5.5.4 Experiments in 3D . . . . . 169

5.6 Conclusions . . . . . 176

6 Conclusions 179

6.1 Conclusions . . . . . 179

6.2 Contributions . . . . . 180

6.3 Future Work... . . . . . 181

6.4 ...and Beyond . . . . . 181

References 183

List of figures

4.1 The CGAPoseNet+GCAN architecture. The output of the InceptionV3 network is reshaped to obtain a set of motor coefficient proposals. Motors are objects in the 1D-Up Algebra \(\mathcal{G}_{4,0}\) with scalar, bivector and quadrivector parts, giving a total of 8 real coefficients. These coefficients are used to build motors \(\in \mathcal{G}_{4,0}\) as input to the Geometric Clifford Algebra Network (GCAN). A motor represents a rotation and a translation, and it is hence a suitable representation for camera poses. The GCAN works in \(\mathcal{G}_{4,0}\) space and has weights, biases and outputs that are also motors, and hence interpretable as poses. The GCAN narrows down the proposals to a single motor through a geometric understanding of the scene. . . . . 89

4.2 The original CGAPoseNet pipeline, which is geometry-agnostic and adds a significant number of parameters compared to baseline PoseNet. . . . 91

4.3 Train (solid line) and validation (dashed line) losses for selected datasets. CGAPoseNet+GCAN attains the lowest loss profile due to its Geometric Clifford Algebra layers. . . . . 93

4.4 Ground truth and predicted translation component of the pose over the test set for selected datasets. . . . . 95

4.5 Ground truth and predicted rotation component of the pose over the test set for selected datasets. . . . . 96

4.6 Translation error over the test set for selected datasets. . . . . 97

4.7 Rotation error over the test set for selected datasets. . . . . 97

4.8 Two views of input and output poses of the GCAN layers for a test image in the Old Hospital dataset (relu activation). . . . . 98

4.9 Two views of input and output poses of the GCAN layers for a test image in the Old Hospital dataset (tanh activation). . . . . 99

4.10 Average input and output poses of the GCAN layers for a test image in the Old Hospital dataset (tanh activation) with 3 (left) and 4 (right) sandwich product layers in the GCAN. . . . . 100

4.11 Effects of \(\lambda\) in \(G(4, 0, 0)\) on pose estimation. . . . . 104

4.12 Define, Refine, Align (DRA): The DRA pipeline to estimate the pose \(M\) that aligns the source line bundle \(\mathcal{L}_S\) to the target line bundle \(\mathcal{L}_T\) via \(\mathcal{L}_T = M\mathcal{L}_S\tilde{M}\) . The Define module is an attention-based feature extractor to identify putative line matches in the two frames \(\mathcal{L}_{d,1}, \mathcal{L}_{d,2}\) . The Refine module consists of an Equivariant Module \(\phi\) and a Rotational Module \(\rho\) , that map lines to poses in \(G(4, 0, 0)\) . Lastly, the Align module averages poses \(M_{\rho,1}, M_{\rho,2}, M_{\phi,1}, M_{\phi,2}\) and normalizes the output to obtain \(\tilde{M}\) , aligning the two bundles. . . . . 105

4.13 The Define block, adapted from the Feature Extractor of [116]. The \(\ell_0\) and \(\ell_0\) components of lines \(\ell\) in \(\mathcal{L}_S, \mathcal{L}_T\) are processed separately through cascaded MLPs. Self- and cross-attention module are used to extract features from a graph representation of the line bundles. Such features are reshaped and normalised into line bundles \(\mathcal{L}_{d,1}, \mathcal{L}_{d,2}\) , containing fewer, matching lines. . . . . 107

4.14 The Equivariant Module \(\phi\) in the Refine block. It consists of a cascade of MV linear and geometric product layers, yielding poses \(M_\phi\) invariant to transformations \(\Gamma(\cdot)\) . . . . . 107

4.15 The Equivariant Module \(\rho\) in the Refine block. Lines are transformed into proposals of poses, which are rotated, translated and downsampled into \(M_\rho\) . . . . . 108

4.16 Two alignment examples (top and bottom) for the test set of Structured3D. (a) The source and target line bundles \(\mathcal{L}_S, \mathcal{L}_T\) ; (b)-(c)-(d) \(\mathcal{L}_S\) rotated and translated with ground truth motor \(M\) (blue) and with estimated motor \(\tilde{M}\) in (b) high curvature space (small \(\lambda\) ); (c) low curvature space (large \(\lambda\) ); (d) combined curvatures (red). . . . . 111

4.17 Two alignment examples (top and bottom) for the test set of Semantic3D. (a) The source and target line bundles \(\mathcal{L}_S, \mathcal{L}_T\) ; (b)-(c)-(d) \(\mathcal{L}_S\) rotated and translated with ground truth motor \(M\) (blue) and with estimated motor \(\tilde{M}\) in (b) high curvature space (small \(\lambda\) ); (c) low curvature space (large \(\lambda\) ); (d) combined curvatures (red). . . . . 115

5.1 Vector part of the (a) coefficient-wise ReLU \(\phi(\mathbf{x})\) (b) phase-dependent ReLU \(f(\mathbf{x})\) and (c) GA-ReLU \(\psi(\mathbf{x})\) . The scalar part, not reported, is identical for all three functions and equal to \(\phi(s)\) . . . . . 124

5.2 MSE error with error bars versus number of training data for Clifford ResNet with ReLU and GA-ReLU activation functions. . . . . 125

5.14 Intermediate layers outputs from \(P\) to \(Q_P\) within the 3D Clifford Geometry block for a test sample in the ShapeNet Car dataset. . . . . 141

5.15 Ablation on the number of hidden channels \(C\) . . . . . 145

5.16 Ablation on the number of FNO blocks \(F\) . . . . . 146

5.17 Ablation on the number of Fourier modes \(m\) . . . . . 147

5.18 Number of parameters (M) and model size (MB) as a function of \(C, R, m\) . 148

5.19 Computational complexity comparison of Fengbo and Transolver [188] for the 3D case. . . . . 150

5.20 Relative \(L_2\) norm promotion versus computational complexity for Fengbo (red arrow and ellipse) and Transolver’s (black arrow and ellipse) ablations. 151

5.21 3D GA approach: solving Maxwell’s PDEs through Clifford ResNet, a ResNet-inspired network working in \(\mathcal{G}_{3,0,0}\) . . . . . 156

5.22 STA approach: solving Maxwell’s PDEs through STAResNet, our ResNet-inspired network working in \(\mathcal{G}_{1,3,0}\) . . . . . 158

5.23 Training and validation losses versus number of epochs for 2D Maxwell’s PDEs for instances sampled at (a) 25s, (b) 50s, (c) 75s, (d) 100s. . . . . 160

5.24 (a) Mean squared error and (b) correlation between estimated and ground truth EM fields in the test set for varying \(\Delta t\) . We included real-valued ResNet for a more comprehensive comparison. . . . . 161

5.25 Ground truth (left column) shown against \(\mathbf{F}^2\) (right column) and \(F^2\) (centre column) over the test set for varying \(\Delta t\) . Values are shown on a \(32 \times 32\) spatial grid. . . . . 163

5.26 The 5 different obstacle configurations. The 3 unseen geometries are highlighted. . . . . 164

5.27 Training and validation losses versus number of epochs for the PDE in the presence of obstacles, for a total of 5 different obstacles configuration during the training phase. . . . . 165

5.28 (a) Mean squared error and (b) correlation between estimated and ground truth fields over the two test sets with seen and unseen obstacle configurations. . . . . 165

5.29 Ground truth (left column) shown against \(\mathbf{F}^2\) (right column) and \(F^2\) (centre column) over the test set with seen obstacles configurations (number 5 of Fig. 5.26) shown in top 2 rows and unseen obstacles configurations (number 7 of Fig. 5.26), shown in bottom 2 rows. Values are shown on a \(48 \times 48\) spatial grid for one chosen timepoint. . . . . 166

5.30 Test error over the estimated EM fields in the presence of seen and unseen obstacle geometries versus the number of trainable parameters. 167

5.31 (a) Mean squared error and (b) correlation between estimated and ground truth EM fields over test set versus rollout steps \(m\) for the 2D case. Different colors refer to datasets with different \(\Delta t\) in between frames. 168

5.32 2D GT \(F^2, \mathbf{F}^2\) , estimated \(F_{est} = \{\hat{F}^2, \hat{\mathbf{F}}^2\}\) and differences \(\Delta^2 = \{|F^2 - \hat{F}^2|, |\mathbf{F}^2 - \hat{\mathbf{F}}^2|\}\) for STAResNet and Clifford ResNet at different rollout steps: (a) \(m = 1\) , (b) \(m = 2\) , (c) \(m = 8\) and (d) \(m = 10\) . Here \(\Delta t = 5s\) . Red circles highlight regions of significant difference. Images are shown on a \(48 \times 48\) spatial grid. . . . . 170

5.33 2D GT \(F^2, \mathbf{F}^2\) , estimated \(F_{est} = \{\hat{F}^2, \hat{\mathbf{F}}^2\}\) and differences \(\Delta^2 = \{|F^2 - \hat{F}^2|, |\mathbf{F}^2 - \hat{\mathbf{F}}^2|\}\) for STAResNet and Clifford ResNet at different rollout steps: (a) \(m = 1\) , (b) \(m = 5\) , (c) \(m = 8\) and (d) \(m = 10\) . Here \(\Delta t = 8s\) . Red circles highlight regions of significant difference. Images are shown on a \(48 \times 48\) spatial grid. . . . . 171

5.34 Training and validation losses versus number of epochs for 2D Maxwell's PDEs for instances sampled at (a) 5s, (b) 8s, (c) 10s, (d) 15s. . . . . 172

5.35 Mean squared error between estimated and ground truth EM fields over test set versus rollout steps \(m\) for the 3D case. Different colors refer to datasets with different \(\Delta t\) in between frames. . . . . 173

5.36 Slices of the 3D GT \(F^2, \mathbf{F}^2\) , estimated \(\hat{F}^2, \hat{\mathbf{F}}^2\) and differences \(|F^2 - \hat{F}^2|, |\mathbf{F}^2 - \hat{\mathbf{F}}^2|\) for STAResNet and Clifford ResNet at different heights \(z\) of the volume at different rollout steps: (a) \(m = 3\) , (b) \(m = 5\) . \(\Delta t = 5s\) . Images are shown on a \(28 \times 28\) spatial grid. . . . . 174

5.37 Slices of the 3D GT \(F^2, \mathbf{F}^2\) , estimated \(\hat{F}^2, \hat{\mathbf{F}}^2\) and differences \(|F^2 - \hat{F}^2|, |\mathbf{F}^2 - \hat{\mathbf{F}}^2|\) for STAResNet and Clifford ResNet at different heights \(z\) of the volume at different rollout steps: (a) \(m = 8\) , (b) \(m = 10\) . \(\Delta t = 5s\) . Images are shown on a \(28 \times 28\) spatial grid. . . . . 175

List of tables

3.6 GDT_HA scores over the PSICOV150 dataset for different 3D projection strategies. . . . . 80

4.1 Number of trainable parameters for the three approaches. . . . . 91

4.2 Median translation and rotation errors over the test set for the 7 approaches. 94

4.3 Ablation study with different backbones for selected dataset. Results superior to the best PoseNet strategy are in bold. . . . . 100

4.4 Number of Parameters for Different Methods . . . . . 109

4.5 Results for the correspondence-free and correspondence-based line registration problem. . . . . 112

4.6 Relative performance of DRA compared to other methods. Positive values (in dark green) indicate a promotion (i.e. DRA’s error is lower), while negative values (in red) indicate a demotion (i.e. DRA’s error is higher). . . . . 113

5.1 The three steps of the Fengbo pipeline. . . . . 135

5.2 Training and testing errors in pressure and velocity prediction on the ShapeNet Car dataset. Fengbo results have been obtained with \(\alpha, \beta = \{5, 1\}\) for pressure and with \(\alpha, \beta = \{1, 50\}\) for velocity (see Section 5.4.6).137

5.3 Training and testing errors in pressure prediction on the Ahmed Body dataset. . . . . 137

5.4 Comparison of different models. \(d\) is the maximum degree of the graph, \(D\) is the feature space dimensionality. *: See notes on computational complexity. . . . . 140

5.5 Comparison of Fengbo’s 3D and 2D configurations. . . . . 142

5.6 Data representations and their processing steps with the Fengbo pipeline.143

5.7 Ablation on the impact of \(M\) . Reporting the relative \(L_2\) norm. . . . . 144

5.8 Ablation on the impact of \(\alpha, \beta\) . Reporting the relative \(L_2\) norm. . . . . 148

5.9 Complexity comparison for different datasets . . . . . 151

Nomenclature

Number Sets

Mathematical Symbols

Acronyms / Abbreviations

Chapter 1

Introduction and Background

I thought I had something so simple
to say. Something useful to
everybody.

— Marcello Mastroianni as Guido,
\(8\frac{1}{2}\) , dir. Federico Fellini (1963)

1.1 Introduction

Geometric Algebra (GA) is a mathematical framework based on the algebra of William Kingdon Clifford. GA is a powerful tool as it offers an intuitive manipulation of geometrical objects (scalars, vectors, oriented surfaces and volumes, etc.) and a simple way to perform geometric transformations by any amount and in any dimension.

Historically, GA has been more successful in some fields than in others. Researchers and engineers in theoretical physics or computer graphics, for example, have been making use of GA for several years. Experiments in GA for Machine Learning (ML) and Deep Learning (DL), on the other hand, have only started to appear in the literature in the last couple of years, despite early notions of a GA-based ML being proposed as early as the nineties. This is surprising, since many ML problems deal with intrinsically geometric data, and such problems are often solved with engineering-intense approaches, which, albeit effective, are very complex and not at all intuitive.

The goal of this dissertation is to bridge the gap between the GA and the ML worlds, and to show how GA can improve on and simplify different aspects of a typical ML pipeline. These aspects include, among others, representing geometrical objects as data in an intuitive, meaningful way; understanding intermediate outputs of a neural

network and boosting its accuracy at a reduced computational cost. These challenges are recurrent problems in the broader Deep Learning community that can be addressed elegantly and effectively using GA.

While it could be argued that GA is merely one alternative among several, more widely recognized approaches in the literature, the versatility of GA is undeniable. This thesis aims to demonstrate that with a few simple GA tools, a wide variety of problems, ranging from bioinformatics to computer vision and engineering, can be tackled using the same methods and the same underlying philosophy.

At the outset of the PhD project, the goal was to show how GA could simplify the handling of proteins. However, it quickly became clear that the same methods applicable to proteins could be extended to a broad range of other geometric objects. This insight revealed that GA could serve as a valuable tool not only in bioinformatics, but also within the larger field of geometric deep learning.

This thesis constitutes a summary of works published starting from September 2021.

• The first half of Chapter 2 has been published as:

  • Pepe, Alberto, Joan Lasenby, and Pablo Chacón. "Learning rotations." Mathematical Methods in the Applied Sciences 47.3 (2024): 1204–1217 [143].
  • The second half is a joint work with Scuola Normale Superiore in Pisa, Italy, and has been presented as: Pepe, Alberto, Federico Lazzari, Vincenzo Barone, Joan Lasenby. "Molecular Geometry Optimization through Rotor-based Evolutionary Algorithm", 13th International Conference on Clifford Algebras and Their Applications in Mathematical Physics (ICCA13), Holon, Israel (2023).

• Elements of Chapter 3 have been published as:

  • Pepe, Alberto, Joan Lasenby, and Pablo Chacon. "Geometric Algebra Models of Proteins for Three-Dimensional Structure Prediction: A Detailed Analysis." International Conference on Advanced Computational Applications of Geometric Algebra (ICACGA). Cham: Springer Nature Switzerland, 2022 [143, 144].
  • Pepe, Alberto, and Joan Lasenby. "Modeling orientational features via Geometric Algebra for 3D protein coordinates prediction." Mathematical Methods in the Applied Sciences 47.14 (2024): 11385–11404 [141].

• Pepe, Alberto, Joan Lasenby, and Pablo Chacon. "Using a graph transformer network to predict 3D coordinates of proteins via Geometric Algebra modelling." International Workshop on Empowering Novel Geometric Algebra for Graphics and Engineering. Cham: Springer Nature Switzerland, 2022 [145].

• Pepe, Alberto, Sven Buchholz, and Joan Lasenby. "Clifford Group Equivariant Neural Network Layers for Protein Structure Prediction." Northern Lights Deep Learning Conference. PMLR, 2024 [137].

• Chapter 4 has been published as:

  • Pepe, Alberto, and Joan Lasenby. "CGAPoseNet: Camera pose regression via a 1D-up approach to Conformal Geometric Algebra." arXiv preprint arXiv:2302.05211 (2023) [140].
  • Pepe, Alberto, Joan Lasenby, and Sven Buchholz. "CGAPoseNet+ GCAN: A Geometric Clifford Algebra Network for Geometry-aware Camera Pose Regression." Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, 2024 [142].
  • Pepe, Alberto, Yuxin Yao, and Joan Lasenby. "Define, Refine, Align: Correspondence-free 3D Line Alignment with Attentional, Equivariant and Rotational Layers." Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), Perception Beyond the Visible Spectrum (PBVS), 2025.

• Elements of Chapter 5 have been published as:

  • Pepe, Alberto, Sven Buchholz, and Joan Lasenby. "GA-ReLU: an activation function for Geometric Algebra Networks applied to 2D Navier–Stokes PDEs." ICLR 2024 Workshop on AI4DifferentialEquations In Science, 2024 [138].
  • Pepe, Alberto, Sven Buchholz, and Joan Lasenby. "STAResNet: a Network in Spacetime Algebra to solve Maxwell's PDEs." arXiv preprint arXiv:2408.13619 (2024) [139].
  • Pepe, Alberto, Mattia Montanari, and Joan Lasenby. "Fengbo: A Clifford Neural Operator pipeline for 3D PDEs in Computational Fluid Dynamics." ICLR 2025.

1.2 Fundamentals of GA

The branch of mathematics called Geometric Algebra (GA) was developed in the second half of the 19th century by William Kingdon Clifford [33]. GA was revived in 1966 thanks to David Hestenes, who reiterated the unifying power of GA and its relevance in physics due to the close ties with geometry [68]. GA provides an intuitive and unified framework for describing geometric entities and transformations acting on them, and it has found application in computer graphics, computer vision, robotics, physics and more. In the rest of this section we will introduce the fundamental operators of GA and the subalgebras that we employed throughout the thesis. For a more rigorous derivation we refer the reader to [103, 43]. All plots in this section have been generated using the Clifford package [61].

1.2.1 Defining a subalgebra

A Geometric Algebra of size \(m\) can be defined over a scalar field and a set of \(m\) independent, orthonormal basis vectors \(\{e_i\}_{i=1,\dots,m}\) . We indicate a generic closed subalgebra with \(\mathcal{G}_{p,q,r}\) or alternatively \(G(p,q,r)\) . A closed subalgebra \(\mathcal{G}_{p,q,r}\) has \(p\) basis vectors that square to 1, \(q\) basis vectors that square to -1 and \(r\) basis vectors that square to 0.

1.2.2 The geometric product

Elements in a GA are called multivectors. Elements of any type can be added or multiplied together. Each element has a grade associated with it. By grade we define the dimension of the hyperplane an object specifies. e.g. scalars are grade 0, vectors are grade 1, bivectors are grade 2, etc. Geometric algebra gets its name from the geometric product. The geometric product between two vectors is given by

\[ab = a \cdot b + a \wedge b \quad (1.1)\]

in which the scalar (or inner) product \(a \cdot b\) is the usual scalar product of linear algebra, which equal to the cosine of the angle between \(a\) and \(b\) when \(a, b\) are unit vectors, while the wedge (or outer) product \(a \wedge b\) produces a bivector (e.g. an oriented plane). The geometric product is hence the sum of a scalar and a bivector, that have different grade. Any multivector of a unique grade \(r\) that can be defined as

\[A = a_1 \wedge a_2 \wedge \dots \wedge a_r \quad (1.2)\]

is called a blade.

The reversion operator for a multivector is given by \(\tilde{\cdot}\) . The reverse of a scalar is equal to the scalar itself and the reverse of a vector is equal to the vector itself. For a multivector, we have that

\[\begin{aligned}(AB)^{\tilde{\cdot}} &= \tilde{B}\tilde{A} \\ (A + B)^{\tilde{\cdot}} &= \tilde{A} + \tilde{B}\end{aligned}\quad (1.3)\]

The general rule to reverse an \(r\) -blade is given by

\[\tilde{A}_r = (-1)^{\frac{r(r+1)}{2}} A_r \quad (1.4)\]

The grade projector operator is denoted by \(\langle A \rangle_r\) , where \(r\) is the grade we want to extract from \(A\) . This comes from the fact that a multivector in a \(m\) dimensional algebra can be written as

\[A = \sum_{i=0}^n \langle A \rangle_i \quad (1.5)\]

The geometric product between a multivector and its reverse gives the squared magnitude: \(|A|^2 = \langle A\tilde{A} \rangle_0\) . The reversion operator can be used to define the inverse of a multivector as

\[A^{-1} = \frac{\tilde{A}}{|A|^2} \quad (1.6)\]

It can be easily shown that \(A^{-1}A = 1\) . The dual of a multivector is defined as

\[A^* = AI_n^{-1} \quad (1.7)\]

where \(I_n\) is called the pseudoscalar, defined as \(I_n = e_1 \wedge e_2 \wedge \dots \wedge e_n\) . The pseudoscalar is the highest grade element in a GA. The product of a grade- \(n\) pseudoscalar \(I_n\) and a grade- \(r\) multivector \(A_r\) is a grade- \((n-r)\) multivector, and is called a duality transformation. The pseudoscalar interchanges inner and outer products:

\[A_r \cdot (B_s I) = \langle A_r B_s I \rangle_{r-(n-s)} = \langle A_r B_s I \rangle_{n-(r+s)} = \langle A_r B_s \rangle_{r+s} I = A_r \wedge B_s I \quad (1.8)\]

1.2.3 Geometric Algebra of the plane and of the space

The GA of the Euclidean plane is denoted by \(\mathcal{G}_{2,0,0}\) , with two basis vectors \(e_1, e_2\) . Being an \(n = 2\) dimensional GA (since \(p + q + r = 2 + 0 + 0 = 2 = n\) ), it is spanned by \(2^2 = 4\) elements, namely a scalar, two vectors \(e_1, e_2\) and the bivector \(e_1 e_2 = e_1 \wedge e_2\) .

The diagram illustrates the elements of the Geometric Algebra (GA) of Euclidean space. It starts with a scalar (represented by a dot), followed by three vectors: \(e_1\) (vertical), \(e_2\) (horizontal), and \(e_3\) (diagonal). These vectors combine to form three bivectors: \(e_{12}\) (a square in the \(e_1\) - \(e_2\) plane), \(e_{23}\) (a square in the \(e_2\) - \(e_3\) plane), and \(e_{13}\) (a square in the \(e_1\) - \(e_3\) plane). Finally, the three vectors combine to form a trivector \(e_{123}\) (a cube).

Figure 1.1: Elements of the Geometric Algebra (GA) of Euclidean space. The diagram shows the progression from a scalar (a dot) to three vectors (e1, e2, e3) to three bivectors (e12, e23, e13) and finally to a trivector (e123).

Fig. 1.1 Elements of \(\mathcal{G}_{3,0,0}\) , the GA of Euclidean space: one scalar, three vectors \(e_1, e_2, e_3\) , three bivectors \(e_{12}, e_{23}, e_{13}\) and one trivector \(e_{123}\) .

\(\mathcal{G}_{2,0,0}\) includes the concept of complex numbers, since the pseudoscalar of this algebra \(I = e_1 \wedge e_2 = e_1 e_2 \equiv e_{12}\) squares to \(-1\) , since \(I^2 = e_{12}^2 = e_{12} e_{12} = (e_1 e_2)(e_1 e_2) = -(e_1 e_1)(e_2 e_2) = -1\) . A scalar plus a bivector can be seen as a representation of a complex number, since \(Z = a + Ib \equiv a + \iota b\) , where \(\iota\) is the imaginary unit. Similarly, if \(X = ae_1 + be_2 = e_1(a + bI) = e_1 Z\) .

Adding a third basis vector \(e_3\) we form \(\mathcal{G}_{3,0,0}\) , the GA of Euclidean space (see Fig. 1.1). It has \(2^3 = 8\) elements, a scalar, three vectors, three bivectors ( \(e_{12}, e_{23}, e_{13}\) ) and one trivector ( \(e_{123} = e_1 \wedge e_2 \wedge e_3 = I_3\) , the pseudoscalar). The GA of space includes the quaternion algebra, since a quaternion \(\mathbf{q} = w + a\mathbf{i} + b\mathbf{j} + c\mathbf{k}\) can be represented as a multivector \(A = w + ae_{12} + be_{13} + ce_{23}\) .

Throughout this thesis, the GA notation will be employed for the majority of mathematical quantities and geometric objects, unless explicitly stated otherwise. For example, lowercase, non-bold letters will be used to represent vectors, uppercase letters will denote bivectors and so on.

1.2.4 Geometric transformations

The concepts below assume a three-dimensional GA, but can be extended to any GA of arbitrary dimension.

Reflections

The reflection of a vector \(a\) in the plane orthogonal to the unit vector \(n\) is given by \(a' = a_{\perp} - a_{\parallel}\) . Hence, we first resolve \(a\) into its parallel and perpendicular components to \(n\) . This is done as \(a = n^2 a = n(n \cdot a + n \wedge a) = a_{\parallel} + a_{\perp}\) , in which \(a_{\parallel} = (a \cdot n)n\) and \(a_{\perp} = n(n \wedge a)\) . This is equivalent to

\[a' = n(n \wedge a) - n(a \cdot n) = -n(a \cdot n + a \wedge n) = -nan \quad (1.9)\]

A 3D diagram illustrating reflection in Geometric Algebra. A magenta plane labeled Π is shown. A cyan plane labeled Π' is its reflection. A unit normal vector n is shown as a black arrow pointing from the magenta plane towards the cyan plane. The reflection is performed across a plane with unit normal vector n.

Fig. 1.2 Reflection in Geometric Algebra: a multivector \(\Pi\) (e.g., a plane represented as a bivector) is reflected across a plane with unit normal vector \(n\) using the formula \(\Pi' = -n\Pi n\) . This reflection preserves the grade of \(\Pi\) .

The formula \(a' = -nan\) works for spaces in any dimension and it leaves lengths and angles unchanged. For bivectors \(B = a \wedge b\) , we have that

\[B' = (-nan) \wedge (-nbn) = \frac{1}{2}(nannbn - nbnnan) = \frac{1}{2}n(ab - ba)n = nBn \quad (1.10)\]

Sandwiching a multivector between a vector always preserves the grade of the multivector (see Fig. 1.2). The dual vector \(IB\) is subject to the same transformation \(I(nBn) = n(IB)n\) .

Rotations

A rotation in the plane is generated by successive reflections in the planes perpendicular to two unit vectors \(m\) and \(n\) . Let's consider two reflections: \(b = -mam\) , and then \(c = -nbn = -n(-mam)n = nmamn\) . We define \(R = nm\) , so to write \(c = Ra\tilde{R}\) , representing a rotation. \(R\) is called a rotor. A rotor is the geometric product of two unit vectors (see Fig. 1.3). We can rewrite it as

\[R = nm = n \cdot m + n \wedge m = \cos(\theta) + n \wedge m \quad (1.11)\]

Diagram illustrating rotation in Geometric Algebra. A multivector Π (red plane) is rotated to Π'' (green plane) by applying two successive reflections across planes with unit normals n and m. The rotation is given by Π'' = mnΠnm = RΠR̃, where R = mn is the rotor. The diagram shows the planes Π, Π', and Π'' in a 3D coordinate system, with the rotor R indicated by a curved arrow.

Fig. 1.3 Rotation in Geometric Algebra: a multivector \(\Pi\) is rotated by applying two successive reflections across planes with unit normals \(n\) and \(m\) . The rotation is given by \(\Pi'' = mn\Pi nm = R\Pi\tilde{R}\) , and preserves the grade of \(\Pi\) . The geometric product \(R = mn\) is called a rotor.

The magnitude of the bivector part is given by \(\sin^2(\theta)\) . We define the unit bivector \(B\) in the \(m \wedge n\) plane as

\[B = \frac{m \wedge n}{\sin(\theta)} \quad (1.12)\]

with \(B^2 = -1\) .

Recalling the complex number duality, the rotor can be written as \(R = \cos(\theta) - B \sin(\theta)\) , or similarly \(R = \exp(-B\theta)\) , in which the exponential function is defined by its Taylor expansion where the multiplications are replaced by the geometric product. For a rotation of an angle \(\theta\) , we have \(R = \exp(-B\theta/2)\) . Summarizing, the rotation operation:

\[a' = Ra\tilde{R} = e^{-B\theta/2}ae^{B\theta/2} \quad (1.13)\]

where \(R\tilde{R} = 1\) . Rotations preserve lengths and angles:

\[a' \cdot b' = \frac{1}{2}(Ra\tilde{R}Rb\tilde{R} + Rb\tilde{R}Ra\tilde{R}) = \frac{1}{2}R(ab + ba)\tilde{R} = a \cdot b R\tilde{R} = a \cdot b \quad (1.14)\]

The reverse of the rotation is given by \(a = \tilde{R}a'R\) .

The single \(R\) acting by itself on a vector has no geometric significance: the formula \(b = Ra\) only works for vectors in the plane of rotation. For this reason, we say that GA has a double-cover representation of the rotation group: since it is double sided, \(R\) and

A 3D visualization showing a 90-degree rotation of geometric objects. A blue vertical plane and a blue circle are shown, along with an orange horizontal plane and an orange circle. A blue vector points upwards, and an orange vector points to the left. Two points, one blue and one orange, are also shown. The rotation is indicated by a curved arrow around the vertical axis.

Fig. 1.4 Visualization of some \(90^\circ\) rotations in 3D Geometric Algebra using sandwich products. The orange objects—vector, circle, plane, and point—are rotated to the blue counterparts by the sandwich product \(X \mapsto RX\tilde{R}\) , where \(\tilde{R}\) is the reverse of \(R\) . Rotation via the sandwich product is independent of the type of geometric object \(X\) , applying uniformly to vectors, bivectors, and points alike.

\(-R\) yield the same rotation. Multivectors can be rotated with the same principle (see Fig. 1.4). For the group combination rule, if \(b = R_1 a \tilde{R}_1\) , \(c = R_2 b \tilde{R}_2\) , then \(c = R a \tilde{R}\) , where \(R = R_1 R_2\) .

In the rest of this thesis we will represent rotations in the Euclidean space in terms of rotors and bivectors, and employing two kinds of rotor-to-bivector maps, namely (i) the exponential map, for which

\[B = -2 \log R \quad (1.15)\]

and (ii) the Cayley transform [68], for which

\[B = \frac{1 - R}{1 + R} \quad (1.16)\]

1.2.5 Conformal Geometric Algebra (CGA)

Let \(x\) be a vector in a Geometric Algebra \(\mathcal{G}_{p,q,r}\) of dimension \(n = p + q + r\) . CGA extends the Geometric Algebra to \(\mathcal{G}_{p+1,q+1,r}\) by introducing two basis vectors, \(e\) and \(\bar{e}\) , with \(e^2 = +1\) and \(\bar{e}^2 = -1\) . Since \(e_i \cdot e = e_i \cdot \bar{e} = 0\) , then also \(x \cdot e = x \cdot \bar{e} = 0\) . Having introduced \(e\) and \(\bar{e}\) , we can create the null vectors following Hestenes' notation:

\[\begin{aligned} n &= e + \bar{e} \\ \bar{n} &= e - \bar{e} \end{aligned} \quad (1.17)\]

An alternative notation, more common today, is given by

\[\begin{aligned} n_\infty &= e + \bar{e} \\ n_0 &= \frac{1}{2}(\bar{e} - e) \end{aligned} \quad (1.18)\]

which represent the point at infinity and the origin, respectively. They are null since \(n_\infty^2 = (e + \bar{e})(e + \bar{e}) = e^2 + 2e \cdot \bar{e} + \bar{e}^2 = 1 + 0 - 1 = 0\) . Moreover, we can show that \(n_\infty \cdot n_0 = e^2 - \bar{e}^2 = 2\) and \(x \cdot n_\infty = x \cdot n_0 = 0\) . CGA defines a mapping of the kind

\[x \in \mathcal{G}_{p,q,r} \longrightarrow F(x) \in \mathcal{G}_{p+1,q+1,r} \quad (1.19)\]

in which \(F(x)\) is defined, using Hestenes' notation, as

\[F(x) = x^2 n + 2x - \bar{n} \quad (1.20)\]

or, alternatively, as

\[\begin{aligned} F(x) &= -\frac{1}{2}(x - e)n_\infty(x - e) \\ F(x) &= \frac{1}{2}(x^2 n_\infty + 2x - n_0) \end{aligned} \quad (1.21)\]

The factor of \(1/2\) is chosen so that \(F(x) \cdot n_\infty = -1\) . \(F(x)\) is always a null vector for any \(x\) . A single point \(x \in \mathcal{G}_{p,q,r}\) corresponds to a family of null vectors \(X = \lambda(x^2 n_\infty + 2x - n_0) \in \mathcal{G}_{p+1,q+1,r}\) . If \(x = 0\) , then \(F(x) = -\frac{1}{2}n_0\) , i.e. the origin. Its inverse is \(n_\infty\) , i.e. mapping points at infinity. If we are dealing with a 3D space, the equivalent CGA will be \(\mathcal{G}_{4,1,0}\) , i.e. a 5D space.

1.2.6 Transformations in CGA

Rotations

Recall that in \(\mathcal{G}_{p,q,r}\) a rotation is expressed by \(x \rightarrow Rx\tilde{R}\) . It can be shown that to rotate a conformal object \(F(x) \rightarrow RF(x)\tilde{R}\) holds:

\[RF(x)\tilde{R} = \frac{1}{2}R(x^2n_\infty + 2x - n_0)\tilde{R} = \frac{1}{2}[x^2Rn_\infty\tilde{R} + 2Rx\tilde{R} - Rn_0\tilde{R}] \quad (1.22)\]

A rotor contains only even blades, so it commutes with \(n_\infty\) . Moreover, \(Rn_0\tilde{R} = n_0\) . We hence have:

\[RF(x)\tilde{R} = \frac{1}{2}(\hat{x}^2n_\infty + 2\hat{x} - n_0) \quad (1.23)\]

where \(\hat{x} = Rx\tilde{R}\) .

Inversions

By inversion we mean the mapping \(x \rightarrow x/x^2\) . It can be shown that reflecting \(F(x)\) in \(e\) gives:

\[-eF(x)e = -\frac{1}{2}[x^2en_\infty e + 2exe - en_0e] = \frac{1}{2}x^2[\frac{1}{x^2}n_\infty + 2\frac{x}{x^2} - n_0] = x^2F\left(\frac{x}{x^2}\right) \quad (1.24)\]

The scale and the sign are irrelevant, so we can safely discard the minus sign and define the inversion as \(e(\cdot)e\) (see Fig. 1.5).

The figure shows a 2D Cartesian coordinate system with a grid. The x-axis and y-axis both range from -1.0 to 3.0 with major grid lines every 0.5 units. A blue vector, labeled \(a = e_1 + 2e_2\) , originates at (0,0) and terminates at (1,2). A red vector, labeled \(F(a_{inv}) = \frac{a}{a^2} = eF(a)e\) , originates at (0,0) and terminates at (0.2, 0.4). The red vector is the inversion of the blue vector.

A 2D coordinate system showing a blue vector 'a = e1 + 2e2' starting at the origin (0,0) and pointing to the point (1,2). A red vector, labeled 'F(a_inv) = a/a^2 = eF(a)e', starts at the origin and points to the point (0.2, 0.4), which is the inversion of vector 'a'.

Fig. 1.5 Inversion in Conformal Geometric Algebra: Euclidean vector \(a\) (blue) is transformed via CGA using sandwich products. Its inversion (red) is given by \(eF(a)e\) .

Translations

A translation takes place if \(x \rightarrow x + a\) . We consider the rotor \(R = T_a = \exp(na/2)\) , and use the Taylor expansion on the exponential. Since \(n\) is a null vector and \(an = -na\) , all higher terms are zero, we can rewrite it as:

\[R = T_a = \exp\left(\frac{na}{2}\right) = 1 + \frac{na}{2} + \frac{1}{2}\left(\frac{na}{2}\right)^2 + \dots \simeq 1 + \frac{na}{2} \quad (1.25)\]

Some useful properties:

• \(Rn_\infty\tilde{R} = n_\infty\) since \(n_\infty\) represents a point at infinity
• \(Rn_0\tilde{R} = n_0 - 2a - a^2n_\infty \rightarrow R(-n_0)\tilde{R} = F(a)\) since \(-n_0\) represents the origin
• \(Rx\tilde{R} = x + n_\infty(a \cdot x)\) by relation of dot and geometric product.

The above properties allow us to show that

\[RF(x)\tilde{R} = \frac{1}{2}((x+a)^2 + 2(x+a) - n_0) = F(x+a) \quad (1.26)\]

An example is shown in Fig. 1.6.

A 3D diagram illustrating a translation in Conformal Geometric Algebra. A pink plane labeled Π is shown at the bottom, and a green plane labeled Π' is shown above it. An arrow labeled T_a points from the pink plane to the green plane, representing the translation vector. The planes are shown within a 3D coordinate system with gray axes.

Fig. 1.6 Translations in Conformal Geometric Algebra: the green plane is the result of translating the original pink plane by the Euclidean vector \(a\) . The translation is performed via the sandwich product \(\Pi' = T_a \Pi \tilde{T}_a\) , with \(T_a = 1 - \frac{1}{2}a\bar{e}\) , \(\tilde{T}_a = 1 + \frac{1}{2}a\bar{e}\) .

Dilations

A dilation in GA is given by \(x \rightarrow \alpha x\) . We consider the rotor \(R = D_\alpha = \exp(\frac{\alpha}{2}e\bar{e})\) . We have that

\[\begin{aligned} D_\alpha F(x) \tilde{D}_\alpha &= \exp\left(\frac{\alpha}{2}e\bar{e}\right) \frac{1}{2}(x^2 n_\infty + 2x - n_0) \exp\left(-\frac{\alpha}{2}e\bar{e}\right) = \\ &= \frac{1}{2}(x^2 \exp(\alpha e\bar{e}) n_\infty + 2x - \exp(\alpha e\bar{e}) n_0) = \\ &= \exp(\alpha) \frac{1}{2}(\hat{x}^2 n_\infty + 2\hat{x} - n_0) \end{aligned} \quad (1.27)\]

with \(\hat{x} = \exp(-\alpha)x\) . Dilations in GA are represented by a rotor \(R = D_\alpha\) in conformal space (see Fig. 1.7). A summary of the transformations in CGA is given in Table 1.1.

Table 1.1 Summary of common geometric transformations in GA and CGA

A 3D visualization showing a dilation in Conformal Geometric Algebra. A small yellow sphere labeled Σ is centered on a light gray plane. A larger teal sphere labeled Σ' is also centered on the same plane, representing the result of dilating Σ by a scale factor α = 6. The dilation is performed via the sandwich product Σ' = D_α Σ D̃_α, with D_α = exp(α/2 e ē).

Fig. 1.7 Dilations in Conformal Geometric Algebra: The teal sphere \(\Sigma'\) is the result of dilating the original yellow sphere \(\Sigma\) by the scale factor \(\alpha = 6\) . The dilation is performed via the sandwich product \(\Sigma' = D_\alpha \Sigma \tilde{D}_\alpha\) , with \(D_\alpha = \exp(\frac{\alpha}{2}e\bar{e})\) .

1.2.7 Objects and incidence

Point pairs, lines, planes, circles and spheres are all represented by blades in the 5D CGA. They can be transformed with rotors. The dual of these objects are used to extract meaningful features of these objects. The simplest object in CGA is represented by a blade of the kind \(A \wedge B\) , that represents a point pair.

Lines

If a line \(L\) passes through two points \(a, b\) with CGA representation \(A\) and \(B\) , then the line can be represented by \(L = A \wedge B \wedge n_\infty\) . The representation \(X\) of any point on the line will satisfy \(X \wedge L = 0\) . This is an incidence relation, i.e. \(X \wedge Y \wedge \dots \wedge Z = 0\) , and it is invariant under rotor/reflection operations. Let's consider two points \(x_1, x_2\) on a line \(L\) with corresponding representation \(X_1, X_2\) in \(\mathcal{G}(p+1, q+1, r)\) . Without loss of generality we assume that the line passes through the origin in the direction of \(e_1\) . Hence \(X_i\) only contains the vectors \(n_\infty, n_0, e_1\) , as any point on the line must have form \(x = \lambda e_1\) . If \(X\) is any point on the line, then the incidence relation is given by

\[X \wedge X_1 \wedge X_2 \wedge n_\infty = X \wedge L = 0 \quad (1.28)\]

Planes

Without loss of generality, consider the plane spanned by \(e_1, e_2\) . A point \(x\) lies on the plane if \(x = \lambda e_1 + \mu e_2\) . Its conformal representation will only contain \(n_\infty, n_0, e_1, e_2\) as \(X = \frac{1}{2}(x^2 n_\infty + 2(\lambda e_1 + \mu e_2) - n_0)\) . Let \(\Pi = X_1 \wedge X_2 \wedge X_3 \wedge n_\infty = \beta(n_\infty \wedge n_0 \wedge e_1 \wedge e_2)\) , so if a point \(X\) lie on the plane we have that

\[X \wedge X_1 \wedge X_2 \wedge X_3 \wedge n_\infty = X \wedge \Pi = 0 \quad (1.29)\]

Circles

Consider \(\mathcal{G}_{2,0,0}\) i.e. the usual Euclidean plane with basis \(e_1, e_2, e_1^2 = 1, e_2^2 = 1\) . Let \(L\) be the line \(x = 1\) and \(a\) a point with coordinates \((x, y)\) . Assume we want to invert the points on this line: this gives the set \(L \rightarrow \left(\frac{1}{1+y^2}, \frac{y}{1+y^2}\right) = (x', y')\) . It can be shown that \([x' - 1/2]^2 + y'^2 = (1/2)^2\) is the equation of a circle centred at \((1/2, 0)\) and radius \(1/2\) . We mapped a line into a circle. If \(X\) is a point on a line, then the inverted point is on the circle is \(X'\) . Since the incidence relation is invariant to inversion, then a point \(X'\) belongs on a circle if

\[\begin{aligned} X' \wedge X'_1 \wedge X'_2 \wedge X'_3 &= e(X \wedge X_1 \wedge X_2 \wedge X_3)e = eXe \wedge eX_1e \wedge eX_2e \wedge eX_3e = \\ &= X'_1 \wedge \mathcal{C} = 0 \end{aligned} \quad (1.30)\]

We know that a line can be represented as a trivector of the form \(L = X_1 \wedge X_2 \wedge n_\infty\) , then any circle as \(\mathcal{C} = X_1 \wedge X_2 \wedge X_3\) , where no \(X_i\) is a multiple of \(n_\infty\) : a line is a circle passing through a point at infinity. Examples of lines and circles in CGA are shown in Fig. 1.8.

A 3D perspective view of a horizontal plane. A straight line, labeled \(L = X_1 \wedge X_2 \wedge n_\infty\) , passes through the plane. A circle, labeled \(C = X_1 \wedge X_2 \wedge X_3\) , is also shown on the plane, appearing as a blue ellipse due to the perspective. The line and the circle intersect at two points.

Figure 1.8: A 3D diagram showing a line L and a circle C on a plane. The line L is labeled L = X_1 \wedge X_2 \wedge n_\infty and the circle C is labeled C = X_1 \wedge X_2 \wedge X_3.

Fig. 1.8 Lines and circles in Conformal Geometric Algebra are trivectors of the type \(X_1 \wedge X_2 \wedge X_3\) . A line can be interpreted as a circle with infinite radius, i.e. containing the point at infinity \(X_3 = n_\infty\) .

Spheres

Given any 4 points with 5D representation \(X_1, X_2, X_3, X_4\) the sphere through those points is given by the 4-vector \(\Sigma = X_1 \wedge X_2 \wedge X_3 \wedge X_4\) . If \(X\) lies on \(\Sigma\) , then \(X \wedge \Sigma = 0\) , that is equivalent to \(X \cdot (\Sigma I) = 0 \rightarrow X \cdot \Sigma^* = 0\) , where \(\Sigma^*\) is the dual of \(\Sigma\) and hence is a vector. Again, since a plane can be represented as a 4-vector of the form \(\Pi = X_1 \wedge X_2 \wedge X_3 \wedge n_\infty\) , it can be seen as a sphere passing through a point at infinity. Examples of planes and spheres in CGA are shown in Fig. 1.9. A summary of CGA objects as blades is given in Table 1.2.

A 3D perspective view showing a sphere, labeled \(\Sigma = X_1 \wedge X_2 \wedge X_3 \wedge X_4\) , which is a semi-transparent red sphere. Below it is a horizontal plane, labeled \(\Pi = X_1 \wedge X_2 \wedge X_3 \wedge n_\infty\) , which is a semi-transparent grey plane. The sphere appears to be resting on or intersecting the plane.

Figure 1.9: A 3D diagram showing a sphere Sigma and a plane Pi. The sphere Sigma is labeled Sigma = X_1 \wedge X_2 \wedge X_3 \wedge X_4 and the plane Pi is labeled Pi = X_1 \wedge X_2 \wedge X_3 \wedge n_\infty.

Fig. 1.9 Planes and spheres in Conformal Geometric Algebra are quadrivectors of the type \(X_1 \wedge X_2 \wedge X_3 \wedge X_4\) . A plane can be interpreted as a sphere with infinite radius, i.e. containing the point at infinity \(X_4 = n_\infty\) .

Table 1.2 Objects in CGA

Intersections

The intersection of any two objects \(W_r, W_s\) is equivalent to taking the dual of the \(2n - r - s\) grade part of the product of \(W_r\) and \(W_s\) . For this purpose, we can define the meet operator as

\[W_r \vee W_s = (W_r^* \wedge W_s^*)^* \quad (1.31)\]

Intuitively, the intersection of a sphere and a circle results in a pair of points, while the intersection of a sphere and a plane yields a circle, provided that the join of the two entities encompasses the entire space. Examples of the use of the meet operator are given in Fig. 1.10-1.11

The figure consists of three separate 3D perspective drawings. The left drawing shows two spheres, one light blue and one light green, overlapping. Their intersection is highlighted as a vertical circle. Above it is the label \(\Sigma_1 \vee \Sigma_2 = \mathcal{C}_1\) . The middle drawing shows a light green sphere resting on a horizontal yellow plane. The intersection is a circle on the sphere's surface. Above it is the label \(\Sigma_1 \vee \Pi = \mathcal{C}_2\) . The right drawing shows two planes, one light blue and one light green, intersecting along a single vertical line. Below the drawing is the label \(\Pi_1 \vee \Pi_2 = L\) .

Figure 1.10: Three diagrams illustrating the meet operator between two grade-4 objects. The first diagram shows two spheres, Σ1 and Σ2, intersecting to form a circle C1. The second diagram shows a sphere Σ1 and a plane Π intersecting to form a circle C2. The third diagram shows two planes, Π1 and Π2, intersecting to form a line L.

Fig. 1.10 The meet operator between two grade-4 objects (two spheres, two planes, or a sphere and a plane), yields a grade-3 object (circles or lines) unless degenerate.

The figure consists of three separate diagrams, each illustrating a geometric meet operation in 3D space:

• Left diagram:** A sphere (grade-4 object) and a circle (grade-3 object) intersect. The circle passes through the sphere, and their intersection is a point pair, labeled A and B. Above the diagram is the equation \(\Sigma \vee \mathcal{C} = A \wedge B\) .
• Middle diagram:** A sphere (grade-4 object) and a line (grade-3 object) intersect. The line passes through the sphere, and their intersection is a point pair, labeled C and D. Below the diagram is the equation \(\Sigma \vee L = C \wedge D\) .
• Right diagram:** A plane (grade-4 object) and a circle (grade-3 object) intersect. The circle lies in the plane, and their intersection is a point pair, labeled E and F. Above the diagram is the equation \(\Pi \vee \mathcal{C} = E \wedge F\) .

Figure 1.11: Three diagrams illustrating the meet operator between grade-4 objects (sphere or plane) and a grade-3 object (line or circle).

Fig. 1.11 The meet operator between a grade-4 objects (sphere or plane) and a grade-3 object (line or circle), yields a grade-2 object (a point pair) unless degenerate.

1.2.8 Other algebras

1d-Up CGA

When dealing with transformations in Euclidean geometry in CGA, the point at infinity \(n_\infty\) is kept constant. However, we can work with non-Euclidean geometries if we keep different quantities constant. For example by keeping \(e\) constant, it can be shown that we are left with a hyperbolic geometry. Similarly, when \(\bar{e}\) is kept constant we are left with a spherical geometry. Note how, by keeping either one of the bases \(e\) or \(\bar{e}\) constant, we have only one additional basis vector compared to 3D GA (hence the name “1D-Up”) instead of two as in CGA (which is a “2D-Up” space compared to 3D GA) [95–97].

We will work with the latter case in Chapter 4, in which \(\bar{e}\) is kept constant and \(e\) is our origin. The main advantages of the 1D-Up approach are: (1) a lower dimensionality of the space compared to CGA, (2) that the Euclidean nature of the space (i.e. that all basis vectors square to \(+1\) ) allows us to construct a simple loss function that is invariant under rigid body transformations and (3) that both translations and rotations can be expressed as rotors, just as in the 2d-Up space. The curvature of the space is controlled by a free parameter \(\lambda\) (see Fig. 1.12), which we will discuss in further detail in Chapter 4.

Spacetime Algebra (STA)

The 3D Spacetime Algebra (STA), i.e. \(\mathcal{G}_{1,3,0}\) , is spanned by four basis vectors \(\{\gamma_\mu\}\) , with \(\gamma_0^2 = -\gamma_k^2 = 1\) and \(k = 1, 2, 3\) . The pseudoscalar is defined as \(I = \gamma_0\gamma_1\gamma_2\gamma_3\) . We

Figure 1.12: Two diagrams illustrating the effect of the parameter lambda on geometry in 1d-Up CGA. The left diagram shows a 3D coordinate system with a plane and a line. As lambda increases, the line transforms into a circle and the plane into a sphere. The right diagram shows a 3D coordinate system with a plane and a sphere. As lambda increases, the sphere transforms into a plane and the plane into a sphere.

Fig. 1.12 The parameter \(\lambda\) controls the curvature of the space in 1d-Up CGA. As \(\lambda\) increases, lines in \(G(4, 0, 0)\) are transformed into circles and planes are transformed into spheres.

say that \(\gamma_0\) is a timelike vector, while \(\gamma_1, \gamma_2, \gamma_3\) are spacelike vectors. A breakdown of the blades in STA is given in Table 1.3.

The even grade elements of the 3D STA form a subalgebra \(\mathcal{G}_{3,0,0}\) which is equivalent to the Pauli algebra. In fact, defining:

\[\sigma_k = \gamma_k \wedge \gamma_0 \quad (1.32)\]

we have that \(\sigma_k^2 = 1\) and that \(\sigma_1 \wedge \sigma_2 = I\sigma_3\) , \(\sigma_2 \wedge \sigma_3 = I\sigma_1\) and \(\sigma_3 \wedge \sigma_1 = I\sigma_2\) . We will work with STA in Chapter 5.

Table 1.3 Blades in Spacetime Algebra (STA) \(\mathcal{G}_{1,3,0}\) . Elements colored in red square to \(+1\) , while elements colored in blue square to \(-1\) .

1.3 Learning with GA

The idea of using GA as a computational substrate for machine learning dates back to early works such as [135], and was further developed in the early 2000s through the contributions of [22, 23, 21, 12]. These studies introduced the notion of Clifford neurons, exploring the use of GA to design learning units and architectures that operate directly on multivectors. While visionary, these efforts remained mostly theoretical and were largely overlooked for over a decade, as the machine learning community increasingly consolidated around real-valued vector spaces and conventional deep learning frameworks.

Interest in GA for machine learning began to re-emerge in the 2020s. This revival was catalyzed by the broader development of geometric deep learning, a subfield of deep learning focused on building architectures that respect the symmetries and structure of geometric data. It leverages tools from group theory and differential geometry to enforce invariance and equivariance on data such as graphs, manifolds, and point clouds.

A notable contribution to this renewed interest came from Melnyk et al. [128], who introduced the Multilayer Geometric Perceptron (MLGP), an architecture composed of geometric neurons that achieves rotation- and translation-equivariant shape classification. This work represented one of the first functioning implementations of a neural network operating directly within a GA framework.

Following this, Brandstetter et al. [17] showed that Clifford neural networks can outperform standard real-valued architectures on partial differential equation (PDE) solving tasks. Their success was attributed to the natural ability of GA to represent vector fields and encode geometric constraints. In parallel, Ruhe et al. [158, 159, 19] extended these foundations by constructing architectures that are equivariant under the full Clifford group and introduced the idea of learning geometric templates and primitives. These developments reinforced the role of GA as a powerful tool for geometric computation in machine learning.

How does this thesis relate to the broader field of geometric deep learning?

Objects and operations introduced in Section 1.2 serve as the mathematical backbone throughout this thesis, as we will exclusively model problems in Bioinformatics (Chapter 3), Computer Vision (Chapter 4), and Physics Modelling (Chapter 5) in terms of multivectors. We will employ multivectors in two ways (see Fig. 1.13):

1. Modelling: we simplify, model, and represent data in terms of GA objects and operations, including lines, planes, circles and rotors.
2. Computing: we design learning units capable of processing multivector-valued data through specialized, GA-native operations, including geometric and sandwich products.

The diagram illustrates the relationship between Machine Learning with Geometric Algebra and Geometric Deep Learning. It is divided into two main sections: 'Machine Learning with Geometric Algebra' (left, yellow box) and 'Geometric Deep Learning' (right, blue box). The 'Machine Learning' section is further divided into 'Modelling' (Chapters 2, 3) and 'Computing' (Chapters 3, 4, 5). 'Modelling' includes molecular geometry optimisation, learning rotation, and cost maps, leading to CGAPoseNet. 'Computing' includes STAResNet, Fengbo, GA-ReLU, CGAPoseNet+GCAN, Define, Refine, Align, and CGENN 3D Projector.

Diagram illustrating the relationship between Machine Learning with Geometric Algebra and Geometric Deep Learning. The diagram is split into two main sections: 'Machine Learning with Geometric Algebra' (left, yellow box) and 'Geometric Deep Learning' (right, blue box). The 'Machine Learning' section is further divided into 'Modelling' (Chapters 2, 3) and 'Computing' (Chapters 3, 4, 5). 'Modelling' includes molecular geometry optimisation, learning rotation, and cost maps, leading to CGAPoseNet. 'Computing' includes STAResNet, Fengbo, GA-ReLU, CGAPoseNet+GCAN, Define, Refine, Align, and CGENN 3D Projector.

Fig. 1.13 Multivectors for Modelling, Understanding, and Computing. This thesis aims to demonstrate the power of GA across the entire machine learning pipeline. Early works, such as those on rotation estimation and molecular geometry optimization (Chapter 2), use multivectors as a modelling tool. Later works like CGAPoseNet+GCAN (Chapter 4) and STAResNet (Chapter 5) exemplify true geometric deep learning, where multivectors serve not only as a modelling tool, but also as computational units within neural architectures. Both approaches serve the goal of understanding how geometric structure informs and constrains the learning process.

Earlier chapters (Chapter 2 and the first part of Chapter 3) concentrate on non-geometric settings: GA is used primarily as a modelling and representational framework, rather than as a means of enforcing geometric constraints or equivariances within architectures. Later chapters (the second half of Chapter 3, and Chapters 4 and 5) explore equivariant and geometric architectures relying on multivector-valued learnable parameters. This second half of the manuscript fits within geometric deep learning.

This division reflects both the natural progression of the PhD research and an intentional choice to highlight the broader capabilities of GA, as captured in the subtitle of the thesis: Multivectors for Modelling, Understanding, and Computing.

1.3.1 GA: the 6^th “G” of Geometric Deep Learning?

A popular taxonomy of the field of geometric deep learning is the one introduced in [20]. This framework identifies five key domains, referred to as the 5 G’s, that capture the structural priors and symmetries commonly found in deep learning problems:

• Grids: Data defined on regular Euclidean domains, such as images, audio signals, or video frames. Convolutional neural networks (CNNs) are the canonical models here, exploiting the translation symmetry of the grid [104].
• Groups: Symmetries formalized via group theory, such as rotations, reflections, or permutations. Learning architectures that are invariant or equivariant under group actions fall into this category, improving generalization by embedding prior symmetry knowledge [34, 35, 40].
• Graphs: Data with relational structure, where nodes and edges represent entities and their interactions. Graph neural networks (GNNs) operate in this regime, enabling learning on irregular, non-Euclidean domains such as social networks or molecular graphs [179, 190, 178].
• Geodesics: Data that lies on smooth manifolds, such as spherical signals or curved surfaces. Models must account for curvature and intrinsic geometry, often using tools from differential geometry to define distances and gradients [125].
• Gauges: The most abstract setting, involving local symmetries and transformations in fiber bundles. Gauge-equivariant architectures are particularly relevant in physics-informed learning and scenarios where data is expressed in locally varying frames [39].

How does GA fit within this taxonomy? GA provides a unified algebraic framework that naturally encompasses each of these domains:

• Grids: GA offers coordinate-free representations for regular grid structures. In 2D or 3D grids (e.g., images or volumetric data), multivectors can efficiently encode geometric primitives such as edges, planes, and volumes. The geometric product captures local spatial relationships naturally, serving as a generalization of convolutional operations.
• Groups: GA excels in representing symmetry groups central to geometric deep learning: Rotations via rotors in the even subalgebra; Translations and scalings

via versors in CGA; Conformal transformations unified in the algebra \(\mathcal{G}_{4,1}\) ; Lorentz transformations for spacetime models in physics. GA makes group actions both mathematically elegant and computationally tractable.

• Graphs: GA can enrich graph neural networks by using multivectors to represent edges and node features. Unlike scalar weights, multivector-valued edge attributes can encode directionality, orientation, and higher-order geometric relationships (e.g., areas and volumes), providing a more expressive relational framework.
• Geodesics: GA is particularly powerful for reasoning on manifolds. In CGA, geodesics (e.g., great circles on spheres or hyperbolic lines) can be expressed algebraically as intersections of geometric primitives. This facilitates learning tasks that depend on intrinsic distances, such as manifold embeddings or molecular conformer prediction.
• Gauges: GA naturally handles basis changes and frame transformations. Its invariance under orthogonal transformations enables the construction of gauge-equivariant networks, useful for modelling systems with local symmetries, as seen in physics-informed and gauge-theoretic machine learning.

We therefore argue that GA functions as a meta-framework: whereas each of the 5 G's typically requires distinct mathematical tools, such as tensors, graphs, Lie groups, or Riemannian manifolds, GA unifies them within a single consistent formalism. This enables end-to-end learning systems in which representation, transformation, and reasoning are all performed natively in the same geometric space.

In summary, GA is not limited to any single domain among the 5 G's. Rather, it subsumes all five within a compact and expressive language, making it a compelling foundation for modelling, understanding, and computing in geometric deep learning (see Table 1.4).

1.3.2 Learning with GA in practice

In Table 1.4, we used the expression “embedding in a GA”. In this section, we clarify what that means in practice. Working with networks in GA means extending traditional models to operate directly on multivector-valued data, enabling more expressive representations of geometry, symmetry, and structure. This approach can be realized through the following three-step workflow:

Table 1.4 Embedding an architecture in a GA \(\mathcal{G}\) enables geometric deep learning on any of the 5 G's of [20]. Left column: architectures on which geometric deep learning is performed; right column: example of possible embedding of such architectures in GA. In this thesis, we provide examples of GA-based learning for grids (Chapter 5), groups (Chapter 3-4), geodesics (Chapter 4).

1. Choose a suitable GA: Select an algebra \(\mathcal{G}_{p,q,r}\) that reflects the geometry of the problem domain. For instance, \(\mathcal{G}_{2,0}\) is appropriate for modelling 2D physical systems such as 2D Navier–Stokes equations, while \(\mathcal{G}_{4,1}\) , a conformal model, is well-suited for tasks involving 3D rigid-body motions.
2. Embed your model in the algebra: Lift all data and learnable parameters into the chosen GA by introducing an additional tensor dimension of size \(2^n\) , where \(n = p + q + r\) is the number of basis vectors. Each element along this new axis corresponds to a basis blade (e.g., scalar, vector, bivector), allowing the network to represent and manipulate the full geometric content of the input.
3. Build your architecture: Adapt standard neural architectures (e.g., CNNs, GNNs, transformers) to perform GA-native operations, such as the geometric product, outer product, and grade projections. These operations replace traditional scalar or matrix multiplications, enabling the model to reason over geometric structure in a principled way.

Note how this workflow is model-agnostic: it can be applied to a wide range of architectures, systematically incorporating geometric structure and symmetry through GA. The rest of this thesis assumes this workflow for building GA networks.

Chapter 2

Rotations as Rotors

There are no answers, only choices.

— Solaris, Stanislaw Lem (1961)

2.1 Introduction

In this chapter we focus on the theme of representation. What defines a good choice when representing rotations? Are rotors a “good” mathematical representation for rotations? Do rotors lend themselves to learning problems? Is GA an optimal space in which to minimise a given objective function? We address these questions through two examples that will lay the foundations for the following chapters and justify our use of rotors in them.

Firstly, we discuss the issue of discontinuities in the representation of rotations, how it affects learning and how rotors (and bivectors) can overcome it in a more compact way compared to other methods in the literature through three toy problems in computer vision.

We then move onto the problem of molecular geometry optimisation in computational chemistry, and show how rotors can more effectively guide an unstable molecular conformation to its energetic minimum as opposed to other representations.

2.2 Learning rotations

2.2.1 Background

In problems such as pose prediction and generation, motion capture and inverse kinematics, neural networks (NNs) are generally trained to perform a regression on rotations [150, 126]. For these tasks, it is common to describe 3D rotations through 3 or 4D representations, including axis-angle representation [79], Euler angles [49] and quaternions [134].

However, it has been shown that 3 and 4D representations present some limitations. Grassia et al. [58], for example, demonstrated how an exponential parametrisation of rotations is more suitable for differentiation and integration as opposed to Euler angles or quaternions. Zhou et al. [201] found the source of large regression errors in the discontinuity of these representations, in which a representation is said to be discontinuous when the mapping \(g : SO(3) \rightarrow \mathbb{R}^D\) from the \(3 \times 3\) rotation matrix in \(SO(3)\) onto a given representation space is a discontinuous function. Saxena et al. [161] also highlighted representation discontinuity as a critical factor.

According to the universal approximation theorem, on the other hand, a NN with one hidden layer can approximate any continuous function for inputs that fall within a specific range [76]. Moreover, the process of training requires extensive computations of derivatives. Hence, trying to learn rotations when parametrised with a discontinuous representation might lead to large errors.

2.2.2 The 6D representation

A possible solution to overcome the discontinuity issue is to increase the number of degrees of freedom (DoF). In [201], 5 and 6D continuous representations have been proposed specifically for deep learning approaches to computer vision in order to represent rotations exclusively (i.e. without translation component). This 6D representation can be visualized as a \(3 \times 2\) matrix with columns \(a_1, a_2\) , which, after a Gram-Schmidt-like orthogonalization procedure, maps to the original rotation matrix in \(SO(3)\) with orthonormal columns \(b_1, b_2, b_3\) (Eq. 2.1).

\[g^{-1} \left( \begin{bmatrix} | & | \\ a_1 & a_2 \\ | & | \end{bmatrix} \right) \mapsto \begin{bmatrix} | & | & | \\ b_1 & b_2 & b_3 \\ | & | & | \end{bmatrix} \in SO(3) \quad (2.1)\]

in which \(g^{-1} : \mathbb{R}^D \rightarrow SO(3)\) is the mapping from the representation space to the the \(3 \times 3\) rotation matrix in \(SO(3)\) . We believe that this representation has some fundamental limitations:

• The 6D representation proposed has more DoF than commonly employed representations, meaning more parameters to be learned.
• It is fully learned. This guarantees continuity, but it also means that no closed form expression exists to convert a rotation from \(SO(3)\) into its 6D representation.
• It is tied to an ad hoc loss function to be minimised, namely the \(L_2\) norm between the original and predicted rotation matrices after the aforementioned orthogonalization procedure, that has to be performed at every iteration. This loss function is less intuitive and more computationally expensive than those normally used in regression tasks.
• No conditions are imposed to avoid \(a_1 = b_1; a_2 = b_2\) i.e. that the 6D representation is learned to be equal to the first two columns of the rotation matrix. This identity mapping is easier to "learn" and indeed a possible scenario, but makes the comparison with other representations biased.

We address these points by reproducing the experiments in [201] employing a full GA formulation.

2.2.3 Learning mappings between representation spaces

The first experiment, which we will refer to as ‘sanity check’, is to use a multi-layer perceptron (MLP) to learn the mapping from rotation matrix \(M \in SO(3)\) to a given representation \(\mathcal{R} \in \mathbb{R}^D\) , in which \(D\) is the number of DoF associated with \(\mathcal{R}\) , or its dimensionality. Then, the representation \(\mathcal{R}\) is converted back to an estimate \(M'\) and the closeness between \(M'\) and \(M\) is measured. We first tested Zhou’s 6D representation along with quaternions, axis-angle, Euler angles and bivectors.

The network has four dense layers with 128 neurons each and a Leaky ReLU activation function between them. The optimiser has been kept as Adam, the learning rate to the default value of 0.001, the number of epochs to \(E = 100\) and the batch size as 64. As a loss function, we adopted either the mean squared error (MSE) or the mean absolute error (MAE) between the target and input representation, defined as:

\[MSE = \frac{1}{D} \sum_{i=1}^D (y_i - y'_i)^2 \quad (2.2)\]

\[MAE = \frac{1}{D} \sum_{i=1}^D \|y_i - y'_i\| \quad (2.3)\]

where \(D\) , \(y\) and \(y'\) are the representation dimensionality, original and predicted representations, respectively. In the case in which \(y\) is a bivector, for example, then the \(y_i\) 's will be the coefficients \(b_{12}, b_{23}, b_{13}\) of the bivectors in \(\mathcal{G}_{3,0,0}\) ; if \(y\) are Euler angles, then the \(y_i\) 's will correspond to the angles \(\alpha, \beta, \gamma\) . We adopted the MSE as it is simple and commonly used in regression tasks. Only for the 6D representation did we employ the \(L_2\) distance between rotation matrix \(M\) and predicted matrix \(M'\) , defined as:

\[L_2 = \sqrt{\sum_{i=1}^3 \sum_{j=1}^3 (y_{ij} - y'_{ij})^2} \quad (2.4)\]

where \(y \equiv M\) , \(y' \equiv M'\) are the original and predicted rotation matrices, respectively. Note that, according to its definition, the 6D representation can only have the \(L_2\) norm between \(M, M'\) as a loss. This is due to two reasons:

• The main feature of the 6D representation, i.e. that the rotation matrix is obtained through orthogonalization of the 6D representation, is implemented within the loss function itself.
• There is no closed form expression to represent a rotation matrix as its 6D counterpart (in other words, the mapping \(g : SO(3) \rightarrow \mathbb{R}^6\) is fully learned). This means that the ground truth, when learning the 6D representation, cannot be expressed in 6D, but only as the original rotation matrix.

It is worth mentioning that in [201] all the experiments share the \(L_2\) norm as a loss function. We believe that this is not only a more complex approach, as the mapping \(g^{-1} : \mathbb{R}^D \rightarrow SO(3)\) has to be evaluated at each iteration, but it is also penalizing discontinuous representations (i.e. those representations for which \(g\) is discontinuous): we want to demonstrate that these representations can be learned successfully when this mapping is omitted from the learning strategy.

The results from [201] have been reproduced and validated on a set of \(10^5\) rotation matrices with a 67-33 train-test split and 30% of the training set employed for validation. The maximum, mean and standard deviation of the geodesic error measured on the testing set are shown in Table 2.1, along with the loss function employed during training and the kind of error employed. The geodesic error is defined as the minimal angular distance between two rotations, where the concept of geodesic is defined in metric geometry as the curve of shortest distance between two points on a manifold. While

different formulations exist in the literature, we will stick to the definition given in [201], expressed as

\[\mathcal{L} = \arccos\left(\frac{\text{tr}(M'') - 1}{2}\right) \quad (2.5)\]

Where \(M'' = MM^T\) . Table 2.1 shows that the 6D representation indeed overcomes the discontinuity issue and outperforms significantly the other representations, including the bivector representation that we tested in addition to the others. Note how, when using MSE between quaternions, the matrix-to-quaternion/rotor learning does not converge to a minimum, most likely due to the antipodal symmetry of the representation, meaning that quaternions \(q\) and \(-q\) represent the same rotation [109]. Our results, however, are in agreement with Zhou’s when employing the same \(L_2\) loss.

Table 2.1 Average Geodesic Error ( \(^\circ\) ), Sanity Check, \(M \mapsto \mathcal{R}\)

GA-based Sanity Check

We believe that the errors shown in Table 2.1 are specific to the kind of mapping learned rather than intrinsic to the representation itself. The maximum \(180^\circ\) error, said to come from the discontinuity of the representation, in fact, is present only when adopting Eq.2.5 as a metric and when the representation is learned starting from the rotation matrix. We now proceed to show that the maximum geodesic error is much smaller when keeping the rotations in GA.

We trained the same network to learn the bivector \(B\) from the rotor \(R\) , either as its Cayley transform or as its logarithm, and then converted back to estimated rotor \(R'\) , where \(R' = e^{-B/2}\) or \(R' = (1 - B)/(1 + B)\) . We measured the regression error via Eq. 2.5 after converting rotors to their corresponding matrix representation.

Results are shown in Table 2.2. Note how comparable results between the 6D representation and the bivector are obtained when expressing the bivector as the Cayley transform of the rotor, with a mean geodesic error of \(1.04^\circ\) for the rotor-to-bivector mapping and of \(1.06^\circ\) for the matrix-to-6D mapping. The bivector, however, only requires 3 DoF to be learned instead of 6. Note also how the same bivector yields a mean geodesic error about three times smaller when learned from the rotor (in Table 2.2) than when learned from the rotation matrix (in Table 2.1)

Table 2.2 Geodesic Error ( \(^\circ\) ), Sanity Check, \(R \mapsto B\)

2.2.4 Pose estimation of 3D point clouds

For this task we trained a network that takes as input a reference and a target point cloud \(P_r, P_t \in \mathbb{R}^{N \times 3}\) respectively, where \(N\) is the number of points per point cloud, to predict a representation \(\mathcal{R} \in \mathbb{R}^D\) of the rotation between the two, where \(D\) is the dimensionality of the representation. These include the rotation matrix ( \(D = 9\) ), the 6D representation ( \(D = 6\) ), Euler angles ( \(D = 3\) ), the axis-angle representation ( \(D = 3\) ), the quaternion/rotor ( \(D = 4\) ) and the bivector ( \(D = 3\) ).

The network is a slightly modified version of PointNet [152], see Figure 2.1. We assume that the registration between the point clouds is known, i.e. we know which points in \(P_r\) correspond to which in \(P_t\) . The input of this network is then the concatenated point clouds \([P_r, P_t]\) of size \(N \times 6\) , passed through five layers of 1D convolutions of size \(6 \times 64 \times 128 \times 256 \times 1024\) , and then through three more dense layers of size \(512 \times 512 \times D\) , where \(D\) is the dimensionality of \(\mathcal{R}\) .

The optimiser has been kept as Adam, the learning rate to the default value of 0.001, the batch size as 32 and the loss to MAE (for this specific task, we empirically verified that it performed better than MSE). A Leaky ReLU activation function has been employed between layers, with the exception of the last layer that uses a linear activation function. The dataset is composed of 726 airplane point clouds taken from the ModelNet10 dataset [191], with 626 training and 100 testing samples. Of the 626 training samples, 30% was employed for validation. Each point cloud is composed of \(N\)

The diagram illustrates the network architecture for pose estimation. The input is the concatenation of reference point cloud \(P_r\) and target point cloud \(P_t\) , each represented by three coordinates \(x_{r1}, x_{r2}, x_{r3}\) and \(x_{t1}, x_{t2}, x_{t3}\) . These are processed through a series of 1D convolutional layers: 1D Conv 6, 1D Conv 64, 1D Conv 128, 1D Conv 256, and 1D Conv 1024. The output of the final convolutional layer is then processed by two dense layers, Dense 512 and Dense 512, leading to an output layer that produces three parameters \(y_1, y_2, y_3\) .

Diagram of a neural network architecture for pose estimation from 3D point clouds. The input consists of concatenated reference (P_r) and target (P_t) point clouds, shown as two 3D plots of an airplane. These are processed through five 1D convolutional layers with 6, 64, 128, 256, and 1024 filters respectively. This is followed by two dense layers with 512 units each, and finally an output layer that produces three parameters y1, y2, and y3.

Fig. 2.1 The network for the pose estimation from 3D point clouds problem. Input: concatenated reference and target point clouds \(P_r, P_t \in \mathbb{R}^{N \times 3}\) . Output: rotation from \(P_r\) to \(P_t\) with \(D\) parameters.

= 3000 points. To each of the 726 point clouds is associated a random rotation matrix from which we obtained the 8 different representations tested.

Results are summarized in Table 2.3, ordered by increasing mean geodesic error along with the number of DoF of the representation. The column "Epochs" specifies the number of training epochs for each representation. Firstly, it can be noticed that using MAE instead of the \(L_2\) norm between \(M\) and \(M'\) , yields a lower error for the matrix representation compared to the 6D representation, as opposed to findings in [201]. Secondly, that the \(180^\circ\) maximum error is present only in some of the measurements, and not consistently, for all the representations analyzed. This leads to a maximum geodesic error averaged over five experiments always below \(180^\circ\) , also for those representations defined as discontinuous by Zhou et al. (e.g. axis-angle, Euler angles, etc.), for which we would expect an error around \(180^\circ\) . We believe that the difference between our results and Zhou's results rely in the different learning strategy and mapping learned, not in the type of representation.

The distribution of the geodesic error over the testing set along with its percentile plot are given in Fig. 2.2, confirming the superiority of our GA approach. Note how the 75th percentile is achieved at \(15^\circ\) for the bivector (expressed as the logarithm of the rotor) compared to \(21^\circ\) and \(25^\circ\) for the 6D and the matrix representations, respectively.

These results show that a pure geometric algebra description with 3 DoF is as good as a representation with 6 DoF, with a simpler loss function employed and an

Table 2.3 Geodesic Error ( \(^\circ\) ), Pose Estimation from 3D Point Clouds.

The figure consists of two side-by-side plots. The left plot is a histogram showing the 'Normalized Count' of geodesic error (in degrees) for eight different rotation representations. The x-axis ranges from 0 to 100 degrees, and the y-axis ranges from 0.00 to 0.05. The Matrix representation (green line) has a very sharp, high peak at approximately 10 degrees. The 6D (cyan), Rotor (blue), Quaternion (red), Bivector (Exp) (orange), Bivector (Cayley) (purple), Euler (black), and Axis Angle (yellow) representations show broader, lower peaks. The right plot is a cumulative distribution function (CDF) showing the 'Probability' of the geodesic error (in degrees). The x-axis ranges from 0 to 100 degrees, and the y-axis ranges from 0.0 to 1.0. The Matrix representation (green line) reaches a probability of 1.0 very quickly, around 20 degrees. The other representations reach 1.0 more gradually, with the Axis Angle (yellow) and Euler (black) representations being the slowest to reach 1.0.

Figure 2.2: Distribution of the geodesic error over the testing set. The left plot shows the Normalized Count of error (0 to 100 degrees) for various representations. The right plot shows the Probability (0.0 to 1.0) of error (0 to 100 degrees). Both plots compare Matrix, 6D, Rotor, Quaternion, Bivector (Exp), Bivector (Cayley), Euler, and Axis Angle representations. The Matrix representation shows a sharp peak at low error, while the Axis Angle representation shows a broader distribution.

Fig. 2.2 Distribution of the geodesic error over the testing set for the proposed representations (left) and related percentile plot (right) for the pose estimation problem.

equivalent geodesic error formulation measured exclusively through rotors, without the need to employ rotation matrices.

Figure 2.3: Two line graphs comparing Bivector and 6D representations. The left graph shows Mean Geodesic Error (°) vs Noise Standard Deviation (0 to 5). The right graph shows Mean Euclidean Distance vs Noise Standard Deviation (0 to 5). In both cases, the Bivector representation (red line) consistently shows lower error and distance than the 6D representation (blue line) as noise increases.

Fig. 2.3 Average geodesic error (left) and Euclidean distance (right) as a function of the standard deviation of the noise applied on \(P_r, P_t\)

Noisy Point Clouds

We then studied the case in which Gaussian noise is added to the point clouds. It is safe to assume imperfect acquisition of samples from a scene (e.g., due to a shaky camera) in a realistic scenario. We restricted the analysis to the two representations with lowest error in the noiseless case, namely the (exponential) bivector and the 6D. We added noise to the reference and target point clouds \(P_r, P_t\) , respectively. Noise has been modelled as additive white Gaussian noise (AWGN) \(\sim \mathcal{N}(0, \sigma^2)\) with variable \(\sigma = \{0.1; 1; 2; 5\}\) . The geodesic error and Euclidean distance measured for the two representations on the testing set are given in Figure 2.3.

It can be seen how, in high noise scenarios, the network likely fits to the noise when trying to predict more parameters of the 6D representation instead of the 3 DoF of the bivector. The average geodesic errors are noticeably higher for the 6D representation than for the bivector case as noise is increased. This is also mirrored in the Euclidean distance between original and predicted point clouds.

The 6D representation is very accurate with little-to-no noise, but it fails when the standard deviation of the noise increases. The bivector representation, on the other hand, is more tolerant to high noise conditions, showing a consistently lower geodesic error at a fraction of the required DoF and with a simpler loss function employed during training.

2.2.5 Inverse kinematics

The goal of this experiment is to estimate the rotation from the T-pose to any arbitrary pose of a skeleton. By T-pose we define the skeleton in standing position and arms stretched out that is commonly used in computer graphics for calibration purposes. The input of the network takes the 3D positions of the \(N = 31\) joints of a skeleton as \(P = (p_1, p_2, p_3, \dots, p_N)\) , where \(p_i = (x_i, y_i, z_i)^T\) . The outputs of the network are the rotations that a T-pose skeleton undergoes to reach the input position. The output is expressed as \(R = (r_1, r_2, r_3, \dots, r_N)\) , with \(r_i \in \mathbb{R}^D\) , where \(D\) is the dimension of the representation. The rotations are hierarchical, in the sense that the rotation of each joint depends on the previous one and on how they are linked to each other. The adopted NN is a four-layer MLP with 1024 neurons per layer (see Figure 2.4). The problem is formulated as a supervised learning problem, that allows us to use MSE loss consistently with the only exception being the 6D representation.

We used 760 clips from the CMU Motion Capture Database (MoCap) [72]. From these clips we selected 10000 frames, with a 67-33 train-test split. 30% of the training set was employed for validation. We centered the position of the root joint (defined as the midpoint between the two hips) at the origin so as to predict solely the rotational component of the joints and not the translation of the skeleton. The network was trained with batch size 64. The average geodesic error between ground truth (from the testing set) and predicted rotation and Euclidean distance of the predicted poses from the ground truth are displayed in ascending order in Tables 2.4 - 2.5, respectively.

Diagram of the neural network architecture for inverse kinematics. On the left, a 3D plot shows a red stick figure (skeleton) in a T-pose. An arrow points from this figure to a neural network diagram. The network consists of an input layer with three nodes labeled x1, x2, and x3. This is followed by four 'Dense' layers, each containing 1024 nodes (represented by blue circles). The final output layer contains three nodes labeled y1, y2, and y3. All nodes in adjacent layers are fully connected.

Fig. 2.4 The network for the inverse kinematic problem. Input: frames with \(31 \times 3\) spatial coordinates. Output: 31 rotations with \(D\) parameters.

Figure 2.5 consists of two plots. The left plot is a normalized count histogram of geodesic error in degrees, ranging from 0 to 160. The y-axis is 'Normalized Count' from 0.000 to 0.014. The right plot is a cumulative probability plot of geodesic error in degrees, ranging from 0 to 200. The y-axis is 'Probability' from 0.0 to 1.0. Both plots compare six representations: Matrix (green), 6D (cyan), Quaternion (yellow), Bivector (Exp) (blue), Bivector (Cayley) (red), Euler (magenta), and Axis Angle (black). In the histogram, Bivector (Exp) has the highest peak at approximately 60 degrees, while Matrix has the lowest. In the cumulative plot, all representations show similar performance, reaching a probability of 1.0 around 150 degrees.

Figure 2.5: Distribution of the geodesic error over the testing set for the proposed representations (left) and related percentile plot (right) for the inverse kinematic problem.

Fig. 2.5 Distribution of the geodesic error over the testing set for the proposed representations (left) and related percentile plot (right) for the inverse kinematic problem.

Figure 2.6 consists of two plots. The left plot is a normalized count histogram of Euclidean distance in cm, ranging from 2 to 10. The y-axis is 'Normalized Count' from 0.0 to 0.7. The right plot is a cumulative probability plot of Euclidean distance in cm, ranging from 2 to 10. The y-axis is 'Probability' from 0.0 to 1.0. Both plots compare the same six representations as in Figure 2.5. In the histogram, the Matrix representation (green) has the highest peak at approximately 4.5 cm, while Bivector (Exp) (blue) has the lowest. In the cumulative plot, Matrix reaches a probability of 1.0 around 6 cm, while Bivector (Exp) reaches it around 8 cm.

Figure 2.6: Distribution of the Euclidean distance between predicted pose and T-pose over the testing set for the proposed representations (left) and related percentile plot (right), inverse kinematic problem.

Fig. 2.6 Distribution of the Euclidean distance between predicted pose and T-pose over the testing set for the proposed representations (left) and related percentile plot (right), inverse kinematic problem.

Table 2.4 Geodesic Error ( \(^\circ\) ), Inverse Kinematics.

Table 2.5 Euclidean Distance (cm), Inverse Kinematics.

It can be noticed that there is no perfect correspondence between a low geodesic error and an accurate motion estimation. Moreover, the difference in Euclidean distance between different representations is small and at most 1.5 cm. The small difference between the representations is probably due to the supervised formulation of the problem. The bivector expressed as the Cayley transform of the rotor allowed the second best pose estimation after the rotation matrix in terms of mean euclidean distance. The rotation matrix has been proven to be consistently the best representation for all three experiments, as expected due to the higher number of DoF.

The distributions over the testing set of the geodesic error and Euclidean distance along with their percentile plots are given in Fig. 2.5 - 2.6, respectively. Note how in Fig. 2.5 the 75th percentile is achieved at about \(55^\circ\) and \(75^\circ\) for the (Cayley) bivector and the 6D and (log) bivector, respectively. Similarly, in Fig. 2.6, the 75th percentile is at about 4.5 cm for the matrix representation, 5 cm for the (Cayley) bivector and above 6 cm for the 6D representation.

To guarantee fairness of comparison between different representations in the pose estimation and inverse kinematics experiments, we regularized the training of the architectures by implementing early stopping with patience of \(P = 10\) epochs. The validation loss curves for the two experiments are shown in Fig. 2.7. They have been normalized to show convergence rather than their absolute value. For the sanity check, instead, we fixed \(E = 100\) as no overfitting was observed.

Figure 2.7: Two line plots showing Normalized Validation Loss versus Epochs for different molecular representations. The left plot is for the pose estimation problem (Epochs 0-70) and the right plot is for the inverse kinematic problem (Epochs 0-120). Both plots compare eight representations: Matrix (blue), 6D (orange), Quaternion / Rotor (green), Bivector (Cayley) (cyan), Bivector (Exp) (magenta), Axis Angle (yellow), and Euler (black). In both cases, the loss for all representations decreases rapidly and stabilizes at low values, with the Matrix and 6D representations generally showing the lowest final loss.

Fig. 2.7 Validation loss with early stopping of each of the 8 representations in the pose estimation problem (left) and inverse kinematic problem (right).

2.3 Molecular geometry optimisation through rotors

2.3.1 Background

Optimising the geometry of a molecule means finding the arrangement in 3D space of its atoms for which the overall energy of the molecule is minimised. The stability of molecular conformations is particularly relevant when designing drugs and studying their effects: a molecule might exist in nature with multiple geometries, or conformers. Different conformers of a same molecule, for example, might or might not be toxic to humans. Hence the study of molecular conformations is of particular relevance.

We will focus our attention on proteins. A protein can be thought of as a chain of amino acids. Each amino acid is characterized by a side chain, which differentiates them, bonded to a \(C_\alpha\) atom. The \(C_\alpha\) is itself bonded to an \(N\) atom and to a \(C\) atom. The repetition of \(N - C_\alpha - C\) triplets in each amino acid in the protein chain defines the protein backbone.

Generally, the 3D structure of a protein is modelled through dihedral angles \(\phi, \psi, \omega\) (i.e. internal coordinates) for the backbone and \(\chi_i\) for the side chains. Each dihedral angle defines a rotation around a bond. Performing a geometry optimisation by

exploring the space of dihedral angles is time consuming and implies different axes of rotation for different dihedral angles.

An example of this approach can be found in [123], in which the fitness function of the optimisation algorithm is a function of the dihedral angles. This strategy has two drawbacks: the first is that the number of degrees of freedom, even for a single amino acid with few atoms, is large. For a single amino acid such as aspartic acid, for example, six dihedral angles were employed \(\{\phi, \psi, \omega, \chi_1, \chi_2, \chi_3\}\) . In fact, Mancini et al. [123] do not report results on dipeptides (i.e. chains of two amino acids, or peptides). The second is that each dihedral angle specifies a rotation with respect to different axes with different pivot points, making rotations difficult to deal with.

In this Section we implement an evolutionary algorithm (EA) to search for stable conformers of molecules by performing an exploration over the rotor space. We associate an amino acid pair with a rotor that describes their relative orientation. Rotors are then employed as a metaheuristic to guide the minimisation of the overall energy of the molecular structure. The proposal of rotors drives the ‘evolution’ of the molecule towards its more stable conformation. We believe that the mathematical formalism of rotors is particularly suitable to describe molecular structures and in particular proteins, which are inherently regular.

We then compare the rotor-guided algorithm with different representations and discuss whether rotors are an optimal representation for the exploration of molecular conformers. We verified this approach on two simple molecules well understood in the literature, the glycine-glycine dipeptide and the glycine-phenylalanine dipeptide.

2.3.2 Molecules as rotors

We model a molecular structure as follows: we associate a rotor \(R\) to each pair of amino acids. A molecule of \(p + 1\) amino acids then characterised by a set of \(p\) rotors \(\mathbf{R} = \{R_0, \dots, R_{p-1}\}\) . Each rotor \(R_i\) applies a uniform rotation to each atom from the \(i + 1\) -th amino acid onwards.

A rotor \(R\) is described by real 4 parameters \(\alpha, \beta_{12}, \beta_{13}, \beta_{23}\) (that could actually be reduced to 3, since \(\alpha^2 = \beta_{12}^2 + \beta_{13}^2 + \beta_{23}^2 = 1\) ). 4 parameters to describe a dipeptide are a significant decrease in the numbers of degrees of freedom of the system compared to the 6 parameters of [123] for a single peptide. We do not include any rotation of the side chain as we decided to simplify as much as possible the treatment by considering only the backbone of the molecule and letting the side chain rotate with the backbone. We optimise the side chain position in 3D space through local geometry optimisation of the software xTB [59] as a refinement step before the fitness evaluation.

In the case of small molecules, in which clashes between planes of amino acids are unlikely to happen, it is possible to leave the rotor unconstrained and to employ completely random rotors, so to maximize the exploration space.

However, to make sure the generated molecular conformations are sound from a chemical point of view, we can further constrain the rotor space by limiting the amino acids to two possible movements, i.e. a rotation along the \(C_i - N_{i+1}\) axis, specified by the vector \(x_{C-N}\) , corresponding to a twist of the amino acid, or a rotation along the \(C_{\alpha,i} - C_i\) axis, specified by the vector \(x_{C_{\alpha}-C}\) , corresponding to a fold of the back bone chain (see Fig. 2.8). The corresponding rotors can be obtained as

\[R_\theta = e^{Ix_{C-N}\theta} = e^{B_{C-N}\theta} \quad (2.6)\]

\[R_\phi = e^{Ix_{C_{\alpha}-C}\phi} = e^{B_{C_{\alpha}-C}\phi} \quad (2.7)\]

in which \(I = e_{123}\) is the pseudoscalar of \(\mathcal{G}_{3,0}\) , \(B_{C-N}\) is the bivector specifying the plane of rotation perpendicular to the \(C - N\) axis, \(B_{C_{\alpha}-C}\) is the bivector specifying the plane of rotation perpendicular to the \(C_{\alpha} - C\) axis, \(\theta\) is the angle of twisting and \(\phi\) the angle of folding.

In this case, for a molecule of \(p + 1\) amino acids, we will then have \(2p\) rotors \(\mathbf{R} = \{R_{\theta,0}, R_{\phi,0}, \dots, R_{\theta,p-1}, R_{\phi,p-1}\}\) , further reducing the numbers of degrees of freedom.

(a) Twisting rotation along the \(x_{C-N}\) axis,
or equivalently in the \(B_{C-N}\) plane.

(b) Folding rotation along the \(x_{C_{\alpha}-C}\) axis,
or equivalently in the \(B_{C_{\alpha}-C}\) plane.

Figure 2.8 shows two diagrams of a tripeptide molecule illustrating the two types of rotors. (a) Twisting rotation along the x_{C-N} axis, or equivalently in the B_{C-N} plane. The diagram shows the molecule with a blue plane B_{C-N} and a red dashed line representing the x_{C-N} axis. (b) Folding rotation along the x_{C_{\alpha}-C} axis, or equivalently in the B_{C_{\alpha}-C} plane. The diagram shows the molecule with a blue plane B_{C_{\alpha}-C} and a red dashed line representing the x_{C_{\alpha}-C} axis.

Fig. 2.8 The two types of rotors employed to model molecular structures.

2.3.3 The evolutionary algorithm

We implemented a differential evolution (DE) algorithm, which belongs to the class of evolutionary algorithms (EAs). DE is a popular optimisation algorithm first introduced in [167] and employed to minimise real multivariate functions. DE starts with an initial population of candidates and keeps the candidate with the best fitness, i.e. candidates

for which the function minimised reaches a (local) minimum as the parameters are varied. The main advantage of DE is a simple problem formulation which does not require the objective function to be differentiable, i.e. it is a gradient-agnostic optimisation method.

In our case, the individual candidates are the 3D structures of molecules, encoded as .xyz files, containing the atomic symbol and the corresponding position in \(\mathbb{R}^3\) for each atom in the molecule. Borrowing from the jargon of genetic algorithms, each individual can be described by a set of genes, each of which is expressed in terms of its alleles. Our goal is to employ rotors as genes of the DE algorithm.

The objective function to be minimised is defined as

\[f(\mathbf{R}) = \Delta E(\mathbf{R}) = E_0 - E_k(\mathbf{R}) \quad (2.8)\]

where \(E_0\) is the energy of the initial molecular conformation, taken as the reference structure, and \(E_k\) is the energy of the \(k\) -th individual of the population, which is obtained applying a set of rotors \(\mathbf{R} = \{R_0, R_1, \dots, R_{p-1}\}\) , with \(p\) being the total number of genes, on the initial reference structure.

The objective function is hence a map of the type \(f : \mathbb{R}^{p \times d} \rightarrow \mathbb{R}\) , in which \(d\) is the dimensionality of the alleles, which depends on their parametrisation. For example, our gene \(R = \alpha + \beta_{12}e_{12} + \beta_{13}e_{13} + \beta_{23}e_{23}\) between each pair of \(N - C_\alpha - C\) triplets can be parametrised by its coefficients. In that case the DE algorithm will explore the space spanned by the coefficients \(\{\alpha, \beta_{12}, \beta_{13}, \beta_{23}\}\) , meaning that \(f : \mathbb{R}^{(p-1) \times 4} \rightarrow \mathbb{R}\) , in which \(p\) is the number of triplets, or amino acids, in the molecule. Conversely, if we fix the axes of rotation to two possible movements, as shown in Fig. 2.8, we only have one angle per rotor as degree of freedom, meaning that the DE algorithm will only explore the space spanned by \(\{\theta, \phi\}\) . This implies that in this case \(f : \mathbb{R}^{2(p-1)} \rightarrow \mathbb{R}\) , in which the factor of 2 comes from the fact that now we have 2 rotors per each pair of \(N - C_\alpha - C\) triplets.

The energy \(E_k(\mathbf{R})\) is evaluated through the single-point xTB method [59], meaning it only requires the structure geometry and avoids element-pair-wise parametrisation. The absolute energy, measured in hartrees (symbol: \(E_h\) ), is not insightful per se, which is why defined our function in terms of deviation from a reference instead. A negative \(\Delta E(\mathbf{R})\) means that the molecular conformation of the \(k\) -th individual, with energy \(E_k\) , is more energy favourable than the starting conformation with energy \(E_0\) .

The rest of the algorithm works as follows and it is summarized in Algorithm 1: for each individual, i.e. candidate molecular structure, with initial genes \(\mathbf{R}_x = \{R_{x,0}, \dots, R_{x,n}\} \in \mathbb{R}^n\) , with \(n = \{(p-1) \times 4, 2(p-1)\}\) based on the different representation:

• Three different structures are picked, obtained by rotating the reference conformation with genes \(\mathbf{R}_a, \mathbf{R}_b, \mathbf{R}_c\) such that \(\mathbf{R}_a \neq \mathbf{R}_b \neq \mathbf{R}_c \neq \mathbf{R}_x\) .
• Select a random index \(P \in \{0, \dots, n\}\) corresponding to a random bond in the chain.
• For each \(i \in \{0, \dots, n\}\) pick a uniformly distributed random number \(r_i \sim U(0, 1)\)
• If \(i = P\) or \(r_i < CP\) , in which \(CP\) is the crossover probability of the algorithm, then select \(R_{y,i} = R_{a,i} + F(R_{b,i} - R_{c,i})\) , in which \(F\) is the differential weight parameter of the algorithm.
• Build \(\mathbf{R}_y = \{R_{y,0}, \dots, R_{y,n}\}\)
• if \(f(\mathbf{R}_y) \leq f(\mathbf{R}_x)\) , then replace the conformation described by the set of rotors \(\mathbf{R}_x\) with that described by the set of rotors \(\mathbf{R}_y\) .
• Repeat until the number of maximum iterations has exceeded or convergence has been reached.

Algorithm 1 Differential Evolution for Molecular Structure Optimization

Require: Initial genes \(\mathbf{R}_x = \{R_{x,0}, \dots, R_{x,n}\} \in \mathbb{R}^n\) , where \(n = \{(p-1) \times 4, 2(p-1)\}\)

Require: Crossover probability \(CP\) , differential weight \(F\)

1: while stopping criteria not met do
2:   Select three different structures  $\mathbf{R}_a, \mathbf{R}_b, \mathbf{R}_c$  such that  $\mathbf{R}_a \neq \mathbf{R}_b \neq \mathbf{R}_c \neq \mathbf{R}_x$ 
3:   Select a random index  $P \in \{0, \dots, n\}$ 
4:   for  $i \in \{0, \dots, n\}$  do
5:     Sample  $r_i \sim U(0, 1)$ 
6:     if  $i = P$  or  $r_i < CP$  then
7:       Set  $R_{y,i} = R_{a,i} + F(R_{b,i} - R_{c,i})$ 
8:     else
9:       Set  $R_{y,i} = R_{x,i}$ 
10:    end if
11:  end for
12:  Construct  $\mathbf{R}_y = \{R_{y,0}, \dots, R_{y,n}\}$ 
13:  if  $f(\mathbf{R}_y) \leq f(\mathbf{R}_x)$  then
14:    Replace  $\mathbf{R}_x$  with  $\mathbf{R}_y$ 
15:  end if
16: end while

2.3.4 Examples

The Gly - Gly dipeptide

A ball-and-stick model of the Gly-Gly dipeptide in its initial conformation. The backbone atoms (carbon, oxygen, nitrogen) are colored green, red, and blue respectively. The side chains (hydrogens) are white. The structure is in a high-energy, energetically unfavourable state.

Initial conformation of Gly-Gly dipeptide

(a) The initial conformation of Gly-Gly. Both oxygen atoms (in red) and nitrogen atoms (in blue) have been manually placed on the same side of the structure, which results into an energetically unfavourable structure.

A ball-and-stick model of the Gly-Gly dipeptide in its final conformation, the

\(C_7\)

conformer. The backbone atoms are colored green, red, and blue. The structure is more compact and represents the most stable conformer of the molecule.

Conformation reached at the end of the EA

(b) The conformation reached at the end of the EA. It corresponds to the \(C_7\) conformer, which is known to be the most stable conformer of the Gly-Gly molecule.

Fig. 2.9 Starting and ending points of the DE algorithm for the Gly-Gly molecule.

Glycylglycine is the simplest dipeptide, composed of two glycine amino acids joined by a peptide bond. Glycine itself is the simplest amino acid as it has a single hydrogen atom as its side chain. We chose the glycine dipeptide because it is a molecule of particular interest: it has been found in comets by the NASA’s Stardust spacecraft [160], corroborating the hypothesis of an exogenous origin of life on Earth and the hypothesis of the existence of other forms of life in the Universe [83, 1]. Moreover, the geometry of its two most stable conformers, known as \(C_5\) and \(C_7\) , has been thoroughly characterized both theoretically and experimentally [151].

The EA for the glycine dipeptide has one gene, corresponding to one rotor, which rotates in 3D space the \(N - C_\alpha - C\) triplet of the second peptide plane and all the atoms bonded to it, while keeping the first fixed. We tested both cases with an unconstrained rotor per amino acid pair (i.e. a gene with 4 alleles, the rotors coefficients) and two constrained rotors per amino acid pairs (i.e. a gene with 2 alleles, the folding and twisting angles).

We run the EA with an initial population size of \(5 \times 4 = 20\) individuals for the unconstrained case and \(5 \times 2 = 10\) individuals for the constrained case, randomly generated with Latin hypercube sampling [67] to maximize the coverage of the specified range, and we let it run for 15 iterations. The initial structure for the dipeptide has been arbitrarily picked to have high conformational energy (see Fig. 2.9a).

The rotor spaces explored for two different allele parametrisations are shown in Fig. 2.10. It is not surprising to see how, by limiting the axes of rotation, the space searched is constrained to two unit circles on \(\Sigma\) , helping the algorithm to reach convergence. In the specific case of the Gly-Gly dipeptide, however, as few clashes happen, both approaches yield identical results, converging to the \(C_7\) conformation, shown in Fig. 2.9b.

(a) Alleles parametrised as \(\{\alpha, \beta_{12}, \beta_{13}, \beta_{23}\}\) . \(\mathbf{R} = \{R_0\}\) is unconstrained. The blue marker indicates the minimum.

(b) Alleles parametrised as \(\{\theta, \phi\}\) . \(\mathbf{R} = \{R_{\theta,0}, R_{\phi,0}\}\) have fixed rotation axes. Green: \(R_\theta\) , magenta: \(R_\phi\) . Larger markers indicate the minima.

Figure 2.10: Bivector coefficients of R on the unit sphere Σ for the Gly-Gly dipeptide. (a) shows a dense grid of red and blue points on a sphere, with a single blue marker at the bottom indicating the minimum. (b) shows a sphere with two distinct circles of points (green and magenta) and larger markers at the bottom indicating minima.

Fig. 2.10 Bivector coefficients of \(\mathbf{R}\) on the unit sphere \(\Sigma\) for the Gly-Gly dipeptide for each iteration.

The Gly - Phe dipeptide

We then studied the glycine-phenylalanine dipeptide. It has been chosen to verify whether the presence of a larger side chain, such as the benzene ring of the phenylalanine, could impact the quality of the prediction and if our pipeline could break the chirality of the molecule, i.e. its handedness. Just as in the previous case, the initial conformation has been arbitrarily picked to be energetically unfavourable (see Fig. 2.11a), so that the evolution wouldn't be biased by the initial condition. We run the EA with an initial population size of 20 individuals, randomly generated with Latin hypercube sampling to maximize the coverage of the specified range, and we let it run for 15 iterations.

Bivector coefficients are shown in 2.12. As in the previous case, whenever the proposed rotor breaks the molecular geometry, or makes two amino acid clash with

A ball-and-stick model of the Gly-Phe dipeptide in its initial conformation. The molecule is shown in a relatively flat, planar arrangement. Carbon atoms are green, hydrogen atoms are white, nitrogen atoms are blue, and oxygen atoms are red.

Initial conformation of Gly-Phe molecule

(a) The initial conformation of Gly-Phe. The structure is planar.

A ball-and-stick model of the Gly-Phe dipeptide at the end of the Evolutionary Algorithm (EA). The molecule's conformation has changed significantly from the initial planar state, now adopting a more compact, folded shape that resembles the

\(C_7\)

conformation. The chirality of the phenylalanine (Phe) residue's side chain is maintained.

Conformation reached at the end of the EA

(b) The conformation reached at the end of the EA. Gly closes to resemble \(C_7\) , and the chirality of the benzene ring is preserved.

Fig. 2.11 Starting and ending points of the DE algorithm for the Gly-Phe molecule.

A 3D plot showing the bivector coefficients of the rotation matrix

\(\mathbf{R}\)

on a unit sphere

\(\Sigma\)

. The sphere is represented by a blue wireframe grid. Two great circles are highlighted: a green one representing the

\(R_{\theta,0}\)

rotation axis and a magenta one representing the

\(R_{\phi,0}\)

rotation axis. Several markers of varying sizes are placed on the sphere's surface, with larger markers indicating the minima of the objective function.

Bivector coefficients on the unit sphere

Fig. 2.12 Bivector coefficients of \(\mathbf{R}\) on the unit sphere \(\Sigma\) for the Gly-Phe dipeptide for each iteration. Alleles are parametrised as \(\{\theta, \phi\}\) . \(\mathbf{R} = \{R_{\theta,0}, R_{\phi,0}\}\) has fixed rotation axes. Green: \(R_{\theta,0}\) , magenta: \(R_{\phi,0}\) . Larger markers indicate the minima.

each other, or breaks a bond, the resulting 3D structure is not fed into xTB and the value of \(f(\mathbf{R})\) is set to an arbitrarily large value, in our case \(10^4\) .

As in the previous case, despite the larger number of iterations, the algorithm converges to a \(C_7\) conformation, which has the lowest energy out of the \(\sim 3000\) conformations explored and hence it is energetically more favourable (see Fig. 2.11b).

Table 2.6 Comparison of rotation representations and convergence speed, expressed as number of iterations of the EA, for Gly-Gly and Gly-Phe dipeptides.

Table 2.6 summarizes the comparison. Rotor-based representations offer a competitive alternative to dihedral-angle parametrizations, which involve a higher number of degrees of freedom and require sequential rotations about distinct local axes. Although dimensionality can be reduced through constrained dihedral formulations, rotor-based approaches benefit from an efficient implementation via the sandwich product, resulting in faster convergence and reduced computational overhead during rotational updates.

2.4 Conclusions

In this chapter we provided two examples in which rotors are used to parametrise rotations in two different learning problems in an effective and compact way.

In Section 2.2 we have proven that certain representations of rotations might be more suitable than others according to the specific task to be solved. We have also showed how a fully GA description of the rotation allows us to bypass the continuity issue, yielding results which are comparable, if not superior, to the continuous 6D representation without needing a novel representation in all three toy problems presented. We believe that the errors introduced when learning rotations are due to the specific rotation matrix-to-representation mapping, and not an intrinsic limitation of the representation itself: for this reason, we will make rotors and bivectors our representation of choice, since they are readily defined, more compact and show enhanced robustness to noise.

In Section 2.3 we have demonstrated how rotors can be employed as a metaheuristic for global geometry optimisation of small molecules. We showed how to model peptides in 3D as a collections of rotors, that have fewer degrees of freedom compared to dihedral angles, and how they can be used to drive an evolutionary algorithm to reach stable conformers in fewer iterations.

Chapter 3

Protein Structure Prediction (PSP)

How do you know that the side you
are used to is better than the one to
come?

— The Forty Rules of Love, Elif
Shafak (2009)

3.1 Introduction

In this chapter we focus on the modelling of proteins and on the estimation of their 3D structure. The 3D structure of a protein - known as its tertiary structure - is the arrangement in space of its amino acid chain. It determines the protein functionality, which is important in drug design and discovery, disease understanding and treatment and more. Determining the structure experimentally, however, is expensive and time consuming, and researchers have been employing computational methods instead, hence the field of protein structure prediction (PSP). Here we explore the advantage of using GA in PSP-related problems.

We first use GA as a modelling tool, to generate a novel set of features that capture information about the proteins local folding (its secondary structure) more effectively than commonly used dihedral angles.

We then use GA as a computational tool, to construct an equivariant neural network in Clifford Algebra that improves on the the accuracy of PSP pipelines over linear layers.

3.2 Protein structures as 3D CGA planes

As anticipated in Section 2.3, a protein can be simplified into a backbone chain and several side chains. The backbone is mostly responsible for the 3D shape of the protein, and it is composed of a sequence of carbon, nitrogen, and oxygen atoms. The \(\alpha\) -carbons are the main feature of the backbone, to which the side chains that differentiate each amino acid are bonded. Each \(\alpha\) -carbon is preceded by a nitrogen atom and followed by a carbon atom. Hence, to each amino acid \(i\) we can associate a triplet of atoms \(\{N, C_\alpha, C\}_i\) .

Each \(\{N, C_\alpha, C\}\) triplet lies on a plane, constraining the protein folding (see Figure 3.1). We can conveniently model a protein backbone in CGA so any three \(\{N, C_\alpha, C\}\) atoms will lie on a plane (not too dissimilar to the residue gas of [81]): let \(A_i, B_i\) and \(C_i\) be the Euclidean coordinates expressed in Conformal space of the atoms \(\{N, C_\alpha, C\}_i\) , respectively. The plane associated with residue \(i\) can be expressed as the 4-blade:

\[\Pi_i = A_i \wedge B_i \wedge C_i \wedge n_\infty \quad (3.1)\]

Given two planes \(\Pi_i, \Pi_j\) corresponding to the amino acids \(i, j\) , we can compute the rotor that brings one to the other as described in [102]:

\[R_{ij} = \frac{1}{\sqrt{\langle K \rangle_0}} (1 - \Pi_i \Pi_j) \quad (3.2)\]

where \(K = 2 - (\Pi_i \Pi_j + \Pi_j \Pi_i)\) and \(\langle \cdot \rangle\) is the grade projector operator. We now use the cost function \(C_\lambda(R)\) that measures how much the rotor \(R\) varies from the identity, as defined in [46]. \(C_\lambda(R)\) is a weighted sum of a translational and a rotational term:

\[C_{\lambda_1 \lambda_2}(R) = \lambda_1 \langle R_{\parallel} \tilde{R}_{\parallel} \rangle_0 + \lambda_2 \langle (R_{\perp} - 1)(\tilde{R}_{\perp} - 1) \rangle_0 \quad (3.3)\]

in which the translational error is represented by \(R_{\parallel} = R \cdot e\) , the rotational error by \(\langle (R_{\perp} - 1)(\tilde{R}_{\perp} - 1) \rangle_0 = \langle (R - 1)(\tilde{R} - 1) \rangle_0\) and \(\lambda_1, \lambda_2 \in \mathbb{R}\) are two weighting scalars. As we are interested in a orientational features, we will focus exclusively on the rotational part (case \(\lambda_1 = 0, \lambda_2 = 1\) ).

Figure 3.1: HIV-1 protease (ID: 1dmp) modelled as a collection of N - C_alpha - C planes. The left panel shows the full protein structure as a collection of red planes in a 3D coordinate system (x, y, z). The right panel is a close-up view of the first 20 residues, showing two specific planes, Pi_1 and Pi_j, and the rotation R_ij between them.

Fig. 3.1 HIV-1 protease (ID: 1dmp) modelled as a collection of \(N - C_\alpha - C\) planes (left) and close up on the first 20 residues (right). Note \(\Pi_j = R_{ij}\Pi_i\tilde{R}_{ij}\) .

3.2.1 Geometrical interpretation of the cost function

A rigorous interpretation of the cost function is derived in [46]-[184]. A rotor in 3D conformal space can be represented via an exponential:

\[R = e^B \quad (3.4)\]

in which \(B\) is a bivector. Assuming \(B = \phi P + tn\) , it is derived in [185] that

\[R = e^B = [\cos(\phi) + P \sin(\phi)][1 + t_\perp n_\infty] + \text{sinc}(\phi)t_\parallel n_\infty \quad (3.5)\]

where \(P\) is a bivector, representing the normalized rotation plane, such that \(P^2 = -1\) , \(2\phi\) is the rotation in plane \(P\) and \(t \in \mathbb{R}^3\) is the translation vector. By grouping the components including and not including \(n_\infty\) , we can rewrite Eq. 3.5 as the sum of a perpendicular and a parallel component of the rotor:

\[\begin{aligned} R &= R_\perp + R_\parallel n_\infty \\ R_\perp &= \cos(\phi) + P \sin(\phi) \\ R_\parallel &= \cos(\phi)t_\perp + \text{sinc}(\phi)t_\parallel + \sin(\phi)t_\perp \end{aligned} \quad (3.6)\]

In the given decomposition, \(R_\perp\) represents the rotational component of the rotor, as it consists of terms independent of \(n_\infty\) and describes a pure rotation in the conformal model. In contrast, \(R_\parallel\) captures the translational effects, as it includes terms proportional to \(n_\infty\) , which in CGA encode translations in space.

We can then rewrite the rotational component of Eq. 3.3 as:

\[\begin{aligned} \langle (R_{\perp} - 1)(\tilde{R}_{\perp} - 1) \rangle_0 &= \langle R_{\perp} \tilde{R}_{\perp} - R_{\perp} - \tilde{R}_{\perp} + 1 \rangle_0 = \\ \langle 1 - \cos(\phi) - P \sin(\phi) - \cos(\phi) + P \sin(\phi) - \cos(\phi) + 1 \rangle_0 &= \quad (3.7) \\ 2(1 - \cos(\phi)) &= 4 \sin^2(\phi/2) \end{aligned}\]

The expression of the rotational error is hence a convex function exclusively dependent on the rotation angle \(\phi\) between the two planes, with a maximum of 2 for a rotation of \(2\phi = \pm\pi\) and a minimum at \(\phi = 0\) . This shows that the cost function of Eq. 3.3 can be used as a measure of the relative orientation of amino acids within a protein chain.

3.3 Cost maps: a geometric description of amino acids

Inter-residue interactions are commonly represented as matrices - also called maps. A contact map C of a protein consisting of \(M\) residues, for example, is a binary \(M \times M\) matrix of the type:

\[C_{ij} = \begin{cases} 1 & \text{if } d_{ij} < 15 \text{ \AA} \\ 0 & \text{otherwise} \end{cases} \quad (3.8)\]

where \(d_{ij}\) is the distance between residues \(i, j\) expressed in \(\text{\AA}\) measured as the Euclidean distance between the \(C_{\alpha}\) coordinates of residues \(i\) and \(j\) . The threshold of \(15 \text{ \AA}\) (in which \(1 \text{ Angstrom} = 10^{-10} \text{ m}\) ) has been picked as an arbitrary reference to define two amino acids in contact with each other; a threshold of 10 or \(8 \text{ \AA}\) would yield similar results. A contact map can hence be interpreted as: two residues are in contact if they are within a certain distance from each other. A more informative metric, usually real-valued, is given by distance maps, which are similarly defined as:

\[D_{ij} = d_{ij} \quad (3.9)\]

From either or both contact and distance maps it is possible to obtain accurate 3D shape estimation. However, when contact or distance maps are predicted and not exact, errors are introduced into the 3D reconstruction step. Having an additional map capturing the orientation between residues can help to further constrain the search space for the protein folding. We can hence employ our cost function to produce a cost map which contains orientational information as follows:

\[M_{ij} = \begin{cases} C_{\lambda_1 \lambda_2}(R_{ij}) & \text{if } d_{ij} < 15 \text{ \AA} \\ 0 & \text{otherwise} \end{cases} \quad (3.10)\]

A comparison between contact map C, distance map D and cost map M is given in Figure 3.2 for an example protein. Note how the range of our cost maps is \([0, 2]\) , which is consistent with Eq. 3.7.

Figure 3.2: Comparison of contact, distance, and cost maps for protein chain 2HC5A. (a) Contact map showing binary values (0 or 1) with a color scale from 0.0 to 1.0. (b) Distance map showing continuous values with a color scale from 0 to 40. (c) Cost map showing continuous values with a color scale from 0.00 to 1.75. All maps are 100x100 matrices.

Fig. 3.2 (a) contact, (b) distance and (c) cost map for protein chain 2HC5A.

3.3.1 Cost maps and secondary structures

It is possible to establish a relationship between patterns in cost maps and the protein secondary structure. By secondary structure we refer to the local folding of a segment of a protein, e.g. \(\alpha\) -helices, \(\beta\) -sheets or turns. Secondary structure information is a common feature in PSP pipelines and one of the most important in predicting distance and contact maps, as shown in [38, 2].

By assigning a colour to each secondary structure, it is possible visualize the secondary structure of each amino acid pair. We arbitrarily assigned red to \(\alpha\) -helices, green to \(\beta\) -sheets, blue to turns and white to all the others. Any combination of these four colours gives the possible secondary structures of the pair, for a total of 10 different colour combinations. As shown in Figure 3.3, it is possible to find a clear correspondence between secondary structures and patterns in the cost maps; for example, \(\alpha\) -helices patterns include parallel lines, equally spaced with low intensity, while \(\beta\) -sheets present high intensity pixels packed tightly together. These patterns are identical and consistent for different proteins, as shown in Figure 3.4. Note also how patterns between residues which are in contact with each other but belonging to

(a)

(b)

Figure 3.3: (a) Cost map and (b) secondary structures for protein chain 4JZK. (a) shows a cost map with a color scale from 0 to 175 on both axes. (b) shows secondary structures with a color scale from 0 to 200 on both axes.

Fig. 3.3 (a) Cost map and (b) secondary structures for protein chain 4JZK. Visualizing the first 200 residues. Red: \(\alpha\) -helices, green: \(\beta\) -sheets, blue: hydrogen bond turns, white: all the others.

3B43

6J48

1G03

2HC5

1K20

1K20

Figure 3.4: Examples of (zoomed-in) cost maps for different proteins along with their secondary structures. The figure shows six pairs of images for proteins 3B43, 6J48, 1G03, 2HC5, 1K20, and 1K20. Each pair consists of a cost map on the left and a secondary structure map on the right.

Fig. 3.4 Examples of (zoomed-in) cost maps for different proteins along with their secondary structures. Patterns associated to a same secondary structure are consistent for different protein chains.

different secondary structures are consistent, e.g. helix-sheet interactions (in red-green colour), in the off-diagonal areas, present typical braided patterns.

To the best of our knowledge, this is the first example of an orientational map that also encodes the secondary structure of the protein so clearly: the three angle maps presented in [197] obtained via dihedral angles of the protein chains, which are also an orientational measure, do not present patterns as evident as those of our cost maps.

3.3.2 Cost maps and different protein conformations

Proteins structures are not static. A protein can display different conformations, or “states”, based on their activity - an example of this is given by enzymes, i.e. proteins that catalyze chemical reactions and can change their conformation when they bind to different substrates. We will then have an “open” state, when the enzyme is ready to bind to the substrate, and the “closed” state, when the enzyme is inactive.

For example, the Lysine-, Arginine-, Ornithine-binding protein has an open (PDB ID: 2LAO) state and a closed (PDB ID: 1LAF) state [164]. Chains 1LAF and 2LAO are shown in Figure 3.5. The corresponding cost maps are given in Figure 3.6. In Figure 3.6c we see that that a higher activity is associated with the second half of the protein, approximately from residues 80 to 200: this is mirrored in the 3D structure, for which the cyan and magenta chains diverge more for said residues (see Figure 3.5). Note also the red regions in Figure 3.6c: this means that the chain of Figure 3.6a is folded in such a way that residues in the region 50-150 are in contact (i.e. at a distance below 15 Å), differently from the chain of Figure 3.6b.

A 3D ribbon diagram showing two protein conformations, 1LAF (cyan) and 2LAO (magenta), overlaid. The structures show significant conformational differences, particularly in the open/closed states. The cyan chain (1LAF) is more compact and folded, while the magenta chain (2LAO) is more extended and open. The divergence is most prominent in the middle region of the protein chains.

3D ribbon diagram showing two protein conformations, 1LAF (cyan) and 2LAO (magenta), overlaid. The structures show significant conformational differences, particularly in the open/closed states.

Fig. 3.5 3D structures for chains 1LAF (cyan) and 2LAO (magenta).

Another example of different protein conformation is given by the spike proteins of viruses: the open (PDB ID: 6VYB) and closed (PDB ID: 6VXX) states of the

Figure 3.6: Cost maps of the first 300 residues of protein chains (a) 1LAF, (b) 2LAO and (c) their difference. The figure consists of three heatmaps. (a) and (b) show cost maps for chains 1LAF and 2LAO respectively, with axes from 0 to 200. (c) shows the difference map with a color scale from -0.05 (blue) to 0.05 (red).

Fig. 3.6 Cost maps of the first 300 residues of protein chains (a) 1LAF, (b) 2LAO and (c) their difference.

SARS-CoV-2 S glycoprotein are shown in Figure 3.7 [121]. The corresponding cost maps are given in Figure 3.8, with Figure 3.8c showing the difference between the two. Note how the two structures are consistent with each other, with highest deviation around residues 300-400, which correspond to the folding coloured in orange of chain 6VYB in the bottom-left part of Figure 3.7, which is absent in chain 6VXX. Note also how these differences would be very hard to spot by considering Figures 3.8a-b individually.

Figure 3.7: 3D structures for chains 6VXX (blue) and 6VYB (orange). The image shows two ribbon diagrams of protein structures, one in blue and one in orange, overlaid. The structures are complex, with many alpha-helices and beta-sheets. The orange structure (6VYB) has a distinct folding region in the bottom-left that is absent in the blue structure (6VXX).

Fig. 3.7 3D structures for chains 6VXX (blue) and 6VYB (orange).

Figure 3.8: Cost maps of the first 500 residues of protein chains (a) 6VXX, (b) 6VYB and (c) their difference. The figure consists of three heatmaps labeled (a), (b), and (c). Each heatmap has x and y axes ranging from 0 to 400. A color bar on the right indicates values from -0.4 (blue) to 0.4 (red). Heatmaps (a) and (b) show a dark blue background with a magenta/purple diagonal band and some scattered points. Heatmap (c) shows a light pink background with scattered magenta/purple points.

Fig. 3.8 Cost maps of the first 500 residues of protein chains (a) 6VXX, (b) 6VYB and (c) their difference.

3.3.3 Cost maps and NMR ensembles

Figure 3.9: Three different models within the NMR ensemble of protein chain 1BVM with the complete ensemble. The figure shows four ribbon diagrams of a protein structure. The first three are individual models, and the fourth is the complete ensemble. The fourth diagram has two green boxes highlighting specific regions of the protein structure.

Fig. 3.9 Three different models within the NMR ensemble of protein chain 1BVM with the complete ensemble.

Nuclear magnetic resonance (NMR) is a technique in analytical chemistry to determine the purity or the structure of a sample. By NMR ensemble we refer to a collection of possible foldings of a protein, which, being flexible, does not possess a unique conformation. An NMR ensemble is hence a collection of different foldings of the same protein, with some small local differences. An average difference of the cost maps associated with multiple structures resolved via NMR can be obtained. The average difference is obtained in 3 steps:

• Compute the cost map for each structure in the ensemble, for a total of \(K\) structures and hence \(K\) cost maps.

• Take the difference between consecutive pairs of cost maps, obtaining \((K - 1)\) differences of cost maps. This approach is taken because consecutive cost maps exhibit smaller divergences from each other.

• Take the overall average by summing the differences of cost maps and divide the resulting plot by \((K - 1)\) pixel-wise.

Regions with higher intensity will appear in the final plot, meaning that those regions were the most different across the different 3D structures. This highlights residues or regions associated with a higher resolution error or with higher flexibility.

We verified this with an NMR ensemble of 20 chains of protein 1BVM (see Figure 3.9: averaging the pairwise differences for each of the 20 chains, we end up with the plot in Figure 3.10. The higher activity around residues 10-50 (corresponding to the green box in Figures 3.9-3.10) or 50-60 (corresponding to the cyan box in Figures 3.9-3.10) signifies that those residues are associated with a higher uncertainty.

Note how, on the contrary, the helix structures are well resolved, which is clear both from Figure 3.9 and from Figure 3.10d, in which residues 30-60 or 90-110, corresponding to helices, have an average cost map difference close to 0. Note also how it is possible to spot the helix structures immediately from Figure 3.10a-b-c.

Figure 3.10 consists of four panels labeled (a), (b), (c), and (d). Panels (a), (b), and (c) are cost maps for protein 1BVM, each showing a 2D grid of residues. Each map has a dark background with bright, diagonal, and block-like patterns. Three colored boxes (green, cyan, and yellow) are overlaid on each map, highlighting specific regions of interest. Panel (d) is the average difference of the 20 cost maps. It shows a similar 2D grid but with a color scale on the right ranging from -0.20 (dark blue) to 0.20 (yellow). The map is predominantly dark blue, indicating low average difference, with some bright spots corresponding to the regions highlighted in the other panels.

Figure 3.10: Four cost maps (a, b, c, d) for protein 1BVM. (a), (b), and (c) show individual cost maps with colored boxes highlighting specific regions. (d) shows the average difference of the 20 cost maps with a color scale from -0.20 to 0.20.

Fig. 3.10 (a)-(c): Three cost maps associated with 3 resolved structures of the NMR ensemble of chain 1BVM. Note the differences in pattern among the three maps highlighted by the coloured boxes. (d): Average difference of the 20 cost maps. The higher activity in the averaged map can be associated to higher resolution uncertainty or variability in the structure.

3.4 Alternative representations

We are now going to introduce two more maps describing the orientation between pairs of amino acids.

3.4.1 Dihedral angles

Dihedral angles are the most common way to describe the orientation of atoms in the protein backbone. They are the angles between planes formed by different atom triplets in the backbone. Commonly, each amino acid \(i\) can be described with three dihedral angles \(\{\phi_i, \psi_i, \omega_i\}\) .

Ball-and-stick model of Insulin (ID: 3i40) on the left and a close-up of the backbone in position i on the right. The close-up shows the backbone atoms: C_alpha,i+1, N_i+1, C_i, C_alpha,i, C_beta,i, N_i, and C_i-1. The dihedral angles are highlighted: omega_i (green) between N_i+1-C_i-C_alpha,i and C_i-N_i+1-C_alpha,i+1; phi_i (pink) between C_beta,i-C_alpha,i-C_i and N_i-C_alpha,i-C_i; and psi_i (blue) between C_i-1-N_i-C_alpha,i and N_i-C_alpha,i-C_i.

Fig. 3.11 Ball-and-stick model of Insulin (ID: 3i40, left) and close up on the backbone in position \(i\) with dihedral angles \(\omega_i, \phi_i, \psi_i\) highlighted (right).

Following the convention in [27], the angle \(\psi_i\) is defined as the dihedral angle between the \(C_{i-1} - N_i - C_{\alpha,i}\) and the \(N_i - C_{\alpha,i} - C_i\) planes, the angle \(\phi_i\) as the dihedral angle between the \(C_{\beta,i} - C_{\alpha,i} - C_i\) and the \(N_i - C_{\alpha,i} - C_i\) planes, while \(\omega_i\) is the dihedral angle between the \(N_{i+1} - C_i - C_{\alpha,i}\) and the \(C_i - N_{i+1} - C_{\alpha,i+1}\) planes (see Fig. 3.11). A total of 5 planes, all relative to residue \(i\) , are involved to specify the triplet of dihedral angles \(\{\phi_i, \psi_i, \omega_i\}\) .

In order to build angle maps, we need dihedral angles between pairs of residues. We will hence define \(\{\phi_{ij}, \psi_{ij}, \omega_{ij}\}\) as the dihedral angles between the same two pair planes mentioned above, but with the first plane relative to residue \(i\) and the second relative to residue \(j\) . With this assumption we have that \(\phi_i = \phi_{ii}, \psi_i = \psi_{ii}\) and \(\omega_i = \omega_{ii}\) . We can hence define the angle maps as follows:

\[\Phi_{ij} = \begin{cases} \phi_{ij} & \text{if } d_{ij} < 15 \text{ \AA} \\ 0 & \text{otherwise} \end{cases} \quad (3.11)\]

\[\Psi_{ij} = \begin{cases} \psi_{ij} & \text{if } d_{ij} < 15 \text{ \AA} \\ 0 & \text{otherwise} \end{cases} \quad (3.12)\]

(a) \(\Phi\) (b) \(\Psi\) (c) \(\Omega\)

Figure 3.12: Angle maps for protein 12asA. (a) Φ map, (b) Ψ map, (c) Ω map. Each map is a 300x300 heatmap with a color scale from 0.0 (dark blue) to 3.0 (yellow). The Ψ map shows a distinct asymmetry compared to the others.

Fig. 3.12 Angle maps for protein 12asA.

\[\Omega_{ij} = \begin{cases} \omega_{ij} & \text{if } d_{ij} < 15 \text{ \AA} \\ 0 & \text{otherwise} \end{cases} \quad (3.13)\]

where \(d_{ij}\) is the Euclidean distance between the \(C_\alpha\) 's of residues \(i, j\) measured in Angstroms ( \(1 \text{ \AA} = 10^{-10} \text{ m}\) ). The cut-off distance is required to discriminate which amino acids are in contact with each other, and \(15 \text{ \AA}\) is a recurrent value in the literature; a value of 10 or \(8 \text{ \AA}\) could have also been picked, yielding slightly thinner profiles in the maps. Examples of angle maps are given in Figure 3.12. Note how the patterns in \(\Psi\) are asymmetric.

3.4.2 Dot product between oriented points

Oriented points were first introduced in [44] with the name of “tangent bivectors”. In [69], an oriented point \(Q\) in CGA is a trivector, i.e., a circle, with radius \(r = 0\) , defined as

\[Q = I_q \wedge q + \left[ \frac{1}{2} q^2 I_q - q(q \cdot I_q) \right] n_\infty + I_q n_0 + (I_q \cdot q) E. \quad (3.14)\]

In the equation above \(q \in \mathbb{R}^3\) is the 3D position vector defining the centre of \(Q\) , \(I_q\) is the (unit) bivector corresponding to the oriented plane in which \(Q\) lies (orthogonal to the normal \(n_q\) of the plane) and \(E = n_\infty \wedge n_0\) is the origin-infinity bivector.

Given a pair of oriented points \(P, Q\) , the dot product between them was first studied in [71], and it was derived to be equal to

\[P \cdot Q = d^2 \left[ -\frac{1}{2} \cos \alpha_{pq} + \cos \Theta_q \cos \Theta_p \right], \quad (3.15)\]

in which \(d\) is the Euclidean distance between \(p, q \in \mathbb{R}^3\) , \(\alpha_{pq}\) is the dihedral angle between the two planes, \(\cos \alpha_{pq} = n_q \cdot n_q\) , \(\cos \Theta_p = d \cdot n_p\) and \(\cos \Theta_q = d \cdot n_q\) .

Equation 3.15 suggests that, up to a scale factor, the dot product between two oriented points is a function of three angles, hence it encodes the orientation between them. Moreover, similarly to the cost function of Eq. 3.10, Equation 3.15 allows us to obtain a scalar measure from an initial GA description.

As each atom in the backbone can be associated with a set of 3D coordinates and it lies on a plane specified by the bonds to the previous and following atoms in the backbone, it is easy to associate an atom with an oriented point. We will call \(Q_{C_\alpha}\) the oriented point centered on \(C_\alpha\) in the \(N - C_\alpha - C\) plane (see Fig. 3.13), \(Q_{C_\beta}\) the oriented point centered on \(C_\beta\) in the \(C_\alpha - C_\beta - C\) plane and so on.

The figure consists of two panels. The left panel shows a 3D perspective view of a protein structure, Haemoglobin (ID: 1a3n), represented as a dense collection of blue oriented points. The points are clustered together, forming the shape of the protein. The right panel is a close-up view of five residues, showing the oriented points as blue ellipses on a plane. Two specific oriented points are labeled: \(Q_{C_\alpha,i}\) (in red) and \(Q_{C_\alpha,j}\) (in green). Normal vectors \(n_{q,i}\) (red arrow) and \(n_{q,j}\) (green arrow) are shown originating from these points. A dashed line connects the two points, representing the vector \(d\) between them. The axes are labeled x, y, and z.

Figure 3.13: Haemoglobin (ID: 1a3n) modelled as a collection of oriented points centered at C_alpha (left) and close up on 5 residues with labelled oriented point (right).

Fig. 3.13 Haemoglobin (ID: 1a3n) modelled as a collection of oriented points centered at \(C_\alpha\) (left) and close up on 5 residues with labelled oriented point (right).

We can again build an \(M \times M\) map for each protein chain of length \(M\) , which we will call dot product maps, defined as:

\[N_{X,ij} = \begin{cases} \frac{1}{d_{ij}^2} (Q_{X,i} \cdot Q_{X,j}) & \text{if } d_{ij} < 15 \text{ \AA} \\ 0 & \text{otherwise} \end{cases} \quad (3.16)\]

In the equation above, \(X\) is any atom among \(C_\alpha, C_\beta, C\) or \(N\) , and \(d_{ij}\) is the Euclidean distance between \(X_i\) and \(X_j\) . Note the normalization factor as we are only interested in the orientation between atoms. Examples of dot product maps are given in Figure 3.14.

(a) \(\mathbf{N}_{C_\alpha}\) (b) \(\mathbf{N}_{C_\beta}\) (c) \(\mathbf{N}_N\)

Figure 3.14: Dot product maps for protein 3i41A. (a) N_{C_\alpha} map showing a diagonal band of high values. (b) N_{C_\beta} map showing a more complex pattern of high values. (c) N_N map showing a pattern of high values. Each map has a color bar on the right indicating values from -0.4 to 0.4.

Fig. 3.14 Dot product maps for protein 3i41A.

3.5 Dataset: PDNET

We employed an expanded version of PDNET [2] to train our models in Section 3.7 and Section 3.8. The original PDNET includes 1000 proteins for training and 150 proteins for testing, of which the top seven features in contact and distance prediction problems are picked. These features are sequence profiles, secondary structures, solvent accessibility, coevolutionary signals, FreeContact [84], contact potentials alignments and Shannon entropy, all predicted starting from the amino acid sequence. The seven features are encoded in PDNET as a stack of \(57 M \times M\) channels for each of its protein chains of length \(M\) . Of the 57 channels, 3 of them correspond to pairwise features, i.e., FreeContact, coevolutionary signals and contact potentials, while the remaining 54 are 27 individual features repeated twice as a matrix and their corresponding transpose (i.e., \(Y\) and \(Y^T\) ).

We expanded PDNET by adding distance maps (defined as \(\mathbf{D}_{ij} = d_{ij}\) , where \(d_{ij} = \|T_i - T_j\|_2\) , with \(T \in \mathbb{R}^{N \times 3}\) being the ground truth coordinates of the \(C_\alpha\) atoms of the protein) and eight possible combinations of the cost maps (namely \(\mathbf{M}_{C_\alpha}, \mathbf{M}_{C_\beta}, \mathbf{N}_{C_\alpha}, \mathbf{N}_{C_\beta}, \mathbf{N}_N, \Omega, \Phi\) and \(\Psi\) ), for a total of nine different cases (see Table 3.1).

We then recast the dataset in the form of heterogeneous graphs so that it can be input into the Graph Transformer (GrT, see Section 5.2.1). By heterogeneous graph \(\mathcal{G}(V, E)\) with \(V\) and \(E\) being its set of nodes and edges, respectively, we denote a graph with different types of nodes and edges. If \(|V| = M\) is the total number of nodes, the graph can be described as a set of adjacency matrices for each of the \(K\) edge types. Adjacency matrices are encoded as \(\{A_k\}_{k=1}^K\) , where \(A_k \in \mathbb{R}^{M \times M}\) , or in tensor form \(\mathbf{A} \in \mathbb{R}^{M \times M \times K}\) : we say that \(A_{k;i,j}\) is non-zero when there exists an edge of type \(k\)

Table 3.1 Combinations of orientational features. The column “Planes” specifies how many planes are required to build the corresponding set of features. The column “ \(K\) ” indicates the total number of adjacency matrices of the graph.

between nodes \(i, j\) . Along with \(\mathbf{A}\) , we can also define a feature matrix \(X \in \mathbb{R}^{M \times D}\) , where \(D\) is the dimensionality of the features, or equivalently we can say there are \(D\) node types.

Ignoring the transposed matrices, we have \(D = 27\) channels which correspond to features relative to individual amino acids which can be manipulated and arranged in a feature matrix \(X \in \mathbb{R}^{M \times D=27}\) . On the other hand, the pairwise features range from a minimum of \(K = 4\) in case (a) up to a maximum of \(K = 7\) in cases (g)-(i) which correspond to the edges of the protein graph, i.e., the adjacency matrices \(\mathbf{A} \in \mathbb{R}^{M \times M \times K}\) (see Table 3.1).

The input to the Graph Transformer is then given by the pair of tensors \(\{\mathbf{A}, X\}^{(i)}\) for each protein \(i\) in the dataset.

3.6 Estimating orientational features

Diagram illustrating the architecture for predicting orientational features. The input consists of feature channels of size M x M x L. This is processed through a sequence of layers: Batch Norm, ReLU, Conv 3x3 (64 channels), Dropout, ReLU, and Conv 3x3 (L channels). A skip connection bypasses these layers and adds the input. The output is then processed by a final ReLU / tanh layer to produce the predicted orientational feature.

Fig. 3.15 Predicting orientational features. We employed the residual neural network of [2] to predict (a) cost maps \(\mathbf{M}_\alpha\) , (b) dot product maps \(\mathbf{N}_\alpha\) and (c) angle maps \(\Phi, \Psi, \Omega\) . While all (true) orientational features provide similar improvements when predicting coordinates, not all of them are as readily predicted.

Results presented in Section 3.8 onwards assume that the orientational features are built starting from ground truth coordinates of atoms in the protein. In a realistic

scenario, however, the coordinates are the end goal of the pipeline, not the starting point. We would then have to employ predicted orientational features, which are less clear than the ground truth features shown so far and must themselves be estimated from other features or multiple sequence alignments (MSA) of the amino acid sequence.

To this end, we employed the PDNET pipeline of [2], originally designed to predict distance maps, to compare how easily different orientational maps are predicted starting from the same set of features and the same pipeline. We attempted to predict cost maps \(\mathbf{M}_{C_\alpha}\) , dot product maps \(\mathbf{N}_{C_\alpha}\) and the three angle maps. The prediction pipeline is shown in Figure 3.15.

Example of predictions are given in Figures 3.16-3.17: cost maps \(\mathbf{M}\) were consistently more accurately predictable. This is likely because they are closely related to secondary structure information, which is implicitly included in other input features. Less accurate are the predictions of angle maps, which don't present easily recognizable patterns or are asymmetric as in \(\Psi\) . A similar consideration can be made for dot product maps.

While all true orientational maps increase the coordinate prediction accuracy in a similar way, not all orientational maps are readily predicted, meaning that in a realistic PSP pipeline we will prefer features that are easier to estimate for the same added value in terms of GDT scores.

Figure 3.16: Ground truth and predicted orientational maps for protein 1a3aA. The figure is a 2x5 grid of heatmaps. The columns are labeled with mathematical symbols: φ, Ψ, Ω, N, and M. The rows are labeled 'Ground Truth' and 'Prediction'. The heatmaps show patterns of varying complexity and symmetry. Ground truth maps for φ, Ω, and M show clear, symmetric patterns. Ground truth for Ψ is asymmetric. Ground truth for N is more complex and noisy. Predicted maps for φ, Ω, and M are very similar to their ground truth counterparts. Predicted Ψ is more symmetric than ground truth Ψ. Predicted N is significantly less accurate, showing a diagonal band of high values instead of the complex ground truth pattern. Predicted M remains relatively accurate.

Fig. 3.16 Ground truth and predicted orientational maps for protein 1a3aA.

3.6.1 Training details

The loss of the PDNET pipeline has been chosen to be

\[\mathcal{L}_{\mathbf{X}} = \log(\cosh(\mathbf{X}_P - \mathbf{X}_T)) \quad (3.17)\]

Figure 3.17: Ground truth and predicted orientational maps for protein 1a70A. The figure shows two rows of five heatmaps labeled with symbols: Φ, Ψ, Ω, N, and M. The top row is labeled 'Ground Truth' and the bottom row is labeled 'Prediction'. Each heatmap shows a complex pattern of colors (blue, orange, yellow) representing orientational information across a protein structure.

Fig. 3.17 Ground truth and predicted orientational maps for protein 1a70A.

where \(\mathbf{X}_P, \mathbf{X}_T\) are the predicted and true orientational maps in the training set, respectively, with \(\mathbf{X} = \{\mathbf{M}, \mathbf{N}, \Phi, \Psi, \Omega\}\) . The only difference between the approaches is the activation function of the last layer: when predicting \(\mathbf{N}\) , which ranges between \([-0.5, 0.5]\) , we employed a tanh activation instead of a ReLU.

The training features are identical to those of PDNET, namely a stack of images of the type \(\{\mathbf{Y}^{(i)}\}_{i=1}^N\) , with \(L = 57\) and \(\mathbf{Y}^{(i)} \in \mathbb{R}^{M \times M}\) , in which \(M\) is the length of the protein sequence. The loss is evaluated per pixel. The training set has been kept to 1000 proteins from the DEEPCOV dataset, and the testing set to 150 proteins from the PSICOV dataset, as in the original PDNET pipeline.

3.7 PSP with GA features

3.7.1 Architecture: Graph Transformer + 3D projector

The end-to-end architecture, derived from [36], is composed of two parts: (1) a Graph Transformer and (2) a 3D projector. A summary of the architecture is shown in Fig. 3.18. We omitted the MSA Transformer of [36] as the employed dataset allows us to directly perform node and edge embedding on its features.

Figure 3.18: The employed architecture diagram. It shows the flow from input features to 3D coordinates. Input: 'feature channels' (represented by three stacked rectangles) and 'orientational features' (represented by a stack of M x M matrices). These are combined via '+' and then processed through 'node embedding' (X, M x D) and 'edge embedding' (A, M x M x K). The embeddings are fed into a 'Graph Transformer' block. The output is a 'node representation' (Z, M x D), which is then processed by a '3D Projector' block to produce '3D coordinates' (M x 3).

Fig. 3.18 The employed architecture. Features are recast in graph form to predict 3D backbone coordinates via a GrT + 3D projector architecture

The Graph Transformer

The GrT has been implemented as described in [199]. The goal of a GrT is to learn informative meta-paths within the graph, i.e., an ordered sequence of node types and edge types. The GrT also implements an attention mechanism, which is a function of the neighbourhood connectivity for each node in the graph [45]. The output \(Z \in \mathbb{R}^{M \times D}\) of the \(l\) -th layer of a GrT with \(C\) attention heads is a node representation with the same dimensionality as \(X\) , which can be written as

\[Z^{(l)} = \bigoplus_{i=1}^C \sigma(\tilde{\Delta}_i^{-1} \tilde{A}_i^{(l)} X W) \quad (3.18)\]

where:

• \(\bigoplus\) is the concatenation operator,
• \(\sigma(\cdot)\) is the sigmoid function,
• \(\tilde{\Delta}_i\) is the degree matrix of \(\tilde{A}_i^{(l)}\) , defined as \(\Delta_m = \sum_n A_{mn}\) ,
• \(X\) is the feature matrix,
• \(W \in \mathbb{R}^{D \times D}\) is a trainable weight matrix, and
• \(\tilde{A}_i^{(l)} = A_i^{(l)} + I\) , in which \(A_i^{(l)}\) is the adjacency matrix from the \(i\) -th channel of the metapath tensor \(\mathbf{A}^{(l)} \in \mathbb{R}^{M \times M \times C}\) .

\(\mathbf{A}^{(l)}\) is evaluated as \(\mathbf{A}^{(l)} = \Delta^{-1} \mathbf{Q}_1 \mathbf{Q}_2\) , where \(\mathbf{Q}_1\) and \(\mathbf{Q}_2\) , both \(\in \mathbb{R}^{M \times M \times C}\) , are two adjacency tensors selected according to:

\[\mathbf{Q} = \varphi[\mathbf{A}; \zeta(\mathbf{W}_\varphi)] \quad (3.19)\]

where \(\mathbf{A} \in \mathbb{R}^{M \times M \times K}\) is the adjacency tensor, \(\varphi(\cdot)\) is the convolution operator, \(\zeta(\cdot)\) is the softmax function and \(\mathbf{W}_\varphi \in \mathbb{R}^{C \times C \times K}\) are the weights of \(\varphi\) . \(Z\) contains the node representations from \(C\) different meta-path graphs.

The 3D Projector

The 3D projector is a simple fully connected layer obeying

\[P = Z^{(L)} W_P, \quad (3.20)\]

where \(Z^{(L)}\) is the output of the \(L\) -th layer of the GrT, \(W_P \in \mathbb{R}^{D \times 3}\) is the weight matrix of the projector and \(P \in \mathbb{R}^{M \times 3}\) are the 3D coordinates of the \(M C_\alpha\) atoms in the protein chain. To train the model, a distance map is evaluated for each protein from the predicted coordinates \(P\) as \(\tilde{D}_{ij} = d_{ij}\) , where \(d_{ij} = \|P_i - P_j\|_2\) is the Euclidean distance between the 3D coordinates of the \(i\) -th and \(j\) -th amino acid in \(P\) .

3.7.2 Training details

The total loss to minimize is equal to

\[\mathcal{L} = \mathcal{L}_1 + \mathcal{L}_2 \quad (3.21)\]

in which the first term minimizes the \(L_1\) norm between \(\mathbf{D}\) (the ground truth distance map) and \(\tilde{\mathbf{D}}\) , as

\[\mathcal{L}_1 = \frac{1}{N^2} \sum_i^N \sum_j^N \|\tilde{D}_{ij} - D_{ij}\|_1 \quad (3.22)\]

The second term maximizes the structural similarity index (SSIM) between \(\mathbf{D}\) and \(\tilde{\mathbf{D}}\) weighted by an arbitrary coefficient \(\alpha = 10\) to make \(\mathcal{L}_2\) of the same order of magnitude as \(\mathcal{L}_1\) :

\[\mathcal{L}_2 = \alpha \left( 1 - SSIM\{\mathbf{D}, \tilde{\mathbf{D}}\} \right) = \alpha \left( 1 - \frac{(2\mu_{\tilde{\mathbf{D}}}\mu_{\mathbf{D}} + c_1)(2\sigma_{\tilde{\mathbf{D}}\mathbf{D}} + c_2)}{(\mu_{\tilde{\mathbf{D}}}^2 + \mu_{\mathbf{D}}^2 + c_1)(\sigma_{\tilde{\mathbf{D}}}^2 + \sigma_{\mathbf{D}}^2 + c_2)} \right) \quad (3.23)\]

where \(\mu_{\mathbf{D}}\) is the mean of \(\mathbf{D}\) , \(\mu_{\tilde{\mathbf{D}}}\) the mean of \(\tilde{\mathbf{D}}\) , \(\sigma_{\tilde{\mathbf{D}}\mathbf{D}}\) the covariance of \(\tilde{\mathbf{D}}\) and \(\mathbf{D}\) , \(\sigma_{\tilde{\mathbf{D}}}^2\) the variance of \(\tilde{\mathbf{D}}\) , \(\sigma_{\mathbf{D}}^2\) the variance of \(\mathbf{D}\) , \(c_1 = (k_1 L)^2\) , \(c_2 = (k_2 L)^2\) with \(k_1 = 0.01\) , \(k_2 = 0.03\) and \(L\) the dynamic range is set to 255.

Note how the loss is measured over distance maps and not over 3D coordinates as 3D coordinates depend on a reference frame, while distances are rotationally and translationally invariant. It is possible to include an orientational term in the loss (e.g., MAE between original and predicted cost maps or angle maps), but this would require more coordinates to be predicted rather than just the \(C_\alpha\) coordinates, which are enough for distance maps.

The model consists of 108813 trainable parameters, of which 108648 are from the GrT and 165 are from the projector. The training and testing sets are subsets of PDNET. The first is composed of 200 proteins, while for testing the accuracy of the model we tested it on 5 test sets of 150 proteins each (labelled D1, D2, etc.). The optimizer has been set to Adam with exponentially decaying learning rate, with initial

learning rate \(\eta_0 = 1 \times 10^{-2}\) and decay rate per epoch \(\gamma = 0.9\) . The GrT has \(C = 4\) attention heads and 3 layers. The batch size has been fixed to \(B = 1\) and the network has been trained for \(E = 5\) epochs, for a total of 1000 training iterations.

Combinations of \(\eta \in \{1 \times 10^{-1}, 1 \times 10^{-2}, 1 \times 10^{-3}, 3 \times 10^{-4}\}\) , \(E \in \{3, 5, 10\}\) , \(B \in \{1, 50, 100\}\) , \(L \in \{3, 6, 10\}\) , \(C \in \{1, 4, 5\}\) , as well as \(\mathcal{L} = \mathcal{L}_1\) and additional layers in the 3D projector have also been implemented and tested, but the hyperparameters above were found to be optimal for our problem.

3.7.3 Alignment of predicted coordinates to ground truth

The predicted coordinates \(P\) are relative to a different reference frame compared to \(T\) , the ground truth coordinates from the PDB database. \(P\) and \(T\) must be aligned before evaluating the GDT scores. The goal of alignment is to find \(\mathbf{R} \in \mathbb{R}^{3 \times 3}\) , \(t \in \mathbb{R}^3\) such that \(P = \mathbf{R}T + t\) . In the GA case, we aim at finding the rotor \(R\) such that \(P = RT\tilde{R}\) . We do so by initially centering \(P, T\) at the origin (i.e by placing their centre of mass at the origin) and align the two point clouds via characteristic multivector technique [100].

Results in Table 3.4 are obtained after aligning \(P, T\) . We picked the best GDT scores out of two different alignments procedures, namely (i) singular value decomposition (SVD) and (ii) GA-based alignment via characteristic multivectors. Results in Tables 3.2-3.3, as they refer to 750 proteins, have been obtained only through (i) as it is faster computationally. The two algorithms are summarized in Algorithms 2-3.

3.7.4 Results

For each protein of length \(M\) we obtained an \(M \times 3\) point cloud \(P\) of predicted 3D coordinates. We then aligned the predicted coordinates to the ground truth coordinates \(T\) (obtained from the Protein Data Bank (PDB) [24]) via singular value decomposition (SVD) (see Appendix A) and performed the GDT, and evaluated the GDT_TS (total score) and GDT_HA (half size) between \(P\) and \(T\) as follows:

\[\text{GDT\_TS} = \frac{p_{<1\text{\AA}} + p_{<2\text{\AA}} + p_{<4\text{\AA}} + p_{<8\text{\AA}}}{4} \quad (3.24)\]

\[\text{GDT\_HA} = \frac{p_{<0.5\text{\AA}} + p_{<1\text{\AA}} + p_{<2\text{\AA}} + p_{<4\text{\AA}}}{4} \quad (3.25)\]

where \(p_{<n\text{\AA}}\) indicates the percentage of an amino acid's coordinates in \(P\) whose distance from the corresponding amino acid's coordinates in \(T\) is below \(n \text{ \AA}\) . Results for 5 test sets D1 to D5 are summarized in Tables 3.2-3.3, consistently with prior art [3, 36]. Note

Algorithm 2 SVD Alignment

i ← 0
M ← 103
while i ≤ M do

     $C_P \leftarrow \frac{1}{N} \sum_i^N P^i$                                      ▷ compute centroids
     $C_T \leftarrow \frac{1}{N} \sum_i^N T^i$ 

     $H \leftarrow (P - C_P)(T - C_T)^T$                                 ▷  $H \in \mathbb{R}^{3 \times 3}$ 
     $U, \Sigma, V^T \leftarrow \text{SVD}(H)$ 
                                                    ▷ perform SVD
     $\mathbf{R} \leftarrow VU^T$ 
     $t \leftarrow C_P - \mathbf{R}C_T$ 

     $T \leftarrow \mathbf{R}T + t$                                             ▷ rotate and translate T
    i ← i + 1
end while

that these results are intended to provide a relative comparison of different orientational features on a common, simple dataset like PDNET. They are not meant to compete with current state-of-the-art methods such as AlphaFold, which typically achieve GDT scores \(\geq 80\) [81].

The highest GDT score has been highlighted in bold and the second highest has been underlined.

Table 3.2 GDT_TS scores over the five datasets.

There are several things to note here: firstly, by adding orientational information, an improvement on the protein coordinates quality corresponding to at least \(\sim 4\%\) compared to the approach without orientational information can be measured, as found in [197].

Secondly, for all the analyzed cases (b)-(i) the relative improvement is generally \(< 2\%\) , which implies that no clear superior approach to modeling the amino acid

Algorithm 3 GA Alignment

i ← 0
M ← 103
while i ≤ M do
  CP ←  $\frac{1}{N} \sum_i^N P^i$  ▷ compute centroids
  CT ←  $\frac{1}{N} \sum_i^N T^i$ 

  F ← (P − CP)(T − CT)T
  G ← (T − CT)(T − CT)T ▷ F, G ∈ ℝ3×3

  f1, f2, f3 ← F:,1, F:,2, F:,3 ▷ extract columns of F, G
  g1, g2, g3 ← G:,1, G:,2, G:,3

  f1 ← (f2 ∧ f3)/(f1 ∧ f2 ∧ f3) ▷ reciprocal frames
  f2 ← (f1 ∧ f3)/(f1 ∧ f2 ∧ f3)
  f3 ← (f1 ∧ f2)/(f1 ∧ f2 ∧ f3)

  X ← 1 + [(f1g1 + f2g2 + f3g3)] + [(f2 ∧ f1)(g1 ∧ g2) + (f3 ∧ f2)(g2 ∧ g3) + (f3 ∧ f1)(g1 ∧ g3)] + (f3 ∧ f2 ∧ f1)(g1 ∧ g2 ∧ g3)
  α ← X $\bar{X}$ 
   $\bar{R}$  ← X/√α

  T ← R $\bar{R}$ 
  i ← i + 1
end while

Table 3.3 GDT_HA scores over the five datasets.

orientation exists. This means that, as long as orientational information is added, an improvement is going to be seen regardless of how this information is encoded.

However, the orientational information encoded through our GA-based metrics generally require much less information about the protein backbone than encoding orientation as angle maps. Compare, for example, case (d) with (i) in Table 3.3: in case (d) we have a single dot product map \(\mathbf{N}_{C_\alpha}\) ( \(K = 5\) ), constructed over oriented

Table 3.4 GDT_TS score for 10 example proteins (GDT_HA score in parenthesis)

points at position \(C_\alpha\) in the plane specified by the \(N - C_\alpha - C\) triplet, that yields comparable results to case (i), in which angle maps are three additional features ( \(K = 7\) ) constructed over a total of 5 different planes.

This indicates two things: GA-based metrics can (1) distill more information about the relative orientation of amino acids from fewer geometrical objects, and (2) that this information can be condensed in fewer features, which are equivalent to the 3 angle maps approach.

Examples of GDT scores for 10 proteins are given in Table 3.4. Again, it can be seen that no approach is consistently superior for randomly selected proteins. Examples of predicted coordinates for selected cases are given in Figures 3.19-3.20-3.21.

Geometrical meaning of features

The fact that our GA approach yields comparable results to angle maps by employing fewer features obtained from less geometrical information can be explained with how the information is packed in GA based features: it is possible to establish a relationship between the secondary structure and the patterns in the cost maps. By secondary structure we refer to local folding patterns of a protein, most commonly \(\alpha\) -helices and \(\beta\) -sheets.

We illustrate this relationship by assigning an arbitrary colour to each secondary structure: red to \(\alpha\) -helices, green to \(\beta\) -sheets, blue to turns and white to all the others. In Fig. 3.22 we see how there is almost one-to-one correspondence between colour patches of secondary structures and patterns in cost maps. This shows how cost maps, despite being built starting from a single plane in the backbone ( \(N - C_\alpha - C\) ), contain information about the protein folding that is among the most relevant in PSP pipelines, namely secondary structures.

To the best of our knowledge, this is the first example of a single orientational map that clearly matches the secondary structures.

Figure 3.19 displays the results for protein 2gomA. The original 3D protein model is shown in green, with the original distance map \(D\) below it. The ground truth and predicted coordinates \(T, P\) are given on the top row for selected cases in red and blue, respectively, with their corresponding GDT scores. Below, the distance map \(\tilde{D}\) built from \(P\) and the mean absolute error (MAE) and structural similarity index measure (SSIM) measured with respect to \(D\) .

Case (f):

• GDT_TS = 25.00
• GDT_HA = 9.84
• MAE: 5.83
• SSIM: 0.689

Case (h):

• GDT_TS = 27.46
• GDT_HA = 9.02
• MAE: 3.07
• SSIM: 0.768

Case (l):

• GDT_TS = 25.00
• GDT_HA = 11.07
• MAE: 3.56
• SSIM: 0.695

Figure 3.19: Results for protein 2gomA. The figure displays the original 3D protein model (green ribbon) and its distance map (D) on the left. To the right, three cases are shown: 'case (f)' (red/blue lines), 'case (h)' (red/blue lines), and 'case (l)' (red/blue lines). Each case includes its GDT_TS and GDT_HA scores, a distance map (D-tilde), and its MAE and SSIM scores.

Fig. 3.19 Results for protein 2gomA. The original 3D protein model is shown in green with the original distance map \(D\) below. The ground truth and predicted coordinates \(T, P\) are given on the top row for selected cases in red and blue, respectively, with their corresponding GDT scores. Below, the distance map \(\tilde{D}\) built from \(P\) and the mean absolute error (MAE) and structural similarity index measure (SSIM) measured with respect to \(D\) .

Figure 3.20 displays results for protein 1yqhA. The top row shows the original 3D protein model (red) and the predicted coordinates (blue) for selected cases (f, h, and i). The bottom row shows the distance map \(D\) built from \(P\) and the MAE and SSIM values measured with respect to \(D\) .

Legend for coordinate plots:
— original (blue)
— predicted - case (f) (red)
— predicted - case (h) (red)
— predicted - case (i) (red)

Case (f) scores:
GDT_TS = 23.80
GDT_HA = 6.49
MAE: 3.64
SSIM: 0.746

Case (h) scores:
GDT_TS = 21.15
GDT_HA = 6.49
MAE: 3.93
SSIM: 0.685

Case (i) scores:
GDT_TS = 16.83
GDT_HA = 4.09
MAE: 4.72
SSIM: 0.663

Figure 3.20: Results for protein 1yqhA. The figure displays four panels. The first panel shows the original 3D protein model in red. The subsequent three panels show the original (blue) and predicted (red) coordinates for selected cases (f, h, and i), along with their corresponding GDT scores. Below each coordinate plot is a distance map D (heatmap) and the MAE and SSIM values.

Fig. 3.20 Results for protein 1yqhA. The original 3D protein model is shown in red with the original distance map \(D\) below. The ground truth and predicted coordinates \(T, P\) are given on the top row for selected cases in red and blue, respectively, with their corresponding GDT scores. Below, the distance map \(D\) built from \(P\) and the MAE and SSIM measured with respect to \(D\) .

Figure 3.21 displays results for protein 1z0jB. The top row shows the original 3D protein model (yellow) and the original distance map \(\mathbf{D}\) (bottom left). The middle row shows the ground truth ( \(T\) , blue) and predicted ( \(P\) , red) coordinates for three cases (g, h, i), along with their corresponding GDT scores.