# Three-Phase Infeasibility Analysis for Distribution Grid Studies

Elizabeth Foster, Amritanshu Pandey, and Larry Pileggi  Dept. of Electrical
and Computer Engineering
Carnegie Mellon University
Pittsburgh, PA USA
{emfoster, amritanp<EMAIL_ADDRESS>
###### Abstract

With the increase of distributed energy resources in the distribution grid,
planning to ensure sufficient infrastructure and resources becomes critical.
Planning at the distribution level is limited by the complexities of
optimizing unbalanced systems. In this paper we develop a three-phase
infeasibility analysis that identifies weak locations in a distribution
network. This optimization is formulated by adding slack current sources at
nodes in the system and minimizing their norm subject to distribution power
flow constraints. Through this analysis we solve instances of power flow that
would otherwise be infeasible and diverge. Under conditions when power flow is
feasible, our approach is equivalent to standard three-phase power flow;
however, for cases where power flow fails, the nonzero slack injection
currents compensate for missing power to make the grid feasible. Since an
uncountable number of injected currents can provide feasibility, we further
explore the optimization formulation that best fits the solution objective
through use of both a least squares and an L1 norm objective. Our L1 norm
formulation localizes power deficient locations through its inherent sparsity.
We show the efficacy of this approach on realistic unbalanced testcases up to
8500 nodes and for a scenario with a high penetration of electric vehicles.

###### Index Terms:

equivalent circuit approach, distribution system, L1-norm sparse solution,
nonlinear optimization, three- phase unbalanced networks

Submitted to the 22nd Power Systems Computation Conference (PSCC 2022).

##  I Introduction

An increasingly complex distribution grid with more distributed energy
resources (DERs) and inverter-based resources necessitates more proactive
planning to ensure there is adequate infrastructure and resources for changing
demand [1]. Planning at the distribution level is often reactive because there
is a limited offering of tools capable of performing three-phase steady-state
optimizations and utilities themselves may lack high fidelity data at the
distribution level [2]. As a greater number of sensors and smart meters are
introduced on the distribution grid and as utilities collect more detailed
information on the distribution grid layout, one of the impediments to more
robust and proactive distribution analysis is the lack of three-phase
optimization methods for large-scale networks.

Figure 1: Representation of the IEEE 4 node case with infeasibility current
sources at each node

When considering planning for future load growth at the grid-edge, simulating
three-phase power flow is not enough to provide advanced planning insight.
Power flow simulations can fail because the testcase misrepresents a real
physical system or because no feasible solution exists for the given set of
operating conditions. It may also fail due to choice of poor initial
conditions or unavailability of good initial conditions [3]. System planners
require clear indications regarding why power flow simulations fail since
different causes correspond to different remedies as it relates to planning
capabilities. We lack methods that can not only identify an infeasible or
collapsed network, but also provide useful insight on the failure.

There is some prior work in identifying infeasible grid settings in both
positive sequence networks [4] [5] [6] and distribution networks [7]. [4]
introduced the concept of missing power at the transmission level based on
power mismatch equations augmented with slack variables. [7] introduced a
similar concept for distribution networks to identify problematic optimal
power flow constraints. More recently, we developed a more robust approach
using equivalent circuit formulation (ECF) for transmission systems in [5].
Our work in [6] built on [5] by adding an L1 norm regularization to induce
greater sparsity in the solution vector. Nonetheless, none of these methods
are robust, scalable, and localized for three-phase distribution networks.
Unbalance between phases is common in distribution networks, so positive
sequence methods in [4] [5] [6] are not directly applicable because they lack
the necessary but complicated phase relationships. [4][7] do not demonstrate
scalability beyond a few hundred nodes. [4] [5] [7] do not attempt to localize
the sources of infeasibility. Difficulties in three-phase distribution systems
further arise from a low X-R ratio; significant number of direct, non-
controllable customers; and significant levels of unbalance [8]. It is an
imminent need in distribution planning studies to have large scale three-phase
analysis capable of detecting localized sources of infeasibility, and none of
the existing works in this space can fill the void.

Therefore, to address this research gap, we propose an ECF-based three-phase
optimization method called the three-phase infeasibility analysis (TPIA) that
can quantify and localize weak locations in the distribution grid. TPIA can
inform system planners where new or updated infrastructure may be needed and
it helps grid operators locate areas in the system where reliability issues
may occur under abnormal or extreme conditions. Section II details relevant
prior work.

In TPIA, we first construct aggregated equivalent circuits to represent three-
phase distribution grids. To solve otherwise infeasible networks, we add
variable slack current sources at each node in the constructed circuit. These
sources then compensate for the mismatch in the Kirchhoff’s Current Law (KCL)
equations. Due to the discrepancy between the number of unknown variables and
constraints, the solution of this overall circuit will have infinitely many
solutions. Because of this mismatch, we minimize the norm of the slack sources
in TPIA through an optimization approach.

An illustration of the slack current source injection is shown for the IEEE 4
node testcase in Fig. 1. The slack current sources, which we call
“infeasibility currents” are the red independent current sources labeled
$i_{f}$. Within the optimization, these new infeasibility currents are
minimized and constrained by KCL at each node. In this paper, we formulate
TPIA with two different objective functions. We first use the common least
squares objective and then we develop TPIA using the least absolute value
objective (L1 norm). We use the L1 norm to induce greater sparsity in results
in comparison to the least squares approach and to localize sources of a power
flow simulation failure. Methodology for the two formulations is shown in
Section III. Implementation of TPIA is shown in Section IV.

Section V describes experimental set-ups and results for three-phase networks
up to 8500 nodes and discusses the empirical differences between the two TPIA
formulations. We show that both the least squares and L1 norm formulations of
TPIA converge when three-phase power flow simulations in GridLAB-D and SUGAR-D
(from our work in [3]) otherwise fail. We show that the L1 formulation
localizes grid infeasibilities whereas the least squares formulation does not.
Finally, we provide a use case with the L1 TPIA formulation for evaluating
electric vehicle (EVs) penetration and validate our results by constructing a
feasible network based on information from the L1 TPIA results.

##  II Prior Work

TPIA solves an optimization problem to locate weak-spots in a three-phase
unbalanced distribution network especially when power flow solution does not
exist. No standard tool today (e.g. OpenDSS, GridLAB-D) can identify these
weak spots for infeasible large-scale realistic three-phase networks. We
develop the TPIA algorithm building on state-of-the-art research in three-
phase distribution grid analysis. We briefly discuss the work here.

###  II-A Steady-state Analysis of Distribution Grid

Distribution grid analysis can be fundamentally categorized as either a
simulation (or power flow) or an optimization problem.

####  II-A1 Three-Phase Power Flow

Three-phase power flow has been formulated using a number of techniques like
the Backward-Forward Sweep Method (BFS) [9], Current Injection Method (CIM)
[10], and Equivalent Circuit Formulation (ECF) [3]. BFS works well for weakly-
meshed or radial systems but is known to suffer otherwise [9]. CIM improves
solution robustness by formulating the problem as current-balance rather than
power-balance [10]. However, CIM is known to suffer when there is a high
penetration of generation models (PV nodes), which is increasingly common in
distribution grids. ECF uses circuit theory to model the electric grid as a
three-phase circuit and uses circuit-theoretic heuristics to obtain the
solution to the three-phase power flow problem. We build the TPIA framework
with ECF as the basis.

In ECF, the AC network constraints are written using KCL, which requires that
the sum of currents at each node add up to zero. Mathematically, the
formulation is shown in (1)-(4), where $\odot$ and $\oslash$ represent
element-wise multiplication (Hadamard Product) and element-wise division
(Hadamard Division) respectively. Variables $V_{R}$ and $V_{I}$ are vectors of
$n$ unknown variables in rectangular components; $\Omega$ represents the
phases $(\Omega={A,B,C})$; $G$ and $B$ are respectively the admittance and
susceptance matrices of size $n\times n$; and $I^{NL}_{R}$ and $I^{NL}_{I}$
are $n\times 1$ vectors representing the nonlinear real and imaginary load or
generator currents. $I^{NL}_{R}$ and $I^{NL}_{I}$ are functions of unknown
variables and are shown in (3) and (4) where $P^{\Omega}$ represents real
power and $Q^{\Omega}$ represents reactive power. For a PV bus, $Q^{\Omega}$
is an unknown variable. At the PQ bus, $Q^{\Omega}$ is known so the current
equations can be simplified. While the ECF approach can handle PV buses, for
brevity, (1)-(4) do not include additional PV model terms and constraints.

$\displaystyle GV^{\Omega}_{R}-BV^{\Omega}_{I}+I^{NL}_{R}(V^{\Omega}_{R},$
$\displaystyle V^{\Omega}_{I},Q^{\Omega})=0$ (1) $\displaystyle
BV^{\Omega}_{R}+GV^{\Omega}_{I}+I^{NL}_{I}(V^{\Omega}_{R},$ $\displaystyle
V^{\Omega}_{I},Q^{\Omega})=0$ (2) $\displaystyle
I^{NL}_{R}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})=P^{\Omega}\odot$
$\displaystyle V_{R}^{\Omega}+Q^{\Omega}\odot V_{I}^{\Omega}$
$\displaystyle\oslash((V_{R}^{\Omega})^{2}+(V_{I}^{\Omega})^{2})$ (3)
$\displaystyle
I^{NL}_{I}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})=P^{\Omega}\odot$
$\displaystyle V_{I}^{\Omega}-Q^{\Omega}\odot V_{R}^{\Omega}$
$\displaystyle\oslash((V_{R}^{\Omega})^{2}+(V_{I}^{\Omega})^{2})$ (4)

####  II-A2 Optimizations in Three-Phase Distribution Grids

While many earlier works on distribution grid optimization relied on
linearized DC constraints that poorly represent the system [8], there has been
significant advances in AC-network constrained optimization methods for
distribution grids. For instance, [7] maps the AC network constraints using
power balance equations and [11] uses CIM to map the AC network constraints
with KCL. Both of these approaches can work with general objectives. However,
because these existing distribution grid optimization formulations do not
demonstrate robustness and scalability to solve large-scale optimization
problems ($>$5000 nodes), we extend the ECF framework to include AC-network
constrained optimization for use with TPIA.

###  II-B Optimization Techniques for Identifying Infeasible Power Grids

