# Monotone projection lower bounds from extended formulation lower bounds

Joshua A. Grochow111Department of Computer Science, University of Colorado at
Boulder, Boulder, CO, USA and Santa Fe Institute, Santa Fe, NM, USA,
<EMAIL_ADDRESS>
###### Abstract

In this short note, we show that the Hamiltonian Cycle polynomial is not a
monotone sub-exponential-size projection of the permanent; this both rules out
a natural attempt at a monotone lower bound on the Boolean permanent, and
shows that the permanent is _not_ complete for non-negative polynomials in
$\mathsf{VNP}_{\mathbb{R}}$ under monotone p-projections. We also show that
the cut polynomials, $\sum_{A\subseteq[n]}\prod_{i\in A,j\notin
A}x_{ij}^{q-1}$, and the perfect matching polynomial (or “unsigned Pfaffian”)
$\frac{1}{2^{n}n!}\sum_{\pi\in S_{2n}}\prod_{i=1}^{n}x_{\pi(2i-1),\pi(2i)}$
are not monotone p-projections of the permanent. The latter can be interpreted
as saying that there is no monotone projection reduction from counting perfect
matchings in general graphs to counting perfect matchings in bipartite graphs,
putting at least one theorem behind the well-established intuition. As the
permanent is universal for monotone formulas, these results also imply
exponential lower bounds on the monotone formula size and monotone circuit
size of these polynomials. To prove these results we introduce a new
connection between monotone projections of polynomials and extended
formulations of linear programs that may have further applications.

##  1 Introduction

The permanent

$\operatorname{perm}_{n}(X)=\sum_{\pi\in S_{n}}x_{1,\pi(1)}x_{2,\pi(2)}\dotsb
x_{n,\pi(n)}$

