Erik Jansson
WASP and KAW postdoctoral fellow in the Cambridge Image Analysis Group and CRA at Wolfson College.
Numerical analysis: the art of being wrong in the mathematically correct way
Concerned with the approximate solution of continuous problems
Geometry: the study of space: distance, shape, size, angles...
Geometric numerics: How to approximate solutions to geometric problems
Shape analysis and matching problems
Structure-preserving numerics for geometric mechanics
Elliptic SPDEs on manifolds
for sampling of random fields
GEOMETRIC
NUMERICAL
METHODS
Paper I and II research questions
Definition: Stochastic process indexed by spatial domain: Assigns random value to each point on domain
We consider Gaussian random fields on surfaces generated by coloring white noise with linear operators that are functions of differential operators.
Classical model inspired by the commonly used Matérn kernel. Whittle–Matérn fields.
Novel, more flexible model proposed by us
Contains Whittle–Matérn fields.
Statistics phrasing: How to sample the fields?
Numerics phrasing: How to approximately solve for \(\mathcal Z\)
Two approaches
Paper I
Paper II
Main scientific question: Convergence of method
Both fundamentally based on surface FEM methods
Main scientific question: Convergence of method
paper IV and VI research questions
But what is a shape?
Point cloud on torus
Function on \([0,1]^2\)
A protein backbone
Shapes are in a (metric) space \(V\) acted upon by a Lie group of deformations \(G\).
\(V = (\mathbb{T}^2)^N\),
\(G = \operatorname{Diff}(\mathbb{T}^2)\)
\(V = C^\infty([0,1]^2)\),
\(G = \operatorname{Diff}([0,1]^2)\)
\(V = \mathbb{R}^{3N}\),
\(G = \operatorname{SO}(3)^N\)
Deformations \(g\) is endpoint of curve \(\gamma \colon [0,1] \to G \)
Right-invariance of metric
What if \(A\) and \(B\) are not in the same space?
Application: Indirect observation of B?
Include a forward operator \(\mathcal{F}\colon V \to W\)
Initial shape: \(A \in V\), target shape: \(B \in W\)
"Reconstructs" the \(\gamma(1).A\) that best would map to \(B\)
Paper IV
Paper VI
Constrain landmarks to move only in certain directions
The reconstruction of protein confirmations is difficult.
Applications cases of (indirect) shape matching
Main scientific question:
Can constrained landmark matching be seen as a neural net?
Main scientific question:
Is indirect shape matching a viable approach for protein conformation reconstruction?
Paper III and VII research questions
Focus here: Lie–Poisson systems: Related to the evolution of observables in Hamiltonian mechanics
Evolves on dual of Lie algebra \(\mathfrak{g}^*\)
Complexified Euler equations:
Rigid body equations:
LP system has a rigid geometric structure:
they evolve on coadjoint orbits
In addition: conserved quantites
They look like onions, sometimes
Discretizations in time for ODEs and time and space for PDEs with LP structure should respect LP structure and conserved quantities
Paper III
Paper VII
Main scientific question: Development, analysis and convergence of method
Main scientific question: How to detect blow-up with geometric numerics?
Paper I: Generate random fields of well-known type on the sphere by SPDE and SFEM approach.
Paper II: Generate random fields of new type by coloring approximated white noise
Paper III: Development of structure-preserving integration of stochastic Lie–Poisson systems
Paper VII: Trustworthy structure preserving numerics for blow-up detection
Paper IV: Sub–Riemannian landmark matching and connections to neural networks
Paper V: Convergence analysis of a gradient flow
Paper VI: Indirect shape matching for protein conformation reconstruction
Paper I: Generate random fields of well-known type on the sphere by SPDE and SFEM approach.
Paper II: Generate random fields of new type by coloring approximated white noise
Paper IV: Sub–Riemannian landmark matching and connections to neural networks
Paper V: Convergence analysis of a gradient flow
Paper VI: Indirect shape matching for protein conformation reconstruction
Paper III: Development of structure-preserving integration of stochastic Lie–Poisson systems
Paper VII: Trustworthy structure preserving numerics for blow-up detection
Paper II: Strong error rates for the algorithms proven using functional calculus approach
Paper III: IsoMP is made stochastic. Method is analyzed by using the geometric structure of LP systems.
Paper VI: Shape matching is adapted to the protein conformation reconstruction setting and is showcased as a promising approach
Paper II results
Our approach is best understood by spectral expansions.
Eigenpairs of \( \mathcal{L} \): \((\lambda_i,e_i)\)
Field approximation using eigenpairs \((\Lambda_{i,h},E_{i,h})\) of \(\mathsf L_h\)
Geometry
approximation
In practice:
Key benefit:
Find a mesh and you can simulate GRFs on it! No need to know the geodesic distance
Research question: Strong error?
\((Z_i)\) are Gaussian with covariance matrix determined by \(\zeta\) and FEM matrices
Key method:
Paper VI results
in 30 seconds (by a maths person)
Electron
microscope
Low dose - low SnR
Many images
Many copies of the same protein but at different orientations and conformations
Proteins in 30 seconds (by a mathematician)
In this work: forget about everything but the \(C_\alpha\)s
Relative positions
\((\mathbb R^3)^N\)
Mathematical model for proteins
Idea: Take a known conformation of the protein as a "prior"
Deform it until it matches target
How to deform? And what energy to use?
Optimization problem!
Shape space: space of relative positions \(V = \mathbb{R}^{3N}\)
Data space - 2D images, \(L^2(\mathbb{R}^2)\) (very noisy!)
We want rigid deformations, so \(G = \operatorname{SO}(3)^N\)
Forward model
Minimization of \(\mathcal E(\nu) \) reduces to:
+ shooting
Gradient is available (even easy to compute)
Deformation
What we want
What we see
Where we start
Research question: It's looking viable
Paper III results
LP systems on \(\mathfrak{g}^*\), stochastic or not, are reductions of Hamiltonian systems on \(T^*G \cong G \times \mathfrak{g}^*\)
There is a mapping \(\mu\) that takes LP system to a left-invariant Hamiltonian system (reconstruction)
upstairs
downstairs
Stochastic Hamiltonian system:
know how to numerically integrate.
Properties of method known
Stochastic Lie–Poisson system:
Don't know how to integrate
Idea: make sure the "stairs" \(\mu\) has nice properties and takes integrator upstairs to integrator downstairs
Intuitive picture:
\(\Phi_h \colon T^*G \to T^* G\): \(G\)-equivariant symplectic integrator \(\rightarrow\)
Reduces to a Lie–Poisson integrator by \(\mu\)
Use Implicit midpoint
Exploit the reconstruction!
Check the literature for error analysis for implicit midpoint
Exploit the reconstruction!
Scientific question is resolved
Idea I: Extend to spatio-temporal fields
Idea II: Compute weak errors
Idea III: Stationarity on general manifolds
Idea IV: Regularization by noise in structure-preserving simulations
Idea V: Blow-up for axisymmetric 3D Euler with Zeitlin
Idea VI: GNNs for better fidelity terms in protein matching for more realistic behavior
Idea VII: Side chain inclusion
Idea VIII: Indirect shape matching via gradient flows
Paper I: Surface finite element approximation of spherical Whittle-Matérn Gaussian random fields. With Annika Lang and Mihály Kovács
Paper II: Non-stationary Gaussian random fields on hypersurfaces: Sampling and strong error analysis. With Annika Lang and Mike Pereira
Paper III: An exponential map free implicit midpoint method for stochastic Lie-Poisson systems. Joint with Annika Lang, Sagy Ephrati and Erwin Luesink
Paper VII: On spectral scaling laws for averaged turbulence on the sphere. Joint with Klas Modin
Paper IV: Sub-Riemannian landmark matching and its interpretation as residual neural networks. With Klas Modin
Paper V: Convergence of the vertical gradient flow for the Gaussian Monge problem. With Klas Modin
Paper VI: Geometric shape matching for recovering protein conformations from single-particle Cryo-EM data. With Ozan Öktem, Jonathan Krook and Klas Modin
Extra slides: Paper I
Details of methods and proof
Error
Main difference: the fractional part is approximated by a Dunford–Taylor integral sinc quadrature!
In general
The error from the geometry discretization is the error of the previous problem in the recursion!
Extra slides: stationarity of RFs on surfaces
G-homogeneous spaces
G-iostropic fields
Geodesic stationarity
Counterexample
\(G\) is a Lie group. Let \(M\) be a G-homogeneous space.
\(G\) acts transitively on \(M\):
For all \(x,y \in M\) there is a \(g \in G\) so that \(g\cdot x = y \)
Examples:
Stationarity = probabilistic invariance under action of group.
Consequence: covariance function depends only on geodesic distance, mean function is constant
Hints at a definition of "geodesic stationarity"
A field is geodesically stationary if its mean function is constant and its covariance functions depends only on the geodesic distance
However- seems to not work!
Feeling: even on a manifold \((M,g)\) that has trivial isometry group, a field generated by \(\mathcal Z = \zeta(\Delta_M) \mathcal W\)
However: covariance function
Insert \(\zeta = \exp(-\Delta_M/2)\)
Heat kernel!
In other words: No hope of only dependence on geodesic distance.
However: asymptotic expansions exist for kernels. There might be some hope to consider only very local expansions to hopefully obtain expressions for the "local covariance function" => maybe the correct notion is one of "local stationarity"?
Well-studied object in Euclidean case.
First step for generalization is not to go to manifolds directly, but to homogenous spaces
(an actual thesis in the strict etymological sense of the word)
Extra slides: surface finite elements
On Euclidean spaces
Just triangulate the domain!
Put the pieces back together, get the original domain!
On manifolds!
Step 1: Triangulate the domain
Issue: Approximate solutions live on \(\mathcal{M}_h\), not \(\mathcal{M}\)!
Put the pieces back together, don't get the same domain!
Step 2: FEM space \(S_h \subset H^1(\mathcal{M}_h)\) of p.w., continuous, linear functions
On manifolds!
Step 1: Triangulate the domain
Issue: Approximate solutions live on \(\mathcal{M}_h\), not \(\mathcal{M}\)!
Put the pieces back together, don't get the same domain!
Step 2: FEM space \(S_h \subset H^1(\mathcal{M}_h)\) of p.w., continuous, linear functions
Step 3: Key tool in surface finite elements: the lift
Takeway: FEM error similar to flat case, up to a "geometry error" term
Extra slides: some notes on kernel functions
To simulate, specify mean and covariance
The sampling question: resolved in isotropic, Euclidean case since the beginning
To simulate, specify mean and covariance
The sampling question: resolved in isotropic, Euclidean case since the beginning
Solution: Hand-craft covariance functions?
(Whittle, 1953) Other, more geometric idea: Matérn field on \(\R^d\) solves:
Write down equivalent SPDE on manifold
The Laplace–Beltrami operator is a geometric object, so this works for all manifolds
(terms and conditions apply)
Sampling by solving the PDE approximately (numerics has entered!)
Extra slides: paper II
Concrete examples
The operator + model details
The FEM operators
Why not compute e.f.s
Actual matrices used
Details on the
(Whittle-) Matérn random fields given by SPDE \((\kappa^2-\Delta_{\mathbb{S}^2})^\beta \mathcal Z = \mathcal{W}\)!
\(\rho_1\) small, \(\rho_2\) large: field is elongated tangentially along level sets of \(f\)
\(\rho_1\) large, \(\rho_2\) small: field is elongated orthogonally along level sets of \(f\)
Eigenpairs of \( \mathcal{L} \): \((\lambda_i,e_i)\)
+ Conditions
\(\mathsf A_\mathcal{M}\) defines elliptic operator \(\mathcal L\)
Discretization of operators
Original
Eigenpairs
Discrete 1
Eigenpairs
Discrete 2
Eigenpairs
Error analysis:
Problem:
Generally: Approximating eigenfunctions are hard!
Eigenvalues are fine, however!
Generally: Approximating eigenfunctions are hard! Multiplicity!
From D. Boffi, Acta Num. (2010)
Chebyshev quadrature approximation of \(\gamma\)
Various white noise approximations, various approximations of \(\mathcal Z\)
sd
sd
sd
??
Important result: errors on the form
\(\|\gamma(\mathcal L)f -\gamma(\mathcal L_h) f\|_{L^2}\)
Brief note on proof:
Extra slides: paper III
Geometry of LP systems
Isospectrality
geometric structure of stochastic LP
why not Itô?
proof of convergence, some more details
Important geometric structure:
LP Systems evolve on coadjoint orbits
(symplectic manifolds that foliate space!)
+ Several preserved quantities, Hamiltonian, Casimirs etc.
Coadjoint orbits
Dynamics on orbit
\(G \subset \operatorname{GL}(n)\) is compact, simply connected and \(J\)-quadratic, i.e., \(A \in \mathfrak{g} \iff AJ + JA^* = 0 \)
Examples: \(\mathfrak{so}(N), \mathfrak{su}(N), \mathfrak{sp}(N)\).
In this setting, LP flow is isospectral:
Why Strato?
Why these coefficients?
We need a chain rule
Same bracket-type coefficient: same geometric structure
Stochastic LP Systems evolve on coadjoint orbits and have Casimirs!
What to assume to prove the existence of solutions?
E.g. Lipschitz is not realistic to assume
But! Some smoothness is sufficient for global existence
System remains on compact level sets of Casimirs
\(\rightarrow\) truncation argument to prove existence
Stochastic LP Systems evolve on coadjoint orbits and have Casimirs!
Coadjoint orbits
With transport noise
With non-transport noise
What happens if we apply an off the shelf integrator?
Non-physical behavior!
What happens with structure-preserving integrators?
Intuitive picture: implicit midpoint in disguise
unresolved unscientific question: if implicit midpoint could wear a disguise, how would it wear it?
Note: No a priori guarantee that this system remains on \(T^*G\)!
Embedded into \((R^{n \times n})^2\) but remains on \(T^*G\)
\(X_t = \mu(Q_t,P_t) \in \mathfrak{g}^*\) satisfies
Stochastic isospectral midpoint:
Take \(Q_n,P_n\) such that \(\mu(Q_n,P_n) = X_n\)
Take \(X_{n+1} = \mu(\Phi_h(Q_n,P_n))\)
1) Show that the implicit midpoint method does not travel too far:
\(\sup_{h \geq 0} \sup_{n \geq 0} \|Q_n,P_n\| \leq R(Q_0,P_0). \)
2) Truncate the CHS system on \(T^*G\)
Upstairs
Downstairs
3) Apply error analysis from literature
Reduction by \(\mu\)
4) Convergence of method for truncated LP system
5) Truncated LP system = LP system
Recipe:
Geometric structure of equations and its preservation is central. Used to prove
Without structure preservation, no convergence guarantees
Extra slides: IsoMP vs other methods
Extra slides: generalities on shape matching
Idea: Use a distance \(d_V\colon V \times V \to \mathbb R\) on \(V\), look for \(g \in G\) that solves
How to match \(A \in V \) with \(B \in V\)?
Problem: In general, no \(g\) takes \(A\) to \(B\) \(\implies\)
more and more complex warps of \(A\)
(the action is non-transitive).
Add a regularization term!
Penalize strange deformations
How to compute \(d_G(e,g)\)?
Pricey! An optimization problem in itself!
The simpler way: Let \(d_G\) be the distance function of a right-invariant Riemannian metric and generate \(g\) as a curve! (Geometry!)
\({\color{#a64d79}{\nu}}\colon [0,1] \to \mathfrak{g}\) is a curve in the vector space \(\mathfrak{g}\)!
Let right-invariance strike!
An optimization problem over curves in a vector space!
\({\color{#a64d79}{\nu}}\colon [0,1] \to \mathfrak{g}\) is a curve in the vector space \(\mathfrak{g}\)!
Let right-invariance strike!
Flexible choice
Right-invariant metric (not so flexible)
Extra slides: paper IV
In usual LM, all vector fields are available
Idea: Constrain set of vector fields
Sub-Riemannian landmark matching!
Landmark matching:
\(V = (\mathbb{T}^2)^N\),
\(G = \operatorname{Diff}(\mathbb{T}^2)\)
\(\mathcal{S}\)
\(v = F(u)(x)\)
\(\mathfrak{X}(M)\)
\(F(u)\)
\(\mathcal{U}\)
\(u\)
\(\mathcal{S}\)
\(v = F(u)(x)\)
\(\mathfrak{X}(M)\)
\(F(u)\)
\(\mathcal{U}\)
\(u\)
Dynamics of \(u\) are determined from the dynamics of \(v\)
Instead of optimizing over vector fields, optimize over \(\mathcal{U}\)
\(\mathcal{U}\) is finite dimensional \(\implies\) discretization
How to find dynamics of \(u\)?
Idea: Apply chain rule to \(\dot m = \operatorname{ad}_v^T m\)
\(\mathcal{S}\)
\(v = F(u)(x)\)
\(\mathfrak{X}(M)\)
\(F(u)\)
\(\mathcal{U}\)
\(u\)
Infinite number of layers \(\implies\)
Time-continuous optimal control problem
Skip connection \(\implies \) Explicit Euler
When matching landmarks, points are moved by a vector field parametrized by some control variable.
A neural networks moves points along a vector field determined by weights and biases.
Extra slides: paper V
Geometric structure
In brief: Gradient flow horizontally or vertically
(More info: Modin, 2017 and references therein)
Gaussian Monge Problem:
\(\mu_0\) and \(\mu_1\) are both (zero-mean) normal distributions on \(\mathbb{R}^n\).
Normal distributions \(\cong\) \(P(n)\), positive-definite symmetric matrices
In brief: Gradient flow horizontally or vertically
How to prove convergence?
Idea: Show \(\frac{\mathrm d} {\mathrm d t} J \to 0\), and that this means we hit polar cone
Extra slides: paper VI
Shooting,
Gradient
Quantitative results
Matching by shooting
Recall example II from intro: reduction to dynamic formulation is made easy by right-invariance
Gradient is available
x
y
z
Deformed
Target
Template
Increase noise
Increase #projections
Quantative measure: Procrustes score
Extra slides: paper VII
Euler equations
Zeitlin
Complexified Euler equations
Simulations for Blow-up
The \(L^2\)-geodesic equation on \(\operatorname{Diff}_{\mu}(\mathbb{S}^2)\):
Euler eqs: Lie-Poisson system on \(C^\infty_0(\mathbb{S}^2) \cong (\mathfrak{X}_\mu(\mathbb{S}^2))^*\)
=> Sphere is Kähler
Initially for the torus, but we work on the sphere!
Goal: Find a mapping \(T_N:C_0^\infty(\mathbb{S^2}) \to \mathfrak{su}(N)\) such that:
Hoppe, 1989:
How to find \(T_{l,m}^N\)?
Coordinate function representation of \(\mathfrak{so}(3)\) in \((C^\infty_0(\mathbb{S}^2),\{\cdot,\cdot\})\).
Spin \(s = (N-1)/2\)-representation of \(\mathfrak{so}(3)\) in \(\mathfrak{u}(N)\)
Hoppe–Yau Laplacian:
Hoppe–Yau (1998):
The first \(N^2\) eigenvalues of \(\Delta\) coincides with those of \(\Delta_N\)
\(T_{l,m}^N\) are eigenfunctions to \(\Delta_N\)
Quantized system on \(\mathfrak{su}(N)\):
Continuous system on \(C^\infty_0(\mathbb{S}^2)\):
Canonical expression (EP-equations):
Relationship between \(\omega\) and \(\psi \)?
Consider complexified Euler equations as LP system on complexified functions:
Relationship between \(\omega\) and \(\psi \)?
Determined by Hamiltonian!
Math
Proven to blow up in finite time for toroidal case. Conjecture same is true for numerics
How to detect in numerics?
Grows very rapidly
How to be sure that the observed behavior is not a numerical artifact?
Growth of norm!
Check spectral deviation
Check reversibility
Check mean time to blow up wrt step size -stabilizes
Check mean time to blow up wrt matrix size
Linear growth after transient phase?
Linear rate vs matrix size - linear relationship
Extra slides: notes on Zeitlin
Why not the torus?
Generator of isometries:
Generator of isometries:
By Erik Jansson
WASP and KAW postdoctoral fellow in the Cambridge Image Analysis Group and CRA at Wolfson College.