There are some approaches that can identify infeasible positive-sequence
transmission networks [4], [5], [6] as well as distribution networks [7]. For
transmission positive-sequence networks, [4], [5], [6] have all developed
techniques to detect grid infeasibilities. [4]’s approach is based on power-
mismatch equations. This approach is known to struggle with three-phase
distribution networks due to the use of power balance equations [3]. Another
distinct approach for identifying grid infeasibility was introduced for
transmission positive sequence networks using a circuit-theoretic ECF
foundation [5]. However, both [4] and [5] cannot localize the sources of
infeasibility in the grid. Our follow-on work in [6] localized grid
infeasibilities in transmission networks by formulating the L2 norm objective
function with an L1 norm regularization term. This formulation encourages some
sparsity in the optimization solution, which can be used to find regions of
grid infeasibility. However, none of these works apply directly to three-phase
distribution networks due to its unbalanced and complex nature.

Research in identifying infeasible three-phase distribution networks is sparse
but exists. For instance, [7] formulates a three-phase optimization problem
using power balance equations to identify infeasible distribution grids. It
adds and minimizes slack variables to accommodate constraint violations and
identify infeasible networks. However, many challenges persist with this
approach. The slack variables do not have a physical meaning, so their values
provide no grid insight outside of a constraint violation; the method has not
been shown to scale; and the approach is unable to localize sources of
infeasibilities. With TPIA, we will address all these challenges and develop
an infeasibility analysis framework that provides practical, interpretable
causes of distribution grid infeasibility.

##  III Three-Phase Infeasibility Analysis

The goal of TPIA is to identify and localize weak sections of three-phase
networks when the three-phase power flow (TPF) problem has no solution. Thus
we begin by constructing a three-phase equivalent circuit of the evaluated
distribution network (compare to [3]) and transforming the TPF simulation
problem into a TPIA optimization problem.

To illustrate further, consider the 4-bus network in Fig. 1. Without the
presence of current sources in red, the network can be solved directly using
the ECF formulation in [3] to obtain a TPF solution. However, if there is no
physical solution for the problem, the traditional TPF technique would diverge
even with the best heuristics. Instead, in TPIA, we add variable current
sources at each phase of each node, namely the infeasibility sources
$i_{f,k}^{\Omega}$ (shown in red for node $k$), to ensure solvability of the
network. The goal of the variable infeasibility sources is to supply any
missing current where the KCL equations cannot be satisfied within the
confines of the original network. While the infeasibility sources ensure
feasibility of the circuit (and correspondingly the TPF problem) for all
possible operating conditions, infinitely many solutions exist as the system
is underdetermined with these additional variable $i_{f}^{\Omega}$ sources. To
ensure that TPIA returns either the TPF solution (when the problem is
feasible) or the least number of infeasibility currents (when the problem is
infeasible), we minimize the norm of $i_{f}^{\Omega}$ subject to TPF AC
network constraints.

###  III-A General Formulation

Mathematically, TPIA takes the form in (III-A), where $i_{f}^{\Omega}$
represents a vector of infeasibility currents; $\mathcal{G}(i_{f}^{\Omega})$
is a generic objective function; and $\mathcal{H}(X)$ represents three phase
AC network constraints from CIM or ECF where $X$ is the vector of unknown
variables in TPF. Unlike positive-sequence power flow, $\mathcal{H}(X)$ models
each phase of $\Omega$ in the distribution network and can capture both
balanced and unbalanced networks.

—s— X, i_f^ΩG(i_f^Ω) H(X) - i_f^Ω = 0

In this paper, we derive and solve two distinct forms of the objective
functions $\mathcal{G}(i_{f}^{\Omega})$ in (III-A). The first minimizes the
infeasibility currents using least squares and the second minimizes using the
L1 norm.

###  III-B Location of Infeasibility Sources

Infeasibility current sources can be added at any defined set of nodes. The
ability to optimize over a predefined set of nodes is beneficial because it
allows a system planner to identify which locations in the system they are
capable of modifying before conducting the optimization. The L1 norm approach
inherently returns sparse solutions which effectively localize weak spots in
the distribution grid whereas the least squares approach spreads out
infeasibility currents. Therefore, defining a subset of notes is especially
beneficial in the least squares approach. Depending on available resources,
system planners may prefer either to make larger adjustments at a few nodes
based on the localized results from L1 norm approach or to make a greater
number of smaller adjustments via the least squares solution. While the
inherent sparsity of the L1 norm is an advantage, it is non-differentiable.
The standard least squares approach is empirically much easier and faster to
solve in comparison. We explore both formulations in the following subsections
because each has distinct advantages.

###  III-C Least Squares Formulation

As a first approach to TPIA, we use least squares to represent the objective
function in (III-A). The least squares objective function is well studied
within optimization and adds no inequality constraints to the Lagrangian
formulation as it has a continuous derivative. Reformulating (III-A) with a
least-squares objective function and using the TPF equations in (1)-(4), we
obtain the formulation in (5a) - (5c).

$\displaystyle min$
$\displaystyle\IEEEeqnarraymulticol{9}{l}{\frac{1}{2}||i_{f,R}^{\Omega}||^{2}+\frac{1}{2}||i_{f,I}^{\Omega}||^{2}}$
(5a) $\displaystyle s.t.$ $\displaystyle
GV^{\Omega}_{R}-BV^{\Omega}_{I}+I^{NL}_{R}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})-i_{f,R}^{\Omega}=0$
(5b) $\displaystyle
BV^{\Omega}_{R}+GV^{\Omega}_{I}+I^{NL}_{I}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})-i_{f,I}^{\Omega}=0$
(5c)

To find an optimal solution to this problem, we formulate the Lagrangian of
the problem in (6) where $\lambda_{R}^{\Omega}$ and $\lambda_{I}^{\Omega}$ are
the vector of dual variables or KKT multipliers. In least squares based TPIA,
we find that these dual variables are closely tied to the infeasibility
currents, such that non-zero magnitudes for dual variables
$(|\lambda_{k}^{\Omega}|>0)$ at node $k$ exist only when the corresponding
node $k$ has non-zero infeasibility currents, $|i_{f,k}^{\Omega}|>0$.

$\displaystyle\mathcal{L}(i_{f,R}^{\Omega},$ $\displaystyle
i_{f,I}^{\Omega},V_{R}^{\Omega},V_{I}^{\Omega},Q^{\Omega},\lambda_{R}^{\Omega},\lambda_{I}^{\Omega})=$
(6)
$\displaystyle\frac{1}{2}||i_{f,R}^{\Omega}||^{2}+\frac{1}{2}||i_{f,I}^{\Omega}||^{2}$
$\displaystyle+$
$\displaystyle\lambda_{R}^{\Omega,T}\left(GV^{\Omega}_{R}-BV^{\Omega}_{I}+I^{NL}_{R}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})-i_{f,R}^{\Omega}\right)$
$\displaystyle+$
$\displaystyle\lambda_{I}^{\Omega,T}\left(BV^{\Omega}_{R}+GV^{\Omega}_{I}+I^{NL}_{I}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})-i_{f,I}^{\Omega}\right)$

With the Lagrangian formulated, we derive the first-order optimality
conditions for the KKT conditions. The KKT conditions for the primal variables
are shown in Equations (7a) - (7e) and are solved in our optimization routine
using Newton’s Method.

$\displaystyle\frac{\partial\mathcal{L}}{\partial i_{f,R}^{\Omega}}=$
$\displaystyle i_{f,R}^{\Omega}-\lambda_{R}^{\Omega}=0$ (7a)
$\displaystyle\frac{\partial\mathcal{L}}{\partial i_{f,I}^{\Omega}}=$
$\displaystyle i_{f,I}^{\Omega}-\lambda_{I}^{\Omega}=0$ (7b)
$\displaystyle\frac{\partial\mathcal{L}}{\partial Q^{\Omega}}=$
$\displaystyle\frac{\partial I^{NL}_{R}}{\partial
Q^{\Omega}}\lambda_{R}^{\Omega}+\frac{\partial I^{NL}_{I}}{\partial
Q^{\Omega}}\lambda_{I}^{\Omega}=0$ (7c)
$\displaystyle\frac{\partial\mathcal{L}}{\partial V^{\Omega}_{R}}=$
$\displaystyle(G^{T}+\frac{\partial I^{NL}_{R}}{\partial
V^{\Omega}_{R}})\lambda_{R}^{\Omega}+(B^{T}+\frac{\partial
I^{NL}_{I}}{\partial V^{\Omega}_{R}})\lambda_{I}^{\Omega}=0$ (7d)
$\displaystyle\frac{\partial\mathcal{L}}{\partial V^{\Omega}_{I}}=$
$\displaystyle(-B^{T}+\frac{\partial I^{NL}_{R}}{\partial
V^{\Omega}_{I}})\lambda_{R}^{\Omega}+(G^{T}+\frac{\partial
I^{NL}_{I}}{\partial V^{\Omega}_{I}})\lambda_{I}^{\Omega}=0$ (7e)

In (7a) and (7b), we can see that per KKT conditions,
$i_{f,R}^{\Omega}={\lambda}_{R}^{\Omega}$ and
$i_{f,I}^{\Omega}=\lambda_{I}^{\Omega}$, which validates our earlier
observation about the tight coupling between dual variables and infeasibility
currents. Equations (7c) - (7e) have nonlinear components and will be
linearized using a Taylor Series approximation which allows us to use Newton’s
Method to solve this optimization. The primal feasibility constraints are the
TPF equations shown in Equations (5b) - (5c).

The disadvantage of the least squares approach is that it does not return a
sparse solution vector for infeasibility currents $i_{f}^{\Omega}$. While a
system planner would be able to identify an infeasible network with the least
squares formulation of TPIA, they would not be able to localize the regions of
the network that caused the infeasibility because the infeasibility currents
would be spread throughout the system. We show this key characteristic in Fig.
3A which demonstrates that infeasibility currents are present at nearly every
node in the system during an infeasible operating scenario. Planners are
unable to use this information to perform corrective actions. To gain any
realistic insight on physical remedies, the planner would need to define a
smaller subset of nodes where infeasibility current sources should be added
rather than across the whole system. Therefore, to get a more sparse solution
vector such that it provides more immediately localized information to the
planner, we consider a second approach to TPIA using an L1 norm objective
function.

###  III-D L1 Norm Formulation

We next consider the L1 norm based objective, which can localize the nonzero
injection currents to the weakest sections of the grid because of its inherent
property to provide sparse solutions.

With this mechanism, when a network file incorrectly represents the real
physical system and TPF fails, system planners can isolate where the modeling
error is based on looking at the location of the nonzero injection currents.
This will prevent significant manual overhead that would otherwise be required
to find the error in the system topology (e.g., typo in the input file). In
situations where TPF fails due to inherent lack of resources, the L1
formulation of TPIA can be used as a planning tool to ascertain where new
infrastructure may be needed in the grid to accommodate growing or changing
electricity demand.

Mathematically, the L1 norm formulation of TPIA is formulated as:

$\displaystyle min$
$\displaystyle\IEEEeqnarraymulticol{9}{l}{|i_{f,R}^{\Omega}|+|i_{f,I}^{\Omega}|}$
(8a) $\displaystyle s.t.$ $\displaystyle
GV^{\Omega}_{R}-BV^{\Omega}_{I}+I^{NL}_{R}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})-i_{f,R}^{\Omega}=0$
(8b) $\displaystyle
BV^{\Omega}_{R}+GV^{\Omega}_{I}+I^{NL}_{I}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})-i_{f,I}^{\Omega}=0$
(8c)

Unlike the least squares formulation in (5a) - (5c), the L1 norm formulation
contains non-differentiable absolute value terms. Due to these non-
differentiable terms, we cannot directly apply Newton’s Method to find
stationary points from the optimization problem.

To solve this problem, we can model infeasibility sources as two parallel
current sources where one is flowing into the node and the other is flowing
out of the node, as shown in Fig. 2B, instead of the injection used in least
squares shown in Fig. 2A.

(a) A

(b) B

Figure 2: [A] shows how infeasibility sources are modeled in the least squares
formulation’s circuit. [B] shows how infeasibility sources are modeled in the
L1 norm formulation. In [B] the infeasibility source is represented as both an
incoming and outgoing current source.

Mathematically we can represent these parallel sources by setting
$i_{f,R}^{\Omega}=i_{f,R+}^{\Omega}-i_{f,R-}^{\Omega}$ and
$i_{f,I}^{\Omega}=i_{f,I+}^{\Omega}-i_{f,I-}^{\Omega}$ and then constraining
$i_{f,R+}^{\Omega},i_{f,R-}^{\Omega},i_{f,I+}^{\Omega}$ and
$i_{f,I-}^{\Omega}$ to be positive values as is done in [12]. Our new L1
objective function then becomes
$|i_{f,R+}^{\Omega}|+|i_{f,R-}^{\Omega}|+|i_{f,I+}^{\Omega}|+|i_{f,I-}^{\Omega}|$
which can be simplified to
$i_{f,R+}^{\Omega}+i_{f,R-}^{\Omega}+i_{f,I+}^{\Omega}+i_{f,I-}^{\Omega}$
since we constrained each of these variables to be positive. Now each variable
is differentiable and can be used within Newton’s Method. This new formulation
in shown in Equations (9a) - (9d).

$\displaystyle\min\quad$ $\displaystyle
i_{f,R+}^{\Omega}+i_{f,R-}^{\Omega}+i_{f,I+}^{\Omega}+i_{f,I-}^{\Omega}$ (9a)
s.t. $\displaystyle
GV^{\Omega}_{R}-BV^{\Omega}_{I}+I^{NL}_{R}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad-(i_{f,R+}^{\Omega}-i_{f,R-}^{\Omega})=0$
(9b) $\displaystyle
BV^{\Omega}_{R}+GV^{\Omega}_{I}+I^{NL}_{I}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad-(i_{f,I+}^{\Omega}-i_{f,I-}^{\Omega})=0$
(9c) $\displaystyle
i_{f,R+}^{\Omega},i_{f,R-}^{\Omega},i_{f,I+}^{\Omega},i_{f,I-}^{\Omega}\geq 0$
(9d)

We then formulate the Lagrangian of the L1 norm formulation in (10). Since we
now have inequality constraints, we add new variables
$\mu_{R+}^{\Omega},\mu_{R-}^{\Omega},\mu_{I+}^{\Omega},$ and
$\mu_{I-}^{\Omega}$ as the inequality dual variable vectors.

$\displaystyle\mathcal{L}(i_{f,R+}^{\Omega},$ $\displaystyle
i_{f,R-}^{\Omega},i_{f,I+}^{\Omega},i_{f,I-}^{\Omega},V_{R}^{\Omega},V_{I}^{\Omega},Q^{\Omega},\lambda_{R}^{\Omega},\lambda_{I}^{\Omega},$
$\displaystyle\mu_{R+}^{\Omega},$
$\displaystyle\mu_{R-}^{\Omega},\mu_{I+}^{\Omega},\mu_{I-}^{\Omega})=$
$\displaystyle
i_{f,R+}^{\Omega}+i_{f,R-}^{\Omega}+i_{f,I+}^{\Omega}+i_{f,I-}^{\Omega}$
$\displaystyle+$
$\displaystyle(\lambda_{R}^{\Omega})^{T}\left(GV^{\Omega}_{R}-BV^{\Omega}_{I}+I^{NL}_{R}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})\right)$
$\displaystyle+$
$\displaystyle(\lambda_{I}^{\Omega})^{T}\left(BV^{\Omega}_{R}+GV^{\Omega}_{I}+I^{NL}_{I}(V^{\Omega}_{R},V^{\Omega}_{I},Q^{\Omega})\right)$
$\displaystyle-$
$\displaystyle(\lambda_{R}^{\Omega})^{T}(i_{f,R+}^{\Omega}-i_{f,R-}^{\Omega})-(\lambda_{I}^{\Omega})^{T}(i_{f,I+}^{\Omega}-i_{f,I-}^{\Omega})$
$\displaystyle+$
$\displaystyle(\mu_{R+}^{\Omega})^{T}i_{f,R+}^{\Omega}+(\mu_{R-}^{\Omega})^{T}i_{f,R-}^{\Omega}$
$\displaystyle+$
$\displaystyle(\mu_{I+}^{\Omega})^{T}i_{f,I+}^{\Omega}+(\mu_{I-}^{\Omega})^{T}i_{f,I-}^{\Omega}$
(10)

Having formulated the Lagrangian, we then derive the first-order optimality
conditions for the KKT conditions which are shown in (11a) - (11d), (12) -
(14). The stationarity constraints for $Q^{\Omega},V^{\Omega}_{R},$ and
$V^{\Omega}_{I}$ remain the same as in (7c) - (7e). The new stationarity
constraints for the $i_{f}^{\Omega}$ variables are shown in (11a) - (11d).

$\displaystyle\frac{\partial\mathcal{L}}{\partial i_{f,R+}^{\Omega}}$
$\displaystyle=$ $\displaystyle 1-\lambda_{R}^{\Omega}+\mu_{R+}^{\Omega}=0$
(11a) $\displaystyle\frac{\partial\mathcal{L}}{\partial i_{f,R-}^{\Omega}}$
$\displaystyle=$ $\displaystyle 1+\lambda_{R}^{\Omega}+\mu_{R-}^{\Omega}=0$
(11b) $\displaystyle\frac{\partial\mathcal{L}}{\partial i_{f,I+}^{\Omega}}$
$\displaystyle=$ $\displaystyle 1-\lambda_{I}^{\Omega}+\mu_{I+}^{\Omega}=0$
(11c) $\displaystyle\frac{\partial\mathcal{L}}{\partial i_{f,I-}^{\Omega}}$
$\displaystyle=$ $\displaystyle 1+\lambda_{I}^{\Omega}+\mu_{I-}^{\Omega}=0$
(11d)

The equality primal feasibility equations in (5b) and (5c) remain the same
with just the appropriate substitution for $i^{f}_{R}$ and $i^{f}_{I}$. We now
also have the inequality primal feasibility equations shown in (12) and the
dual feasibility equations shown in (13).

$\displaystyle
i_{f,R+}^{\Omega},i_{f,R-}^{\Omega},i_{f,I+}^{\Omega},i_{f,I-}^{\Omega}\geq 0$
(12)
$\displaystyle\mu_{R+}^{\Omega},\mu_{R-}^{\Omega},\mu_{I+}^{\Omega},\mu_{I-}^{\Omega}\geq
0$ (13)

Because we introduced inequality constraints, we also add the complimentary
slackness equations shown in (14). The new variable $\epsilon$ in (14)
represents a small valued number that stems from the perturbation of KKT
equations (see primal-dual interior point algorithm in Boyd [13]).

$\displaystyle(\mu_{R+}^{\Omega})^{T}i_{f,R+}^{\Omega}$
$\displaystyle+(\mu_{R-}^{\Omega})^{T}i_{f,R-}^{\Omega}$ $\displaystyle+$
$\displaystyle(\mu_{I+}^{\Omega})^{T}i_{f,I+}^{\Omega}+(\mu_{I-}^{\Omega})^{T}i_{f,I-}^{\Omega}=\epsilon$
(14)

The sparsity achieved with the L1 norm formulation is discussed in Section V
and is shown in Fig. 3B.

##  IV Implementation of TPIA

We developed this framework within our tool called SUGAR-distribution
(SUGAR-D) [14], a version of which was used in the TPF framework in [3]. We
built our own tool rather than using a commercial solver in order to best
exploit the physics of the problem; ensure the real, rather than relaxed,
problem is solved; and to incorporate circuit heuristics for improved
convergence [3]. SUGAR-D uses ECF in conjunction with Newton’s Method to solve
TPF simulations. Algorithm 1 succinctly describes the TPIA methodology.

The TPIA implementation includes but is not limited to the following electric
power grid elements: capacitors; voltage regulators; transformers; switches;
overhead, underground, and triplex transmission lines; inverter-based
distributed generators (IBDGs); fuses; shunts; and loads. ECF element three-
phase models draw on the physical nature of the component, so physical
properties like different combinations of delta and wye connections, unbalance
amongst phases, or the inclusion of the neutral wire can be incorporated
without loss of generality.

With the L1 norm formulation we used a special limiting heuristic called diode
limiting [15] to ensure convergence. Diode limiting is used to mitigate large
step sizes between Newton’s Method iterations amongst inequality dual
variables, ensuring that the dual variables remain within their feasible
space. To accommodate the sensitivity of the L1 norm formulation to initial
conditions, the formulation was also “warm started” for certain testcases with
the TPF solution to the testcase at a less stressed set of conditions. Since
the solution to the L1 norm formulation will be the TPF solution if there are
no infeasibility currents, the warm start enables the L1 norm formulation to
have a reasonable initial starting point.

procedure:

input distribution network models, warm start (TPIA L1 norm only) or initial
conditions, solver settings, Newton’s Method tolerance ($\epsilon)$

add dual variables $(\mu,\lambda)$ and current sources $(i_{f})$

initialize primal and dual variables $(\mu,\lambda)$, infeasibility currents
$(i_{f})$

create sparse matrix structure

loop through linear elements to create the linear admittance matrix $Y_{Lin}$

while iteration count $k$ $<$ max iterations do

loop through nonlinear elements and update values

solve equations for unknown values at $k+1$ using Newton’s Method

limit dual variables (TPIA L1 norm only only)

calculate the error between iterations $k$ and $k+1$

if $error>tolerance$ then

update variables dependent on $k$ to $k+1$ values

else

break the loop and return converged states

end if

end while