(where $S_{n}$ denotes the symmetric group of all permutations of
$\\{1,\dotsc,n\\}$) has long fascinated combinatorialists [13, 25, 14], more
recently physicists [26, 2], and since Valiant’s seminal paper [22], has also
been a key object of study in computational complexity. Despite its beauty,
the permanent has some computational quirks: in particular, although the
permanent of integer matrices is $\mathsf{\\#P}$-complete and the permanent is
$\mathsf{VNP}$-complete in characteristic zero, the permanent _mod 2_ is the
same as the determinant, and hence can easily be computed. In fact, computing
the permanent mod $2^{k}$ is easy for any $k$ [22], though the proof is more
involved. Modulo any other number $n$, the permanent of integer matrices is
$\mathsf{Mod}_{n}\mathsf{P}$-complete.

In contrast, the seemingly similar Hamiltonian Cycle polynomial,

$HC_{n}(X)=\sum_{\text{$n$-cycles }\sigma}x_{1,\sigma(1)}x_{2,\sigma(2)}\dotsb
x_{n,\sigma(n)},$

where the sum is only over $n$-cycles rather than over all permutations, does
not have these quirks: The Hamiltonian Cycle polynomial is
$\mathsf{VNP}$-complete over any ring $R$ [21] and
$\mathsf{Mod}_{n}\mathsf{P}$-complete for all $n$ (that is, counting
Hamiltonian cycles is complete for these Boolean counting classes).

Jukna [11] observed that, over the Boolean semi-ring, if the Hamiltonian Cycle
polynomial were a monotone p-projection of the permanent, there would be a
$2^{n^{\Omega(1)}}$ lower bound on monotone circuits computing the permanent,
a lower bound that still remains open (the current record is still Razborov’s
$n^{\Omega(\log n)}$ [15]). Even over the real numbers, such a monotone
p-projection would give an alternative proof of a $2^{n^{\Omega(1)}}$ lower
bound on the permanent (Jerrum and Snir [9] already showed the permanent
requires monotone circuits of size $2^{\Omega(n)}$ over $\mathbb{R}$ and over
the tropical $(\min,+)$ semi-ring). Here we show that no such monotone
reduction exists—over $\mathbb{R}$, nor over the tropical semi-ring, nor over
the Boolean semi-ring—by connecting monotone p-projections to extended
formulations of linear programs.

We use the same technique to show that the perfect matching polynomial or
“unsigned Pfaffian”

$\frac{1}{2^{n}n!}\sum_{\pi\in
S_{2n}}\prod_{i=1}^{n}x_{\pi(2i-1),\pi(2i)}=\sum_{\begin{subarray}{c}\pi\in
S_{2n}\\\ \pi(1)<\pi(3)<\dotsb<\pi(2n-1)\\\ \pi(2k-1)<\pi(2k)\;\forall
k\end{subarray}}\prod_{i=1}^{n}x_{\pi(2i-1),\pi(2i)}$

is not a monotone p-projection of the permanent. As the perfect matching
polynomial counts perfect matchings in a general graph, and the permanent
counts perfect matchings in a bipartite graph, we have:

######  Corollary 1.1.

Any efficient projection reduction from counting perfect matchings in general
graphs to counting perfect matchings in bipartite graphs must be non-monotone
(and, therefore, seemingly quite unintuitive!).

######  Remark 1.2 (On the Boolean semi-ring).

Our results also hold for _formal polynomials_ over the Boolean semi-ring
$\mathbb{B}=(\\{0,1\\},\vee,\wedge)$. Over the Boolean semi-ring, the
permanent is the indicator function of the existence of a perfect matching in
a bipartite graph, and the unsigned Pfaffian is the indicator function of a
perfect matching in a general graph. However, over $\mathbb{B}$, each function
is represented by more than one formal polynomial, and we do not yet know how
to extend our results to the setting of _functions_ over $\mathbb{B}$. See
Section 5 for details and specific questions.

We also use the same technique to show that the cut polynomials
$\operatorname{Cut}^{q}=\sum_{A\subseteq[n]}\prod_{i\in A,j\notin
A}x_{ij}^{q-1}$ are not monotone p-projections of the permanent. Perhaps the
main complexity-theoretic interest in the cut polynomials is that
$\operatorname{Cut}^{q}$ _over the finite field $\mathbb{F}_{q}$_ was (until
recently [12]) the only known example of a natural polynomial that is neither
in $\mathsf{VP}_{\mathbb{F}_{q}}$ nor $\mathsf{VNP}_{\mathbb{F}_{q}}$-complete
under a standard complexity-theoretic assumption (that $\mathsf{PH}$ doesn’t
collapse) [3]; there it was also shown that if
$\mathsf{VP}_{\mathbb{F}_{q}}\neq\mathsf{VNP}_{\mathbb{F}_{q}}$ then such
polynomials of intermediate complexity must exist. In that paper, it was asked
whether the cut polynomials, considered as polynomials _over the rationals_ ,
were $\mathsf{VNP}_{\mathbb{Q}}$-complete.222$\operatorname{Cut}^{2}$ was
subsequently shown to be $\mathsf{VNP}_{\mathbb{Q}}$-complete under circuit
reductions [5]; its completeness under projections remains open. Although our
results don’t touch on this question, these previous results motivate the
study of these polynomials over $\mathbb{Q}$.

Because the permanent is universal for monotone formulas, our lower bounds
also imply exponential lower bounds on the monotone algebraic formula
size—and, by balancing algebraic circuits, monotone algebraic circuit size—of
these polynomials; see Section 4.2.

Finally, we note that our results shed a little more light on the
complicatedness of the known $\mathsf{VNP}$-completeness proofs for the
permanent [22, 1]. Namely, prior to our result, the fact that the permanent is
not hard modulo $2$ already implied that any completeness result must use $2$
in a “bad” way: for example, dividing by 2 somewhere. This is indeed true of
both Valiant’s original proof [22] and Aaronson’s independent quantum linear-
optics proof [1]. One might hope for a classical analogue of Aaronson’s
quantum proof, using the characterization of $\mathsf{BPP}$ in terms of
stochastic matrices as a replacement for the characterization of
$\mathsf{BQP}$ using unitary matrices. However, our result says that any
completeness proof for the permanent must use non-monotone reductions, so such
a classical analogue is not possible:

######  Corollary 1.3.

Aaronson’s quantum linear optics proof [1] that the permanent is
$\mathsf{\\#P}$-hard cannot be replaced by one using classical randomized
algorithms in place of quantum algorithms.

Our results also imply the necessity of the use of negative numbers in
Valiant’s $4\times 4$ gadget [22, p. 195]. In light of these results,
Valiant’s $4\times 4$ gadget may perhaps seem less mysterious than the fact
that such a gadget exists that is _only_ $4\times 4$!

To prove these results, we show that a monotone projection between non-
negative polynomials essentially implies that the Newton polytope of one
polynomial is an extension of the Newton polytope of the other (Lemma 3.1),
and then apply known lower bounds on the extension complexity of certain
polytopes. We hope that the connection between Newton polytopes, monotone
projections, and extended formulations finds further use.

##  2 Preliminaries

A polynomial $f(x_{1},\dotsc,x_{n})$ is a _(simple) projection_ of a
polynomial $g(y_{1},\dotsc,y_{m})$ if $f$ can be constructed from $g$ by
replacing each $y_{i}$ with a constant or with some $x_{j}$. The polynomial
$f$ is an _affine projection_ of $g$ if $f$ can be constructed from $g$ by
replacing each $y_{i}$ with an affine linear function $\pi_{i}(\vec{x})$. When
we say “projection” we mean simple projection. Given two families of
polynomials $(f_{n})$, $(g_{n})$, if there is a function $p(n)$ such that
$f_{n}$ is a projection of $g_{p(n)}$ for all sufficiently large $n$, then we
say that $(f_{n})$ is a projection of $(g_{n})$ with _blow-up_ $p(n)$. If
$(f_{n})$ is a projection of $(g_{n})$ with polynomial blow-up, we say that
$(f_{n})$ is a _p-projection_ of $(g_{n})$.

Over any subring of $\mathbb{R}$—or more generally any totally ordered semi-
ring (see below)—a _monotone projection_ is a projection in which all
constants appearing in the projection are non-negative. Monotone p-projection
is defined analogously.

To each monomial $x_{1}^{e_{1}}\dotsb x_{n}^{e_{n}}$ we associate its exponent
vector $(e_{1},\dotsc,e_{n})$, as a point in
$\mathbb{N}^{n}\subseteq\mathbb{R}^{n}$. We then have:

######  Definition 2.1.

The _Newton polytope_ of a polynomial $f(x_{1},\dotsc,x_{n})$, denoted
$\operatorname{Newt}(f)$, is the convex hull in $\mathbb{R}^{n}$ of the
exponent vectors of all monomials appearing in $f$ with non-zero coefficient.

A polytope is _integral_ if all its vertices have integer coordinates; note
that Newton polytopes are always integral. A _face_ of a polytope $P$ is the
intersection of $P$ with a linear space $L$ such that none of the interior
points of $P$ lie in $L$.

For a polytope $P$, let $c(P)$ denote the “complexity” of $P$, as measured by
the minimal number of linear inequalities needed to define $P$ (equivalently,
the number of faces of $P$ of dimension $\dim P-1$). A polytope
$Q\subseteq\mathbb{R}^{m}$ is an _extension_ of $P\subseteq\mathbb{R}^{n}$ if
there is an affine linear map $\ell\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ such
that $\ell(Q)=P$. The _extension complexity_ of $P$, denoted $xc(P)$, is the
minimum complexity of any extension of $P$ (of any dimension):
$xc(P)=\min\\{c(Q)|Q\text{ is an extension of }P\\}$.

The $m$-th _cycle cover polytope_ (also known as the bipartite perfect
matching polytope) is the convex hull in $\mathbb{R}^{m^{2}}$ of the
$\\{0,1\\}$-indicator functions of the directed cycle covers of the complete
directed graph with self-loops on $m$ vertices. The cycle cover polytope is
the Newton polytope of the permanent, as each monomial in the permanent
corresponds to such a cycle cover.

A totally ordered semi-ring (we only consider commutative ones here) is a
totally ordered set together with two operations, denoted
$(R,\leq,\times,+,1,0)$ such that $(R,\times,1)$ and $(R,+,0)$ are both
commutative monoids, $\times$ distributes over $+$, $0\times a=0$ for all $a$,
$a+c\leq b+c$ whenever $a\leq b$, and $ac\leq bc$ whenever $a\leq b$ and
$0\leq c$. An element $c$ of a totally ordered semi-ring is non-negative if
$0\leq c$, and is positive if furthermore $c\neq 0$. We will restrict our
attention to _non-zero_ totally ordered semi-rings; equivalently, we assume
$1\neq 0$. (There is only one polynomial over the zero semi-ring—the zero
polynomial—so that case was of no interest anyways.)

Note that (non-zero) totally ordered semi-rings always have “characteristic
zero:” $1+\dotsb+1\neq 0$ (any number of times). In a non-zero totally ordered
semi-ring, either $0<1$ or $1<0$; we handle the first case, the second being
analogous. If $0<1$, then by adding $1$ to both sides $k$ times we get that
$\sum_{i\in[k]}1<\sum_{i\in[k+1]}1$. If the ring had nonzero characteristic,
then for some $k>1$, the right-hand side here would become zero, and we would
get $0<1<1+1<\dotsb<0$, a contradiction.

The following totally ordered semi-rings are of particular interest:

  * •

The real numbers with its usual ordering and algebraic operations
$(\mathbb{R},\leq,\times,+)$

  * •

The so-called “tropical semi-ring” $(\mathbb{R},\leq,+,\min)$, which is the
real numbers with its usual ordering, where the product is taken to be real
addition and the addition operation is taken to be the minimum.

  * •

The Boolean and-or semi-ring $\mathbb{B}=(\\{0,1\\},\leq,\wedge,\vee)$, where
$0\leq 1$.

To get a feel for the latter two semi-rings, note that polynomials over the
tropical semi-ring generally compute some optimization problem, and over
$\mathbb{B}$ generally compute a decision problem. For example, the
Hamiltonian Cycle polynomial over the tropical semi-ring computes the
Traveling Salesperson Problem, and over $\mathbb{B}$ is the indicator function
of the existence of a Hamiltonian cycle. Note that over $\mathbb{R}$, if two
formal polynomials compute the same function then they must be identical, but
this is not true over the tropical nor Boolean semi-rings.

##  3 Main Lemma

######  Lemma 3.1.

Let $R$ be a totally ordered semi-ring, and let $f(x_{1},\dotsc,x_{n})$ and
$g(y_{1},\dotsc,y_{m})$ be polynomials over $R$ with non-negative
coefficients. If $f$ is a monotone projection of $g$, then some face of
$\operatorname{Newt}(g)$ is an extension of $\operatorname{Newt}(f)$. In
particular, $xc(\operatorname{Newt}(f))\leq c(\operatorname{Newt}(g))$.

###### Proof.

Under simple projections, a monomial in the $y$’s maps to some scalar multiple
of a monomial in the $x$’s (possibly the empty monomial, resulting in a
constant term, or possibly the zero multiple, resulting in zero). Let $\pi$ be
a monotone projection map, defined on the variables $y_{i}$, and extended
naturally to monomials and terms in the $y$’s. (Recall that a _term_ of a
polynomial is a monomial together with its coefficient.) Since each term $t$
of $g$ is a monomial multiplied by a positive coefficient, and since $\pi$ is
non-negative, $\pi(t)$ is either zero or a single monomial in the $x$’s with
nonzero coefficient. The former situation can happen only if $t$ contains some
variable $y_{i}$ such that $\pi(y_{i})=0$. Let $\ker(\pi)$ denote the set
$\\{y_{i}|\pi(y_{i})=0\\}$. Thus, for every term $t$ of $g$ that is disjoint
from $\ker(\pi)$, $\pi(t)$ actually appears in $f$—possibly with a different
coefficient, but still non-zero—since no two terms can cancel under projection
by $\pi$.

Let $e_{1},\dotsc,e_{m}$ be the coordinates on $\mathbb{R}^{m}$, the ambient
space of $\operatorname{Newt}(g)$. Let $K$ denote the subspace of
$\mathbb{R}^{m}$ defined by the equations $e_{i}=0$ for each $i$ such that
$y_{i}\in\ker(\pi)$. Let $P$ be the intersection of $\operatorname{Newt}(g)$
with $K$, considered as a polytope in $K$; since all vertices of
$\operatorname{Newt}(g)$ are non-negative, intersecting
$\operatorname{Newt}(g)$ with a coordinate hyperplane, $e_{i}=0$, results in a
face of $\operatorname{Newt}(g)$, and thus $P$ is a face of
$\operatorname{Newt}(g)$. Note that $P$ is exactly the convex hull of the
exponent vectors of monomials in $g$ that are disjoint from $\ker(\pi)$. In
particular, since $\pi$ is multiplicative on monomials, it induces a _linear_
map $\ell_{\pi}$ from $K$ to $\mathbb{R}^{n}$ (the ambient space of
$\operatorname{Newt}(f)$). By the previous paragraph, the exponent vectors of
$f$ are exactly $\ell_{\pi}$ applied to the exponent vectors of monomials in
$g$ that are disjoint from $\ker(\pi)$. By the linearity of $\ell_{\pi}$ and
the convexity of $P$ and $\operatorname{Newt}(f)$, we have that
$\operatorname{Newt}(f)=\ell_{\pi}(P)$, so $P$ is an extension of
$\operatorname{Newt}(f)$. Since $P$ is defined by intersecting
$\operatorname{Newt}(g)$ with additional linear equations, the lemma follows.
∎

Several partial converses to our Main Lemma also hold. Perhaps the most
natural and interesting of these is:

######  Observation 3.2.

Let $R$ be a totally ordered semi-ring. Given any sequence of integral
polytopes $(P_{n}\subseteq\mathbb{R}^{n})$ such that the $\text{poly}(n)$-th
cycle cover polytope is an extension of $P_{n}$ along an affine linear map
$\ell_{n}\colon\mathbb{R}^{\text{poly}(n)}\to\mathbb{R}^{n}$ with integer
coefficients of polynomial bit-length, there is a sequence of polynomials
$(f_{n})\in\mathsf{VNP}_{R}$ such that $\operatorname{Newt}(f_{n})=P_{n}$ and
$f$ is a monotone p-projection of the permanent.

###### Proof.

Let $C_{m}$ denote the $m$-th cycle cover polytope, let $m(n)$ be a polynomial
such that $C_{m(n)}$ is an extended formulation of $P_{n}$, and let $b(n)$ be
a polynomial upper bound on the bit-length of the coefficients of $\ell_{n}$.
Let $V_{m}$ denote the vertex set of the cycle cover polytope, i.e. the
incidence vectors of cycle covers. Define $f_{n}$ as $\sum_{\vec{e}\in
V_{m}}\vec{y}^{\ell_{n}(\vec{e})}$, where $\vec{y}=(y_{1},\dotsc,y_{n})$, and
the vector notation $\vec{y}^{\vec{e}^{\prime}}$ is defined as
$y_{1}^{e_{1}^{\prime}}y_{2}^{e_{2}^{\prime}}\dotsb y_{n}^{e_{n}^{\prime}}$.
Note that $\ell_{n}$ has only integer coefficients by assumption, and each
$\vec{e}\in V_{m}$ is integral, so the vector $\ell_{n}(\vec{e})$ is also
integral, and the above expression is well-defined. By construction, every
exponent vector of $f_{n}$ is in $\ell_{n}(C_{m(n)})=P_{n}$. Conversely, every
vertex of $P_{n}$ is an exponent vector of $f_{n}$, for any non-zero totally
ordered semi-ring has characteristic zero (see Section 2). (Without noting
this, it would be possible that $k$ distinct vertices in $V_{m}$ would get
mapped to the same point under $\ell_{n}$, for $k$ a multiple of the
characteristic, and then the corresponding monomials
$\vec{y}^{\ell_{n}(\vec{e})}$ would add up to $0$ in $f_{n}$.) Thus
$\operatorname{Newt}(f_{n})=P_{n}$. Furthermore, $f_{n}$ is a monotone
_nonlinear_ projection of the permanent using the map
$x_{ij}\mapsto\vec{y}^{\ell((0,0,\dotsc,1,\dotsc,0))}$, where the $1$ is in
the $(i,j)$ position. Using the universality of the permanent and repeated
squaring, this can easily be turned into a monotone _simple_ projection of the
permanent of size $\text{poly}(m(n),b(n))$. ∎

This can be generalized from the cycle cover polytopes and the permanent to
arbitrary integral polytopes and the natural associated polynomial (the sum
over all monomials whose exponent vectors are vertices of the polytope), but
at the price of using “monomial projections”—in which each variable is
replaced by a monomial—rather than simple projections. There ought to be a
version of this observation over sufficiently large fields and allowing
rational coefficients in $\ell$, using Strassen’s division trick [20], but the
only such versions the author could come up with had so many hypotheses as to
seem uninteresting.

##  4 Applications

###  4.1 Projection Lower Bounds

######  Remark 4.1.

The following theorems hold over any totally ordered semi-ring, including the
Boolean and-or semi-ring, the non-negative real numbers under multiplication
and addition, and the tropical semi-ring of real numbers under addition and
$\min$. To see that this introduces no additional difficulty, note that over
any totally ordered semi-ring $R$, the Newton polytope of a polynomial over
$R$ is still a polytope in a vector space over the real numbers, so standard
results on polytopes and the cited results on extension complexity still
apply.

######  Theorem 4.2.

Over any totally ordered semi-ring, the Hamiltonian Cycle polynomial is not a
monotone affine p-projection of the permanent; in fact, any monotone affine
projection from the permanent to the Hamiltonian Cycle polynomial has blow-up
at least $2^{\Omega(n)}$.

###### Proof.

First, recall that if an $n$-variable polynomial is an affine projection of
the $m\times m$ permanent, then it is a simple projection of the
$(n+1)m\times(n+1)m$ permanent. For completeness we recall the brief proof:
Let $\pi_{ij}(\vec{x})$ be the affine linear function corresponding to the
variable $y_{ij}$ of the $m\times m$ permanent, and write
$\pi_{ij}=a_{0}+a_{1}x_{1}+\dotsb+a_{n}x_{n}$. Let $G$ be the complete
directed graph with loops on $m$ vertices and edge weights $y_{ij}$. Replace
the edge $(i,j)$ by $n+1$ parallel edges with weights
$a_{0},a_{1}x_{1},\dotsb,a_{n}x_{n}$. Add a new vertex on each of these
parallel edges, splitting each parallel edge into two. For the edge weighted
$a_{0}$, the two edges have weights $1,a_{0}$, and for the remaining edges the
new edges get weights $a_{i},x_{i}$. It is a simple and instructive exercise
to see that this has the desired effect. Note also that if the original affine
projection $\pi$ was monotone, then so is the constructed simple projection.

Now we show the result for simple projections. If the Hamiltonian Cycle
polynomial were a monotone projection of the permanent, then by the Main
Lemma, some face of $\operatorname{Newt}(\operatorname{perm})$ would be an
extension of $\operatorname{Newt}(HC)$.

The Newton polytope of the permanent is the cycle cover polytope (see Section
2). The cycle cover polytope can easily be described by the $m^{2}$
inequalities saying that all variables $x_{i,j}$ are non-negative, together
with the equalities saying that each vertex has in-degree and out-degree
exactly 1, namely $\sum_{i}x_{i,j}=1$ for all $j$ and $\sum_{j}x_{i,j}=1$ for
all $i$ (it is easy to see that these are necessary; for sufficiency, see, e.
g., [19, Theorem 18.1]). Since equalities do not count towards the complexity
of a polytope, we have $c(\operatorname{Newt}(\operatorname{perm}_{m}))\leq
m^{2}$.

But the Newton polytope of the $n$-th Hamiltonian Cycle polynomial is exactly
the TSP polytope, which by [17, Corollary 2] requires extension complexity
$2^{\Omega(n)}$. ∎

######  Theorem 4.3.

Over any totally ordered semi-ring, the perfect matching polynomial (or
“unsigned Pfaffian”) is not a monotone affine p-projection of the permanent;
in fact, any monotone affine projection from the permanent to the perfect
matching polynomial has blow-up at least $2^{\Omega(n)}$.

###### Proof.

The proof is the same as for the Hamiltonian Cycle polynomial, using [17,
Theorem 1], which gives a lower bound of $2^{\Omega(n)}$ on the extension
complexity of the perfect matching polytope, which is the Newton polytope of
the perfect matching polynomial. ∎

######  Theorem 4.4.

Over any totally ordered semi-ring, for any $q$, the $q$-th cut polynomial is
not a monotone affine p-projection of the permanent; in fact, any monotone
affine projection from the permanent to the $q$-th cut polynomial has blow-up
at least $2^{\Omega(n)}$.

###### Proof.

Use [6, Theorem 7], which says that
$xc(\operatorname{Newt}(\operatorname{Cut}^{2}))\geq 2^{\Omega(n)}$, as
$\operatorname{Newt}(\operatorname{Cut}^{2})$ is the cut polytope. The one
additional observation we need is that
$\operatorname{Newt}(\operatorname{Cut}^{q})$ is just the $(q-1)$-scaled
version of $\operatorname{Newt}(\operatorname{Cut}^{2})$, and this rescaling
does not affect the extension complexity. ∎

###  4.2 Monotone Formula and Circuit Lower Bounds

As pointed out by an anonymous reviewer, the universality of the permanent
also holds in the monotone setting, so lower bounds on monotone projections
from the permanent imply the same lower bounds on monotone formula size, and
therefore quasi-polynomially related lower bounds on monotone circuit size. We
assume circuits only have gates of bounded fan-in; with unbounded fan-in,
rather than losing a factor of a half in the exponent of the exponent, we lose
a factor of a third.

######  Proposition 4.5.

Any polynomial computable by a monotone formula of size $s$ is a monotone
projection of $\operatorname{perm}_{s+1}$.

###### Proof.

The proof of the universality of the permanent given in [4, Proposition 2.16]
works _mutatis mutandis_ in the monotone setting. ∎

As a consequence of this, Theorems 4.2–4.4 are nearly tight, since every
monotone polynomial in $n$ variables of $\text{poly}(n)$ degree can be written
as a monotone formula of size $2^{O(n\log n)}$ (write it as a sum of
monomials).

######  Corollary 4.6.

Over any totally ordered semi-ring, any monotone formula computing the
Hamiltonian Cycle polynomial, the perfect matching polynomial, or the $q$-th
cut polynomial has size at least $2^{\Omega(n)}$. Consequently, any monotone
circuit computing these polynomials has size at least $2^{\Omega(\sqrt{n})}$.

For the cut polynomials, we believe this result to be new. For the other
polynomials, this provides a new proof of (slightly weaker versions of)
previously known lower bounds. Namely, Jerrum and Snir gave a lower bound of
$(n-1)((n-2)2^{n-3}+1)=2^{n+\Omega(\log n)}$ on the monotone circuit size of
$HC$ [9, Section 4.4], and a lower bound of $n(2^{n-1}-1)$ on the monotone
circuit size of the permanent [9, Section 4.3]. As the permanent is a monotone
projection of the perfect matching polynomial—namely, restrict the perfect
matching polynomial to a bipartite graph, e. g., by setting $x_{ij}=0$
whenever $i$ and $j$ have the same parity—the same lower bound holds for the
perfect matching polynomial.

###### Proof.

The first part follows by combining Proposition 4.5 with Theorems 4.2–4.4. The
second part follows from the fact that monotone circuits of size $s$ can be
balanced to have size $\text{poly}(s)$ and depth $O(\log^{2}s)$ (the proof in
[24] works _mutatis mutandis_ in the monotone setting), which can then be
converted to monotone formulas of size $s^{O(\log s)}=2^{O(\log^{2}s)}$ by the
usual conversion from bounded fan-in circuits to formulas. If there is a
monotone circuit of size $s$ computing any of these polynomials, there is thus
a monotone formula of size $2^{O(\log^{2}s)}$, which must be at least
$2^{\Omega(n)}$, so $s\geq 2^{\Omega(n^{1/2})}$. ∎

##  5 Open Questions

Despite the common feeling that Razborov’s super-polynomial lower bound [15]
on monotone Boolean circuits for CLIQUE “finished off” monotone Boolean
circuit lower bounds, several natural and interesting question remain. For
example, does Directed $s$-$t$ Connectivity require monotone Boolean circuits
of size $\Omega(n^{3})$? (A matching upper bound is given by the Bellman–Ford
algorithm.) Is there a monotone Boolean reduction from general perfect
matching to bipartite perfect matching? A positive answer to the following
question would rule out such monotone (projection) reductions.

######  Open Question 5.1.

Extend Theorem 4.3 from formal polynomials over the Boolean semi-ring to
Boolean functions.

However, there are even easier questions, intermediate between the Boolean
function case and the algebraic case considered in this paper; Jukna [10]
discusses the notion of one polynomial “counting” another, which means that
they agree on all $\\{0,1\\}$ inputs.

######  Open Question 5.2.

Prove that no monotone polynomial-size projection of the permanent agrees with
the perfect matching polynomial on all $\\{0,1\\}$ inputs (“counts the perfect
matching polynomial”). Similarly, prove that no monotone polynomial-size
projection of the permanent counts the Hamiltonian cycle polynomial.

S. Jukna points out (personal communication) that projections of the $s$-$t$
connectivity polynomial correspond, even in the Boolean setting, to switching-
and-rectifier networks, so the known lower bounds on monotone switching-and-
rectifier networks (see, e. g., the survey [16]) imply that the Hamiltonian
path polynomial and the permanent are not monotone p-projections of the
$s$-$t$ connectivity polynomial, even over the Boolean semi-ring. This helps
explain why the only known monotone lower bound on the $s$-$t$ connectivity
polynomial that we are aware of [10] goes by a somewhat roundabout proof:
Razborov’s lower bound on CLIQUE [15], followed by Valiant’s reduction from
the clique polynomial to the Hamiltonian path polynomial [21], followed by a
standard reduction from Hamiltonian path to counting $s$-$t$ paths. In the
course of discussing this, we were led to the following question; although the
motivation for the question has since disappeared, it still seems like an
interesting question about polytopes, whose answer may require new methods.

######  Open Question 5.3 (S. Jukna, personal communication).

Is the $m$-th $s$-$t$ path polytope an extension of the $n$-th TSP polytope
(or $n$-th cycle cover polytope) with $m\leq\text{poly}(n)$?

Since the separation problem for the $s$-$t$ path polytope is
$\mathsf{NP}$-hard (see, e. g., [19, §13.1])—and the cycle cover polytope has
low (extension) complexity—answering this question negatively seems to require
more subtle understanding of these polytopes than “simply” an extended
formulation lower bound.

Another example of a natural polytope question with a similar flavor comes
from the cut polynomials. In combination with Bürgisser’s results and
questions on the cut polynomials [3] (discussed in Section 1), we are led to
the following question.

######  Open Question 5.4.

Is the $m$-th cut polytope an extension of the $n$-th TSP polytope, for
$m\leq\text{poly}(n)$?

A negative answer would show that $\operatorname{Cut}^{q}$ is not complete for
non-negative polynomials in $\mathsf{VNP}_{\mathbb{Q}}$ under monotone
p-projections, though as with the example of the permanent, this is not
necessarily an obstacle to being $\mathsf{VNP}$-complete under general
p-projections. Yet even the monotone completeness of the cut polynomials
remains open. In fact, even more basic questions remain open:

######  Open Question 5.5.

Is every non-negative polynomial in $\mathsf{VNP}$ a monotone projection of
the Hamiltonian Cycle polynomial? Is there any polynomial that is “positive
$\mathsf{VNP}$-complete” in this sense?

To relate this to the current proofs of $\mathsf{VNP}$-completeness of
$HC_{n}$, we need to draw a distinction. Let $\mathsf{VP}_{\mathbb{R}}^{\geq
0}$ denote the polynomial families in $\mathsf{VP}_{\mathbb{R}}$ all of whose
coefficients are non-negative, and let $\mathsf{mVP_{\mathbb{R}}}$ (“monotone
$\mathsf{VP}$”) denote the class of families of polynomials with polynomially
many variables, of polynomial degree, and computable by polynomial-size
_monotone_ circuits over $\mathbb{R}$. Similarly, define
$\mathsf{VNP}_{\mathbb{R}}^{\geq 0}$ to be the non-negative polynomials in
$\mathsf{VNP}_{\mathbb{R}}$, and $\mathsf{mVNP}_{\mathbb{R}}$ to be the
function families of the form
$f_{n}=\sum_{\vec{e}\in\\{0,1\\}}^{\text{poly}(n)}g_{m}(\vec{e},\vec{x})$,
where $m\leq\text{poly}(n)$ and $(g_{m})\in\mathsf{mVP}_{\mathbb{R}}$.

Valiant’s original completeness proof for the Hamiltonian Cycle polynomial
[21] is “mostly” monotone: It uses polynomial-size formulas for the
coefficients of the monomials (coming from the definition of $\mathsf{VNP}$),
but otherwise is entirely monotone. In other words, the proof shows that $HC$
is $\mathsf{mVNP}$-hard under monotone projections. However, we note that it’s
not clear whether $HC$ is even _in_ $\mathsf{mVNP}$! Question 5.5 asks whether
$HC$, or indeed any polynomial, is $\mathsf{VNP}^{\geq 0}$-complete under
monotone projections; the question of whether there exist polynomials that are
$\mathsf{mVNP}$-complete under monotone projections also seems potentially
interesting.

Finally, we ask about stronger notions of monotone reduction, which seem to
require a different kind of proof technique. Recall that a _c-reduction_ from
$f$ to $g$ is a family of algebraic circuits for $f$ with oracle gates for
$g$.

######  Open Question 5.6.

Do the analogues of Theorems 4.2–4.4 hold for monotone bounded-depth
c-reductions in place of affine p-projections? What about weakly-skew or even
general monotone c-reductions?

##  6 Subsequent Developments

Since the appearance of the preliminary version of this paper [7], our Main
Lemma 3.1 has been used to prove that several other polynomials of
combinatorial and complexity-theoretic interest are not sub-exponential-size
projections of the permanent.

  1. 1.

The $n$-th _satisfiability polynomial_ over $\mathbb{F}_{q}$ is a polynomial
in $n+8\binom{n}{3}$ variables denoted $X_{1},\dotsc,X_{n}$ and $\\{Y_{c}:c\in
C_{n}\\}$ where $C_{n}$ denote the set of clauses on 3 literals in $n$
variables. It is defined as:

$\operatorname{Sat}^{q}_{n}(X,Y)=\sum_{a\in\\{0,1\\}^{n}}\left(\prod_{i\in[n]}X_{i}^{q-1}\right)\left(\prod_{c\in
C_{n}:c(a)=1}Y_{c}^{q-1}\right)$

  2. 2.

A _clow_ in an $n$-vertex graph is a closed walk of length exactly $n$, in
which the minimum-numbered vertex appears exactly once. The $n$-th _clow
polynomial_ over $\mathbb{F}_{q}$ is a polynomial in $\binom{n}{2}+n$
variables $X_{e}$ for each edge $e$ in the complete undirected graph $K_{n}$
on $n$ vertices and $Y_{v}$ for each $v\in[n]$. It is defined as:

$\operatorname{Clow}^{q}_{n}(X,Y)=\sum_{w:\text{clow of length
$n$}}\left(\prod_{e:\text{edges in
$w$}}X_{e}^{q-1}\right)\left(\prod_{v:\text{distinct vertices in
$w$}}Y_{v}^{q-1}\right),$

or more precisely,

$\operatorname{Clow}^{q}_{n}(X,Y)=\sum_{w=[v_{0},\dotsc,v_{n-1}]}\left(\prod_{i\in[n]}X_{(v_{i-1},v_{i\text{
mod
}n})}^{q-1}\right)\left(\prod_{v\in\\{v_{0},\dotsc,v_{n-1}\\}}Y_{v}^{q-1}\right),$

where the sum is over clows $w$ and $v_{0}$ denotes the minimum-numbered
vertex in $w$.

  3. 3.

The _clique polynomial_ is a polynomial in $\binom{n}{2}$ variables $X_{e}$:

$\operatorname{Clique}_{n}(X)=\sum_{T\subseteq\binom{[n]}{2}:T\text{ is a
clique in $K_{n}$}}\prod_{e\in T}X_{e}.$

######  Theorem (Mahajan and Saurabh [12, Theorems 2 and 6]).

Over any totally ordered semi-ring, any monotone affine projection from the
permanent to $\operatorname{Sat}^{q}_{n}$ or to the clique polynomial requires
blow-up at least $2^{\Omega(\sqrt{n})}$. Any monotone affine projection from
the permanent to $\operatorname{Clow}^{q}_{n}$ requires blow-up at least
$2^{\Omega(n)}$.

As in Section 4.2, we get the following corollary. Again, we note that the
lower bound on the clique polynomial over the Boolean semi-ring only works for
the formal clique polynomial (in contrast to Razborov’s result [15], which
works for any monotone Boolean circuit computing the CLIQUE _function_).

######  Corollary 6.1.

Over any totally ordered semi-ring, any monotone formula computing
$\operatorname{Clow}^{q}_{n}$ has size at least $2^{\Omega(n)}$ and any
monotone circuit computing $\operatorname{Clow}^{q}_{n}$ has size at least
$2^{\Omega(\sqrt{n})}$. Any monotone formula computing
$\operatorname{Sat}^{q}_{n}$ or $\operatorname{Clique}_{n}$ has size at least
$2^{\Omega(\sqrt{n})}$ and any monotone circuit computing these polynomials
has size at least $2^{\Omega(n^{1/4})}$.

For the clow and satisfiability polynomials, we believe this result to be new.
For the clique polynomials, this provides a new proof of (a slightly weaker
version of) the exponential monotone circuit lower bound due to Schnorr
[18].333Schnorr showed a $\binom{n}{k}-1$ lower bound on the monotone circuit
size of the the $k$-th clique polynomial $\operatorname{Clique}^{k}_{n}$, the
sum over all cliques of size $k$, rather than all cliques. For $k=n/2$, this
lower bound is asymptotically equal to $2^{n}/\sqrt{\pi n/2}$. The $k$-th
clique polynomial $\operatorname{Clique}_{n}^{k}$ is the degree $k$
homogeneous component of the clique polynomial $\operatorname{Clique}_{n}$; by
homogenization (implicit in Strassen’s work, explicit in Valiant [23, Lemma
2]), any monotone circuit of size $s$ for $\operatorname{Clique}_{n}$ can be
converted into a monotone circuit of size $s(n/2+1)^{2}$ computing
$\operatorname{Clique}^{n/2}_{n}$. Thus Schnorr’s result implies a lower bound
of $\Omega(2^{n}/n^{5/2})$ on the monotone circuit complexity of
$\operatorname{Clique}_{n}$.

Acknowledgment. We would like to thank Stasys Jukna for the question that
motivated this paper [11], and cstheory.stackexchange.com for providing a
forum for the question. We also thank Stasys for comments on a draft, pointing
out the paper by Schnorr [18], and interesting discussions leading to
Questions 5.2 and 5.3. We thank Leslie Valiant for an interesting conversation
that led to Question 5.5. We thank Ketan Mulmuley and Youming Qiao for
collaborating on [8], which is why the author had Newton polytopes on the
mind. We thank an anonymous reviewer for pointing out the monotone
universality of the permanent and therefore the implications for monotone
formula and circuit size (Section 4.2 and Corollary 6.1). The author was
supported during this work by an Omidyar Fellowship from the Santa Fe
Institute.

## References

  * [1] Scott Aaronson.  A linear-optical proof that the permanent is $\mathsf{\\#P}$-hard.  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 467(2136):3393–3405, 2011.  doi:10.1098/rspa.2011.0232.
  * [2] Scott Aaronson and Alex Arkhipov.  The computational complexity of linear optics.  Theory Comput., 9:143–252, 2013.  Initially appeared in STOC ’11.  doi:10.4086/toc.2013.v009a004.
  * [3] Peter Bürgisser.  On the structure of Valiant’s complexity classes.  Discrete Math. Theor. Comput. Sci., 3(3):73–94, 1999.  Initially appeared in STACS ’98.
  * [4] Peter Bürgisser.  Completeness and reduction in algebraic complexity theory, volume 7 of Algorithms and Computation in Mathematics.  Springer-Verlag, Berlin, 2000.  doi:10.1007/978-3-662-04179-6.
  * [5] Nicolas de Rugy-Altherre.  A dichotomy theorem for homomorphism polynomials.  In 37th International Symposium on Mathematical Foundations of Computer Science, MFCS ’12, pages 308–322, 2012.  doi:10.1007/978-3-642-32589-2\\_29.
  * [6] Samuel Fiorini, Serge Massar, Sebastian Pokutta, Hans Raj Tiwary, and Ronald de Wolf.  Exponential lower bounds for polytopes in combinatorial optimization.  J. ACM, 62(2):Art. 17, 23, 2015.  Initially appeared in STOC ’12.  doi:10.1145/2716307.
  * [7] Joshua A. Grochow.  Monotone projection lower bounds from extended formulation lower bounds.  ECCC Tech. Report TR15-171 and arXiv:1510.08417 [cs.CC], 2015.
  * [8] Joshua A. Grochow, Ketan D. Mulmuley, and Youming Qiao.  Boundaries of $\mathsf{VP}$ and $\mathsf{VNP}$.  In 43rd International Colloquium on Automata, Languages, and Programming, ICALP ’16, 2016.  Preprint of the full version available as arXiv:1605.02815 [cs.CC].
  * [9] Mark Jerrum and Marc Snir.  Some exact complexity results for straight-line computations over semirings.  J. ACM, 29(3):874–897, 1982.  doi:10.1145/322326.322341.
  * [10] Stasys Jukna.  Lower bounds for monotone counting circuits.  ECCC Tech. Report TR14-169, 2014.
  * [11] Stasys Jukna.  “Why is Hamilton Cycle so different from permanent?”.  http://cstheory.stackexchange.com/q/27496/129, 2014\.
  * [12] Meena Mahajan and Nitin Saurabh.  Some complete and intermediate polynomials in algebraic complexity theory.  In Proceedings of 11th Computer Science Symposium in Russia, volume 9691 of CSR ’16, Lecture Notes in Computer Science, pages 251–265, 2016.  Full version available as ECCC Tech. Report TR16-038.  doi:10.1007/978-3-319-34171-2\\_18.
  * [13] Henryk Minc.  Permanents.  Cambridge University Press, 1984.
  * [14] Thomas Muir and William H. Metzler.  A treatise on the theory of determinants.  Dover, New York, 1960.
  * [15] Alexander A. Razborov.  Lower bounds of monotone complexity of the logical permanent function.  Mat. Zametki, 37(6):887–900, 942, 1985.
  * [16] Alexander A. Razborov.  Lower bounds for deterministic and nondeterministic branching programs.  In L. Budach, editor, Fundamentals of Computation Theory, volume 529 of Lecture Notes in Computer Science, pages 47–60. Springer Berlin Heidelberg, 1991.  URL: http://people.cs.uchicago.edu/~razborov/files/fct.pdf, doi:10.1007/3-540-54458-5\\_49.
  * [17] Thomas Rothvoß.  The matching polytope has exponential extension complexity.  In Proceedings of the 46th Annual ACM Symposium on Theory of Computing, STOC ’14, pages 263–272, New York, NY, USA, 2014. ACM.  Preprint available as arXiv:1311.2369 [cs.CC].  doi:10.1145/2591796.2591834.
  * [18] Claus-Peter Schnorr.  A lower bound on the number of additions in monotone computations.  Theoret. Comput. Sci., 2(3):305–315, 1976.  doi:10.1016/0304-3975(76)90083-9.
  * [19] Alexander Schrijver.  Combinatorial optimization. Polyhedra and efficiency. Vol. A, volume 24 of Algorithms and Combinatorics.  Springer-Verlag, Berlin, 2003.  Paths, flows, matchings, Chapters 1–38.
  * [20] Volker Strassen.  Vermeidung von Divisionen.  J. Reine Angew. Math., 264:184–202, 1973.
  * [21] Leslie G. Valiant.  Completeness classes in algebra.  In Proceedings of the Eleventh Annual ACM Symposium on Theory of Computing, STOC ’79, pages 249–261, New York, NY, USA, 1979. ACM.  doi:10.1145/800135.804419.
  * [22] Leslie G. Valiant.  The complexity of computing the permanent.  Theoret. Comput. Sci., 8(2):189–201, 1979.  doi:10.1016/0304-3975(79)90044-6.
  * [23] Leslie G. Valiant.  Negation can be exponentially powerful.  Theoret. Comput. Sci., 12(3):303–314, 1980.  doi:10.1016/0304-3975(80)90060-2.
  * [24] Leslie G. Valiant, Sven Skyum, Stuart Berkowitz, and Charles Rackoff.  Fast parallel computation of polynomials using few processors.  SIAM J. Comput., 12(4):641–644, 1983.  doi:10.1137/0212043.
  * [25] Jacobus H. van Lint and Richard M. Wilson.  A course in combinatorics.  Cambridge University Press, 2001.
  * [26] Tzu-Chieh Wei and Simone Severini.  Matrix permanent and quantum entanglement of permutation invariant states.  J. Math. Phys., 51(9), 2010.  doi:10.1063/1.3464263.


# SfePy - Write Your Own FE Application

Robert Cimrman∗† * Corresponding author<EMAIL_ADDRESS>New Technologies
Research Centre, University of West Bohemia, Plzeň, Czech RepublicCopyright ©
2014 Robert Cimrman. This is an open-access article distributed under the
terms of the Creative Commons Attribution License, which permits unrestricted
use, distribution, and reproduction in any medium, provided the original
author and source are credited. http://creativecommons.org/licenses/by/3.0/

###### Abstract

SfePy (Simple Finite Elements in Python) is a framework for solving various
kinds of problems (mechanics, physics, biology, …) described by partial
differential equations in two or three space dimensions by the finite element
method. The paper illustrates its use in an interactive environment or as a
framework for building custom finite-element based solvers.

###### Index Terms:

partial differential equations, finite element method, Python

##  1 Introduction

SfePy (Simple Finite Elements in Python) is a multi-platform (Linux, Mac OS X,
Windows) software released under the New BSD license, see http://sfepy.org. It
implements one of the standard ways of discretizing partial differential
equations (PDEs) which is the finite element method (FEM) [R1].

So far, SfePy has been employed for modelling in materials science, including,
for example, multiscale biomechanical modelling (bone, muscle tissue with
blood perfusion) [R2], [R3], [R4], [R5], [R6], [R7], computation of acoustic
transmission coefficients across an interface of arbitrary microstructure
geometry [R8], computation of phononic band gaps [R9], [R10], finite element
formulation of Schroedinger equation [R11], and other applications, see also
Figure 3. The software can be used as

  * •

a collection of modules (a library) for building custom or domain-specific
applications,

  * •

a highly configurable "black box" PDE solver with problem description files in
Python.

In this paper we focus on illustrating the former use by using a particular
example. All examples presented below were tested to work with the version
2013.3 of SfePy.

##  2 Development

The SfePy project uses Git for source code management and GitHub web site for
the source code hosting and developer interaction, similarly to many other
scientific python tools. The version 2013.3 has been released in September 18.
2013, and its git-controlled sources contained 761 total files, 138191 total
lines (725935 added, 587744 removed) and 3857 commits done by 15 authors.
However, about 90% of all commits were done by the main developer as most
authors contributed only a few commits. The project seeks and encourages more
regular contributors, see http://sfepy.org/doc-devel/development.html.

##  3 Short Description

The code is written mostly in Python. For speed in general, it relies on fast
vectorized operations provided by NumPy [R12] arrays, with heavy use of
advanced broadcasting and "index tricks" features. C and Cython [R13] are used
in places where vectorization is not possible, or is too difficult/unreadable.
Other components of the scientific Python software stack are used as well,
among others: SciPy [R14] solvers and algorithms, Matplotlib [R15] for 2D
plots, Mayavi [R16] for 3D plots and simple postprocessing GUI, IPython [R17]
for a customized shell, SymPy [R18] for symbolic operations/code generation
etc.

The basic structure of the code allows a flexible definition of various
problems. The problems are defined using components directly corresponding to
mathematical counterparts present in a weak formulation of a problem in the
finite element setting: a solution domain and its sub-domains (regions),
variables belonging to suitable discrete function spaces, equations as sums of
terms (weak form integrals), various kinds of boundary conditions,
material/constitutive parameters etc.

The key notion in SfePy is a _term_ , which is the smallest unit that can be
used to build _equations_ (linear combinations of terms). It corresponds to a
weak formulation integral and takes usually several arguments: (optional)
material parameters, a single virtual (or test) function variable and zero or
more state (or unknown) variables. The available terms (currently 105) are
listed at our web site (http://sfepy.org/doc-devel/terms_overview.html). The
already existing terms allow to solve problems from many scientific domains,
see Figure 3. Those terms cover many common PDEs in continuum mechanics,
poromechanics, biomechanics etc. with a notable exception of electromagnetism
(work in progress, see below).

SfePy discretizes PDEs using a continuous Galerkin approximation with finite
elements as defined in [R1]. Discontinuous element-wise constant approximation
is also possible. Currently the code supports the 2D area (triangle,
rectangle) and 3D volume (tetrahedron, hexahedron) elements. Structural
elements like shells, plates, membranes or beams are not supported with a
single exception of a hyperelastic Mooney-Rivlin membrane.

Several kinds of basis or shape functions can be used for the finite element
approximation of the physical fields:

  * •

the classical nodal (Lagrange) basis can be used with all element types;

  * •

the hierarchical (Lobatto) basis can be used with tensor-product elements
(rectangle, hexahedron).

Although the code can provide the basis function polynomials of a high order
(10 or more), orders greater than 4 are not practically usable, especially in
3D. This is caused by the way the code assembles the element contributions
into the global matrix - it relies on NumPy vectorization and evaluates the
element matrices all at once and then adds to the global matrix - this allows
fast term evaluation and assembling but its drawback is a very high memory
consumption for high polynomial orders. Also, the polynomial order has to be
uniform over the whole (sub)domain where a field is defined. Related to
polynomial orders, tables with quadrature points for numerical integration are
available or can be generated as needed.

We are now working on implementing Nédélec and Raviart-Thomas vector bases in
order to support other kinds of PDEs, such as the Maxwell equations of
electromagnetism.

Once the equations are assembled, a number solvers can be used to solve the
problem. SfePy provides a unified interface to many standard codes, for
example UMFPACK [R19], PETSc [R20], Pysparse [R21] as well as the solvers
available in SciPy. Various solver classes are supported: linear, nonlinear,
eigenvalue, optimization, and time stepping.

Besides the external solvers mentioned above, several solvers are implemented
directly in SfePy.

Nonlinear/optimization solvers:

  * •

The Newton solver with a backtracking line-search is the default solver for
all problems. For simplicity we use a unified approach to solve both the
linear and non-linear problems - former (should) converge in a single
nonlinear solver iteration.

  * •

The steepest descent optimization solver with a backtracking line-search can
be used as a fallback optimization solver when more sophisticated solvers
fail.

Time-stepping-related solvers:

  * •

The stationary solver is used to solve time-independent problems.

  * •

The equation sequence solver is a stationary solver that analyzes the
dependencies among equations and solves smaller blocks first. An example would
be the thermoelasticity problem described below, where the elasticity equation
depends on the load given by temperature distribution, but the temperature
distribution (Poisson equation) does not depend on deformation. Then the
temperature distribution can be found first, followed by the elasticity
problem with the already known temperature load. This greatly reduces memory
usage and improves speed of solution.

  * •

The simple implicit time stepping solver is used for (quasistatic) time-
dependent problems, using a fixed time step.

  * •

The adaptive implicit time stepping solver can change the time step according
to a user provided function. The default function adapt_time_step() decreases
the step in case of bad Newton convergence and increases the step (up to a
limit) when the convergence is fast. It is convenient for large deformation
(hyperelasticity) problems.

  * •

The explicit time stepping solver can be used for dynamic problems.

##  4 Thermoelasticity Example

This example involves calculating a temperature distribution in an object
followed by an elastic deformation analysis of the object loaded by the
thermal expansion and boundary displacement constraints. It shows how to use
SfePy in a script/interactively. The actual equations (weak form) are
described below. The entire script consists of the following steps:

Import modules. The SfePy package is organized into several sub-packages. The
example uses:

  * •

sfepy.fem: the finite element method (FEM) modules

  * •

sfepy.terms: the weak formulation terms - equations building blocks

  * •

sfepy.solvers: interfaces to various solvers (SciPy, PETSc, …)

  * •

sfepy.postprocess: post-processing & visualization based on Mayavi

import numpy as npfrom sfepy.fem import (Mesh, Domain, Field, FieldVariable,
Material, Integral, Equation, Equations, ProblemDefinition)from sfepy.terms
import Termfrom sfepy.fem.conditions import Conditions, EssentialBCfrom
sfepy.solvers.ls import ScipyDirectfrom sfepy.solvers.nls import Newtonfrom
sfepy.postprocess import ViewerLoad a mesh file defining the object geometry.

mesh = Mesh.from_file(’meshes/2d/square_tri2.mesh’)domain = Domain(’domain’,
mesh)Define solution and boundary conditions domains, called regions.

omega = domain.create_region(’Omega’, ’all’)left =
domain.create_region(’Left’, ’vertices in x < -0.999’, ’facet’)right =
domain.create_region(’Right’, ’vertices in x > 0.999’, ’facet’)bottom =
domain.create_region(’Bottom’, ’vertices in y < -0.999’, ’facet’)top =
domain.create_region(’Top’, ’vertices in y > 0.999’, ’facet’)Save regions for
visualization.

domain.save_regions_as_groups(’regions.vtk’)Use a quadratic approximation for
temperature field, define unknown $T$ and test $s$ variables.

field_t = Field.from_args(’temperature’, np.float64, ’scalar’, omega, 2)t =
FieldVariable(’t’, ’unknown’, field_t, 1)s = FieldVariable(’s’, ’test’,
field_t, 1, primary_var_name=’t’)Define numerical quadrature for the
approximate integration rule.

integral = Integral(’i’, order=2)Define the Laplace equation governing the
temperature distribution:

$\int_{\Omega}\nabla s\cdot\nabla T=0\;,\quad\forall s\;.$

term = Term.new(’dw_laplace(s, t)’, integral, omega, s=s, t=t)eq =
Equation(’temperature’, term)eqs = Equations([eq])Set boundary conditions for
the temperature: $T=10\mbox{ on }\Gamma_{\rm left}$, $T=30\mbox{ on
}\Gamma_{\rm right}$.

t_left = EssentialBC(’t_left’, left, {’t.0’ : 10.0})t_right =
EssentialBC(’t_right’, right, {’t.0’ : 30.0})Create linear (ScipyDirect) and
nonlinear solvers (Newton).

ls = ScipyDirect({})nls = Newton({}, lin_solver=ls)Combine the equations,
boundary conditions and solvers to form a full problem definition.

pb = ProblemDefinition(’temperature’, equations=eqs, nls=nls,
ls=ls)pb.time_update(ebcs=Conditions([t_left, t_right]))Solve the temperature
distribution problem to get $T$.

temperature = pb.solve()out = temperature.create_output_dict()Use a linear
approximation for displacement field, define unknown $\underline{u}$ and test
$\underline{v}$ variables. The variables are vectors with two components in
any point, as we are solving on a 2D domain.

field_u = Field.from_args(’displacement’, np.float64, ’vector’, omega, 1)u =
FieldVariable(’u’, ’unknown’, field_u, mesh.dim)v = FieldVariable(’v’, ’test’,
field_u, mesh.dim, primary_var_name=’u’)Set Lamé parameters of elasticity
$\lambda$, $\mu$, thermal expansion coefficient $\alpha_{ij}$ and background
temperature $T_{0}$. Constant values are used here. In general, material
parameters can be given as functions of space and time.

lam = 10.0 # Lame parameters.mu = 5.0te = 0.5 # Thermal expansion
coefficient.T0 = 20.0 # Background temperature.eye_sym = np.array([[1], [1],
[0]], dtype=np.float64)m = Material(’m’, lam=lam, mu=mu, alpha=te *
eye_sym)Define and set the temperature load variable to $T-T_{0}$.

t2 = FieldVariable(’t’, ’parameter’, field_t, 1, primary_var_name=’(set-to-
None)’)t2.set_data(t() - T0)Define the thermoelasticity equation governing
structure deformation:

$\int_{\Omega}D_{ijkl}\
e_{ij}(\underline{v})e_{kl}(\underline{u})-\int_{\Omega}(T-T_{0})\
\alpha_{ij}e_{ij}(\underline{v})=0\;,\quad\forall\underline{v}\;,$

where $D_{ijkl}=\mu(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk})+\lambda\
\delta_{ij}\delta_{kl}$ is the homogeneous isotropic elasticity tensor and
$e_{ij}(\underline{u})=\frac{1}{2}(\frac{\partial u_{i}}{\partial
x_{j}}+\frac{\partial u_{j}}{\partial x_{i}})$ is the small strain tensor. The
equations can be built as linear combinations of terms.

term1 = Term.new(’dw_lin_elastic_iso(m.lam, m.mu, v, u)’, integral, omega,
m=m, v=v, u=u)term2 = Term.new(’dw_biot(m.alpha, v, t)’, integral, omega, m=m,
v=v, t=t2)eq = Equation(’temperature’, term1 - term2)eqs = Equations([eq])Set
boundary conditions for the displacements: $\underline{u}=0\mbox{ on
}\Gamma_{\rm bottom}$, $u_{1}=0.0\mbox{ on }\Gamma_{\rm top}$ ($x$
-component).

u_bottom = EssentialBC(’u_bottom’, bottom, {’u.all’ : 0.0})u_top =
EssentialBC(’u_top’, top, {’u.[0]’ : 0.0})Set the thermoelasticity equations
and boundary conditions to the problem definition.

pb.set_equations_instance(eqs,
keep_solvers=True)pb.time_update(ebcs=Conditions([u_bottom, u_top]))Solve the
thermoelasticity problem to get $\underline{u}$.

displacement = pb.solve()out.update(displacement.create_output_dict())Save the
solution of both problems into a single VTK file.

pb.save_state(’thermoelasticity.vtk’, out=out)Display the solution using
Mayavi.

view = Viewer(’thermoelasticity.vtk’)view(vector_mode=’warp_norm’,
rel_scaling=1, is_scalar_bar=True, is_wireframe=True, opacity={’wireframe’ :
0.1})

###  4.1 Results

The above script saves the domain geometry as well as the temperature and
displacement fields into a VTK file called ’thermoelasticity.vtk’ and also
displays the results using Mayavi. The results are shown in Figures 1 and 2.

Figure 1: The temperature distribution.

Figure 2: The deformed mesh showing displacements.

##  5 Alternative Way: Problem Description Files

Problem description files (PDF) are Python modules containing definitions of
the various components (mesh, regions, fields, equations, …) using basic data
types such as dict and tuple. For simple problems, no programming at all is
required. On the other hand, all the power of Python (and supporting SfePy
modules) is available when needed. The definitions are used to construct and
initialize in an automatic way the corresponding objects, similarly to what
was presented in the example above, and the problem is solved. The main script
for running a simulation described in a PDF is called simple.py.

###  5.1 Example: Temperature Distribution

This example defines the problem of temperature distribution on a 2D
rectangular domain. It directly corresponds to the temperature part of the
thermoelasticity example, only for the sake of completeness a definition of a
material coefficient is shown as well.

from sfepy import data_dirfilename_mesh = data_dir +
’/meshes/2d/square_tri2.mesh’materials = { ’coef’ : ({’val’ : 1.0},),}regions
= { ’Omega’ : ’all’, ’Left’ : (’vertices in (x < -0.999)’, ’facet’), ’Right’ :
(’vertices in (x > 0.999)’, ’facet’),}fields = { ’temperature’ : (’real’, 1,
’Omega’, 2),}variables = { ’t’ : (’unknown field’, ’temperature’, 0), ’s’ :
(’test field’, ’temperature’, ’t’),}ebcs = { ’t_left’ : (’Left’, {’t.0’ :
10.0}), ’t_right’ : (’Right’, {’t.0’ : 30.0}),}integrals = { ’i1’ : (’v’,
2),}equations = { ’eq’ : ’dw_laplace.i1.Omega(coef.val, s, t) = 0’}solvers = {
’ls’ : (’ls.scipy_direct’, {}), ’newton’ : (’nls.newton’, {’i_max’ : 1,
’eps_a’ : 1e-10, }),}options = { ’nls’ : ’newton’, ’ls’ : ’ls’,}Many more
examples can be found at http://docs.sfepy.org/gallery/gallery or
http://sfepy.org/doc-devel/examples.html.

Figure 3: Gallery of applications. Perfusion and acoustic images by Vladimír
Lukeš.

##  6 Conclusion

We briefly introduced the open source finite element package SfePy as a tool
for building domain-specific FE-based solvers as well as a black-box PDE
solver.

###  6.1 Support

Work on SfePy is partially supported by the Grant Agency of the Czech
Republic, projects P108/11/0853 and 101/09/1630.

## References

  * [R1] Thomas J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover Publications, 2000.
  * [R2] R.~Cimrman and E.~Rohan. Two-scale modeling of tissue perfusion problem using homogenization of dual porous media. _International Journal for Multiscale Computational Engineering_ , 8(1):81–102, 2010.
  * [R3] E.~Rohan, R.~Cimrman, S.~Naili, and T.~Lemaire. Multiscale modelling of compact bone based on homogenization of double porous medium. In _Computational Plasticity X - Fundamentals and Applications_ , 2009\.
  * [R4] E.~Rohan and R.~Cimrman. Multiscale fe simulation of diffusion-deformation processes in homogenized dual-porous media. _Mathematics and Computers in Simulation_ , 82(10):1744–1772, 2012.
  * [R5] R.~Cimrman and E.~Rohan. On modelling the parallel diffusion flow in deforming porous media. _Mathematics and Computers in Simulation_ , 76(1-3):34–43, 2007.
  * [R6] E.~Rohan and R.~Cimrman. Multiscale fe simulation of diffusion-deformation processes in homogenized dual-porous media. _Mathematics and Computers in Simulation_ , 82(10):1744–1772, 2012.
  * [R7] E.~Rohan, S.~Naili, R.~Cimrman, and T.~Lemaire. Hierarchical homogenization of fluid saturated porous solid with multiple porosity scales. _Comptes Rendus - Mecanique_ , 340(10):688–694, 2012.
  * [R8] E.~Rohan and V.~Lukeš. Homogenization of the acoustic transmission through a perforated layer. _Journal of Computational and Applied Mathematics_ , 234(6):1876–1885, 2010.
  * [R9] E.~Rohan, B.~Miara, and F.~Seifrt. Numerical simulation of acoustic band gaps in homogenized elastic composites. _International Journal of Engineering Science_ , 47(4):573–594, 2009.
  * [R10] E.~Rohan and B.~Miara. Band gaps and vibration of strongly heterogeneous reissner-mindlin elastic plates. _Comptes Rendus Mathematique_ , 349(13-14):777–781, 2011.
  * [R11] R.~Cimrman, J.~Vackář, M.~Novák, O.~Čertík, E.~Rohan, and M.~Tůma. Finite element code in python as a universal and modular tool applied to kohn-sham equations. In _ECCOMAS 2012 - European Congress on Computational Methods in Applied Sciences and Engineering, e-Book Full Papers_ , pages 5212–5221, 2012.
  * [R12] T. E. Oliphant. Python for scientific computing. _Computing in Science & Engineering_, 9(3):10-20, 2007. http://www.numpy.org.
  * [R13] R. Bradshaw, S. Behnel, D. S. Seljebotn, G. Ewing, et al. The Cython compiler. http://cython.org.
  * [R14] E. Jones, T. E. Oliphant, P. Peterson, et al. SciPy: Open source scientific tools for Python, 2001-. http://www.scipy.org.
  * [R15] J. D. Hunter. Matplotlib: A 2d graphics environment. _Computing in Science & Engineering_, 9(3):90-95, 2007. http://matplotlib.org/.
  * [R16] P. Ramachandran and G. Varoquaux. Mayavi: 3d visualization of scientific data. _IEEE Computing in Science & Engineering_, 13(2):40-51, 2011.
  * [R17] F. Pérez and B. E. Granger. IPython: A system for interactive scientific computing. _Computing in Science & Engineering_, 9(3):21-29, 2007. http://ipython.org/.
  * [R18] SymPy Development Team. Sympy: Python library for symbolic mathematics, 2013. http://www.sympy.org.
  * [R19] T. A. Davis. Algorithm 832: UMFPACK, an unsymmetric-pattern multifrontal method. _ACM Transactions on Mathematical Software_ , 30(2):196–199, 2004.
  * [R20] S. Balay, J. Brown, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang. PETSc Web page, 2013. http://www.mcs.anl.gov/petsc.
  * [R21] R. Geus, D. Wheeler, and D. Orban. Pysparse documentation. http://pysparse.sourceforge.net.