Home / geometry calculators and solvers / Fast and Exact (Poisson) Solvers on Symmetric …

Fast and Exact (Poisson) Solvers on Symmetric … - geometry calculators and solvers

Fast and Exact (Poisson) Solvers on Symmetric …-geometry calculators and solvers

Eurographics Symposium on Geometry Processing 2015 Volume 34 (2015), Number 5
Mirela Ben-Chen and Ligang Liu
(Guest Editors)
Fast and Exact (Poisson) Solvers on Symmetric Geometries
M. Kazhdan
Johns Hopkins University, USA
Figure 1: Applications of our solver: (top) spherical image stitching, followed by Laplacian sharpening, (middle) successive
frames of incompressible fluid simulation on a surface of revolution, and (bottom) successive frames of wave propagation on a
surface of revolution with angular boundary.
In computer graphics, numerous geometry processing applications reduce to the solution of a Poisson equation.
When considering geometries with symmetry, a natural question to consider is whether and how the symmetry can
be leveraged to derive an efficient solver for the underlying system of linear equations. In this work we provide
a simple representation-theoretic analysis that demonstrates how symmetries of the geometry translate into block
diagonalization of the linear operators and we show how this results in efficient linear solvers for surfaces of
revolution with and without angular boundaries.
Categories and Subject Descriptors (according to ACM CCS): I.3.5 [Computer Graphics]: Geometric algorithms,
languages, and systems--Fluid Simulation
c 2015 The Author(s)
Computer Graphics Forum c 2015 The Eurographics Association and John
Wiley & Sons Ltd. Published by John Wiley & Sons Ltd.
M.Kazhdan / Fast and Exact (Poisson) Solvers on Symmetric Geometries
1. Introduction In general, approaches for solving the Poisson equation
can be categorized as either iterative or direct depending on
Solving the Poisson equation is an essential step in many whether an approximate or exact solution is returned.
graphics applications. For example, it enables gradient do-
main image processing, it provides the Hodge decomposi-
tion that supports fluid simulation, and it describes the oscil- Iterative Solvers
lation of a membrane governed by the wave equation (Fig-
ure 1). As such, there has been a large body of research fo- This class of approaches is categorized by algorithms that
cused on efficiently finding exact or approximate solutions. iteratively improve an estimate of the solution. For general
formulations of the Poisson equation, greedy techniques like
In this work, we consider the problem of solving Poisson- Jacobi and Gauss-Seidel [Saa03] relaxation have been used.
like systems on symmetric geometries. On the one hand, While these methods are only guaranteed to improve the so-
these geometries contain sufficient regularity to enable the lution, other methods like conjugate gradients [She94] return
design of an efficient exact solver. On the other, they in- the exact solution in a finite number of iterations.
clude surfaces with varying curvature - so an efficient solver
makes it possible to interactively explore how the surface On their own, these methods are often inefficient because
metric interacts with the Laplace-Beltrami operator. they take too many iterations to produce a reasonable so-
lution. However, these approaches can be integrated with a
The key to our approach derives from the observation that hierarchical multigrid solver to produce accurate answers in
when a linear operator commutes with the symmetry group linear time [BHM00]. Though multigrid was initially used
of the geometry, the linear operator becomes block diagonal for planar lattices, where a hierarchical representation is nat-
in the irreducible representations of the group, with (at least) urally obtained by coarsening the grid, extensions to other
one block associated with each class of irreducible represen- regular geometries [KH10], as well as extensions to triangle
tations. Thus, rather than having to solve one large linear meshes [KCVS98,AKS05,SYBF06,CLB09] and graphs in
system of equations, we can obtain the solution to the linear general [RS87, CFH00] have also been proposed.
system by solving a set of smaller ones.
Taking a representation-theoretic perspective allows us Direct Solvers
to generalize earlier work in the area. Considering di-
hedral symmetry, our work explains Hockney's [Hoc65] Leveraging the fact that the Laplacian is a (semi-)definite
FFT/tridiagonal solver for Poisson equations on rectangu- systems makes it possible to apply methods like Cholesky
lar domains with Neumann/Dirichlet boundaries and pro- factorization [GL96] to express the system as the product of
vides an extension to surfaces of revolution. And, using the upper and lower triangular matrices. Then a solution can be
spherical harmonic decomposition, we show how existing obtained using forward and backward substitution. A chal-
approaches for processing surfaces of revolution using the lenge of using this type of approach is that even though the
structure of circulant matrices [Dav79, BB06] extends to 3D discretization of the Laplacian is often sparse, its inverse
volumes with full rotational symmetry. does not need to be, and computing/storing the Cholesky fac-
torization can become prohibitively expensive. For a recent
We begin our discussion with a brief description of earlier survey of direct solvers in geometry-processing applications,
work on Poisson solvers on Euclidean and non-Euclidean we refer the reader to [BBK05].
domains (?2). Next, we review some basic results from
representation theory and show how these imply block- Spectral solvers the fact that the cosine and sine functions
diagonalizability (?3). We then look in detail at the impli- are eigenvectors of the planar Laplace operator. Thus, a solu-
cations for surfaces of revolution (?4) and also consider the tion to the planar Poisson problem can be obtained by com-
case of volumes with rotational symmetry (?5). We moti- puting the Fourier transform, scaling the Fourier coefficients
vate the utility of our approach by considering applications by the reciprocal of their associated frequency, and then ap-
in spherical image processing, as well as stable fluid simu- plying the inverse Fourier transform [SS88]. As the Fourier
lation, and wave propagation on surfaces of revolution (?6) decomposition on a planar grid can be obtained quickly us-
and conclude by summarizing our work and discussing po- ing the FFT [CT65], this algorithm is quite efficient in prac-
tential directions for future research (?7). tice. Although spectral decompositions of the Laplace opera-
tor has also been used for signal processing on more general
geometries [VL08], this is not a practical method for solv-
2. Related Work ing the Poisson system, as computing the spectral decompo-
The ubiquity of the Poisson equation has made it well- sition is significantly more expensive than solving the sys-
studied across numerous communities. While a comprehen- tem (e.g. the shift-invert implementation requires solving the
sive survey of related methods is beyond the scope of the pa- Poisson equation to compute the spectral decomposition).
per, we briefly review some of the key approaches for solv- In [Hoc65], Hockney describes an approach that is a hy-
ing the system of equations. brid of spectral and direct solvers: Computing the Fourier
c 2015 The Author(s)
Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.
M.Kazhdan / Fast and Exact (Poisson) Solvers on Symmetric Geometries
transform of along just the rows of a planar grid, the Pois- 3.2. Block Diagonalization
son system reduces to a set of tridiagonal systems of equa- Lemma 3.1. Given a representation (,V ), if L is a G-
tions - one for each column of Fourier coefficients. This is linear map then L(V ) V for all .
generalized in [Dav79] which describes how the FFT can be
used to diagonalize (block-)circulant matrices and has been Proof Let W V be an irreducible sub-representations.
used for finding the modes of vibration of surfaces of revo- By Schur's Lemma it follows that L(W ) must also be an ir-
lution [BB06]. Taking a representation theoretic perspective, reducible sub-representation. However, if L(W ) V then,
we extend this approach to general symmetry groups. using Mashke's Theorem, we can obtain a different decom-
position of V into irreducible components, contradicting the
uniqueness of the decomposition.
3. Representation Theory and Block Diagonalization Corollary 3.2. If L is G-linear and {v , . . . , v } is a
1 d
We begin by briefly reviewing some basic representation the- basis for V , then L is block diagonal in the basis
ory ([Ser77, FH91]). Using this theory, we show that if a lin- {v , . . . , v }. In particular, solving the system Lx = b
1 d
ear operator commutes with the symmetry group, the linear can be done by solving || systems each of size d ? d ,
operator can be expressed in block-diagonal form. rather than solving a single system of size d ? d .
At first glance, Corollary 3.2 seems counter-intuitive. On
3.1. Representation Theory Review the one hand, it suggests that as the symmetry group gets
more complicated and the dimensions of the irreducible rep-
Definition Given a finite/compact group G and a complex resentations grow, the diagonal blocks become larger and the
vector space V , a representation of G on V is a map from solver becomes less efficient. On the other hand, we would
G into the automorphisms of V satisfying: expect that as the linear system has more symmetries, there
(e) = Id. should be more opportunity for designing an efficient solver.
We show that this is in fact the case.
(g ? h) = (g) ? (h) g, h G Notation For simplicity, we assume that V = V (as a G-
where e is the identity element in G. linear map must send each V back into itself). We write:
Definition Given a representation (,V ) of G, a subspace m
W V is a sub-representation if (g)(W ) W for all g G. V = Vk
Theorem (Maschke). Given a representation (,V ), if k=1
W1 V is a sub-representation, then there exists W2 V , where the Vk are isomorphic irreducible representations. We
with V = W1 W2 such that W2 is also a sub-representation. set n to be the dimension of Vk. And we set k : V1 Vk to
be the isomorphism between the irreducible representations.
Definition A representation (,V ) is irreducible if the only
sub-representations of V are {0} and V itself. Definition We say that a basis {v1, ? ? ? , v1
1 n, ? ? ? , vm, ? ? ? , vm
1 n },
with vk Vk for all 1 i n and 1 k m, is consistent if:
Definition Given two representations (1,V1) and (2,V2),
a linear map L : V1 V2 is said to be G-linear if: vk
i = k(v1
i ).
L 1(g) = 2(g) L g G. Lemma 3.3. Given a G-linear map L : V V and given a
consistent basis for V , the matrix representation of L in this
If L is an isomorphism, then the irreducible representations basis consists of m ? m blocks where each block is a multiple
are said to be isomorphic. of the identity.
Lemma (Schur). If L is a G-linear map between irreducible Proof It suffices to show that if j : V Vj is the projection
representations then L = 0 or L is an isomorphism. onto Vj then in a consistent basis the matrix expression for:
Corollary (Schur). The only G-linear maps from an irre-
ducible representation into itself are multiples of the identity. j L : Vi Vj
Theorem (Canonical Decomposition). Given a representa- is a multiple of the identity.
tion (,V ), V can be decomposed as the direct sum of irre-
ducible representations (with multiplicity): Using Maschke's theorem, it follows that j is G-linear.
Since L and {k} are also G-linear, and since G-linearity is
preserved under composition and inversion, the map:
V = V with V = Vk
k=1 -1 j L i : V1 V1
is isomor- must be G-linear as well.
are irreducible representations and Vk
where Vk
phic to Vk~
~ if and only if = ~ .
Finally, since V1 is irreducible, the corollary to Schur's
Remark While the isomorphism between V and the direct Lemma implies that the map is a multiple of the identity.
is not unique, the decomposition of V into the Thus, j L : Vi Vj is a multiple of the identity matrix
sum of Vk
direct sum of V is. when expressed in the consistent basis.
c 2015 The Author(s)
Computer Graphics Forum c 2015 The Eurographics Association and John Wiley & Sons Ltd.