Calculate residual values from original nonlinear equations and determine if
any KKT conditions are violated

Algorithm 1 Infeasibility Analysis  TABLE I: Results from running varying
testcases using GridLAB-D, SUGAR-D TPF, and the TPIA formulations.

| GridLAB-D | SUGAR-D | TPIA Least Squares | TPIA L1 norm
---|---|---|---|---
TPF
Testcase | Converged | Converged | Matrix | Iterations | Time (s) | Nonzero $i_{f}$ | Matrix | Iterations | Time (s) | Nonzero $i_{f}$
Size | Size
R4-12.47-1 | Yes | Yes | 42614 | 5 | 1.28 | 0 | 75254 | 5 | 2.08 | 0
R1-12.47-3_OV | No | No | 2244 | 11 | 0.23 | 73 | 4032 | $92^{1}$ | $1.84^{1}$ | 3
R2-25.00-1_OV | No | No | 22612 | 20 | 2.5 | 766 | 39364 | $86^{1}$ | $15.43^{1}$ | 1
R3-12.47-3_OV | No | No | 142594 | 43 | 25.7 | 1504 | 251750 | $439^{1}$ | $443.4^{1}$ | 7
R4-12-47-1_EV | No | No | 46294 | 28 | 7.033 | 324 | 83074 | $20^{1}$ | $6.14^{1}$ | 1
8500node_OV | No | No | 116724 | 235 | 118.43 | 4760 | 219576 | $826^{1}$ | $711.25^{1}$ | 4

111

These values were obtained with a warm start. The reported values represent
the sum of the respective quantity from the L1 approach with the quantities it
took to get the warm start values.

##  V Experiments

###  V-A Experimental Set-up

To validate the efficacy of the TPIA approach, we run experiments on realistic
but modified taxonomy feeders and a modified IEEE 8500node testcase. We
simulate these networks using GridLAB-D, SUGAR-D TPF (described in [3]), and
the least squares and L1 norm TPIA formulations. Table I documents the cases
that were used in these experiments. They are available on Github at
https://github.com/emfoster/TPIA-Experiment-Testcases in the Gridlabd file
format .glm.

The taxonomy feeders R4-12.47-1, R1-12.47-3, R2-25.00-1, R3-12.47-3, and
R4-12.47-1, respectively represent a small urban center, moderate suburban
area, heavy suburban area, and a heavy urban area with a rural spur in four
different climate regions [16]. R4-12.47-1 is used in its original form but
the other networks are modified to make them infeasible by increasing their
loading. We identify these cases in Table I with the addition of “OV” to their
titles. The testcase R4-12.47-1_EV has electric vehicles added to it and is
described in Section V.c. For the purpose of this paper, new infeasibility
current sources are added at all nodes in each testcase with the exception of
the slack bus.

All experiments were ran using an IDE on an i7 processor. Table I documents
the results of our experiments on these networks. The matrix size column
represents the size of the inverted square matrix used in the Newton’s Method
after any rows and columns of zeros are removed. Time, iterations, and nonzero
infeasibility currents are reported for both TPIA formulations. Our code was
not optimized for speed and could be improved in future versions.

(a) A

(b) B

Figure 3: A graph representation of R1-12.47-3_OV with a heatmap of the
greatest infeasibility current magnitude amongst the phases at each node
resulting from the least squares [A] and L1 norm [B] formulations. Currents
are depicted on a logarithmic scale.

###  V-B Comparison of results between GridLAB-D, SUGAR-D TPF, and TPIA

Table I shows the results for each testcase obtained using GridLAB-D, SUGAR-D
TPF, and the least squares and L1 norm TPIA formulations. GridLAB-D and
SUGAR-D TPF only converge for the standard taxonomy feeder R4-12.47-1 which
represents a traditional TPF case. On this testcase, we see that both TPIA
approaches returned zero infeasibility currents which is what we would expect
when TPF simulations converge. We also validated that the solution from the
TPIA approaches matches with that from SUGAR-D TPF.

Neither GridLAB-D nor SUGAR-D TPF converged on any of the modified testcases
that were designed to stress the system. Neither TPF simulations provided any
insight into the physical reason for the failure. However, the two TPIA
formulations converged to a solution with nonzero infeasibility currents for
all modified testcases. The presence of non-zero infeasibility currents
indicate that the problem in its original state is undispatchable.

For reporting purposes and to categorize arbitrarily small infeasibility
current as zero, we create an infeasibility current threshold of
$\epsilon_{If}=1\text{e-}3$ where values less than $\epsilon_{If}$ are set to
zero. Table I gives the number of maximal infeasibility current magnitudes
above $\epsilon_{If}$. There is a stark difference in the number of nonzero
infeasibility currents between the least squares and L1 norm formulations. In
Fig. 3, the magnitude of the infeasibility currents at each phase was
calculated and the maximum value amongst these phase infeasibility currents is
shown at its respective node for R1-12.47-3_OV.

The least squares formulation is less sensitive to initial conditions than the
L1 norm formulation and it did not require any additional heuristics to
converge. As shown in Table I, least squares was faster and required fewer
iterations than the L1 norm approach except for $R4-12.47-1\\_EV$. Speed for
both approaches could be improved through using additional circuits
heuristics. For the L1 norm, a more optimal approach to finding a good warm
start would also improve the time.

The infeasibility currents from TPIA least squares are spread throughout the
network as shown in Fig. 3A, which makes it challenging for grid planners or
operators to identify grid weaknesses. We can see from both Table I and Fig.3A
that 73 of 76 nodes in the system have nonzero infeasibility currents for the
least squares formulation. As noted in the Section III.b, infeasibility
current sources can be added at only a subset of nodes, so for least squares,
a better option may be for system planners with knowledge of the physical
system to identify locations with potential for new or updated infrastructure
and to conduct TPIA least squares over only that subset. In Fig. 3B, TPIA L1
norm had only 3 nodes with infeasibility currents. TPIA L1 norm is able to
return sparse vectors of infeasibility currents, which in turn allows grid
planners to isolate grid weakness and take realistic corrective action to make
the network feasible. In a practical world setting, this implies that if the
grid planner had access to the sparse L1 norm solution, they could add
infrastructure to at most 3 locations in the network to make the case
feasible; whereas, no such insight is available with the least squares
solution since infeasibilities are spread across almost 96% of the locations
in Fig. 3.

Figure 4: A graph representation of the testcase R4-12.47-1_EV with a heatmap
of the greatest infeasibility current magnitude amongst the phases at each
node resulting from L1 norm formulation

###  V-C Use-case scenario on electric vehicle penetration and infrastructure
additions

To further illustrate how this research may be used in a planning scenario, we
did a general study on the impact of adding battery electric vehicles (BEVs)
to an overloaded version of testcase R4-12.47-1. From [17, Tab. DP04], 91.3%
of households in the U.S. have at least one vehicle available with 58.6% of
households having access to two or more vehicles. We conduct our study for the
future scenario where internal combustion vehicles are phased out and we
assume similar vehicle availability. We say that 90% of households have at
least one BEV available and that at some point 55% of those vehicles are
charging via a Level 2 charger conservatively at their maximum charge rate. We
model the charging BEVs as simple constant PQ loads and add them at random
triplex nodes in the system. We ran this testcase in GridLAB-D and SUGAR-D TPF
and both failed to converge. We next used the L1 norm TPIA on the BEV network
(R1-12.47-1_EV) to localize areas of infeasibility in the system and found
that there was one infeasible node in the network that was power deficient.
The scenario is shown in Fig. 4 where we highlight the area of grid weakness.

To demonstrate, how a grid planner could use this information to take
corrective actions, we took the infeasibility currents and voltages for each
phase for the node and then calculated the amount of missing power at each
node. Then, to make the case feasible or dispatchable, we added a battery at
the node rated at the calculated missing power. After adding the battery to
the system, we were able to validate that the network is feasible by running
the modified network using SUGAR-D TPF. We additionally ran both TPIA
formulations and found that both return zero infeasibility currents with the
added battery. The experiment also served as the validation of our choice of
infeasibility current threshold $\epsilon_{If}$. This is the value of
infeasibility currents in the solution vector below which they are treated as
inconsequential and as numerical noise.

The BEV use-case is a pilot example of how the TPIA methodology proposed
within this paper can be used in practice. Other more sophisticated
applications may significantly improve their performance with the use of this
proposed methodology in the future.

##  VI Conclusion

We developed a novel TPIA framework that converges for unbalanced networks
that otherwise fail with traditional TPF. For both the least squares and L1
norm approaches, we show scalability of the method up to 8500 nodes. From the
L1 norm formulation of TPIA we also extract critical information on where the
cause of TPF failures might be for a given network. This proposed approach can
be a critical tool for distribution planning studies moving forward. It can be
used to study the transition of today’s distribution grids to ones with a
high-penetration of DERs or inverter based technologies. In future work, we
plan to more rigorously evaluate the impact of new-age technologies on the
distribution grid using this method.

## Acknowledgment

This work was sponsored in part by the National Science Foundation under
contract ECCS-1800812.

