# A Unified Framework for Probabilistic Verification of AI Systems
via Weighted Model Integration

Paolo Morettin DISI
University of Trento
Italy  Andrea Passerini  DISI
University of Trento
Italy  Roberto Sebastiani  DISI
University of Trento
Italy

###### Abstract

The probabilistic formal verification (PFV) of AI systems is in its infancy.
So far, approaches have been limited to ad-hoc algorithms for specific classes
of models and/or properties. We propose a unifying framework for the PFV of AI
systems based on Weighted Model Integration (WMI), which allows to frame the
problem in very general terms. Crucially, this reduction enables the
verification of many properties of interest, like fairness, robustness or
monotonicity, over a wide range of machine learning models, without making
strong distributional assumptions. We support the generality of the approach
by solving multiple verification tasks with a single, off-the-shelf WMI
solver, then discuss the scalability challenges and research directions
related to this promising framework.

##  1 Introduction and motivation

The development of verification techniques for AI systems is nowadays
considered a critical problem for both the scientific community [Seshia et
al., 2016, Amodei et al., 2016] and regulatory bodies [EC, 2021].

In the last three decades, model checking has made huge progress in the
verification of both hardware and software systems, guaranteeing the
deployment of safe systems in many fields [Clarke, 1997]. Qualitative
verification techniques can answer whether a property is satisfied or not, but
they fall short when we need to quantify how likely a property will hold in an
uncertain environment, or how close our system is to fulfilling a formal
requirement.

Due to the large variety of models and their inherent combinatorial
complexity, verification efforts in the machine learning (ML) community
typically target specific classes of predictors [Baluta et al., 2019],
properties [Bastani et al., 2019], or both [Weng et al., 2019]. However, the
complexity and versatility of modern AI systems calls for a unified framework
to assess their trustworthiness, which cannot be captured by a single
evaluation metric or formal property. This papers aims to introduce such a
framework. We show how by leveraging the Weighted Model Integration (WMI)
[Belle et al., 2015] formalism, it is possible to devise a unified formulation
for the probabilistic verification of combinatorial AI systems. Broadly
speaking, WMI is the task of computing probabilities of arbitrary combinations
of logical and algebraic constraints given a structured joint distribution
over both continuous and discrete variables.

The proposed framework satisfies three crucial desiderata for general-purpose
verification of AI systems:

  * D1)

Support for _arbitrary distributions_ over continuous and discrete domains,
including complex models learned from data.

  * D2)

Support for a _broad range of ML models_ via a single representation language,
from decision trees to neural networks (NNs) and ensemble models.

  * D3)

Support for a _broad range of properties_ of interest, ranging from ML-
specific to system-wide properties that are typically considered in PFV.

The contribution of this paper is threefold: 1) it bridges AI verification
with the field of probabilistic inference with algebraic and logical
constraints; 2) it provides theoretical means of studying the relations among
different PFV tasks and 3) it introduces another significant use case for WMI
solvers, whose performance so far have been mostly evaluated on synthetic
problems only.

After discussing the related work and providing the necessary background on
WMI, we present our unified framework. We showcase its flexibility on multiple
verification tasks, concluding with a discussion of the challenges and
opportunities offered by the framework.

##  2 Related Work

Probabilistic model checking [Vardi, 1985, Katoen, 2016] traditionally
addresses the verification of sequential systems, achieving high scalability
with strong distributional assumptions on the priors [Hartmanns and Hermanns,
2015], such as independence or unimodality. For instance, random variables
such as transition delays in timed automata are usually modelled with
independent exponential distributions. These techniques cannot be directly
applied to modern AI systems with ML components solving predictive problems
over multivariate, multi-modal distributions. Adding support for the
verification of sequential systems with ML components and complex priors is a
challenging research goal.

The common thread among the related work in ML verification is targeting
specific subclasses of problems (models and/or properties) in order to achieve
higher scalability. The largest body of work considers deterministic models,
and reduces the decision problem to combinatorial reasoning problems like
Satisfiability Modulo Theories (SMT) [Huang et al., 2017, Katz et al., 2017,
2019, Einziger et al., 2019, Devos et al., 2021], linear (LP) [Bastani et al.,
2016] or mixed integer-linear programming (MILP) [Tjeng et al., 2017, Chen et
al., 2019]. This line of work cannot address quantitative tasks, like
verifying fairness properties or certifying a nondeterministic system.

The probabilistic verification of ML models is in its infancy. Most works
focus on verifying if a specific property, such as robustness [Mangal et al.,
2019, Weng et al., 2019, Mirman et al., 2021] or fairness [Bastani et al.,
2019], holds with probability higher than a given threshold. Alternatively,
quantitative verification tasks can be reduced to a combinatorial counting
problem. This approach was first applied to the evaluation of simple
probabilistic programs with univariate priors on random variables [Chistikov
et al., 2015, Gehr et al., 2016]. This approach was later extended to fairness
verification of decision procedures based on small ML models [Albarghouthi et
al., 2017]. Similarly to our approach, a probabilistic model of the
environment is part of the specification that the learned model must satisfy,
albeit the only property considered in this work is group fairness.
Quantitative robustness verification of binarized NNs was also reduced to a
counting problem, leveraging the advances in approximate counting over
propositional logic theories [Baluta et al., 2019]. This approach is however
limited to discrete systems. The framework presented here is a generalization
of the formalism above to logical/numerical problems.

##  3 Background on Weighted Model Integration

We assume the reader is familiar with the basic syntax, semantics, and results
of propositional and first-order logic. We adopt the notation and definitions
in [Morettin et al., 2019] — including some terminology and concepts from
Satisfiability Modulo Theories (SMT) — which we summarize below. We will
restrict SMT to quantifier-free formulas over linear real arithmetic
($\text{SMT}(\mathcal{LRA})$) which consist in atomic propositions over some
set $\mathbf{A}$ (aka Boolean atoms), linear (in)equalities
$(\sum_{i}a_{i}x_{i}\bowtie b)$ over a set of continuous variables
$\mathbf{x}$ s.t. $a_{i},b$ are rational constants and $\bowtie\
\in\\{{\leq,\geq,>,<,=,\neq}\\}$ (aka $\mathcal{LRA}$ atoms), all of them
combined with the standard Boolean operators
$\\{{\neg,\wedge,\vee,\rightarrow,\leftarrow,\leftrightarrow}\\}$. (Our work
easily extends to the integers [Kolb et al., 2018].) The notions of literal,
clause, partial and total truth assignment extend directly to
$\text{SMT}(\mathcal{LRA})$ atoms. We represent truth assignments $\mu$ as
conjunctions of literals, s.t. positive and negative literals in $\mu$ mean
that the atom is assigned to true and false respectively.

For brevity, we introduce interval formulas
$(x\in[l,u])\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}(l\leq x)\land(x\leq
u)$ and if-then-else terms $ite(.;.;.)$ s.t.
$t=ite(cond;e_{1};e_{2})\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(cond\rightarrow t=e_{1})\land(\neg cond\rightarrow t=e_{2})$.

$\text{SMT}(\mathcal{LRA})$ is concerned with finding a truth assignment $\mu$
to the atoms of a formula $\Delta$ that both tautologically entails $\Delta$
($\mu\models_{\mathbb{B}}\Delta$) and is $\mathcal{LRA}$-satisfiable. We
denote with $\mathcal{T}\mathcal{T}\mathcal{A}(\Delta)$ the set of all the
$\mathcal{LRA}$-satisfiable total truth assignments satisfying $\Delta$. For
instance, if $\Delta\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(x<0)\lor(x>1)$, then $\mu\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(x<0)\wedge(x>1)\models_{\mathbb{B}}\Delta$ but it is not
$\mathcal{LRA}$-satisfiable. Thus, the formula has only two
$\mathcal{LRA}$-satisfiable total truth assignments:
$\mathcal{T}\mathcal{T}\mathcal{A}(\Delta)=\\{(x<0)\land\neg(x>1),\neg(x<0)\land(x>1)\\}$.
(Hereafter, we may drop “$\mathcal{LRA}$-satisfiable” when referring to
assignments in $\mathcal{T}\mathcal{T}\mathcal{A}(\Delta)$ since this is clear
from context.) Notice that every $\mu\in\mathcal{T}\mathcal{T}\mathcal{A}$ is
a convex subregion of $\Delta$. We denote with $\mu^{\mathcal{LRA}}{}$ the
portion on $\mu$ mapping $\mathcal{LRA}$-atoms to truth values.

###### Example 1

The formula $\Delta\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(x\in[0,1])\land(y\in[0,1])\land(A\rightarrow(x+y\leq 1))$ defines
a region over $\mathbf{A}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\\{A\\}$ and $\mathbf{x}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\\{x,y\\}$, depicted in Figure 2 (left). It results in 3 (total)
$\mathcal{LRA}$-satisfiable truth assignments (we omit
$(x\in[0,1])\wedge(y\in[0,1])$ for short):

$\displaystyle\mathcal{T}\mathcal{T}\mathcal{A}(\Delta)=$
$\displaystyle\\{\overbrace{A\land(x+y\leq 1)}^{\mu_{1}},$
$\displaystyle\overbrace{\neg A\land(x+y\leq 1)}^{\mu_{2}},\overbrace{\neg
A\land\neg(x+y\leq 1)}^{\mu_{3}}\\}.\quad\diamond$

$w(x,y,A)$$\times$$x+y\leq 1$$A$$x$$1$$y$$\top$$\bot$$\bot$$\top$

Figure 1:  (Left) The hybrid region defined in Ex. 1. (Center) The weighted
formula in Ex. 2. (Right) An equivalent XADD representation of the weight.

We introduce non-negative weight functions
$w:\mathbf{x}\cup\mathbf{A}\mapsto\mathbb{R}^{+}$, which intuitively defines a
(possibly unnormalized) density function over $\mathbf{A}\cup\mathbf{x}$.
$w(\mathbf{x},\mathbf{A})$ is defined by: (1) a $\text{SMT}(\mathcal{LRA})$
formula $\chi$, called the _support_ of $w$, outside of which $w$ is $0$; (2)
a combination of real functions whose integral over $\mu^{\mathcal{LRA}}$ is
computable, structured by means of nested if-then-elses with
$\mathcal{LRA}$-conditions (we refer to Morettin et al. [2019] for a formal
definition of $w(\mathbf{x},\mathbf{A})$). When (2) is a DAG with polynomial
leaves, as depicted in Fig. 2 (right), it is equivalent to the notion of
_extended algebraic decision diagram_ (XADD) introduced by Sanner et al.
(2012). $w_{[\mu]}$ denotes $w$ restricted to the truth values of $\mu$. The
Weighted Model Integral of a formula $\Delta(\mathbf{x},\mathbf{A})$ and a
weight function $w$($\mathbf{x}$,$\mathbf{A}$) over $\chi$ is defined as:

$\displaystyle\textstyle WMI(\Delta,w)\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\sum_{\mu\in\mathcal{T}\mathcal{T}\mathcal{A}(\Delta\wedge\chi)}\int_{\mu^{\mathcal{LRA}}}w_{[\mu]}(\mathbf{x})\>d\mathbf{x}$
(1)

where $\int_{\Delta}$ denotes integration over the solutions of a formula
$\Delta$. WMI generalizes Weighted Model Counting (WMC) [Sang et al., 2005].
As with WMC, WMI can be interpreted as the total unnormalized probability mass
of the pair $\langle\Delta,w\rangle$, i.e. its _partition function_. The
normalized conditional probability of $\Gamma$ given $\Delta$ can then be
computed as:

$\displaystyle P(\Gamma|\Delta)=WMI(\Gamma\land\Delta,w)/WMI(\Delta,w).$ (2)

###### Example 2

Let $\Delta$ be as in Ex.1, with weights:

$\displaystyle w_{(x+y\leq 1)}(x,y)=x,\quad\quad w_{\neg A}(x,y)=y,$
$\displaystyle w_{\ell}(x,y)=1\quad\quad\forall\ell\notin\\{(x+y\leq 1),\neg
A\\}$

That is, $w$ defines the following piecewise density over the convex
subregions of $\Delta$ (Fig. 2, center):

$\displaystyle w(\mathbf{x},\mathbf{A})\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}ite((x+y\leq 1);x;1)\times ite(A;1;y)$

with $\chi\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}\top$, so that:
$\left\\{\begin{array}[]{ll}w_{[\mu_{1}]}(x,y)&=x\times 1=x\\\
w_{[\mu_{2}]}(x,y)&=x\times y=xy\\\ w_{[\mu_{3}]}(x,y)&=1\times
y=y\end{array}\right\\}$.

Then, the weighted model integral of $\langle\Delta,w\rangle$ is:

$\displaystyle WMI(\Delta,w)$
$\displaystyle=\int_{0}^{1}\\!\\!\\!\int_{0}^{(1-y)}\\!\\!\\!x\>dx\>dy+\int_{0}^{1}\\!\\!\\!\int_{0}^{(1-y)}\\!\\!\\!xy\>dx\>dy$
$\displaystyle+\int_{0}^{1}\\!\\!\\!\int_{(1-y)}^{1}\\!\\!\\!y\>dx\>dy=\frac{1}{6}+\frac{1}{24}+\frac{1}{3}=\frac{13}{24}.$

We can now compute normalized probabilities such as

$P(A|\Delta)=\frac{WMI(A\land\Delta,w)}{WMI(\Delta,w)}=\frac{1/6}{13/24}=\frac{4}{13}\quad\diamond$

Crucially for our purposes, $\Gamma$ and $\Delta$ in (2) can be _arbitrarily
complex_ SMT-LRA formulas. This gives as a very powerful computational tool
for quantifying the probability associated to complex events or properties.

##  4 WMI-based verification of AI systems

In this section we introduce our unifying perspective on the probabilistic
verification of AI systems. We envision two roles, the _developer_ of the
system under verification $\mathcal{S}$ and the _verifier_ (e.g. a regulatory
body) that provides a specification of the requirements $\mathcal{R}$.

The developer is responsible for faithfully modelling $\mathcal{S}$ with a
logical encoding $\Delta_{\mathcal{S}}$ of its deterministic behaviour, such
as the functional relationships between its input and output. In many
practical scenarios, however, $\mathcal{S}$ is not fully deterministic. For
instance, the images captured by a camera might be subject to noise introduced
by the sensor. The developer is then expected to model the uncertainty of
$\mathcal{S}$ with a probabilistic model $P_{\mathcal{S}}$.

Similarly, the requirements $\mathcal{R}$ provided by the verifier include a
logical formula $\Delta_{\mathcal{R}}$, encoding a desirable property that the
system should satisfy with high probability. This probability is computed with
respect to a probabilistic model that encodes the uncertainty of the
environment where $\mathcal{S}$ is expected to be deployed, denoted by
$P_{\mathcal{R}}$.

The goal of verification is computing the probability that system
$\mathcal{S}$ satisfies requirement $\mathcal{R}$ given: (1) a probabilistic
model that jointly represents the uncertainty of $\mathcal{S}$ and of its
surrounding environment and (2) a logical encoding of $\mathcal{S}$ and of the
property described in $\mathcal{R}$. In what follows, we denote with
$\mathbf{x}$ and $\mathbf{y}$ the input and output of $\mathcal{S}$
respectively, and with $\mathbf{e}$, we denote extra environmental variables
that are not observed by $\mathcal{S}$ but might be part of the specification,
such as protected attributes that shouldn’t be accessed by a fair hiring
system (e.g. the gender of a job applicant). The PFV task can be formalized as
computing :

$\displaystyle\\!\\!P(\mathcal{R}\>|\>\mathcal{S})=\frac{\int_{\Delta_{\mathcal{R}}(\mathbf{e},\mathbf{x},\mathbf{y})\land\Delta_{\mathcal{S}}(\mathbf{x},\mathbf{y})}P(\mathbf{e},\mathbf{x},\mathbf{y})\>d\mathbf{e}\>d\mathbf{x}\>d\mathbf{y}}{\int_{\Delta_{\mathcal{S}}(\mathbf{x},\mathbf{y})}P(\mathbf{e},\mathbf{x},\mathbf{y})\>d\mathbf{e}\>d\mathbf{x}\>d\mathbf{y}},$
(3)

where $P(\mathbf{e},\mathbf{x},\mathbf{y})$ models the probabilistic
relationship between inputs, outputs and environment, and factorizes as:

$\displaystyle
P(\mathbf{e},\mathbf{x},\mathbf{y})=P_{\mathcal{S}}(\mathbf{y}\>|\>\mathbf{e},\mathbf{x})\cdot
P_{\mathcal{R}}(\mathbf{e},\mathbf{x}).$

The definition of $P_{\mathcal{S}}(\mathbf{y}\>|\>\mathbf{e},\mathbf{x})$
accommodates uncontrollable sources of uncertainty as well as any deliberate
probabilistic mechanism in $\mathcal{S}$, even those that are not conditioned
on the input, such as randomly sampled latent variables in a generative model.

Figure 2:  An overview of the PFV task in Example 3, showcasing the modularity
of the framework. The requirement combines a prior $P_{\mathcal{R}}$ with a
formal property $\Delta_{\mathcal{R}}$ over the inputs and outputs of
$\mathcal{S}$. The latter is composed of a noisy ($P_{\mathcal{S}}$) camera
followed by two NNs $f$ and $g$, feeding their input into a simple decision
procedure $d$

.

###### Example 3

Consider a system composed by a security camera connected to a gate, meant to
grant access to authorized individuals only. Two convolutional NNs $f$ and $g$
receive the camera signal $c(\mathbf{x})$. The first network predicts
$f(c(\mathbf{x}))\in\\{auth,\neg auth\\}$. The second network predicts the
distance of the individual from the gate $g(c(\mathbf{x}))\in\mathbb{R}$.
Their output is then fed into a decision procedure that controls the lock,
$d(f(c(\mathbf{x})),g(c(\mathbf{x})))=y\in\\{lock,\neg lock\\}$. The system is
mostly deterministic and is logically modelled with $\Delta_{\mathcal{S}}$,
with the exception of the camera noise, which is modelled with
$P_{\mathcal{S}}(c(\mathbf{x})\>|\>\mathbf{x})$. The verifier provides both a
probabilistic model of the environment that the system would observe through
the camera, $P_{\mathcal{R}}(\mathbf{x})$ and a set of properties
$\Delta_{\mathcal{R}}$ over $\mathbf{x},y$ that should hold with high
probability. Fig. 2 summarizes the use case above.

Following (2), the inference problem can be solved by computing the following
weighted model integrals:

$\displaystyle
P(\mathcal{R}\>|\>\mathcal{S})=\frac{WMI(\Delta_{\mathcal{R}}\land\Delta_{\mathcal{S}},w)}{WMI(\Delta_{\mathcal{S}},w)}$
(4)

where $\Delta_{\mathcal{R}}$ and $\Delta_{\mathcal{S}}$ are SMT formulas over
$\mathbf{e},\mathbf{x},\mathbf{y}$ encoding $\mathcal{R}$ and the
deterministic aspect of $\mathcal{S}$ respectively, while the weight function
$w$ encodes the probabilistic aspect of both $\mathcal{S}$ and the
environment:

$\displaystyle
w(\mathbf{e},\mathbf{x},\mathbf{y})=w_{\mathcal{S}}(\mathbf{y}\>|\>\mathbf{e},\mathbf{x})\cdot
w_{\mathcal{R}}(\mathbf{e},\mathbf{x})\propto
P(\mathbf{e},\mathbf{x},\mathbf{y}).$

All the properties described in this manuscript can be expressed in terms of
_algebraic_ and/or _logical_ preconditions and postconditions
($\mathcal{R}_{pre}\rightarrow\mathcal{R}_{post}$). For the purpose of
verification, the relevant quantity (or ratios of these quantities, as shown
in (5) and (6)) to compute is:

$\displaystyle
P(\mathcal{R}_{post}\>|\>\mathcal{R}_{pre},\mathcal{S})=\frac{WMI(\Delta_{\mathcal{R}_{post}}\land\Delta_{\mathcal{R}_{pre}}\land\Delta_{\mathcal{S}},w)}{WMI(\Delta_{\mathcal{R}_{pre}}\land\Delta_{\mathcal{S}},w)}$

Since SMT encodings can be naturally conjoined and the class of weight
functions supported by the formalism is closed under multiplication, it is
trivial to combine arbitrary encodings of $\mathcal{S}$ and $\mathcal{R}$ from
both logical and probabilistic perspectives. Now, we answer the question:
"what kind of systems and properties can be encoded in this framework?".

###  4.1 D1: Probabilistic modelling

The choice of density functions is dictated by the integration procedure
adopted by the WMI solver. Any function can be used as long as it is
integrable, exactly or approximately, in a region defined by a conjunction of
$\mathcal{LRA}$ atoms, i.e. a convex polytope. Prior work on exact WMI has
almost exclusively focused on using polynomials as building blocks for their
structured weight functions. The reason is twofold: 1) polynomials are easy to
work with, being closed under sum, product and integration over polytopes; 2)
they can approximate any density with arbitrary precision. Notably, the
complexity of exact integration grows only polynomially with respect to the
degree [Baldoni et al., 2011], making higher-degree piecewise functions
preferrable over a drastically higher number of lower-degree pieces. When
Monte Carlo-based procedures are used, any non-negative density can be
adopted, including Gaussians [Merrell et al., 2017, Dos Martires et al.,
2019], as long as statistical guarantees on the approximation are provided.

By combining these density functions by means of $\mathcal{LRA}$-conditions,
sums and products, a variety of popular probabilistic models can be modelled,
including (hybrid variants of) graphical models like Bayesian and Markov
networks [Belle et al., 2015, Albarghouthi et al., 2017], probabilistic logic
programs with continuous variables [Dos Martires et al., 2019], and tractable
models like density estimation trees (DETs) [Ram and Gray, 2011] or sum-
product networks (SPNs) [Poon and Domingos, 2011]. For illustrative purposes,
we report the encodings of the latter two models. Typically, the support is an
axis-aligned rectangle inferred from the training set, albeit more complex and
possibly non-convex supports can be used [Morettin et al., 2020].

#### DETs.

The density estimation variant of decision or regression trees shares the same
internal structure with binary decision nodes of the form $(x_{i}\leq k)$ or a
proposition $A_{i}$. Different from its predictive variants, each leaf of a
DET encodes the probability mass corresponding to the subset of the support
induced by the decisions in its path from the root. DETs can be trivially
encoded as a tree of if-then-elses with constant leaves.

#### SPNs.

This class of models combine tractable univariate distributions by means of
mixtures (weighted sums) over the same variables or factorizations (products)
over disjoint sets of variables, resulting in a tractable but expressive joint
probability. The common choices for the univariate distributions, Gaussians
and piecewise polynomials [Molina et al., 2018], and combinations by means of
sums and products, are integrable in $\mathcal{LRA}$ and encodable in the
formalism.

