Klas Modin PRO
Mathematician at Chalmers University of Technology and the University of Gothenburg
Riemannian gradient flows in shape analysis
Klas Modin
What can you do with a metric and a function?
- Newton's equation on M
- Riemannian gradient flow on M
\nabla_{\dot x}\dot x = -\nabla V(x)
∇x˙x˙=−∇V(x)
\dot x = -\nabla V(x)
x˙=−∇V(x)
Formal geometric setting
\dot g = -\nabla^{\mathcal G} E(g)
g˙=−∇GE(g)
- Lie group G acting on Q
- Template (source) q0∈Q
- Target q1∈Q
- Riemannian metric G on G
- Symmetric under isotropy subgroup Gg⋅q(g⋅q˙,g⋅q˙)=Gq(q˙,q˙)∀g∈Gq0
- "Energy" functional E(g)=Fq1(g⋅q0)
Def: momentum map J:T∗Q→g∗ for cotangent lifted action
\langle J(q,p),\xi \rangle = \langle p,\xi\cdot q \rangle
⟨J(q,p),ξ⟩=⟨p,ξ⋅q⟩
Proposition: gradient flow is
\dot g = - u(g\cdot q_0)\cdot g
g˙=−u(g⋅q0)⋅g
u(q) := \mathcal A^{-1}J(q,d F_{q_1}(q))
u(q):=A−1J(q,dFq1(q))
Computing the gradient
Simplest special case: right-invariant metric Ge(ξ,ξ)=⟨Aξ,ξ⟩,A:g→g∗
-\nabla^{\mathcal G} E(g)
−∇GE(g)
G
G
Q
Q
e
e
q_0
q0
g(t)
g(t)
q_1
q1
q(t) = g(t)\cdot q_0
q(t)=g(t)⋅q0
Geometric picture
Proposition: gradient flow on OrbG(q0) is
\dot q = - u(q)\cdot q
q˙=−u(q)⋅q
u(q) := \mathcal A^{-1}J(q,d F_{q_1}(q))
u(q):=A−1J(q,dFq1(q))
Induced gradient flow on q0-orbit
G induces metric Gˉ on OrbG(q0)
\dot q = -\nabla^{\bar{\mathcal G}} F_{q_1}(q)
q˙=−∇GˉFq1(q)
Typical form
F_{q_1}(q) = d(s,s_1) + \varepsilon e(r,r_0)
Fq1(q)=d(s,s1)+εe(r,r0)
Choice of Fq1 and regularization
- Shape space S
- Regularization space R
- Take Q=S×R
- Take q0=(s0,r0)
- Take q1=(s1,r0)
S
S
R
R
F_{q_1}(g\cdot q_0) = d(g\cdot s_0,s_1) + \varepsilon e(g\cdot r_0,r_0)
Fq1(g⋅q0)=d(g⋅s0,s1)+εe(g⋅r0,r0)
distance or divergence
r_0
r0
s_0
s0
s_1
s1
Gradient flow
\dot g = - u(g\cdot q_0)\cdot g
g˙=−u(g⋅q0)⋅g
u(q) = \mathcal A^{-1}J(q,F_{q_*}(q))
u(q)=A−1J(q,Fq∗(q))
Lie-Euler method
q_k = g_k\cdot q_0
qk=gk⋅q0
\xi_k = \mathcal{A}^{-1}J(q_k,d F_{q_*})
ξk=A−1J(qk,dFq∗)
g_{k+1} = \exp(-h\xi_k)g_k
gk+1=exp(−hξk)gk
Generic time discretization
Guides object-oriented design of shape analysis software
List of examples
-
Polar decompositions (in optimal transport)
multivariate Gaussians, smooth, Dirac deltas?, entropic regularization? (Schrödinger bridge) -
Toda flow
connections to KdV, QR-algorithm, etc, etc -
QR decomposition
by "silly" entropy gradient flow -
Brockett's flow
new interpretation as entropy gradient flow -
Density matching
using optimal information transport -
LDDMM
also simple regularization using deformation tensor - Shape reconstruction in tomography
horizontal slice
I
I
R
R
fiber
\pi
π
fiber
I
I
W_1
W1
K
K
P(n)
P(n)
A
A
Q
Q
Example: QR decomposition
\mathrm{Hor}_A = \{ V\in T_A\mathrm{GL}(n) \mid \ell(VA^{-1}) = 0 \}
HorA={V∈TAGL(n)∣ℓ(VA−1)=0}
K = \{ R\in \mathrm{GL}(n)\mid \ell(R)=0, R_{ii}>0 \} \Rightarrow T_R K = \mathrm{Hor}_R
K={R∈GL(n)∣ℓ(R)=0,Rii>0}⇒TRK=HorR
- G=GL(n)
- Q=P(n) (sym. pos. def. matrices)
- A such that:
QR example
Manifold of inverse covariance matrices
P(n) = \{ W\in \mathbb{R}^{n\times n}\mid W=W^\top, W>0 \}
P(n)={W∈Rn×n∣W=W⊤,W>0}
\displaystyle p(x;W^{-1}) = \sqrt{\frac{|W|}{(2\pi)^n}}\exp(-\frac{1}{2}x^\top W x)
p(x;W−1)=(2π)n∣W∣exp(−21x⊤Wx)
QR example
Fisher-Rao metric on P(n)
T_W P(n) = \{ U\in \mathbb{R}^{n\times n}\mid U=U^\top \}
TWP(n)={U∈Rn×n∣U=U⊤}
U
U
W
W
\bar{\mathcal G}_W(U,U) = \frac{1}{2}\mathrm{tr}(W^{-1}UW^{-1}U)
GˉW(U,U)=21tr(W−1UW−1U)
QR example
Fisher-Rao invariance
U
U
W
W
g_W(U,U) = \frac{1}{2}\mathrm{tr}(W^{-1}UW^{-1}U)
gW(U,U)=21tr(W−1UW−1U)
\mathrm{GL}(n)\times P(n) \ni (A,W) \mapsto A^\top W A \in P(n)
GL(n)×P(n)∋(A,W)↦A⊤WA∈P(n)
Right action of GL(n) on P(n) is transitive
\bar{\mathcal G}_{A^\top W A}(A^\top U A,A^\top U A) = \bar{\mathcal G}_W(U,U)
GˉA⊤WA(A⊤UA,A⊤UA)=GˉW(U,U)
QR example
Compatible metric on GL(n)
\displaystyle \mathcal G_A(V,V) = \frac{1}{2}\mathrm{tr}\big(\ell(VA^{-1})^\top\ell(VA^{-1})+\sigma(VA^{-1})\sigma(VA^{-1}) \big)
GA(V,V)=21tr(ℓ(VA−1)⊤ℓ(VA−1)+σ(VA−1)σ(VA−1))
V
V
A
A
QR example
Use relative entropy as "energy"
\displaystyle
F_{W_1}(W) = \frac{n}{2}- \frac{1}{2}\mathrm{tr}(W_1W^{-1}) + \frac{1}{2}\log(\det(W_1W^{-1}))
FW1(W)=2n−21tr(W1W−1)+21log(det(W1W−1))
Notice: no regularization used here
QR example
Final "silly QR" flow
\displaystyle \dot R = \frac{1}{2} R^{-\top}(W_1-R^\top R) + ZR, \qquad Z\in\mathfrak{o}(n)
R˙=21R−⊤(W1−R⊤R)+ZR,Z∈o(n)
QR example
Convergence
Convexity lemma:
-\mathrm{Hess}(E)_R = \mathcal G_R
−Hess(E)R=GR
Corollary:
d^2(R(t),R_\infty) \leq \, \mathrm{e}^{-2 t} d^2(R(0),R_\infty)
d2(R(t),R∞)≤e−2td2(R(0),R∞)
QR example
R = \begin{bmatrix}
3 & -1 \\
0 & 2
\end{bmatrix}
R=[30−12]
W_1 = \pi(R) = R^\top R
W1=π(R)=R⊤R
Example
QR example
Example
-\bar H(R(t))
−Hˉ(R(t))
\bar d(R(t),R_\infty)^2
dˉ(R(t),R∞)2
\text{slope of }\exp(-2t)
slope of exp(−2t)
QR example
Example: Brockett's flow
E(Q) = \mathrm{tr}(NQ^\top M Q), \qquad M\in P(n), \; N \in D(n)
E(Q)=tr(NQ⊤MQ),M∈P(n),N∈D(n)
- G=O(n) with canonical metric (Killing form)
- Brockett's energy functional:
Brockett example
- Gives gradient flow:
\dot Q = - Q(NQ^\top M Q - Q^\top MQN)
Q˙=−Q(NQ⊤MQ−Q⊤MQN)
New interpretation
- Action of O(n) on P(n) not transitive
- Energy functional defined by relative entropy
Brockett example
\displaystyle
F_{N}(W) = \frac{n}{2}- \frac{1}{2}\mathrm{tr}(NW^{-1}) + \frac{1}{2}\log(\det(NW^{-1}))
FN(W)=2n−21tr(NW−1)+21log(det(NW−1))
Proposition: Fisher-Rao gradient flow restricted to orbits is
\dot W = [W,[W^{-1},N]]
W˙=[W,[W−1,N]]
Corollary: Expressed in Σ=W−1 we get
\dot \Sigma = [\Sigma,[\Sigma,N]]
Σ˙=[Σ,[Σ,N]]
Double bracket form of Brockett's flow
Example: density matching
- (M,g) compact Riemannian manifold
- G=Diff(M)
- S=Dens(M)={ρ∈C∞(M)∣ρ>0,∫Mρμ=1}
- R=Dens(M), so Q=Dens(M)×Dens(M)
- Action by pushforward (not composition!)
- Right-invariant metric ⟨v,v⟩=∫Mg(v,Δv)μ+⋯
- Energy functional Fρ1(ρ,σ)=dF2(ρ1,ρ)+ϵdF2(1,σ)
Density example
lots of structure!
(with M. Bauer and S. Joshi)
Riemannian submersion
\mathrm{Diff}(M)
Diff(M)
\mathrm{Prob}(M)
Prob(M)
\mathrm{Id}
Id
\mu
μ
\mu_0
μ0
\pi(\varphi)=\varphi^*\mu
π(φ)=φ∗μ
H1 metric
\mathrm{Hor}
Hor
Fisher-Rao metric = explicit geodesics
Density gradient flow
Density example
- Momentum map J:T∗Dens(M)→X(M)∗J(ρ,p)=ρ∇p
Gradient flow on orbits of Dens(M)×Dens(M)
\dot \rho = - v\cdot \nabla \rho - \rho\, \mathrm{div}(v)
ρ˙=−v⋅∇ρ−ρdiv(v)
\dot \sigma = - v\cdot \nabla \sigma - \sigma\, \mathrm{div}(v)
σ˙=−v⋅∇σ−σdiv(v)
\displaystyle
\Delta v = \rho \nabla \frac{\delta F_{\rho_1}}{\delta\rho}(\rho,\sigma) + \varepsilon\sigma \nabla \frac{\delta F_{\rho_1}}{\delta\sigma}(\rho,\sigma)
Δv=ρ∇δρδFρ1(ρ,σ)+εσ∇δσδFρ1(ρ,σ)
Application: draw samples from non-uniform distribution on M
Problem 1: given μ∈Dens(M) generate N samples from μ
Most cases: use Monte-Carlo based methods
Special case here:
- M low dimensional
- μ very non-uniform
- N very large
transport map approach
might be useful
Density example
Optimal transport problem
Problem 2: given μ∈Dens(M) find φ∈Diff(M) minimizing
under constraint φ∗μ0=μ
Studied case: (Moselhy and Marzouk 2012, Reich 2013, ...)
- dist = L2-Wasserstein distance
- ⇒ optimal mass transport problem
- ⇒ solve Monge-Ampere equation (heavily non-linear PDE)
E(\varphi) = \mathrm{dist}(\mathrm{id},\varphi)^2
E(φ)=dist(id,φ)2
Our notion:
- use optimal information transport
Density example
Simple 2D example
Warp computation time (256*256 gridsize, 100 time-steps): ~1s
Sample computation time (10^7 samples): < 1s
Density example
NOT good for image registration
ρ0
ρ1
Density example
...but works with regularization!
ρ0
ρ1
Density example
Application: CT of breathing lung
Data: breathing cycle of rat, CT imaging
Density example
Results
Regularized density flow
LDDMM
Density example
Example: LDDMM
- M Riemannian manifold
- G=Diff(M)
- S=C∞(M)
- R=Diff(M)
- Action by composition
- Right-invariant H1 metric from Met
- Energy functional FI1(I,η)=∣I1−I∣L22+ϵdA2(id,η)
LDDMM example
Why is LDDMM computationally expensive?
Because δηδdA2(id,⋅) is expensive
Deformation regularization
- (M,g) Riemannian manifold
- G=Diff(M)
- S=C∞(M)
- R=Met(M)
- Action by composition and push-forward
- Right-invariant metric defined by kernel
- Energy functional FI1(I,h)=∣I1−I∣L22+ϵdMet2(g,h)
- That is E(φ)=∣I1−I0∘φ−1∣L22+ϵdMet2(g,φ∗g)
LDDMM example
explicit formula (cf. Peter's talk)
(cf. Joshi, Pennec, and others)
deformation
tensor
Resulting gradient flow
- Momentum map J:T∗(C∞(M)×Met(M))→X(M)∗J(I,h,p,α)=p∇I+∣h∣1/2div(α)
Gradient flow on orbit in C∞(M)×Met(M)
\dot I = - v\cdot \nabla I
I˙=−v⋅∇I
\dot \mathsf h = - \mathcal L_v \mathsf h
\displaystyle
\mathcal A v = (I-I_0)\nabla I + \varepsilon \vert\mathsf h\vert^{1/2}\mathrm{div}(\mathsf h - \mathsf g)
Av=(I−I0)∇I+ε∣h∣1/2div(h−g)
LDDMM example
Example: reconstruction
- (M,g) Riemannian manifold
- G=Diff(M)
- S= some shape space (image or density)
- R=Dens(M) or Met(M)
- Forward model T:S→H (think ray-transform)
- Energy functional
Reconstruction example
(with O. Öktem)
E(\varphi) = \vert T(\varphi_*\rho_0) - \mathbf{d}\vert_{\mathcal H}^2 + \varepsilon \mathrm{dist}^2(\mathsf{g},\varphi_*\mathsf{g})
E(φ)=∣T(φ∗ρ0)−d∣H2+εdist2(g,φ∗g)
Resulting gradient flow
Reconstruction example
Gradient flow on orbit in Dens(M)×Met(M)
\dot \rho = - \mathrm{div}(\rho v)
ρ˙=−div(ρv)
\dot \mathsf h = - \mathcal L_v \mathsf h
\displaystyle
\mathcal A v = \rho\nabla T^*(T\rho-\mathbf d) + \varepsilon \vert\mathsf h\vert^{1/2}\mathrm{div}(\mathsf h - \mathsf g)
Av=ρ∇T∗(Tρ−d)+ε∣h∣1/2div(h−g)
Preliminary numerical results
Reconstruction example
THANKS!
References
- M.
Geometry of Matrix Decompositions Seen Through Optimal Transport and Information Geometry
J. Geom. Mech., 9(3):335-390, 2017
- Bauer, Joshi, M.
Diffeomorphic density matching
by optimal information transport
SIAM J. Imaging Sci., 8(3):1718-1751, 2015
Slides available at: slides.com/kmodin
Riemannian gradient flows in shape analysis Klas Modin
Riemannian gradient flows in shape analysis
By Klas Modin
Riemannian gradient flows in shape analysis
Presentation given 2017-11-13 at the Isaac Newton Institute in Cambridge.
- 2,135