## References

  * [1] National Academy of Sciences, Engineering, and Medicine. The Future of Electric Power in the United States. Washington, D.C: the National Academies Press, 2021.
  * [2] J.S. Homer, D. Narang, Y. Tang, M. Coddington, J.D. Taft, M. Ingram, A.C. Orrell, and A. Hoke, Electric Distribution Systems Planning with DERs - High-level Assessment of Tools and Methods, Pacific Northwest National Lab, Richland, WA, March 2020.
  * [3] A. Pandey, M. Jereminov, M. R. Wagner, D. M. Bromberg, G. Hug and L. Pileggi, “Robust Power Flow and Three-Phase Power Flow Analyses,” in IEEE Transactions on Power Systems, vol. 34, no. 1, pp. 616-626, Jan. 2019.
  * [4] T. J. Overbye, “A power flow measure for unsolvable cases,” in IEEE Transactions on Power Systems, vol. 9, no. 3, pp. 1359-1365, Aug. 1994.
  * [5] M. Jereminov, D. M. Bromberg, A. Pandey, M. R. Wagner and L. Pileggi, “Evaluating Feasibility Within Power Flow,” in IEEE Transactions on Smart Grid, vol. 11, no. 4, pp. 3522-3534, July 2020.
  * [6] S. Li, A. Pandey, A. Agarwal, M. Jereminov and L. Pileggi, “A LASSO-Inspired Approach for Localizing Power System Infeasibility,” 2020 IEEE Power & Energy Society General Meeting (PESGM), 2020, pp. 1-5. Systems, vol. 12, no. 1, pp. 383-388, Feb. 1997.
  * [7] A. S. Zamzam, N. D. Sidiropoulos, and E. Dall’Anese, “Beyond Relaxation and Newton-Raphson: Solving AC-OPF for Multi-Phase Systems with Renewables” in IEEE Transactions on Smart Grid, vol. 9, no. 5, pp. 3966-3975, Sept. 2016.
  * [8] S. Claeys, M. Vanin, F. Geth, and G. Deconinck, “Applications of optimization models for electricity distribution networks” in WIREs Energy and Environment, vol. 19, no. 5, May 2021.
  * [9] D. Shirmohammadi, H. W. Wong, A. Semlyen, and G .X .Luo, “A compensation-based power flow method for weakly meshed distribution and transmission networks” in IEEE Transactions on Power Systems, vol. 3, no. 2, pp.753-762, May 1988.
  * [10] P. A. N. Garcia, J. L. R. Pereira, S. Carneiro, V. M. da Costa and N. Martins, “Three-phase power flow calculations using the current injection method,” in IEEE Transactions on Power Systems, vol. 15, no. 2, pp. 508-514, May 2000.
  * [11] A. O’Connell and A. Keane, “Multi-Period Three-Phase Unbalanced Optimal Power Flow” 2014 5th IEEE PES Innovative Smart Grid Technologies Europe (ISGT Europe), Oct. 12-15.
  * [12] H. Singh, F. L. Alvarado, and E. Liu, “Constrained LAV State Estimation Using Penalty Functions,” in IEEE Transactions on Power Systems, vol. 12, no. 1, pp. 383-388, Feb. 1997.
  * [13] S. Boyd, S. P. Boyd, and L. Vandenberge, Convex Optimization, Cambridge University Press, 2004.
  * [14] A. Pandey, N. Turner-Bandele, and E. Foster, SUGAR-D: A Distribution Systems Analysis and Optimization Tool v2.02, CMU ECE Pileggi Lab Group, Pittsburgh, PA, 2021.
  * [15] M. Jerminov, A. Pandey, and L. Pileggi, “Equivalent Circuit Formulation for Solving AC Optimal Power Flow” IEEE Transactions on Power Systems, vol. 34, no3, pp. 2354-2365.
  * [16] K.P. Schneider, Y. Chen, D. P. Chassin, R. G. Pratt, D. W. Engel, and S. E. Thompson, Modern Grid Initiative Distribution Taxonomy Final Report, Pacific Northwest National Lab, Richland, WA, Nov. 2008.
  * [17] American Community Survey, U.S. Census Bureau, Washington D.C., 2019.


:
85 2018 Machine Learning for Healthcare 1

# A Domain Guided CNN Architecture for Predicting Age from Structural Brain
Images

Pascal Sturmfels<EMAIL_ADDRESS>Department of Electrical Engineering
University of Michigan and Computer Science
Ann Arbor  MI  Saige Rutherford<EMAIL_ADDRESS>Department of Psychiatry
University of Michigan
Ann Arbor  MI  Mike Angstadt<EMAIL_ADDRESS>Department of Psychiatry
University of Michigan
Ann Arbor  MI  Mark Peterson<EMAIL_ADDRESS>Department of Physical
Medicine
University of Michigan and Rehabilitation
Ann Arbor  MI  Chandra Sripada<EMAIL_ADDRESS>Department of Psychiatry
University of Michigan
Ann Arbor  MI  Jenna Wiens<EMAIL_ADDRESS>Department of Electrical
Engineering
University of Michigan and Computer Science
Ann Arbor  MI

###### Abstract

Given the wide success of convolutional neural networks (CNNs) applied to
natural images, researchers have begun to apply them to neuroimaging data. To
date, however, exploration of novel CNN architectures tailored to neuroimaging
data has been limited. Several recent works fail to leverage the 3D structure
of the brain, instead treating the brain as a set of independent 2D slices.
Approaches that do utilize 3D convolutions rely on architectures developed for
object recognition tasks in natural 2D images. Such architectures make
assumptions about the input that may not hold for neuroimages. For example,
existing architectures assume that patterns in the brain exhibit translation
invariance. However, the meaning of a pattern may differ depending on its
location in the brain. There is a need to explore novel architectures that are
tailored to neuroimaging data. We present two simple modifications to existing
CNN architectures to better learn from structural neuroimaging data. Applied
to the task of brain age prediction, our network achieves a mean absolute
error (MAE) of 1.4 years and trains 30% faster than a CNN baseline that
achieves a MAE of 1.6 years. Our results suggest that lessons learned from
developing models on natural images may not directly transfer to neuroimaging
tasks. Instead, there remains a large space of unexplored questions regarding
model development in this area, whose answers may differ from conventional
wisdom.

##  1 Introduction

Convolutional neural networks (CNNs) have achieved state-of-the-art
performance across a variety of tasks involving natural images, including
object and action recognition. Convolutional layers come with a host of
benefits, including a more efficient usage of parameters compared to standard
neural networks and the ability to learn patterns at multiple scales. In
addition, convolutional layers allow the network to learn spatially invariant
relationships between input and output. This is particularly useful in vision
tasks involving natural images, in which an object’s identity is often
independent from its location in the image.

Researchers from the neuroimaging community have recently begun exploring the
utility of CNNs applied to publicly available brain-image datasets, predicting
neurological conditions like Alzheimer’s (Sarraf et al., 2016) and autism
(Anirudh and Thiagarajan, 2017), as well as age (Cole et al., 2017) and
survival time of high-grade glioma patients (Nie et al., 2016). Although
structural brain scans and natural images share similarities, there are
several key differences that make adapting CNNs to brain images non-trivial.

Much of the existing work on CNNs has focused on 2-dimensional (2D) images,
but structural MRI scans of the human brain are 3-dimensional (3D). Some
previous work has treated 3D volumes as a set of independent 2D slices (Sarraf
et al., 2016; Farooq et al., 2017), but doing so fails to fully leverage the
spatial structure of the brain. A recent review (Bernal et al., 2017) suggests
that 3D models outperform 2D models on several neuroimaging tasks. This
increase in performance, however, comes at the cost of increased computation.
The computational complexity involved in applying CNNs to 3D volumes of the
brain makes it difficult to explore architectures even on small datasets. In
Cole et al. (2017), a 10-layer 3D CNN took 83 hours to train on only 2000
images.

Existing applications of 3D convolutions to neuroimaging data (Cole et al.,
2017; Cao et al., 2016) use architectures based on existing 2D CNN
architectures. However, choices like i) the size and number of convolutional
filters, ii) the type and size of pooling, and iii) the relationship between
pooling, convolutional and fully connected layers, are among a myraid choices
that, while extensively studied for natural images, have not been well studied
in the context of 3D structural neuroimages.

Architecture design choices are often inspired by known invariances in the
data. The types of invariances in natural images are typically affine, as
evidenced by the data augmentation techniques, such as translations and
reflections, used in existing applications of CNNs (Krizhevsky et al., 2012).
However, it is unclear which of these invariances arise in neuroimaging data.
Given that structural brain images are aligned, translation invariance may not
be as important in brain images as it is in natural images. In fact, if local
patterns in the brain are location-dependent, then standard convolutions may
be inappropriate.

In this work, we begin to explore CNN architectural choices in the context of
structural neuroimaging data. We propose two modifications to existing CNN
architectures. Our modifications produce a network that is able to better
learn patterns from brain images and trains faster than architectures
developed on natural image tasks. Furthermore, these modifications are
straightforward to implement in Tensorflow (Abadi et al., 2016). By sharing
our implementation111https://gitlab.eecs.umich.edu/mld3/brain_age_prediction,
we hope to inspire other researchers to continue to challenge current
assumptions. We summarize both the technical and clinical significance of our
work below.

#### Technical Significance

Based on the structure of neuroimaging data, we present two modifications to
existing CNN architectures: 1) learning different parameters in different
regions of the brain, and 2) applying more filters in the early stages of the
network. These modifications enable learning distinct patterns in different
brain regions. By designing these modifications to target 3D brain images, we
demonstrate improved performance and training times compared to a 3D CNN
baseline when predicting age from brain scans. These improvements are robust
to changes to the amount of training data and the number of network
parameters. Our work suggests a greater space of improvements can be achieved
by redesigning CNNs for brain images.

#### Clinical Relevance

Neuroimages are a rich form of data from which a wealth of clinically relevant
labels can be predicted. We focus on the task of predicting age from
neuroimages. Age prediction from neuroimaging data allows researchers to
quantify the difference in predicted age and chronological age, which can be a
powerful marker of deviation from expected aging trajectories. There are a
number of psychiatric and neurological conditions that are closely linked to
such age-related deviations, including attention-deficit/hyperactivity
disorder (Shaw et al., 2007) and Alzheimer’s (Stern, 2012). Differences in
predicted and chronological age have already been shown to be reliable and
heritable biomarkers for a variety of neurodegenerative conditions (Cole et
al., 2017). In addition, age prediction from brain scans enables further
investigation into mechanistic factors correlated with accelerated or
decelerated brain aging and associated changes in cognitive function. Our work
takes the first step toward designing deep learning models that capture the
invariances present in structural brain scans. We demonstrate that our ideas
improve age prediction from these scans, but we hypothesize that the proposed
ideas could generalize to other neuroimaging tasks.

####

The rest of this paper is organized as follows. Section 2 presents related
work, both in developing CNNs on neuroimaging data and in generalizing CNNs
for data with different types of invariances. Section 3 discusses the proposed
architectural modifications in more detail. Section 4 details the dataset we
use for evaluation, and our experimental setup. In Section 5, we present an
empirical evaluation of our proposed approach as well as several follow-up
analyses.

##  2 Related Work

There is a broad scope of work in applying CNNs to brain images. The most
common tasks include anomaly detection, which covers tumor and micro-bleed
detection, segmentation, which includes skull stripping and tumor
segmentation, and label prediction. Label prediction is involves disease
prediction, as well as the task we focus on: age prediction.

Predicting age from brain images is an important step towards understanding
typical brain development and further understanding developmental disorders
(Dosenbach et al., 2010). While there has been a lot of progress on the task
of developing machine learning techniques for predicting age from brain images
(Franke et al., 2012), here we focus primarily on those works that utilize
CNNs. Cole et al. (2017) propose a CNN architecture with repeated blocks of
3D, $3\times 3\times 3$ convolutions. Their focus is not on comparing
different CNN architectures, but rather comparing the performance of their
proposed architecture to a Gaussian Process model.