###  4.2 D2: Encoding ML models in $\text{SMT}(\mathcal{LRA})$

A broad range of (deterministic) ML models can be encoded as
$\text{SMT}(\mathcal{LRA})$ formulas. Some representative examples are
reported below.

#### Tree-based predictors.

As shown above for DETs, the structure of tree-based predictors can be encoded
by means of nested if-then-elses with propositional or linear conditions.
Leaves are simply encoded with (a conjunction of) atoms $(y=c_{i})$, mapping
output variables to their respective values. If axis-aligned splits are not
expressive enough, arbitrary $\mathcal{LRA}$-conditions enable more complex
piecewise decompositions of the joint density, such as those employed in
Optimal Classification Trees [Bertsimas and Dunn, 2017].

#### Non-linear predictors.

Support for linear models, which can be trivially encoded in
$\text{SMT}(\mathcal{LRA})$ as $(y=(\mathbf{w}\mathbf{x}+b))$, can be extended
to non-linear cases. For instance, we can model NNs with rectified linear
activations by encoding each unit $i$ with inputs $\mathbf{x}_{i}$, weights
$\mathbf{w}_{i}$ and $b_{i}$ as 111Equipping WMI solvers with specialized
reasoning procedures for ReLUs [Katz et al., 2017, 2019] would be ideal.:

$\displaystyle(h_{i}=(\mathbf{w}_{i}\mathbf{x}_{i}+b_{i}))\land(y_{i}=ite(h_{i}>0,\>h_{i},\>0))$

The full network is then encoded by conjoining the $\text{SMT}(\mathcal{LRA})$
representation of each unit. Other non-linear activation functions can be
approximated with arbitrary precision. Common operations like convolutions and
max/average pooling have $\text{SMT}(\mathcal{LRA})$ representations. Other
non-linear predictors, such as support vector machines with piecewise linear
feature maps [Huang et al., 2013] can be similarly encoded.

#### Complex models.

The compositional nature of SMT can support ensembles of $K$ ML models,

$\displaystyle\textstyle(\mathbf{y}=a(\mathbf{y}_{1},...,\mathbf{y}_{K}))\land\bigwedge_{i=1}^{K}\Delta_{f_{i}}(\mathbf{x})$

as long as the aggregation function $a$ can be encoded, such as the average.
More in general, if every component of an AI system can be modelled with
$\text{SMT}(\mathcal{LRA})$, its behaviour can be verified as a whole. This is
in stark contrast with most approaches in literature, which are only able to
verify ML components in isolation. As verification of "traditional" software
and hardware systems often relies on SMT modelling and solving, using the same
paradigm for the verification of modern AI systems is a promising direction.

###  4.3 D3: Encoding the properties via WMI

So far we demonstrated the flexibility of WMI in modelling and reasoning
probabilistically over a wide range of AI systems, but this would be a
pointless exercise if we couldn’t use it to quantify properties of practical
interest. The ML community has identified a number of important properties
that learned models should satisfy. For instance, in contexts with high socio-
economic stakes, predictions over individuals of a population should be fair.
Many definitions of fairness are probabilistic in nature, being based on the
notion of _population model_ $P(\mathbf{x})$.

#### Set-based properties.

For instance, if $f$ is determining a positive vs. negative outcome for an
individual 222These concepts apply to regression tasks by considering a
suitable notion of distance and a threshold., it is often desirable to
guarantee the _demographic parity_ of a system relative to a protected
subpopulation $\mathcal{M}$ [Calders et al., 2009]:

$\displaystyle
1-\delta\leq\frac{P(f(\mathbf{x})=pos\>|\>\mathbf{x}\in\mathcal{M})}{P(f(\mathbf{x})=pos\>|\>\mathbf{x}\in\overline{\mathcal{M}})}\leq
1/(1-\delta)$ (5)

In short, the ratio of positive outcomes among members of $\mathcal{M}$ and
its complement $\overline{\mathcal{M}}$ should be close. A limitation of this
notion is that it does not take into account whether the candidate individuals
are qualified for the positive outcome in the first place. If the qualified
subpopulation $\mathcal{Q}$ is known to the verifier, _equality of
opportunity_ can be required instead [Hardt et al., 2016]:

$\displaystyle
1-\delta\leq\frac{P(f(\mathbf{x})=pos\>|\>\mathbf{x}\in\mathcal{M}\cap\mathcal{Q})}{P(f(\mathbf{x})=pos\>|\>\mathbf{x}\in\overline{\mathcal{M}}\cap\mathcal{Q})}\leq
1/(1-\delta)$ (6)

We notice that these preconditions can be trivially encoded when the subsets
$\mathcal{M}$ and $\mathcal{Q}$ are characterized by a combination of
categorical features, arguably the most common case in fairness scenarios.
Additionally, numerical constraints can be encoded with arbitrary piecewise
linear sets. For instance, assuming that the output of the model is a binary
decision $y_{hire}$, that
$\Delta_{\mathcal{M}}\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}x_{female}$
and that the qualified group can be defined in terms of GPA and years of work
experience $\Delta_{\mathcal{Q}}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(x_{GPA}>3.5)\lor(x_{experience}>20)$, we can compute the ratio in
(6) as [Albarghouthi et al., 2017]:

$\displaystyle\frac{WMI(y_{hire}\land\Gamma,w)}{WMI(\Gamma,w)}\quad/\quad\frac{WMI(y_{hire}\land\Psi,w)}{WMI(\Psi,w)}$
$\displaystyle\text{where }\Gamma\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\Delta_{\mathcal{M}}\land\Delta_{\mathcal{Q}}\land\Delta_{\mathcal{S}},\Psi\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\neg\Delta_{\mathcal{M}}\land\Delta_{\mathcal{Q}}\land\Delta_{\mathcal{S}}.$

#### Distance-based properties.

The above notions, which are defined on subpopulations or groups, are commonly
referred to as group fairness properties. An alternative notion is _individual
fairness_ [Dwork et al., 2012], stating that similar individuals, according to
a suitable metric $d$ and threshold $\epsilon$, should receive similar
treatment/outcomes:

$\displaystyle
P(f(\mathbf{x})=f(\mathbf{x}^{\prime})\>|\>d(\mathbf{x},\mathbf{x}^{\prime})<\epsilon)>1-\delta$
(7)

Notice that the property above is akin to the notion of _probabilistic
robustness_ [Weng et al., 2019, Mirman et al., 2021]. In contrast with
fairness verification, where $P(\mathbf{x})$ is defined globally, the
distributional assumptions over the perturbations are typically local. For
instance, one might verify the robustness of predictions of an image
classifier when a Gaussian noise is added to the pixels.

Encoding distance-based properties poses additional challenges, requiring to
define the system’s behaviour inside the local neighborhood of infinitely many
points. In turn, this requires multiple instantiations of the same random
variable $\mathbf{x},\mathbf{x}^{\prime}\sim P(\mathbf{x})$ in the
$\text{SMT}(\mathcal{LRA})$ formula. We observe that a simple solution is
“cloning” the weight function and its support:

$\displaystyle cl(\Delta)\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}\Delta\land\Delta^{\prime},\quad\quad
cl(w)\stackrel{{\scriptstyle\text{\tiny def}}}{{=}}w\cdot w^{\prime}$

where $expr^{\prime}$ denotes the expression obtained by substituting every
variable occurrence $v$ with a fresh copy $v^{\prime}$ in $expr$. This
approach, which guarantees that $\mathbf{x}$ and $\mathbf{x}^{\prime}$ are
independently drawn from the same distribution, is called _self-composition_
in program analysis [Barthe et al., 2011]. We can then compute queries
involving both, such as $x^{\prime}\in[x-\epsilon,x+\epsilon]$ (Fig. 3).

Figure 3:  (_Left_) A simple distribution $\Delta=(x\in[0,1])$ and $w(x)=x$.
(_Right_) the distribution obtained by self-composition, enabling the
computation of queries like $P(x^{\prime}\in[x-\epsilon,x+\epsilon])$.

In terms of _distances_ , both $L_{1}$ and $L_{\infty}$ can be encoded by
defining absolute values as auxiliary variables $ite((a\leq
b),(abs_{a,b}=b-a),(abs_{a,b}=a-b))$:

$\displaystyle\textstyle(L_{1}^{\mathbf{x},\mathbf{x}^{\prime}}=\sum_{i=1}^{N}\>abs_{x_{i},x_{i}^{\prime}}),\quad(L_{\infty}^{\mathbf{x},\mathbf{x}^{\prime}}=max_{i=1}^{N}\>abs_{x_{i},x_{i}^{\prime}})$

where $max(.)$ is encoded with nested if-then-elses:

$\displaystyle max(\\{v_{1},v_{2}\\})\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}ite(v_{1}<v_{2};\>$ $\displaystyle v_{2};\>v_{1})$ $\displaystyle
max(\\{v_{1},v_{2}\\}\cup V)\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}ite(v_{1}<v_{2};\>$ $\displaystyle max(\\{v_{2}\\}\cup V);$
$\displaystyle max(\\{v_{1}\\}\cup V))$

Following (7) and assuming w.l.o.g. a binary classifier $y=f(\mathbf{x})$, we
can quantify individual fairness as:

$\displaystyle\frac{WMI((y\leftrightarrow
y^{\prime})\land(L_{\infty}^{\mathbf{x},\mathbf{x}^{\prime}}<\epsilon)\land
cl(\Delta_{\mathcal{S}}),cl(w))}{WMI((L_{\infty}^{\mathbf{x},\mathbf{x}^{\prime}}<\epsilon)\land
cl(\Delta_{\mathcal{S}}),cl(w))}$

Probabilistic robustness can be similarly encoded, with the difference that
the distribution of noisy inputs $\mathbf{x}^{\prime}$ is conditioned on
$\mathbf{x}$ and explicitly provided as part of $\mathcal{R}$:
$w_{\mathcal{R}}(\mathbf{x},\mathbf{x}^{\prime})=w_{\mathcal{R}}(\mathbf{x})\cdot
w_{\mathcal{R}}(\mathbf{x}^{\prime}\>|\>\mathbf{x})$.

#### Other algebraic properties.

In general, with arbitrary combinations of linear inequalities and logical
constraints, we can encode many useful algebraic properties. Beyond fairness
and robustness, there has been considerable work in enforcing and verifying
_monotonic behaviour_ of learned predictors [Sill, 1997, Dvijotham et al.,
2018, Wehenkel and Louppe, 2019, Liu et al., 2020, Sivaraman et al., 2020]. In
a probabilistic setting, this can be quantified as
$P(f(\mathbf{x})<f(\mathbf{x}^{\prime})\>|\>\mathbf{x}<\mathbf{x}^{\prime})$.
Again, this can be computed in WMI by leveraging self-composition:

$\displaystyle\frac{WMI((y<y^{\prime})\land(\mathbf{x}<\mathbf{x}^{\prime})\land
cl(\Delta_{\mathcal{S}}),cl(w))}{WMI((\mathbf{x}<\mathbf{x}^{\prime})\land
cl(\Delta_{\mathcal{S}}),cl(w))}$ (8)

Checking the _equivalence_ of two predictors $f$ and $g$ finds applications in
verifying that a compressed model that has to be deployed in a resource-
constrained setting behaves consistently w.r.t. the original model [Narodytska
et al., 2018]:

$\displaystyle
P(f(\mathbf{x})=g(\mathbf{x}))=\frac{WMI((y_{f}=y_{g})\land\Delta_{\mathcal{S}},w)}{WMI(\Delta_{\mathcal{S}},w)}$
(9)

Monotonicity and equivalence are not the only algebraic properties of interest
for the ML community. In their influential paper, Katz et al. [Katz et al.,
2017] verify a NN-based collision detection system for unmanned aircrafts
against properties such as: “If the intruder is directly ahead and is moving
towards the ownship, the system won’t issue a clear-of-conflict advisory”,
encoded as combinations of inequalities over linear and angular quantities.
Since all these cyber-physical safety constraints are effectively encoded in
$\text{SMT}(\mathcal{LRA})$, they can be encoded in our framework.

###### Example 4

We show how to quantify the robustness of the system described in Ex. 3 to the
camera noise, given a validated prior on $\mathbf{x}$ (e.g. an SPN), whose
unnormalized density is denoted with $w_{\mathcal{R}}$. For simplicity, we
consider a piecewise polynomial approximation of a Gaussian function
$w_{\mathcal{S}}(\mathbf{n};\sigma)\simeq\mathcal{N}(\mathbf{0},\sigma\mathbb{I})$
that conservatively over-estimates the noise level added to each pixel
independently: $\Delta_{c}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(c(\mathbf{x})=\mathbf{x}+\mathbf{n})$. $\Delta_{f}$ and
$\Delta_{g}$ denote the $\text{SMT}(\mathcal{LRA})$ encoding of the CNNs $f$
and $g$, where $(f(c(\mathbf{x}))\geq 0)$ iff an individual is _authorized_.
$d$ unlocks the gate only if an authorized individual is at most 4 meters away
$\Delta_{d}\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(y\leftrightarrow((f(c(\mathbf{x}))\geq
0)\wedge(g(c(\mathbf{x}))\leq 4)))$. We can compute how robust are the
decisions in $\sigma$-noisy ($y$) vs. noiseless ($y^{\prime}$) settings,
applying self-composition to the relevant part of $\mathcal{S}$ and computing
how likely the decision would not change:

$\displaystyle WMI((y\leftrightarrow y^{\prime})\wedge\Gamma,w)\quad/\quad
WMI(\Gamma,w)$ where $\displaystyle\Gamma\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}(c(\mathbf{x})^{\prime}=\mathbf{x})\wedge\Delta_{c}\wedge
cl(\Delta_{f}\wedge\Delta_{g}\wedge\Delta_{d})$ $\displaystyle
w(\mathbf{x},\mathbf{n})\stackrel{{\scriptstyle\text{\tiny
def}}}{{=}}w_{\mathcal{R}}(\mathbf{x})\cdot
w_{\mathcal{S}}(\mathbf{n};\sigma)$

##  5 Experiments

We provide evidence of the flexibility of the proposed approach by
implementing a procedure that is able to jointly verify multiple properties on
a feedforward NN. Thanks to our unifying perspective, we were able to solve
every task with a non-specialized, off-the-shelf WMI solver.

To have full control over the distributions, we implemented a synthetic
generative model of a proof-of-concept binary classification task, where a NN
$f$ has to predict whether an individual’s yearly income is greater than $35K$
dollars or not, based on: the working hours per week ($hpw$), the years of
working experience ($yexp$) and the hourly wage ($hw$). Additionally, the
ground-truth joint distribution includes a binary variable indicating gender.
The predictive task we considered is inspired by the Adult UCI dataset, a
common benchmark in fairness verification, modified in order to jointly
support other verification tasks. In the generative model, we implemented a
toggleable mechanism that statistically penalizes the hourly wage of female
individuals. We generated both biased and unbiased versions of separate
training datasets for the population model $p$ and $f$, with $10000$ instances
each. In our experiment $p$ is the combination of two DETs trained separately
on individuals of different gender, achieving a log-likelihood on a test set
close to $1\cdot 10^{-4}$. We quantified multiple properties of
$y=f(hpw,hw,yexp)$ throughout the training, reducing the PFV tasks to WMI and
leveraging the exact solver SA-WMI-PA [Spallitta et al., 2022]. We implemented
$f$ in PyTorch [Paszke et al., 2019] as a feed-forward ReLU architecture with
two latent layers of size $8$, trained with SGD minimizing the binary cross-
entropy loss. The size of $f$ is comparable with the models used in the
FairSquare benchmarks and was expressive enough to achieve high empirical
accuracy. Notice that we can avoid encoding the last sigmoid, defyining the
output of $f$ on the logit value
$(f(\mathbf{x})=1)\leftrightarrow(logit(\mathbf{x})>0)$.

Using $p$, we monitored the demographic parity (5) and the robustness of $f$
to noise in the input. Assuming that individuals tend to overestimate certain
features, we modelled an additive noise to $yexp$ with a triangular
distribution with mode $0$ and support $[0,15]$ and computed how likely the
output would not change in spite of statistically incorrect input. We
additionally verified the correctness of $f$ by quantifying its monotonicity
(9) with respect to each feature. Given two inputs $\mathbf{x}$ and
$\mathbf{x}^{\prime}$ differing on a single feature $i$, we measured how
likely it is that $f(\mathbf{x})>f(\mathbf{x}^{\prime})$ if
$x_{i}>x_{i}^{\prime}$, using a uniform prior in place of $p$.

Despite the simplicity of the predictive task, our results clearly show that
looking at the empirical accuracy or a single formal property can be deceptive
(Fig. 4). We remark that the probabilities computed with our procedure are
exact: aggregate results over multiple runs would only provide information on
the training stability. The networks $f$ were trained for $1000$ epochs,
converging to an accuracy of $0.9$ around epoch $300$. The plot focuses on the
first epochs for increased readability. The accuracy curves suggest that the
unbiased version seems easier to learn for $f$. Yet, solely looking at the
empirical performance is misleading: the monotonicity of the unbiased model
with respect to $hw$ is very low, with a 20% chance of not respecting the
monotonicity constraint. Moreover, the unbiased network is not as robust to
noise in $yexp$ as the biased one, suggesting that it learned to disregard
$hw$ in its predictions. In the biased setting, the network displayed higher
robustness and correctly learned all the monotonicity constraints. However,
with the increasing predictive performance we can observe also a drastic
decrease in the demographic parity, which manifests in a difference in $hw$
among individuals of different gender.

Figure 4:  Quantitative evaluation while training with unbiased (_left_) and
biased (_right_) datasets. Monotonicity results mostly overlap, indicating
that the NN correctly learned the relation, except for _hours per week_ in the
biased case.

While more specialized techniques for the verification of single properties
exist, this is the first time that a single computational procedure is used to
jointly address fundamentally different tasks and provide a diverse set of
insights. For instance, FairSquare [Albarghouthi et al., 2017] was used to
verify group fairness in similar settings but 1) it is more restricted in
terms of encodable population models, requiring them to be defined by means of
univariate priors; 2) it is not fully quantitative, only approximately
answering whether 5 is greater than a given threshold, 3) it cannot be used to
verify other properties. The same holds for the other properties we
considered. Overall, these experiments confirm the potential of the proposed
solution as a unifying framework for jointly addressing multiple diverse PFV
tasks.

The source code and the material needed to reproduce the experiments are
freely available online333Available in the supplementary material. It will be
made publicly available upon acceptance.. The details of the experimental
setup are available in the supplementary material.

##  6 Conclusion and future work

Our WMI-based perspective on the probabilistic verification of AI systems has
many benefits. WMI enables flexible probabilistic reasoning over a broad range
of probabimlistic models, including distributions over mixed logical/numerical
domains (D1). These probabilistic models can be combined with multiple machine
learning models, using a representation language that is commonly used in PFV.
This aspect shows promise in addressing (D2). These aspects, paired with the
ability to verify many different properties (D3), offer unprecedented
flexibility in a single framework.

Practical PFV algorithms need to trade off generality for scalability.
Nonetheless, we believe that this unifying perspective can substantially
contribute to the systematization of this emerging field, and boost further
research on the topic. In the following, we discuss its challenges and
relevant research directions.

#### Scalability.