Few works exist that explicitly attempt to adapt CNNs to brain image tasks
like age prediction or disease prediction. Meszlényi et al. (2017) and
Kawahara et al. (2017) propose novel architectures to predict mild cognitive
impairment and age, respectively, from functional connectivity networks
extracted from functional MRI. They model the network as a graph and apply
convolutions across the edges. Neither consider applying CNNs directly to
brain scans. Zheng et al. (2017) introduce a method involving a new pooling
operation to predict HIV and ADHD from directly from functional imaging.
However, their method requires an ensemble of CNNs and a separate algorithm to
mix the ensemble, which increases computational complexity during training.
Many of the works that focus on learning from brain data do not discuss
architectural modifications. As described above, those that do either focus on
learning from functional connectivity networks rather than actual images, or
require computationally expensive models.

More generally, there is a large body of work that aims to tailor CNNs to data
with different structural properties. Dieleman et al. (2016) introduce four
new operations to help CNNs learn invariance to rotations. Gens and Domingos
(2014) and Cohen and Welling (2016) attempt to generalize CNNs to arbitrary
groups of affine transformations. Ngiam et al. (2010) modify convolutions in a
way that helps CNNs learn invariances present in the training data. These
works focus on natural images. None consider 3D spatial data. In contrast, we
investigate a different and relatively unexplored class of images. We propose
a novel architecture based on the structure of 3D brain volumes. Our approach
increases performance without sacrificing computational efficiency, and is
straightforward to implement.

##  3 Methods

This section describes two proposed modifications to existing architecture
choices. The first, regional segmentation, is motivated by the idea that a CNN
should learn different parameters in different regions of the brain. We
segment the brain into consecutive regions and treat each region as an input
channel. This encourages the network to learn region-dependent patterns.
However, region-specific information is lost after the initial convolutional
layers. To address this, the second modification applies more filters earlier
in the network and fewer filters later. This is in contrast to the existing
practice of learning fewer filters early on in the network and more in the
later layers. Combined, these two modifications both improve performance on
age prediction and decrease training time.

###  3.1 Regional Segmentation

CNNs are designed to capture spatially invariant patterns (LeCun et al.,
1998). This model works well for natural images, since most objects in images
retain their identity regardless of where they appear in an image. In
contrast, brain images are typically aligned: each brain in a given dataset
occupies the same spatial region. Complete spatial invariance is thus
unnecessary since the network is not required to deduce the location of
objects of interest. Furthermore, across different regions of the brain, the
same pattern may have different meaning. Requiring a network to learn the same
parameters over the entirety of the brain image may force those parameters to
be too general and lose region-specific information.

To better capture region-specific information, a CNN should learn different
parameters in different brain regions. To address this, we propose the
following modification: before applying convolutions to brain images, divide
the brain image into distinct regions and concatenate those regions along the
channel dimension. Learned weights will not be shared across regions, allowing
the network to learn region-specific patterns.

Given that we aim to learn different patterns in different regions, how should
these regions be chosen? It may seem appealing to segment the brain into known
anatomical regions. However, following known anatomical regions is
challenging. Anatomical regions in the brain may vary in size and shape, but
it is computationally desirable to have regions with equal dimensions.
Formatting these regions to share the same dimensions may lead to large
quantities of zero padding and wasted computation.

Figure 1: A depiction of regional segmentation. The operation is simple: it
divides an input volume into separate, consecutive regions and then treats
each region as a channel. Different colors in the figure indicate different
parameters. The convolutional filter does not share parameters across regions,
but applies the separate parameters simultaneously over separate regions.

In light of the difficulty of using expert-defined anatomical regions, we
segment the brain into consecutive, equally-sized cubes. More formally: let
$I$ be a 3D, single-channel volume with dimensions $(X,Y,Z)$. Let $k$ be an
integer regional segmentation rate, a hyperparameter. Regional segmentation
divides $I$ into $k^{3}$ adjacent regions of size
$(\lfloor\frac{X}{k}\rfloor,\lfloor\frac{Y}{k}\rfloor,\lfloor\frac{Z}{k}\rfloor)$.
The regions are then concatenated along the channel dimension. In practice,
each region is given an additional 3-voxel boundary in each dimension to avoid
losing information along the edges of the regions. This approach, depicted in
Figure 1 is easy to implement, requires no prior knowledge, and allows for
different levels of granularity. Setting $k$ to be large corresponds to a
finer segmentation of the brain, with less spatial parameter sharing;
conversely, setting $k$ to be small corresponds to more spatial parameter
sharing. Existing architectures set $k=1$, while fully connected networks set
$k$ to be the input dimension, treating each voxel as independent of the
others.

Separate parameters are learned for each region since regions are no longer
spatially connected. Given that these regions are now treated as separate
input channels, depthwise convolutions may be more appropriate than standard
convolutions, because they peform completely separate convolutions across
channels (Chollet, 2016). However, current deep learning libraries lack
implementations of 3D depthwise convolutions, and non-native implementations
are computationally expensive. Normal convolutions do learn different
parameters over different channels, but combine information across channels
after convolution. We mitigate this problem in two ways. First, we rotate
regions into the orientation that minimizes spatial overlap across the channel
dimension as depicted in Figure 2, ensuring that as much region-specific
information as possible is retained during convolutions. Second, we alter the
number of filters in the network to focus on the earlier layers, as discussed
in the next section.

Figure 2:  An illustration of the overlap in the channel dimension between
regions of a center-aligned image after regional segmentation with $k=2$. The
original image is depicted on the left. The top pathway depicts the default
rotation of the regions, which minimizes spatial overlap. The bottom pathway
depicts the rotation that maximizes spatial overlap. For illustration
purposes, this image depicts a 2D circle, but the same principle applies for
3D brain scans.

###  3.2 Filter Layouts

Existing CNN architectures have a common trend: the number of filters in
convolutional layers starts small and increases with each layer. This models
human intuition: at finer detail, objects share local patterns like edges and
corners. But with increasing abstraction, objects become more and more
distinct. Having a large number of specialized filters in later layers of a
CNN is thus critical for an object recognition task with a large number of
distinct classes (Zeiler and Fergus, 2014).

The input to our network, however, consists of segmented regions of the brain
concatenated across the channel dimension. This poses a problem: the first
convolutional operation across the input image combines information from all
regions. With each convolution, we lose information distinguishing regions. To
mitigate this issue, we propose reversing the traditional filter scheme:
instead of starting with a small number of filters and increasing them, start
with a large number of filters and decrease that number in later layers. We
hypothesize that it is critical to extract as much information as possible
from earlier stages of the image, and less important to learn information in
the later layers, when information from all regions is blended together.

Our approach forces the network to learn more parameters at the lower levels,
since these levels contain more information from distinct regions. This
approach also bottlenecks the representation of these regions at later layers
in the network, which acts as a form of regularization.

Usually, it is difficult to experiment with a large number of filters early on
in a network, because earlier filters are applied to larger images. Loading
more filters on larger images results in reduced training speed and increased
memory demands. However, by first applying regional segmentation and then
reversing the filter layout, we find that the training time is still decreased
relative to a baseline, as discussed in Section 5.

##  4 Experimental Setup & Baselines

In this section, we describe the brain image dataset and task we used to
evaluate our proposed modifications. We detail our exclusion criteria, the
baseline we compare our proposed method to, and the specifics of how our
models are trained. We also describe the architectural details of both the
baseline and our proposed approach.

###  4.1 Dataset & Task

We evaluate our modifications using data from the Philadelphia
Neurodevelopmental Cohort (PNC) study (Satterthwaite et al., 2014). The PNC
study is a population-based sample of children ages 8-21 ($\mu$=15.48,
$\sigma$=3.19) living in the greater Philadelphia area. Neuroimaging data were
collected from 1,445 subjects. Data from 997 of those subjects were made
publicly available. These data were initially downloaded and preprocessed for
the purposes of functional connectivity analysis. Therefore, data quality
inclusion criteria follow functional connectivity guidelines (low motion, at
least 4min usable data). Data from 724 of the publicly available subjects met
these criteria and were included. For our analysis, we use T1-weighted
structural images (voxel size $0.94\times 0.94\times 1$, FOV dimensions
$196\times 256\times 160$). These images were preprocessed using the DARTEL
toolbox in SPM8. Preprocessing steps included: bias field correction, brain
extraction, and spatially normalizing to MNI152 template space. After
preprocessing, the structural images have dimensions $121\times 145\times
121$.

From these images, we aim to predict subject age. To learn a model mapping 3D
images to age, we randomly split the subjects into training (n=524),
validation (n=100) and test (n=100) sets. Validation data were used for model
selection and stopping criteria. All reported performance metrics pertain to
the test set unless otherwise specified.

###  4.2 Models

We compare a standard baseline architecture from Cole et al. (2017) to an
architecture with our proposed modifications. Both the baseline and our
proposed network are made up of repeated blocks of (convolution, convolution,
batch normalization, pooling), as depicted in Figure 3. All architectures
consist of 4 blocks total, followed by a hidden layer with 256 units, and then
a single unit to output age prediction. All convolutional and hidden layers
are followed by the eLU activation function. All models use $3\times 3\times
3$ convolutions of stride length 1. Before convolution, the inputs are padded
so that the dimensions of the outputs are equal to the dimensions of the
inputs, e.g. “same” padding in Tensorflow. Pooling layers are all $2\times
2\times 2$ max pooling with stride length 2.

Figure 3:  Both the baseline and proposed method use repeated copies of the
convolutional block above. The convolutional layers within a block always have
the same number of filters.  Table 1:  Architectures of the baseline and
proposed CNNs, from left to right. Blocks 1 through 4 are copies of the block
above. Both architectures have a single hidden layer with 256 units. The
hidden layer uses the same activation function as the convolutional layers,
namely eLu.

# Filters Approach Segmentation Input Channels Block 1 Block 2 Block 3 Block 4
Hidden Units Output Units Baseline No 1 8 16 32 64 256 1 Proposed Yes $k^{3}$
64 32 16 8 256 1

The Baseline architecture has 8 filters in its first block of convolutions,
and double that number in every successive block, for a total of 64 filters in
the final block. Our Proposed approach has 64 filters in its first block, and
half that every block after, for a total of 8 filters in the final block. A
summary of both architectures is presented in Table 1. Each layer of
convolutional filters can be thought of as a 5-D tensor of shape $C\times
3\times 3\times 3\times M$, where $C$ is the output dimension of the previous
layer, and $M$ is the number of filters (the output dimension of the current
layer). For example, in the first layer of the baseline network, the original
input image is convolved with a 5-D tensor of shape $1\times 3\times 3\times
3\times 8$, while in the first layer of the proposed network, the segmented
input image is convolved with a 5-D tensor of shape $k^{3}\times 3\times
3\times 3\times 64$.

###  4.3 Implementation Details

For all experiments, every model is trained five times with different random
initializations until performance on the validation set no longer improves,
approximately 700 epochs. Weights are initialized using Xavier initialization
(Glorot and Bengio, 2010). For each run, the model weights that achieved the
best performance on the validation set during training were saved. These
weights were restored after each run and the model was evaluated on the test
set. We report the average test performance (and standard deviation) over the
different initializations.