Being a strict generalization of model counting (#SAT), an archetypal
#P-complete problem, exact WMI is intractable in general. Recent efforts have
been devoted to addressing the challenges related to WMI-based PFV, such as
mitigating the combinatorial explosion of integrals [Spallitta et al., 2022],
providing guarantees on the approximation [Abboud et al., 2022] and
approximating the integrals [Spallitta et al., 2024]. We hope that this
concrete application will become a driving force in this novel but vibrant
research area [Morettin et al., 2021].

While our experiments feature the naive encodings introduced in Sec. 4.1, our
formulation can readily leverage many techniques addressing scalability of ML
verification, including _abstraction_ [Elboher et al., 2020, Ashok et al.,
2020] or _preimage approximation_ [Zhang et al., 2023] techniques. Combining
approximations with guarantees and problem encodings is a crucial step towards
the adoption of WMI-based PFV techniques and a promising research direction
for future work. Further investigation on the factors hindering scalability
might also lead to the development of novel models with a favourable trade-off
between empirical accuracy and verifiability.

#### Arithmetic extensions.

Although WMI is in principle not restricted to $\text{SMT}(\mathcal{LRA})$
formulas, most works have focused on this setting due to the complexity of
reasoning over non-linear constraints while providing bounds on the
approximation. Further extensions of WMI to other theories , or the support of
other families of weight functions, could further push the boundaries of WMI-
based verification.

#### Probabilistic neuro-symbolic verification.

All the properties that we presented are defined on the _concrete input space_
of the system under verification. Yet, being able to define and verify
properties using _semantic_ features is deemed crucial step for advancing the
trustworthiness of our AI systems [Seshia et al., 2016]. Since joint logical
and probabilistic reasoning over $\mathcal{S}$ and $\mathcal{R}$ is achieved
via a unified representation language and computational tools, the classes of
ML model that can be encoded in our framework can appear in the definition of
$\mathcal{R}$. For instance, since convolutional NNs with ReLU can be encoded,
then $\mathcal{R}$ can also be defined in terms of these models, effectively
implementing a probabilistic variant of _neuro-symbolic verification_ [Xie et
al., 2022]. Given a neural predicate that detects stop signs in images, we
could quantify the probability of semantic properties like “If a stop sign is
in front of the camera, a deceleration command is issued” on systems operating
at pixel level.

#### Sequential systems.

In this work, we focused on non-sequential systems. Yet, most ML models are
part of larger software and/or hardware systems with memory. Extending
practical PFV tools like PRISM [Kwiatkowska et al., 2006] and STORM [Dehnert
et al., 2017] with the concepts presented here would significantly advance our
ability to verify AI systems. In this regard, our framework offers a less
cumbersome and error-prone integration process with respect to implementing
multiple system- or property-specific approaches and algorithms into the
existing verification software.

The core idea and formalization is a joint work of PM, AP and RS. PM authored
the code and performed the experiments.

###### Acknowledgements.

Funded by the European Union. Views and opinions expressed are however those
of the author(s) only and do not necessarily reflect those of the European
Union or the European Health and Digital Executive Agency (HaDEA). Neither the
European Union nor the granting authority can be held responsible for them.
Grant Agreement no. 101120763 - TANGO. We acknowledge the support of the PNRR
project FAIR - Future AI Research (PE00000013), under the NRRP MUR program
funded by the NextGenerationEU. The work was partially supported by the
project “AI@TN” funded by the Autonomous Province of Trento.

## References

  * Abboud et al. [2022] Ralph Abboud, İsmail İlkan Ceylan, and Radoslav Dimitrov.  Approximate weighted model integration on dnf structures.  _Artificial Intelligence_ , page 103753, 2022.
  * Albarghouthi et al. [2017] Aws Albarghouthi, Loris D’Antoni, Samuel Drews, and Aditya V Nori.  Fairsquare: probabilistic verification of program fairness.  _Proceedings of the ACM on Programming Languages_ , 1(OOPSLA):1–30, 2017.
  * Amodei et al. [2016] Dario Amodei, Chris Olah, Jacob Steinhardt, Paul Christiano, John Schulman, and Dan Mané.  Concrete problems in ai safety.  _arXiv preprint arXiv:1606.06565_ , 2016.
  * Ashok et al. [2020] Pranav Ashok, Vahid Hashemi, Jan Křetínskỳ, and Stefanie Mohr.  Deepabstract: neural network abstraction for accelerating verification.  In _International Symposium on Automated Technology for Verification and Analysis_ , pages 92–107. Springer, 2020.
  * Baldoni et al. [2011] Velleda Baldoni, Nicole Berline, Jesus De Loera, Matthias Köppe, and Michèle Vergne.  How to integrate a polynomial over a simplex.  _Mathematics of Computation_ , 80(273):297–325, 2011.
  * Baluta et al. [2019] Teodora Baluta, Shiqi Shen, Shweta Shinde, Kuldeep S Meel, and Prateek Saxena.  Quantitative verification of neural networks and its security applications.  In _Proceedings of the 2019 ACM SIGSAC Conference on Computer and Communications Security_ , pages 1249–1264, 2019.
  * Barthe et al. [2011] Gilles Barthe, Pedro R D’argenio, and Tamara Rezk.  Secure information flow by self-composition.  _Mathematical Structures in Computer Science_ , 21(6):1207–1252, 2011.
  * Bastani et al. [2016] Osbert Bastani, Yani Ioannou, Leonidas Lampropoulos, Dimitrios Vytiniotis, Aditya Nori, and Antonio Criminisi.  Measuring neural net robustness with constraints.  _Advances in neural information processing systems_ , 29, 2016.
  * Bastani et al. [2019] Osbert Bastani, Xin Zhang, and Armando Solar-Lezama.  Probabilistic verification of fairness properties via concentration.  _Proceedings of the ACM on Programming Languages_ , 3(OOPSLA):1–27, 2019.
  * Belle et al. [2015] Vaishak Belle, Andrea Passerini, and Guy Van den Broeck.  Probabilistic inference in hybrid domains by weighted model integration.  In _Twenty-Fourth International Joint Conference on Artificial Intelligence_ , 2015.
  * Bertsimas and Dunn [2017] Dimitris Bertsimas and Jack Dunn.  Optimal classification trees.  _Machine Learning_ , 106:1039–1082, 2017.
  * Calders et al. [2009] Toon Calders, Faisal Kamiran, and Mykola Pechenizkiy.  Building classifiers with independency constraints.  In _2009 IEEE international conference on data mining workshops_ , pages 13–18. IEEE, 2009.
  * Chen et al. [2019] Hongge Chen, Huan Zhang, Si Si, Yang Li, Duane Boning, and Cho-Jui Hsieh.  Robustness verification of tree-based models.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Chistikov et al. [2015] Dmitry Chistikov, Rayna Dimitrova, and Rupak Majumdar.  Approximate counting in smt and value estimation for probabilistic programs.  In _International Conference on Tools and Algorithms for the Construction and Analysis of Systems_ , pages 320–334. Springer, 2015.
  * Clarke [1997] Edmund M Clarke.  Model checking.  In _International Conference on Foundations of Software Technology and Theoretical Computer Science_ , pages 54–56. Springer, 1997.
  * Dehnert et al. [2017] Christian Dehnert, Sebastian Junges, Joost-Pieter Katoen, and Matthias Volk.  A storm is coming: A modern probabilistic model checker.  In _Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part II 30_ , pages 592–600. Springer, 2017.
  * Devos et al. [2021] Laurens Devos, Wannes Meert, and Jesse Davis.  Verifying tree ensembles by reasoning about potential instances.  In _Proceedings of the 2021 SIAM International Conference on Data Mining (SDM)_ , pages 450–458. SIAM, 2021.
  * Dos Martires et al. [2019] Pedro Zuidberg Dos Martires, Anton Dries, and Luc De Raedt.  Exact and approximate weighted model integration with probability density functions using knowledge compilation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 7825–7833, 2019.
  * Dvijotham et al. [2018] Krishnamurthy Dvijotham, Marta Garnelo, Alhussein Fawzi, and Pushmeet Kohli.  Verification of deep probabilistic models.  _arXiv preprint arXiv:1812.02795_ , 2018.
  * Dwork et al. [2012] Cynthia Dwork, Moritz Hardt, Toniann Pitassi, Omer Reingold, and Richard Zemel.  Fairness through awareness.  In _Proceedings of the 3rd innovations in theoretical computer science conference_ , pages 214–226, 2012.
  * EC [2021] EC.  Proposal for a regulation of the european parliament and of the council laying down harmonised rules on artificial intelligence (artificial intelligence act) and amending certain union legislative acts.  2021\.  URL https://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:52021PC0206.
  * Einziger et al. [2019] Gil Einziger, Maayan Goldstein, Yaniv Sa’ar, and Itai Segall.  Verifying robustness of gradient boosted models.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pages 2446–2453, 2019.
  * Elboher et al. [2020] Yizhak Yisrael Elboher, Justin Gottschlich, and Guy Katz.  An abstraction-based framework for neural network verification.  In _Computer Aided Verification: 32nd International Conference, CAV 2020, Los Angeles, CA, USA, July 21–24, 2020, Proceedings, Part I 32_ , pages 43–65. Springer, 2020.
  * Gehr et al. [2016] Timon Gehr, Sasa Misailovic, and Martin Vechev.  Psi: Exact symbolic inference for probabilistic programs.  In _Computer Aided Verification: 28th International Conference, CAV 2016, Toronto, ON, Canada, July 17-23, 2016, Proceedings, Part I 28_ , pages 62–83. Springer, 2016.
  * Hardt et al. [2016] Moritz Hardt, Eric Price, and Nati Srebro.  Equality of opportunity in supervised learning.  _Advances in neural information processing systems_ , 29, 2016.
  * Hartmanns and Hermanns [2015] Arnd Hartmanns and Holger Hermanns.  In the quantitative automata zoo.  _Science of Computer Programming_ , 112:3–23, 2015.
  * Huang et al. [2013] Xiaolin Huang, Siamak Mehrkanoon, and Johan AK Suykens.  Support vector machines with piecewise linear feature mapping.  _Neurocomputing_ , 117:118–127, 2013.
  * Huang et al. [2017] Xiaowei Huang, Marta Kwiatkowska, Sen Wang, and Min Wu.  Safety verification of deep neural networks.  In _Computer Aided Verification: 29th International Conference, CAV 2017, Heidelberg, Germany, July 24-28, 2017, Proceedings, Part I 30_ , pages 3–29. Springer, 2017.
  * Katoen [2016] Joost-Pieter Katoen.  The probabilistic model checking landscape.  In _Proceedings of the 31st Annual ACM/IEEE Symposium on Logic in Computer Science_ , pages 31–45, 2016.
  * Katz et al. [2017] Guy Katz, Clark Barrett, David L Dill, Kyle Julian, and Mykel J Kochenderfer.  Reluplex: An efficient smt solver for verifying deep neural networks.  In _International conference on computer aided verification_ , pages 97–117. Springer, 2017.
  * Katz et al. [2019] Guy Katz, Derek A Huang, Duligur Ibeling, Kyle Julian, Christopher Lazarus, Rachel Lim, Parth Shah, Shantanu Thakoor, Haoze Wu, Aleksandar Zeljić, et al.  The marabou framework for verification and analysis of deep neural networks.  In _International Conference on Computer Aided Verification_ , pages 443–452. Springer, 2019.
  * Kolb et al. [2018] Samuel Kolb, Martin Mladenov, Scott Sanner, Vaishak Belle, and Kristian Kersting.  Efficient symbolic integration for probabilistic inference.  In _IJCAI_ , pages 5031–5037, 2018.
  * Kwiatkowska et al. [2006] Marta Kwiatkowska, Gethin Norman, and David Parker.  Quantitative analysis with the probabilistic model checker prism.  _Electronic Notes in Theoretical Computer Science_ , 153(2):5–31, 2006.
  * Liu et al. [2020] Xingchao Liu, Xing Han, Na Zhang, and Qiang Liu.  Certified monotonic neural networks.  _Advances in Neural Information Processing Systems_ , 33:15427–15438, 2020.
  * Mangal et al. [2019] Ravi Mangal, Aditya V Nori, and Alessandro Orso.  Robustness of neural networks: A probabilistic and practical approach.  In _2019 IEEE/ACM 41st International Conference on Software Engineering: New Ideas and Emerging Results (ICSE-NIER)_ , pages 93–96. IEEE, 2019\.
  * Merrell et al. [2017] David Merrell, Aws Albarghouthi, and Loris D’Antoni.  Weighted model integration with orthogonal transformations.  In _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence_ , 2017.
  * Mirman et al. [2021] Matthew Mirman, Alexander Hägele, Pavol Bielik, Timon Gehr, and Martin Vechev.  Robustness certification with generative models.  In _Proceedings of the 42nd ACM SIGPLAN International Conference on Programming Language Design and Implementation_ , pages 1141–1154, 2021.
  * Molina et al. [2018] Alejandro Molina, Antonio Vergari, Nicola Di Mauro, Sriraam Natarajan, Floriana Esposito, and Kristian Kersting.  Mixed sum-product networks: A deep architecture for hybrid domains.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018.
  * Morettin et al. [2019] Paolo Morettin, Andrea Passerini, and Roberto Sebastiani.  Advanced smt techniques for weighted model integration.  _Artificial Intelligence_ , 275:1–27, 2019.
  * Morettin et al. [2020] Paolo Morettin, Samuel Kolb, Stefano Teso, and Andrea Passerini.  Learning weighted model integration distributions.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pages 5224–5231, 2020.
  * Morettin et al. [2021] Paolo Morettin, Pedro Zuidberg Dos Martires, Samuel Kolb, and Andrea Passerini.  Hybrid probabilistic inference with logical and algebraic constraints: a survey.  In _Proceedings of the 30th International Joint Conference on Artificial Intelligence_ , 2021.
  * Narodytska et al. [2018] Nina Narodytska, Shiva Kasiviswanathan, Leonid Ryzhyk, Mooly Sagiv, and Toby Walsh.  Verifying properties of binarized deep neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Poon and Domingos [2011] Hoifung Poon and Pedro Domingos.  Sum-product networks: A new deep architecture.  In _2011 IEEE International Conference on Computer Vision Workshops (ICCV Workshops)_ , pages 689–690. IEEE, 2011.
  * Ram and Gray [2011] Parikshit Ram and Alexander G Gray.  Density estimation trees.  In _Proceedings of the 17th ACM SIGKDD international conference on Knowledge discovery and data mining_ , pages 627–635, 2011.
  * Sang et al. [2005] Tian Sang, Paul Beame, and Henry A Kautz.  Performing bayesian inference by weighted model counting.  In _AAAI_ , volume 5, pages 475–481, 2005.
  * Seshia et al. [2016] Sanjit A Seshia, Dorsa Sadigh, and S Shankar Sastry.  Towards verified artificial intelligence.  _arXiv preprint arXiv:1606.08514_ , 2016.
  * Sill [1997] Joseph Sill.  Monotonic networks.  _Advances in neural information processing systems_ , 10, 1997.
  * Sivaraman et al. [2020] Aishwarya Sivaraman, Golnoosh Farnadi, Todd Millstein, and Guy Van den Broeck.  Counterexample-guided learning of monotonic neural networks.  _Advances in Neural Information Processing Systems_ , 33:11936–11948, 2020.
  * Spallitta et al. [2022] Giuseppe Spallitta, Gabriele Masina, Paolo Morettin, Andrea Passerini, and Roberto Sebastiani.  Smt-based weighted model integration with structure awareness.  In _The 38th Conference on Uncertainty in Artificial Intelligence_ , 2022.
  * Spallitta et al. [2024] Giuseppe Spallitta, Gabriele Masina, Paolo Morettin, Andrea Passerini, and Roberto Sebastiani.  Enhancing smt-based weighted model integration by structure awareness.  _Artificial Intelligence_ , page 104067, 2024.
  * Tjeng et al. [2017] Vincent Tjeng, Kai Xiao, and Russ Tedrake.  Evaluating robustness of neural networks with mixed integer programming.  _arXiv preprint arXiv:1711.07356_ , 2017.
  * Vardi [1985] Moshe Y Vardi.  Automatic verification of probabilistic concurrent finite state programs.  In _26th Annual Symposium on Foundations of Computer Science (SFCS 1985)_ , pages 327–338. IEEE, 1985.
  * Wehenkel and Louppe [2019] Antoine Wehenkel and Gilles Louppe.  Unconstrained monotonic neural networks.  _Advances in neural information processing systems_ , 32, 2019.
  * Weng et al. [2019] Lily Weng, Pin-Yu Chen, Lam Nguyen, Mark Squillante, Akhilan Boopathy, Ivan Oseledets, and Luca Daniel.  Proven: Verifying robustness of neural networks with a probabilistic approach.  In _International Conference on Machine Learning_ , pages 6727–6736. PMLR, 2019.
  * Xie et al. [2022] Xuan Xie, Kristian Kersting, and Daniel Neider.  Neuro-symbolic verification of deep neural networks.  2022\.
  * Zhang et al. [2023] Xiyue Zhang, Benjie Wang, and Marta Kwiatkowska.  On preimage approximation for neural networks.  _arXiv preprint arXiv:2305.03686_ , 2023.

A Unified Framework for Probabilistic Verification of AI Systems
via Weighted Model Integration
(Supplementary Material)

##  Appendix A Experimental setup

###  A.1 Ground truth distribution

$yexp$$hpw$$hw$$female$$income$$y$ Figure 5: Graphical model of the generating distribution. The observed variables in the sampled data are circled. The dependency $female\rightarrow hw$ is only present in the biased dataset. Female gender | $female\sim Bernoulli(0.5)$
---|---
Work hours per week | $hpw=5\cdot max(min(2,\mathcal{N}(7,1)),12)$
Years of working experience | $yexp=max(min(0,\mathcal{N}(20,7)),55)$
Hourly wage | $hw=\mathcal{N}(1,0.02)\cdot\frac{500\cdot yexp+30000}{1920}$
Biased wage | $hw_{b}=hw+hw\cdot N(01,0.02)\cdot(1-2\cdot\llbracket feamale\rrbracket)$
Yearly income | $income=48\cdot hpw\cdot hw$
Binary class | $y=(income>35000)$
Table 1: Full specification of the probabilistic models. Figure 6: Data
distributions used to train the population model.

###  A.2 Models and training

For these experiments, we borrowed the DET implementation and WMI encoding of
Morettin et al. 2020. The population models were trained on $10000$ instances
using hyperparameters $n_{min}=1000$ and $n_{max}=2000$. These hyperparameters
respectively dictate the minimum and maximum number of instances per leaf. The
binary predictors were implemented in PyTorch (Paszke et al., 2019) as fully
connected feedforward networks with two latent layers having $8$ neurons each.
The networks were trained on $10000$ labelled instances minimizing the binary
cross entropy loss with stochastic gradient descent (batch size $200$,
learning rate $0.005$ and no momentum). The predictive performance of the
learned models was assessed on an identically distributed test set with $1000$
instances. The WMI problems resulting from the demographic parity,
monotonicity and robustness queries were solved using the exact solver SA-WMI-
PA (Spallitta et al., 2022). All the experiments were run on a machine with 28
Core (2.20GHz frequency) and 128GB of ram, running Ubuntu Linux 20.04.


# Single-hemisphere dynamos in M-dwarf stars

Benjamin P. Brown Department Astrophysical and Planetary Sciences & LASP,
University of Colorado, Boulder, CO 80309, USA  Jeffrey S. Oishi Department of
Physics and Astronomy, Bates College, Lewiston, ME 04240, USA  Geoffrey M.
Vasil School of Mathematics and Statistics, University of Sydney, Sydney,
Sydney, NSW 2006, AUS  Daniel Lecoanet Princeton Center for Theoretical
Science & Department of Astrophysical Sciences, Princeton University,
Princeton, NJ 08544, USA  Keaton J. Burns Department of Mathematics,
Massachusetts Institute of Technology, Cambridge, MA 02139, USA  Center for
Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA

###### Abstract

M-dwarf stars below a certain mass are convective from their cores to their
photospheres. These fully convective objects are extremely numerous, very
magnetically active, and the likely hosts of many exoplanets. Here we study,
for the first time, dynamo action in simulations of stratified, rotating fully
convective M-dwarf stars. Importantly, we use new techniques to capture the
correct full ball geometry down to the center of the star. We find surprising
dynamo states in these systems, with the global-scale mean fields confined
strongly to a single hemisphere, in contrast to prior stellar dynamo
solutions. These hemispheric-dynamo stars are likely to have profoundly
different interactions with their surroundings, with important implications
for exoplanet habitability and stellar spindown.

Astrophysical fluid dynamics(101), Magnetohydrodynamical simulations (1966),
Stellar magnetic fields (1610), M dwarf stars (982)

††journal: ApJL††software:  Dedalus (commit 7efb884, Burns et al., 2020,
http://dedalus-project.org), Dedalus-Sphere (commit c663ff4, Vasil et al.,
2019; Lecoanet et al., 2019,
https://github.com/DedalusProject/dedalus_sphere), MESA (release 11701, Paxton
et al., 2011, 2013, 2015, 2018, http://mesa.sourceforge.net/).

##  1 Introduction

Tiny M-dwarfs are the most numerous stars in the universe and one of the most
promising locations for finding habitable planets. In particular, they have
many Earth-sized habitable planets, and the nearest habitable exoplanet, when
found, is likely to orbit an M-dwarf (e.g., Dressing & Charbonneau, 2015).
M-dwarfs are intrinsically dim, and observations of them are hard, but our
data set is growing rapidly and now includes long-duration campaigns of
slowly-rotating systems (Newton et al., 2018). The habitability of potential
exo-Earths around other stars is at the mercy of the magnetic environment.
M-dwarfs, in particular, display multi-kilogauss surface fields. Photometric
monitoring infers incredibly large starspots (e.g., Newton et al., 2016) or
many, many small spots (e.g., Jackson & Jeffries, 2012), and the combination
of strong fields and complex topologies leads exceptionally high levels of
magnetic activity (e.g., West et al., 2015). Mega-flares larger than any ever
seen on the Sun (e.g., Schmidt et al., 2019) happen often enough to threaten
the ability of planets to support life, or possibly even to retain atmospheres
(e.g., Shkolnik & Barman, 2014).

The interiors of type-M4 stars or later (masses of $\sim 0.35M_{\odot}$ or
less) are convectively unstable from their nuclear burning cores to the
photosphere. Despite the relative simplicity of their internal structure, we
know very little about the internal processes that drive their magnetic
dynamos. Only a limited number of simulations have studied global-scale dynamo
processes in M-dwarfs (e.g., Browning, 2008; Yadav et al., 2015a, b, 2016).
Worse, all of these have suffered a fundamental limitation: heretofore it has
only been possible to simulate spherical shells, with the coordinate
singularity at $r=0$ avoided by making a “cut-out” in the middle of the
simulation. The impact of this cut-out on the global solution has been
unknown. Using new techniques in the open-source Dedalus pseudospectral
framework, we compute the first spherical solutions that include $r=0$ without
any cut-out. These are the first dynamo models of truly, fully convective
M-dwarf stars.

As in solar-type stars, our simulations of M-dwarfs have wreath-building
convection zone dynamos that create global-scale fields. However, our M-dwarf
simulations find one significant difference: We find an abundance of dynamo
solutions localized to a single hemisphere. Single-hemisphere dynamo states
have profound implications for stellar spin-down, exoplanet habitability, and
the interpretation of observations of strong fields on stellar surfaces.

##  2 Model system

The interiors of M-dwarfs are regions of low-Mach number convection with
significant density stratification. To study convection and dynamo action in
these stars, we solve the anelastic equations using the open-source,
spectrally-accurate Dedalus framework (Burns et al., 2020). This work is
enabled by new approaches to global-spherical domains including $r=0$ (Vasil
et al., 2019; Lecoanet et al., 2019).

The non-dimensional energy-conserving anelastic equations with MHD and
momentum and thermal diffusion are (Brown et al., 2012; Vasil et al., 2013;
Lecoanet et al., 2014):

$\frac{\partial\boldsymbol{u}}{\partial
t}+\boldsymbol{\nabla}\varpi-\frac{\mathrm{Ra}\mathrm{Ek}^{2}}{\mathrm{Pr}}T\nabla
s-\mathrm{Ek}\frac{\left(\nabla^{2}\boldsymbol{u}+\text@frac{1}{3}\nabla(\nabla\cdot\boldsymbol{u})\right)}{\rho}=\\\
-\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{u}-\boldsymbol{\hat{z}}\times\boldsymbol{u}+\frac{\boldsymbol{J}\times\boldsymbol{B}}{\rho},$
(1) $\frac{\partial s}{\partial
t}-\frac{\mathrm{Ek}}{\mathrm{Pr}}\frac{\left(\nabla^{2}s+\boldsymbol{\nabla}\ln
T\cdot\boldsymbol{\nabla}s\right)}{\rho}=\\\
-\boldsymbol{u}\cdot\boldsymbol{\nabla}\boldsymbol{s}+\frac{\mathrm{Ek}}{\mathrm{Pr}}\mathcal{S}+\frac{\mathrm{Pr}}{\mathrm{Ra}\mathrm{Ek}}\frac{\Phi_{\nu}}{\rho
T}+\frac{1}{\mathrm{Pr}\mathrm{Pm}}\frac{1}{\mathrm{Ra}\mathrm{Ek}}\frac{J^{2}}{\rho^{2}T},$
(2) $\frac{\partial\boldsymbol{A}}{\partial
t}+\nabla\phi-\frac{\mathrm{Ek}}{\mathrm{Pm}}\frac{\nabla^{2}\boldsymbol{A}}{\rho}=\boldsymbol{u}\times\boldsymbol{B},$
(3)

together with the Coulomb gauge $\boldsymbol{\nabla}\cdot\boldsymbol{A}=0$ and
the anelastic
$\boldsymbol{\nabla}\cdot\boldsymbol{u}+\boldsymbol{u}\cdot\boldsymbol{\nabla}\ln\rho=0$
constraints. $\mathcal{S}$ is a volume heating term and rotation is
represented as $\boldsymbol{\Omega}=\Omega\boldsymbol{\hat{z}}$. The magnetic
field $\boldsymbol{B}=\nabla\times\boldsymbol{A}$, and the electrostatic
potential $\phi$ and reduced pressure $\varpi$ enforce the constraints. We use
the ideal gas equation of state, and diffusive entropy heat flux
$\boldsymbol{Q}=-\kappa T\boldsymbol{\nabla}s$ (Lecoanet et al., 2014). The
boundaries are constant entropy, impenetrable and stress free for velocity,
and potential for magnetic field. The quantities $\nabla\ln\rho$, $\nabla T$,
$1/\rho$, $\nabla\ln T$ and $\mathcal{S}$ represent the internal M-dwarf
structure.

We non-dimensionalize this set of equations on a characteristic lengthscale
$L=R$ the radius of the domain ($0.85R_{\star}$), rotation timescale
$\tau=1/(2\Omega)$, density $\rho_{c}$, magnetic field
$B^{2}=16\pi\rho_{c}\Omega^{2}R^{2}$, and entropy
$\mathcal{S}_{c}\sigma^{2}R^{2}/(\kappa/\rho_{c})$. The dynamic momentum,
$\mu$, and thermal diffusivities, $\kappa$, are constant. The magnetic
diffusivity $\eta\propto\rho^{-1}$. The global non-dimensional parameters are
the Ekman number $\mathrm{Ek}$, Rayleigh number $\mathrm{Ra}$, Prandtl number
$\mathrm{Pr}$ and magnetic Prandtl number $\mathrm{Pm}$,

$\displaystyle\mathrm{Ek}$ $\displaystyle=$
$\displaystyle\frac{\mu/\rho_{c}}{2\Omega
R^{2}},\quad\mathrm{Ra}=\frac{\mathcal{S}_{c}\sigma^{2}R^{2}T_{c}}{\mu\kappa^{2}/\rho_{c}^{3}},$
(4) $\displaystyle\quad\mathrm{Pr}$ $\displaystyle=$
$\displaystyle\frac{\mu}{\kappa},\quad\ \
\quad\,\mathrm{Pm}=\frac{\mu}{\rho_{c}\eta},$ (5)

where $\sigma$ represents the radial extent of the nuclear-burning region. The
Rayleigh number has this dependance on the heating because mean entropy
contrast across the domain scales with the internal heating to within an
$\mathcal{O}(1)$ constant as $\Delta
s\sim\rho\,\mathcal{S}_{c}\sigma^{2}/\kappa$. The values are specified at the
center $r=0$ (denoted with a subscript $c$), which is the extremal value for
$\mathrm{Ek}$ and $\mathrm{Ra}$. In our simulations, the local Ekman
$\mathrm{Ek}$ and Rayleigh $\mathrm{Ra}$ number depend on $\rho(r)/\rho_{c}$.
The Prandtl numbers $\mathrm{Pr}$ and $\mathrm{Pm}$ are constant.

Figure 1: Time-latitude plots of the global-scale dynamo. Top is the toroidal
mangetic field $B_{\phi}$ at mid-convection zone ($r=0.5$). Bottom is the
vector potential of the poloidal magnetic field $A_{\phi}$ at the top of the
domain ($r=1$). Color gives the sense of polarity for both quantities. The
system goes through many dynamo reversals of the global-scale field during the
period simulated. The mean fields are predominantly in a single hemisphere
(here the Northern). At all times the global poloidal fields have a dipolar
component, with strong contributions from higher-order modes.

We use the convective Rossby number as an input proxy for the degree of
rotational constraint,

$\mathrm{Ro}_{\mathrm{C}}^{2}\equiv\frac{\mathrm{Ra}\mathrm{Ek}^{2}}{\mathrm{Pr}}=\frac{\mathcal{S}_{c}\sigma^{2}T_{c}}{(\kappa/\rho_{c})4\Omega^{2}R^{2}}.$
(6)

This combination has a long history of usefulness for boundary-heated
incompressible convection (Julien et al., 1996). In that context, the
combination is independent of microphysical diffusivities and depends on the
domain size. The situation here is less simple. The internal heating
structure, $\mathcal{S}(r)$, and radiative diffusion $\kappa(r)$, determine
the superadiabatic entropy gradient. Even thought our definition produces a
more complex parameter dependance, we find good scaling correlation between
$\mathrm{Ro}_{\mathrm{C}}$ and the output measured local Rossby number; see
Anders et al. (2019) for more details.

Fully-convective M-dwarf stars are nearly adiabatic, and are well represented
by polytropic solutions to the Lane-Emden equations. We simulate the inner 3
density scale heights ($n_{\rho}=3$) of an adiabatically stratified ($m=1.5$)
Lane-Emden polytrope, corresponding to the inner 85% of the star. Our
background state is in excellent agreement with a stellar structure model
computed with MESA (Paxton et al., 2011) for a $0.3M_{\odot}$ star of age 5
Gyrs, at a solar metallicity.111In contrast to prior work (Yadav et al.,
2015a, b, 2016) which erroneously assumed a polytropic structure with a
gravity profile that grows linearly with radius. For $\mathcal{S}$, we extract
the volumetric heating and radiative cooling from the MESA model, and fit the
radial dependence with

$\mathcal{S}=\frac{\epsilon_{\mathrm{nuc}}}{T}-\frac{1}{\rho
T}\nabla\cdot\boldsymbol{F}_{\mathrm{rad}}\approx\left[\left(\frac{Q_{0}}{Q_{1}}\right)e^{-r^{2}/(2\sigma^{2})}+1\right]Q_{1}.$
(7)

A nonlinear least squares fit gives $\sigma\approx 0.12$ and
$Q_{0}/Q_{1}\approx 11$, with an $L_{2}$ error of about 3%. The simulations
are normalized such that $\mathcal{S}(r=0)\sigma^{2}=1$.

In this letter we report on a solution with $\mathrm{Ro}_{\mathrm{C}}=0.41$,
$\mathrm{Ek}=2.5\times 10^{-5}$ and $\mathrm{Pr}=\mathrm{Pm}=1$. The system
starts in an adiabatic state with small entropy noise perturbations. We
initialize a weak poloidal magnetic field with spherical harmonic mode
$\ell=1,2$, $m=1$. The simulation has $L=127$ and $N=127$, where $L$ is the
maximal degree of the spin-weighted spherical harmonics used in the spectral
representation and $N$ is the maximal degree of the Jacobi polynomials used in
the radial expansion. We computed many such dynamo simulations, but focus our
attention on a single representative case.

The convection is a successful dynamo, building both small- and global-scale
magnetic fields which change in strength and orientation over time. This
simulation ran for $t=6250\,\Omega^{-1}$, or about 0.3 diffusion times at the
core $r=0$ and about 6.3 diffusion times at $r=1$. The internal heating from
$\mathcal{S}$ drives a convectively unstable state. By roughly $t\approx
2000\,\Omega^{-1}$, the potential, kinetic and magnetic energies all
equilibrate; fluctuating around well-defined averages.

Figure 2: Magnetic field topologies through a magnetic reversal. The left side
shown the mean toroidal field $B_{\phi}$ and poloidal vector potential
$A_{\phi}$ for the full domain, temporally averaged over one rotation period.
The right side shows the radial magnetic field $B_{r}$ at the upper boundary
$r=1$. These are taken from periods of peak magnetic field strength during two
successive dynamo cycles; the top row from $t=5950\,\Omega^{-1}$, bottom row
from $t=6025\,\Omega^{-1}$.

After time-averaging over the second half of the simulation, the RMS Reynolds
number $\approx$ 240; the fluctuating Reynolds number $\approx$ 120; the RMS
Rossby number $Ro=|\nabla\times\boldsymbol{u}|/2\Omega\approx 0.3$. In this
letter, mean quantities are the $m=0$ components (denoted with angle
brackets). Fluctuating quantities are those with the mean removed. The
time–volume-averaged total magnetic energy is $\approx$ 25% or the total
kinetic energy. Approximately 25% of the magnetic energy exists in mean
toroidal fields, $\approx$ 68% exists in fluctuating fields, and the balance
is in mean poloidal fields. The differential rotation contains $\approx$ 45%
of the kinetic energy, with $\approx$ 55% in the fluctuating convection.

##  3 Hemispheric dynamos

Figure 1 shows the surprising global-scale mean fields. The fields are
organized in space and in time, and undergo regular periodic reversals of
polarity in both the mean toroidal $\langle B_{\phi}\rangle$ and the mean
poloidal field $\langle A_{\phi}\rangle$. The fields share similarities with
the “wreathy dynamo” states found in dynamo simulations of solar-type stars
(e.g., Brown et al., 2010, 2011; Nelson et al., 2013). In the fully convective
dynamo simulations, the mean fields represent a significant fraction of the
total magnetic field. We measure the global magnetic reversals using Lomb-
Scargle periodograms of the mean magnetic field (e.g., Townsend, 2010;
VanderPlas, 2018). The fields change with a period of $P\approx
190\,\Omega^{-1}$, which corresponds to a polarity reversal roughly every 15
rotation periods. The temporal consistency of these cycles is high. The peak
of the periodogram is clearly present in both $\langle B_{\phi}\rangle$
(normalized peak power $\approx 0.3$, $P\approx 190.1\,\Omega^{-1}$) and
$\langle A_{\phi}\rangle$ (peak power $\approx 0.4$, $P\approx
191.8\,\Omega^{-1}$), while the aliasing peaks are at much lower amplitude
($\lesssim 0.05$ and $\lesssim 0.1$ respectively). The surface-poloidal and
depth-toroidal fields change with nearly the same phase.

It is surprising is that the mean fields of these fully-convective dynamos
generally occupy one hemisphere or the other, but not both. In rapidly
rotating solar-type stars with convection-zone shells fields typically
symmetrically fill both hemispheres (e.g., Brown et al., 2010, 2011; Nelson et
al., 2013; Augustson et al., 2013, 2015). Importantly: all previous M-dwarf
dynamo simulations in deep shells also find two-hemisphere solutions. Those
simulations have either found mean fields symmetrically in both hemispheres
(e.g., Yadav et al., 2015a, 2016) or have built mean fields that don’t undergo
very many reversals (Browning, 2008).

The tendency towards mean fields in a single hemisphere is strong: the initial
dynamo action ($t\lesssim 125\,\Omega^{-1}$) is symmetric in the two
hemispheres, but strong southern-hemisphere fields rapidly replace the initial
transients. At about $t\approx 1500\,\Omega^{-1}$ those strong fields shift to
the northern hemisphere, where they remain for the period we have simulated.
The other hemisphere is not devoid of field, but the mean fields there are
much weaker than in the dominant hemisphere. We speculate that the reversals
would continue on some longer timescale if allowed.

Figure 2 shows the detailed structure of the global magnetic structure during
two phases of reversed polarity late in the simulation. The first happens
around $t\approx 5950\,\Omega^{-1}$ averaged over one rotation period ($\Delta
t\approx 2\pi\,\Omega^{-1}$) when positive (red) polarity dominates the
northern hemisphere. We see that a positive $\langle B_{\phi}\rangle$ wreath
fills the star, reaching from the core to the surface. North of this structure
and closer to the surface lie the weak remnant wreaths of the previous cycle.
A distinctly weaker and non-dipolar poloidal field $\langle A_{\phi}\rangle$
accompanies the dominant toroidal structure. There is some global-dipole
component, but it is much weaker than the higher-order modes. The near-surface
cuts also show the poloidal-field structure; indicating a substantial
quadrupolar component. As the wreath migrates up along the rotation axis, a
new opposite-polarity wreath (blue) replaces it. This reversal is visible in
the surface radial fields.

Figure 3: Properties of flows and heat transport in the fully convective
simulation. Left: a cut across the equator of the fluctuating entropy with the
$m=0$ average subtracted. Center: Luminosity profiles throughout the
simulation. The radial profiles are averaged horizontally and in time. Right:
These flows drive a strong solar-like differential rotation; fast near the
equator and slow near the rotation axis. These flows are shown at time
$t=6025\,\Omega^{-1}$, corresponding to the lower row of Figure 2; the
differential rotation is time-averaged over about 10 rotation periods,
centered on this time.

The convective flows fill the interior, as visible in Figure 3. Here we show
entropy fluctuations in a cut at the equator. Cold, low-entropy plumes (blue)
fall from the near-surface regions. The core comprises rising warm, high-
entropy plumes (red). The center at $r=0$ is a region of vigorous dynamics,
and plumes frequently cross the origin. Convective enthalpy transport
($L_{h}$) dominates the radial luminosity balance, except in the upper
boundary layer where thermal diffusion takes over ($L_{Q}$). Though these
simulations are stratified, the kinetic energy flux ($L_{\mathrm{KE}}$) is
small, similar to rotationally-constrained simulations of solar-type stars
(e.g., Brown et al., 2008). The sum of the luminosities closely matches the
nuclear source derived from MESA ($L_{\mathcal{S}}$) indicating that these
models are in proper thermal equilibrium. This convection drives a robust
differential rotation, shown by the mean angular velocity,
$\langle\Omega\rangle$, with a fast equator and slow polar and core regions.
The gradients of angular velocity are strong in both radius and latitude.
These gradients likely build the mean toroidal field $\langle B_{\phi}\rangle$
via shear induction (a.k.a the ”$\Omega$-effect”). Correlations in the
fluctuating flows and fields maintain the polar field. The latitudinal
contrast in $\langle\Omega\rangle$ is about 10%, and the mean magnetic fields
coexist with this shearing flow. This is similar to results from rapidly
rotating simulations of solar-type stars (Brown et al., 2011), but markedly
different from the results of Browning (2008) or Yadav et al. (2015a) where
the fields strongly quench the differential rotation. Whether differential
rotation and fields can coexist likely depends on the Rossby number regimes of
the dynamos, as discussed in Yadav et al. (2016).

##  4 Implications of Hemispheric Dynamos

This letter has focused on a single dynamo solution that shows striking
hemispheric asymmetries. We find the same preference for single-hemisphere
dynamos in every fully convective M-dwarf case computed to present; including
a sweep of cases at $\mathrm{Ek}=2.5\times 10^{-5}$ and
$\mathrm{Ro}_{\mathrm{C}}=0.38$–$0.57$, corresponding to RMS Reynolds numbers
of $210$–$350$ and Rossby numbers of $0.27$–$0.42$. To check if single-
hemisphere preference happens spuriously from a numerical source, we also
rotated around the $\boldsymbol{\hat{x}}$ rather than the
$\boldsymbol{\hat{z}}$ axis (for $\mathrm{Ro}_{\mathrm{C}}=0.41$,
$\mathrm{Ek}=2.5\times 10^{-5}$). In this strange solution, cycling single-
hemispheric dynamo action happened around the east-west pole, relative to the
numerical grid.

So far, we cannot determine why fully convective simulations so strongly
prefer single-hemisphere dynamos; nor why these have not been found
previously. Possibly, hemispheric dynamos are local to this region of
parameter space. It is also possible that small core cutouts, as used in shell
simulations, significantly change the allowed topologies of the global-scale
fields even when those cutouts are small ($\lesssim 10\%R$). Small solid
objects are a relevant source of vorticity injection. In any case, we find
that hemispheric dynamo states are robust and possibly ubiquitous in fully
convective simulations, and this prompts us to consider their impacts if they
exist in real stars.

Strong global-scale mean fields confined to a single hemisphere will couple
dramatically differently to the outflowing stellar wind than a traditional
dipole. The higher-order moments of field fall of much faster with distance
and this likely produces weaker angular momentum transport and less efficient
spindown in hemispheric-dynamo stars. The presence of mean-field in a single
hemisphere could dramatically affect attempts to observe this field and lead
to significant viewing angle effects; this could explain substantial
differences in observed topologies for otherwise similar stars (e.g., Morin et
al., 2010).

Lastly, we expect that the habitability around a hemispheric-dynamo star might
significantly differ from one with a dipolar field. Hemispheric-dynamo stars
would be less likely to turn nearby exoplanets into lava worlds via inductive
heating (e.g., Kislyakova et al., 2017, 2018), as the stellar field at the
location of the planetary orbit will be much weaker. But the bad news is there
might be much higher levels of magnetic activity, flares and coronal mass
ejections; the mean-field topologies are more complicated than simple dipoles.
Alternatively, though, space weather might show a strong directional
preference. A planet with a low-inclination orbit might escape the worst. But
even then we anticipate the preferred hemisphere to switch often enough to
cause trouble for everyone. Altogether, fully convective stars have a rich
landscape of dynamo states, and exploring this landscape shines new light on
our understanding of the dynamos operating in solar-type stars. It seems many
questions remain wide open.

This work was supported by NASA LWS grant NNX16AC92G, NASA SSW grant
80NSSC19K0026, and NASA HTMS grant 80NSSC20K1280. Computations were conducted
with support by the NASA High End Computing (HEC) Program through the NASA
Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades
with allocation GIDs s1647 and s2114.

## References

  * Anders et al. (2019) Anders, E. H., Manduca, C. M., Brown, B. P., Oishi, J. S., & Vasil, G. M. 2019, ApJ, 872, 138, doi: 10.3847/1538-4357/aaff61
  * Augustson et al. (2015) Augustson, K., Brun, A. S., Miesch, M., & Toomre, J. 2015, ApJ, 809, 149, doi: 10.1088/0004-637X/809/2/149
  * Augustson et al. (2013) Augustson, K. C., Brun, A. S., & Toomre, J. 2013, ApJ, 777, 153, doi: 10.1088/0004-637X/777/2/153
  * Brown et al. (2008) Brown, B. P., Browning, M. K., Brun, A. S., Miesch, M. S., & Toomre, J. 2008, ApJ, 689, 1354, doi: 10.1086/592397
  * Brown et al. (2010) —. 2010, ApJ, 711, 424, doi: 10.1088/0004-637X/711/1/424
  * Brown et al. (2011) Brown, B. P., Miesch, M. S., Browning, M. K., Brun, A. S., & Toomre, J. 2011, ApJ, 731, 69, doi: 10.1088/0004-637X/731/1/69
  * Brown et al. (2012) Brown, B. P., Vasil, G. M., & Zweibel, E. G. 2012, ApJ, 756, 109, doi: 10.1088/0004-637X/756/2/109
  * Browning (2008) Browning, M. K. 2008, ApJ, 676, 1262, doi: 10.1086/527432
  * Burns et al. (2020) Burns, K. J., Vasil, G. M., Oishi, J. S., Lecoanet, D., & Brown, B. P. 2020, Physical Review Research, 2, 023068, doi: 10.1103/PhysRevResearch.2.023068
  * Dressing & Charbonneau (2015) Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45, doi: 10.1088/0004-637X/807/1/45
  * Jackson & Jeffries (2012) Jackson, R. J., & Jeffries, R. D. 2012, MNRAS, 423, 2966, doi: 10.1111/j.1365-2966.2012.21119.x
  * Julien et al. (1996) Julien, K., Legg, S., McWilliams, J., & Werne, J. 1996, Journal of Fluid Mechanics, 322, 243, doi: 10.1017/S0022112096002789
  * Kislyakova et al. (2018) Kislyakova, K. G., Fossati, L., Johnstone, C. P., et al. 2018, ApJ, 858, 105, doi: 10.3847/1538-4357/aabae4
  * Kislyakova et al. (2017) Kislyakova, K. G., Noack, L., Johnstone, C. P., et al. 2017, Nature Astronomy, 1, 878, doi: 10.1038/s41550-017-0284-0
  * Lecoanet et al. (2014) Lecoanet, D., Brown, B. P., Zweibel, E. G., et al. 2014, ApJ, 797, 94, doi: 10.1088/0004-637X/797/2/94
  * Lecoanet et al. (2019) Lecoanet, D., Vasil, G. M., Burns, K. J., Brown, B. P., & Oishi, J. S. 2019, Journal of Computational Physics: X, 100012
  * Morin et al. (2010) Morin, J., Donati, J. F., Petit, P., et al. 2010, MNRAS, 407, 2269, doi: 10.1111/j.1365-2966.2010.17101.x
  * Nelson et al. (2013) Nelson, N. J., Brown, B. P., Brun, A. S., Miesch, M. S., & Toomre, J. 2013, ApJ, 762, 73, doi: 10.1088/0004-637X/762/2/73
  * Newton et al. (2016) Newton, E. R., Irwin, J., Charbonneau, D., et al. 2016, ApJ, 821, 93, doi: 10.3847/0004-637X/821/2/93
  * Newton et al. (2018) Newton, E. R., Mondrik, N., Irwin, J., Winters, J. G., & Charbonneau, D. 2018, AJ, 156, 217, doi: 10.3847/1538-3881/aad73b
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Schmidt et al. (2019) Schmidt, S. J., Shappee, B. J., van Saders, J. L., et al. 2019, ApJ, 876, 115, doi: 10.3847/1538-4357/ab148d
  * Shkolnik & Barman (2014) Shkolnik, E. L., & Barman, T. S. 2014, AJ, 148, 64, doi: 10.1088/0004-6256/148/4/64
  * Townsend (2010) Townsend, R. H. D. 2010, ApJS, 191, 247, doi: 10.1088/0067-0049/191/2/247
  * VanderPlas (2018) VanderPlas, J. T. 2018, ApJS, 236, 16, doi: 10.3847/1538-4365/aab766
  * Vasil et al. (2013) Vasil, G. M., Lecoanet, D., Brown, B. P., Wood, T. S., & Zweibel, E. G. 2013, ApJ, 773, 169, doi: 10.1088/0004-637X/773/2/169
  * Vasil et al. (2019) Vasil, G. M., Lecoanet, D., Burns, K. J., Oishi, J. S., & Brown, B. P. 2019, Journal of Computational Physics: X, 100013
  * West et al. (2015) West, A. A., Weisenburger, K. L., Irwin, J., et al. 2015, ApJ, 812, 3, doi: 10.1088/0004-637X/812/1/3
  * Yadav et al. (2015a) Yadav, R. K., Christensen, U. R., Morin, J., et al. 2015a, ApJ, 813, L31, doi: 10.1088/2041-8205/813/2/L31
  * Yadav et al. (2016) Yadav, R. K., Christensen, U. R., Wolk, S. J., & Poppenhaeger, K. 2016, ApJ, 833, L28, doi: 10.3847/2041-8213/833/2/L28
  * Yadav et al. (2015b) Yadav, R. K., Gastine, T., Christensen, U. R., & Reiners, A. 2015b, A&A, 573, A68, doi: 10.1051/0004-6361/201424589