For computational efficiency, all images are pooled by a factor of 3 in each
dimension before being fed into the network, resulting in images of size
$41\times 49\times 41$. We pool the images by averaging consecutive $3\times
3\times 3$ blocks; the specific aggregation function performed is discussed in
the next section.

We train all models with the Adam optimizer (Kingma and Ba, 2014). We found no
significant performance improvement for any model tested with non-default
parameters, so we keep the defaults
($\alpha=0.001,\beta_{1}=0.9,\beta_{2}=0.999,\epsilon=10^{-8}$). We optimize
with respect to mean squared error (MSE) on the training set. We report both
MSE and mean absolute error (MAE) on the held-out test set. We implemented all
models using Tensorflow (Kingma and Ba, 2014) and ran all experiments on a
GeForce GTX 1080 Ti GPU.

To select the hyperparameter $k$ for our proposed model, we searched in the
range $[1,2,3,4]$, where 1 corresponds to no segmentation to and 4 corresponds
64 regions. $k$ was chosen based on the model that provided the lowest average
performance (MSE) on the validation set.

##  5 Results & Discussion

In this section, we compare the performance of our proposed approach to the
baseline architecture described in Section 4.2. We also compare the
performance of each proposed modification independently. We measure
performance in terms of both error and training time. We present an in-depth
follow-up analysis, in which we investigate performance across a range of
different settings.

###  5.1 Prediction Error

In Table 2 we report the MSE and the MAE for the task of predicting age from
structural brain images. As described in Section 4.3, reported performances
below represent the average and standard deviation of five different
initializations of each model evaluated on the held-out test set.

Table 2: Comparison of model performances for the task of predicting age on
the PNC dataset. The regional segmentation rate $k$, as described in Section
3, refers to how finely the brain is segmented before filters are learned in
each region. We find that 2 is the optimal setting for this task. The filter
layout refers to either increasing the number of filters (8-16-32-64), or
decreasing the number of filters (64-32-16-8), as the depth of the network
continues.

Approach Filter Layout $k$ MSE MAE Training Time (minutes) Baseline 8-16-32-64
- $3.72\pm 0.20$ $1.58\pm 0.06$ 52.5 Proposed 64-32-16-8 2 $\mathbf{3.03\pm
0.15}$ $\mathbf{1.43\pm 0.03}$ 40 Segmentation only 8-16-32-64 2 $3.76\pm
0.08$ $1.59\pm 0.03$ 30 Reverse layout only 64-32-16-8 - $3.64\pm 0.35$
$1.53\pm 0.05$ 117

Figure 4: Performance of proposed model plotted against the regional
segmentation rate, $k$. The points are offset from integer values to better
display both test and validation performance. Each point is the mean of five
runs, and the error bars are the standard deviation across those runs. This
chart demonstrates that for this particular problem, $k=2$ is a locally
optimal hyperparameter.

For the proposed model, we selected a regional segmentation rate of $k=2$
based on validation performance, resulting in eight different regions. Our
experiments show that the combination of segmenting the brain into eight
separate regions and reversing the order of filters results in improved age
prediction from structural brain data. Results indicate that it is the
combination of both modifications that improves performance; neither
modification alone lowers the prediction error from the baseline.

To understand how model performance varies with $k$, we plot both the
validation and test performance of our proposed method for each setting of $k$
in Figure 4. Each point is the average performance of the model across five
initializations, and error bars represent the standard deviation. For our
proposed architecture, segmenting the brain too finely or leaving the brain
whole results in decreased performance. We only used validation performances
to do hyperparameter selection, but the same trend holds for both test and
validation data.

###  5.2 Training Time

Table 2 also displays the training time of the evaluated models. Training time
refers to wall clock time, in minutes, per 700 epochs, the number of epochs
taken to train each model. Each model was trained on the same GPU under the
same load, for fair comparison. The table shows that reversing the layout of
the filters in a CNN is a computationally expensive modification. Doing so
more than doubles the training time from the equivalent baseline. More filters
earlier in the network results in more filters being applied when the image
data is larger and convolutions are more expensive.

Table 2 also shows that applying regional segmentation nearly halves training
time. This is due in part to the reduced number of fully connected parameters
in the proposed network, as discussed in section 5.4, and in part to the size
of the feature maps in the network after the first convolution. A similar
volume of data is input to the first convolution of the baseline and the
segmentation only network. In the former, the images are arranged as
($41\times 49\times 41\times 1$) arrays. In the latter, the images are
arranged as ($23\times 27\times 23\times 8$) arrays222Note that the dimensions
of the segmented image are not exactly equal to
($\lfloor\frac{41}{2}\rfloor\times\lfloor\frac{49}{2}\rfloor\times\lfloor\frac{41}{2}\rfloor\times
8$) due to the padding boundary described in section 3.. After the first
convolution, the number of feature maps in each section of the network is the
same, but the maps are roughly twice as small in each spatial dimension in the
segmentation only network, meaning the segmentation only network performs less
computation in the remaining convolutional layers.

The underlying implementation of the convolution operation also affects
training time. Consider an input image of shape ($N\times X\times Y\times
Z\times C$) where $N$ is the number of images in a batch, $X,Y,\textrm{ and
}Z$ are spatial dimensions, $C$ is the number of channels in the image and a
convolutional filter tensor of shape ($C\times\ell\times\ell\times\ell\times
M$) where $\ell$ is the size of the convolutional filter and $M$ is the number
of filters. The convolution can be thought of as a matrix multiplication
between a filter matrix of shape ($M\times C\ell^{3}$) and an input matrix of
shape333This assumes that the output dimension is made to be equal to the
input dimension using same padding. ($C\ell^{3}\times NXYZ$). This
multiplication is more efficient if the products $C\ell^{3}$ and $NXYZ$ are
similar in magnitude, rather than one being significantly larger than the
other (Chetlur et al., 2014). If the input image is a single-channel, 3D brain
image, and the filter size is small, the product $C\ell^{3}$ will be much
smaller than $NXYZ$. In our baseline network, $C\ell^{3}=27$, while
$NXYZ=4*41*49*41=329476$. Regional segmentation increases $C$ while decreasing
$XYZ$ by a factor of $k^{3}$, making the products $C\ell^{3}$ and $NXYZ$
closer in magnitude and resulting in a more efficient convolution.

To demonstrate how regional segmentation impacts the speed of convolutions, we
simulated the forward pass of the convolution operation on a randomly-
generated $4\times\frac{72}{k}\times\frac{72}{k}\times\frac{72}{k}\times
k^{3}$ image. We convolved the image with a randomly generated $k^{3}\times
3\times 3\times 3\times 8$ filter, similar to the first convolutional layer of
our baseline architecture. We varied the segmentation rate $k$ in the range
$[1,2,3,4,6,8,9,12,18,24,36]$, choosing the divisors of 72 and forgoing the
voxel boundary at the edges of regions to ensure that each convolution was
applied over exactly the same amount of data. Each data-point of Figure 5(a)
displays the result of timing 500 forward-pass convolutions at the indicated
segmentation rate. Figure 5(a) demonstrates that initially as the segmentation
rate increases and the products $C\ell^{3}$ and $NXYZ$ become closer in
magnitude, the convolution speeds up. But, as the image becomes too finely
segmented, $C\ell^{3}$ exceeds $NXYZ$ and the convolution slows down.

###  5.3 Robustness to Pooling Strategy

In both the baseline and proposed approach, we average pool the images before
feeding them into the network to speed up training. Average pooling smooths
out signals within the original images, but other pooling strategies may have
different properties. To better understand how this smoothing affects the
relative performance of to our proposed approach to the baseline, we evaluate
both models on two other pooling strategies. Max pooling corresponds to taking
the maximum value of every consecutive $3\times 3\times 3$ block. Naive
pooling refers to taking the first voxel and discarding the rest in every
consecutive $3\times 3\times 3$ block. While average pooling smooths signals,
max pooling emphasizes strong signals and naive pooling arbitrarily selects
voxels based on order. Given these differences, we hypothesized that relative
performance may vary based on the pooling type. To test this hypothesis, we
train the baseline and our proposed model separately on the training images
pooled with each strategy. We evaluated the MSE of each model on the held-out
test set pooled with each strategy, respectively (Figure 5(b)).

Average pooling is the best strategy for both the baseline model and our
proposed method. Our proposed method, however, is more robust to pooling
strategies than the baseline. Most importantly, our proposed method reports a
lower MSE for all pooling strategies than the baseline for any pooling
strategy.

(a)

(b)

Figure 5: Convolution time vs. segmentation rate on the left, and performance
based on pooling type on the right. (a) Time to compute the forward pass of a
convolution on a 3D, single-channel image batch after regional segmentation.
The plot shows that segmenting the image decreases training time up to a
point, but too finely segmenting the image results in increased training time.
(b) Plot of performance on pooling types. Average and Max refer to average and
max pooling, respectively, while Naive refers to downsampling the image by
keeping the first voxel in every consecutive $3\times 3\times 3$ block.
Average pooling is the best pooling type for both proposed and baseline. The
proposed model outperforms the baseline using any pooling type, and is less
sensitive to how the images are pooled.

###  5.4 Sensitivity to Number of Fully Connected Parameters

Our proposed model and the baseline model have a similar number of
convolutional parameters (233,280 compared to 219,672). However, because the
proposed model has fewer convolutional filters in later layers, it also has
fewer parameters in the fully connected layers, despite having the same number
of hidden units. Specifically, the proposed model has 16,640 weights in its
fully connected layers, while the baseline has 590,080 such weights, despite
both having 256 hidden units. In order to examine how much the number of fully
parameters impacts performance, we varied the number of units in the hidden
layer in our proposed approach in the range $[256,2496,4736,6976,9216]$. We
varied the number of hidden units in our baseline in the range
$[7,69,131,194,256]$. These ranges were chosen because they result in a
similar number of fully connected parameters for each model, within 1% of
difference. Note that $256$ is the default number of hidden units for both
models.

The resulting performance of each model across these different settings is
plotted in Figure 6(a). Having fewer fully connected parameters hurts the
baseline. In contrast, the reverse model remains within the margin of error
regardless of the number of fully connected parameters. These results suggest
that the difference in performance between the proposed approach and the
baseline is not due to the proposed approach having fewer fully connected
parameters.

###  5.5 Varying the Amount of Training Data

Lack of data is a large concern in MRI studies, because gathering research-
grade MRI scans is expensive and time consuming. It is not uncommon for
researchers to work with datasets of fewer than 200 patients (Zheng et al.,
2017; Meszlényi et al., 2017). Compared to other single-site datasets, the PNC
dataset is relatively large. In our setup, we have 524 training examples. To
better understand the impact of small training sizes on our conclusions, we
varied the number of training examples in the range $[100,200,300,400,524]$.
For each number of training examples, we independently and randomly selected
examples without replacement from the original training set of 524.

In Figure 6(b), we plot the number of training examples vs. model performance
in terms of MSE for both the baseline and proposed approach. The line in
purple plots the absolute difference in performance between the two models.
The performance of both approaches declines as we decrease the number of
training examples, as expected. However, the difference between the two models
increases as the number of training examples decreases. This suggests that our
method is more robust to small training sets. Furthermore, our proposed
approach exhibits less variation in performance, indicating that it is
potentially less more robust to weight initialization.

(a)

(b)

Figure 6: Investigating the robustness of observed improvements (a) Plot of
performance vs. number parameters in the fully connected layers. The left
side, 16,640 parameters, is the default for the proposed model. The right
side, 590,080 parameters, is the default for the baseline model. We vary the
number of fully connected parameters by changing the number of units in the
fully connected layer. Reducing the number of parameters for the baseline
model does not help performance: rather, it hurts it. (b) Plot of performance
vs. number of examples in training set size. The full training set, 524
examples, is on the right. The performance gap between our model and the
baseline increases as training examples decrease.

##  6 Conclusions

In this paper, we introduced two novel modifications to existing CNN
architectures, inspired by the structure of brain image data. First, we apply
different learned parameters to different regions in the brain. Second, we
start with a large number of filters in the first layer and decrease the
number of filters after each convolutional block. These modifications
encourage the network to learn region-specific patterns in the brain.
Combined, they are simple and easy to implement, and result in a model that
trains faster than a comparable CNN baseline. In addition, using our proposed
architecture, we consistently observed improved age prediction across multiple
downsampling types and varying amounts of training data.

Our work suggests that there is a larger space of possible architectural
choices to explore when adapting CNNs to brain imaging tasks. Due to
fundamental differences between brain images and natural images, previous
findings about how to build CNNs may not apply to neuroimaging data, or other
types of data outside of the domain of natural images. This work provides a
starting point to challenge assumptions existing architectures make when
applying CNNs to brain images.

There exists a variety of ways to expand this work. We did not focus
increasing the depth of the neural network, or techniques associated with
deeper networks, like inception modules (Szegedy et al., 2015) or skip
connections (He et al., 2016). However, if publicly available brain image
datasets increase by an order of magnitude, such techniques may be of
interest. Our work also did not focus on comparing pooling types within the
network, activation functions, or larger convolution sizes. However, a better
understanding of how to best design architectures for 3D brain volumes demands
investigation in these areas as well. We hope that our work inspires more
questions into how to best design CNNs for neuroimaging learning tasks.
Understanding how to challenge and adapt ideas from existing CNN architectures
will be critical to better predicting a multitude of clinically relevant
labels, from age to Alzheimer’s.

This work was supported by the Exercise and Sports Science Initiative (grant
no. U056408), and the National Science Foundation (NSF award no. IIS-1553146).
The views and conclusions in this document are those of the authors and should
not be interpreted as necessarily representing the official policies, either
expressed or implied, of the Exercise and Sports Science Initiative or the
NSF.

## References

  * Abadi et al. (2016) Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  Tensorflow: A system for large-scale machine learning.  In _OSDI_ , volume 16, pages 265–283, 2016.
  * Anirudh and Thiagarajan (2017) Rushil Anirudh and Jayaraman J Thiagarajan.  Bootstrapping graph convolutional neural networks for autism spectrum disorder classification.  _arXiv preprint arXiv:1704.07487_ , 2017.
  * Bernal et al. (2017) Jose Bernal, Kaisar Kushibar, Daniel S Asfaw, Sergi Valverde, Arnau Oliver, Robert Martí, and Xavier Lladó.  Deep convolutional neural networks for brain image analysis on magnetic resonance imaging: a review.  _arXiv preprint arXiv:1712.03747_ , 2017.
  * Cao et al. (2016) Lijun Cao, Zhi Liu, Xiaofu He, Yankun Cao, and Kening Li.  Mental disease feature extraction with mri by 3d convolutional neural network with multi-channel input.  In _Smart Cloud (SmartCloud), IEEE International Conference on_ , pages 224–227. IEEE, 2016.
  * Chetlur et al. (2014) Sharan Chetlur, Cliff Woolley, Philippe Vandermersch, Jonathan Cohen, John Tran, Bryan Catanzaro, and Evan Shelhamer.  cudnn: Efficient primitives for deep learning.  _arXiv preprint arXiv:1410.0759_ , 2014.
  * Chollet (2016) François Chollet.  Xception: Deep learning with depthwise separable convolutions.  _arXiv preprint_ , 2016.
  * Cohen and Welling (2016) Taco Cohen and Max Welling.  Group equivariant convolutional networks.  In _International Conference on Machine Learning_ , pages 2990–2999, 2016.
  * Cole et al. (2017) James H Cole, Rudra PK Poudel, Dimosthenis Tsagkrasoulis, Matthan WA Caan, Claire Steves, Tim D Spector, and Giovanni Montana.  Predicting brain age with deep learning from raw imaging data results in a reliable and heritable biomarker.  _NeuroImage_ , 163:115–124, 2017.
  * Dieleman et al. (2016) Sander Dieleman, Jeffrey De Fauw, and Koray Kavukcuoglu.  Exploiting cyclic symmetry in convolutional neural networks.  _arXiv preprint arXiv:1602.02660_ , 2016.
  * Dosenbach et al. (2010) Nico UF Dosenbach, Binyam Nardos, Alexander L Cohen, Damien A Fair, Jonathan D Power, Jessica A Church, Steven M Nelson, Gagan S Wig, Alecia C Vogel, Christina N Lessov-Schlaggar, et al.  Prediction of individual brain maturity using fmri.  _Science_ , 329(5997):1358–1361, 2010.
  * Farooq et al. (2017) Ammarah Farooq, SyedMuhammad Anwar, Muhammad Awais, and Saad Rehman.  A deep cnn based multi-class classification of alzheimer’s disease using mri.  In _Imaging Systems and Techniques (IST), 2017 IEEE International Conference on_ , pages 1–6. IEEE, 2017.
  * Franke et al. (2012) Katja Franke, Eileen Luders, Arne May, Marko Wilke, and Christian Gaser.  Brain maturation: predicting individual brainage in children and adolescents using structural mri.  _NeuroImage_ , 63(3):1305–1312, 2012.
  * Gens and Domingos (2014) Robert Gens and Pedro M Domingos.  Deep symmetry networks.  In _Advances in neural information processing systems_ , pages 2537–2545, 2014.
  * Glorot and Bengio (2010) Xavier Glorot and Yoshua Bengio.  Understanding the difficulty of training deep feedforward neural networks.  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ , pages 249–256, 2010.
  * He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * Kawahara et al. (2017) Jeremy Kawahara, Colin J Brown, Steven P Miller, Brian G Booth, Vann Chau, Ruth E Grunau, Jill G Zwicker, and Ghassan Hamarneh.  Brainnetcnn: convolutional neural networks for brain networks; towards predicting neurodevelopment.  _NeuroImage_ , 146:1038–1049, 2017.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Krizhevsky et al. (2012) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton.  Imagenet classification with deep convolutional neural networks.  In _Advances in neural information processing systems_ , pages 1097–1105, 2012.
  * LeCun et al. (1998) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner.  Gradient-based learning applied to document recognition.  _Proceedings of the IEEE_ , 86(11):2278–2324, 1998.
  * Meszlényi et al. (2017) Regina J Meszlényi, Krisztian Buza, and Zoltán Vidnyánszky.  Resting state fmri functional connectivity-based classification using a convolutional neural network architecture.  _Frontiers in neuroinformatics_ , 11:61, 2017.
  * Ngiam et al. (2010) Jiquan Ngiam, Zhenghao Chen, Daniel Chia, Pang W Koh, Quoc V Le, and Andrew Y Ng.  Tiled convolutional neural networks.  In _Advances in neural information processing systems_ , pages 1279–1287, 2010.
  * Nie et al. (2016) Dong Nie, Han Zhang, Ehsan Adeli, Luyan Liu, and Dinggang Shen.  3d deep learning for multi-modal imaging-guided survival time prediction of brain tumor patients.  In _International Conference on Medical Image Computing and Computer-Assisted Intervention_ , pages 212–220. Springer, 2016.
  * Sarraf et al. (2016) Saman Sarraf, Ghassem Tofighi, et al.  Deepad: Alzheimer′ s disease classification via deep convolutional neural networks using mri and fmri.  _bioRxiv_ , page 070441, 2016.
  * Satterthwaite et al. (2014) Theodore D Satterthwaite, Mark A Elliott, Kosha Ruparel, James Loughead, Karthik Prabhakaran, Monica E Calkins, Ryan Hopson, Chad Jackson, Jack Keefe, Marisa Riley, et al.  Neuroimaging of the philadelphia neurodevelopmental cohort.  _Neuroimage_ , 86:544–553, 2014.
  * Shaw et al. (2007) P. Shaw, K. Eckstrand, W. Sharp, J. Blumenthal, J. P. Lerch, D. Greenstein, L. Clasen, A. Evans, J. Giedd, and J. L. Rapoport.  Attention-deficit/hyperactivity disorder is characterized by a delay in cortical maturation.  _Proceedings of the National Academy of Sciences_ , 104(49):19649–19654, 2007.  ISSN 0027-8424.  10.1073/pnas.0707741104.  URL http://www.pnas.org/content/104/49/19649.
  * Stern (2012) Yaakov Stern.  Cognitive reserve in ageing and alzheimer’s disease.  _The Lancet Neurology_ , 11(11):1006–1012, 2012\.
  * Szegedy et al. (2015) Christian Szegedy, Wei Liu, Yangqing Jia, Pierre Sermanet, Scott Reed, Dragomir Anguelov, Dumitru Erhan, Vincent Vanhoucke, Andrew Rabinovich, et al.  Going deeper with convolutions.  CVPR, 2015.
  * Zeiler and Fergus (2014) Matthew D Zeiler and Rob Fergus.  Visualizing and understanding convolutional networks.  In _European conference on computer vision_ , pages 818–833. Springer, 2014.
  * Zheng et al. (2017) Lei Zheng, Jingyuan Zhang, Bokai Cao, Philip S Yu, and Ann Ragin.  A novel ensemble approach on regionalized neural networks for brain disorder prediction.  In _Proceedings of the Symposium on Applied Computing_ , pages 817–823. ACM, 2017.