# Development of Energy Models for Design Space Exploration of Embedded Many-
Core Systems

Christian Klarhorst†, Martin Flasskamp†, Johannes Ax†,
Thorsten Jungeblut†, Wayne Kelly*, Mario Porrmann†, and Ulrich Rückert†
<EMAIL_ADDRESS><EMAIL_ADDRESS>† Cognitronics and Sensor
Systems Group, * Science and Engineering Faculty CITEC, Bielefeld University,
Queensland University of Technology Bielefeld, Germany Brisbane, Australia

(2018)

###### Abstract.

This paper introduces a methodology to develop energy models for the design
space exploration of embedded many-core systems. The design process of such
systems can benefit from sophisticated models. Software and hardware can be
specifically optimized based on comprehensive knowledge about application
scenario and hardware behavior. The contribution of our work is an automated
framework to estimate the energy consumption at an arbitrary abstraction level
without the need to provide further information about the system. We validated
our framework with the configurable many-core system CoreVA-MPSoC. Compared to
a simulation of the CoreVA-MPSoC on gate level in a 28 nm FD-SOI standard cell
technology, our framework shows an average estimation error of about 4 %.

energy modeling, many-core architecture, design space exploration

††copyright: none††conference: High Performance Energy Efficient Embedded
Systems Workshop; January 2018; Manchester, United Kingdom††journalyear: 2018

##  1\. Introduction

Embedded system designers often have to deal with trade-offs. One is to
optimize a system that performs acceptable in the general case instead of
building an application-specific system. A general purpose CPU excels an ASIC
in flexibility but underlies it in performance. However, a high
configurability of hardware and software increases the available design space
tremendously.

An embedded system passes through many different design stages, involving
different groups of developers from the hardware, software and system design
domain. All have their own perspective on the system at different abstraction
levels. For instance a software developer does not know in detail how the
usage of different hardware components affects the performance or energy
consumption of an application. Just as a hardware developer may not have
enough information about specific application requirements in an early design
stage. Each group relies on estimations of specific system metrics that are
gathered by tools to perform design decisions. Building optimized designs for
specific use cases often requires an information flow between all involved
groups in order to find solutions that are holistically optimal, which leads
to a hardware software co-developing process. Therefore, it is appropriate to
postpone trade-off decisions during the design process in order to decide at
the right time, as informed as possible and for the right use case scenario.
However, keeping more degrees of freedom lead to a growing design space on all
different abstraction levels and the available and explored space better has
to be handled automatically.

Our target platform is the embedded MPSoC architecture CoreVA-MPSoC, which is
highly configurable at design time. The CoreVA-MPSoC toolchain has the
objective to create an interface between the different developer groups. It
consists of an MPSoC system specification in XML, parameterized VHDL SoC
components, a graphical user interface to configure MPSoCs, a cycle accurate
system simulator, a software backend based on the LLVM compiler framework and
our CoreVA-MPSoC compiler together with a broad range of single and many-core
applications. The VHDL SoC components are bundled with their API specification
into a component library. The interaction between these components is depicted
in Figure 1.

Figure 1. Interaction of CoreVA-MPSoC toolchain components

The system needs to be accurately modeled during the design process of
hardware and software despite of the high configurability. This can be
achieved by combining partial models of individual components like execution
units or memories to a complete system model. There are many different metrics
to evaluate the system depending on the specific abstraction level. In
(Flasskamp et al., 2016) we presented a model to estimate the throughput of an
MPSoC. In this work we focus on modeling the energy consumption and how these
models influence the whole MPSoC design process. Energy efficiency is a
crucial design aspect for all kinds of processor systems. Obviously, battery
powered embedded systems are limited regarding their run time. Additionally,
the power consumption of high performance systems limits their performance due
to heat dissipation. The power consumption can be addressed at all design
stages. On the hardware development side by choosing or dimensioning
components and on the software side by energy-aware compiler optimization.
Hence, we introduce an energy estimation framework for the CoreVA-MPSoC
toolchain.

Current research typically focuses on modeling effects like cache behavior (Li
et al., 2009) or components like NoC routers (Wolkotte et al., 2005) of a
specific system. In contrast, we want to cover a broad range of systems as for
instance the CoreVA-MPSoC as a highly configurable architecture. A huge design
space implies to trade off the required accuracy against the setup time of the
model. Therefore, our framework generalizes the model definition and automates
its creation. It is able to create an energy model based on the needs of
different user groups that variate in granularity, accuracy and available
time.

This work is organized as follows. The hardware and software parts of the
CoreVA-MPSoC project are introduced in Sections 2.1 and 2.2. An overview of
related work in the domain of modeling and estimating energy consumption is
presented in Section 3. Section 4 describes how energy estimation correlates
with system state transitions and Section 5 goes into details how the model
creation is automatically done by our energy estimation framework. The
analysis of our energy estimation framework is presented in Section 6.

##  2\. The CoreVA-MPSoC Project

As a target platform for our estimation framework presented in this work, we
use the self-developed many-core system CoreVA-MPSoC. The CoreVA-MPSoC is used
for application domains like signal processing, vision processing, and
encryption in embedded and energy-limited systems. In (Sievers et al., 2017)
the CoreVA-MPSoC is presented as a platform for Software-Defined-Radio
applications.

###  2.1. Hardware Architecture

Figure 2. Hierarchical CoreVA-MPSoC architecture

The CoreVA-MPSoC can be classified into the domain of many-core architectures.
Many-core architectures are defined by a very high number of small-sized CPU
cores, which are able to run applications with high resource efficiency.
Kalray’s MPPA-256 (de Dinechin et al., 2013) and Adapteva’s Epiphany are
typical examples for such many-core architectures. In contrast to the
Epiphany, our CoreVA-MPSoC features a hierarchical communication
infrastructure. The hierarchical infrastructure consists of a Network-on-Chip
(NoC) interconnect that couples several CPU clusters (cf. Figure 2). Within
each CPU cluster several VLIW CPU cores are tightly coupled via a bus-based
interconnect.

The hardware architecture of our CoreVA-MPSoC is highly configurable at design
time on all different hierarchy levels.

Basic building block and lowest hierarchy level of our MPSoC is the VLIW CPU
CoreVA that is designed to provide high resource efficiency (Lütkemeier et
al., 2013). According to the application domain, the number of VLIW slots,
ALUs, multiplication, and division units can be configured at design time.
Each CPU has local instruction and data memories, which are used as software-
managed scratchpad memories, without the use of caches.

A CPU cluster is used to tightly couple up to 32 of those CPU cores using a
bus-based interconnect. The number of CPUs and the interconnect fabric can be
configured at design time and supports either AMBA AXI4 or Wishbone standard,
each featuring a crossbar or a shared bus (Sievers et al., 2015a). Within a
cluster, a CPU can access the local data memory of every other CPU in a Non-
Uniform Memory Access (NUMA) fashion. Flexibility and communication bandwidth
can further be improved by integrating a shared data memory into a cluster,
which is tightly coupled to each CPU of the same cluster (Sievers et al.,
2015b).

For realizing MPSoCs with dozens or hundreds of CPU cores a Network on Chip
(NoC) is introduced to the CoreVA-MPSoC. The NoC acts as a second
interconnection hierarchy level and connects multiple CPU clusters. The NoC
features a packet based wormhole switching. Packets are segmented into flits
with 64 bit payload data. The NoC interconnect is built up of routers, each
having a configurable number of ports. This flexibility permits the
implementation of most common network topologies. Each CPU cluster is
connected to the NoC via a network interface (NI) (Ax et al., ress), which
bridges between the address-based communication of the clusters and the flow-
and packet-based communication of the NoC. Therefore, it provides a certain
number of independent uni-directional channels. These channels are
synchronized and the packet data is directly stored to and read from each
CPU’s local data memory or the cluster’s shared L1 data memory. Hence, the
CPUs can benefit from the low access latency of its local memory.

###  2.2. Programming Infrastructure

To program the MPSoC, the inter-CPU communication is encapsulated by our
communication C-library (API), transparently for both cluster and NoC
communication. The CPUs communicate via a block-based synchronization model
and unidirectional communication channels (Kelly et al., 2014). To run C
applications, we have created a C compiler tool chain based on the LLVM
compiler infrastructure. Our custom backend supports VLIW and SIMD
vectorization.

Nevertheless, efficiently programming a complete MPSoC with many CPUs is a
challenge for the programmer. For this reason we have developed a
parallelizing compiler for streaming applications to assist in programming the
CoreVA-MPSoC (Kelly et al., 2014; Flasskamp et al., 2016). This compiler,
called CoreVA-MPSoC compiler, processes applications written in the StreamIt
language (Thies et al., 2002). An application is represented by a structured
data flow graph, which describes the inherent parallelism of its tasks. Our
compiler for streaming applications searches for a valid partition with the
best performance. It utilizes an approach based on simulated annealing to map
the tasks of a program onto the individual CPUs of the MPSoC. A typical
optimization goal is to find a placement for the application’s tasks which
maximizes the throughput or minimizes the latency of an application. During
the partitioning process the compiler exploits three degrees of freedom to
alter an application’s data flow graph. Firstly, the compiler decides on which
processor a task is placed. Secondly, a task can be cloned to exploit data
parallelism. Thirdly, the granularity of work done in each iteration can be
increased to reduce the overhead of communication. These changes are called
mutations and further increase the search space for the partitioning
algorithm.

Every partition is judged regarding the achieved performance, by a Simulation
Based Estimation (SBE) model, which combines a single execution-based
simulation and an analytic approach (Flasskamp et al., 2016). MPSoC
configurations with the same number of CPUs can differ in many ways like
topology, communication infrastructure and memory architecture. All these
characteristics are modeled in our SBE to achieve accurate performance
estimations. In addition, a partition is also checked if any hardware limits
of a specified MPSoC configuration are exceeded. Therefore, our SBE model
estimates the consumption of data memory and enables the partitioning
algorithm to reject partitions exceeding the physically available memory.

Finally, the best partition is used to generate an individual C code file for
each CPU of the MPSoC. With the energy estimation framework presented in this
paper, our SBE model can be extended to estimate the expected energy
consumption of an application and perform energy-aware compilers optimization.

##  3\. Related Work

Energy consumption is still a research topic of much interest. On the one hand
power consumption imposes limits on the system performance due to e.g. heat
dissipation. On the other hand power efficiency often limits the diversity of
target applications. Research focuses on energy models, simulators and
configuration space exploration, often to characterize various system aspects.
Differences are in the purpose, target user group, level of detail and the
method to infer and analyze those models. The purpose might be mainly driven
by the needs of the various development process levels to provide benefits for
the involved user groups. Research was often used to limit the design-space
and point out design inefficiencies as early as possible (Vasilakis, 2015).

In the area of hardware development research focuses on different levels of
detail. On gate level Roselló et al. (Rossello et al., 2005) present an energy
model for domino CMOS gates based on a detailed description of internal
capacitance switching and discharging currents. Many models target specific
processor components for example Dramsim (Rosenfeld et al., 2011) and Cacti
(Wilton and Jouppi, 1996) are well-known tools used for timing and energy
analysis of DRAMS and caches. In (Wolkotte et al., 2005) Wolkotte et al. focus
on NoC routers while the Orion2 (Kahng et al., 2009) network simulator models
complete Network on Chips. Energy simulation tools for different kinds of
processors are provided by the frameworks Watch (Brooks et al., 2000) for CPUs
and GPUWatch (Leng et al., 2013) for GPUs. Just as by the work of Vasilakis et
al. (Vasilakis, 2015) for ARM processors and Jordans et al. (Jordans et al.,
2013) for VLIW processors. An instruction-level energy characterization of the
Adapteva Epiphany SoC was done by Ortiz et al. in (Ortiz et al., 2017). They
measured the energy consumption of the SoC with a digital current meter while
executing a number of microbenchmarks. The McPat framework (Li et al., 2009)
focuses on complete SoCs by describing fundamental components of multi-core
and many-core processor configurations in a flexible XML interface. Specific
properties and implementation details like cache behavior and timing are
described quantitatively. By evaluating the model McPat estimates power, area
and timing. Hesse et al. created a power consumption model of a complete
circuit board of a wireless sensor node (Hesse et al., 2016). The model
containing sensor, microcontroller, wireless transceiver and power supply is
utilized for choosing components at an early design stage. In the compiler
research area the configuration space of compiler options and its impact on
power consumption was explored by Pallister et al. (Pallister et al., 2013).
They performed power measurements on five embedded platforms to ensure all
architectural effects are captured. All user groups focus on different levels
of detail or granularity and admit different error margins.

Besides the target platform or components the aforementioned work further
differs in their methodology. The McPat framework (Li et al., 2009) estimates
energy based on a model specification while other work gathers data from
various sources to fit models like Ortiz et al. did by measuring the MPSoC
(Ortiz et al., 2017). Multiple options to create and gather data for the
second option may be found in research. Microbenchmarks which are often
handwritten are used for the behavior creation in contrast to real
applications or benchmark suits. Furthermore data can be gathered by
simulating those systems or by measuring final ASICs under real world
conditions. The simulation method provides fine granular power analysis, e.g.
it is possible to gather power data for all components of a CPU using gate
level analysis at the cost of an often compute intensive but scalable
simulation Finally, the simulation method is able to provide data before tape
out and it is able to simulate various environment conditions. While the
measurement method provides better accuracy because it incorporates all
physical effects and uses real word conditions.

Figure 3. Exemplary energy estimation of instruction set simulation by model
created at gate level

##  4\. Energy Estimation Methodologies

Energy consumption of digital circuits originates from different physical
effects. The static power consumption is the result of leakage current of the
transistors. While the dynamic power consumption is composed of switching
current, a transistor switches from one logic state to another, and the
switching frequency dependent load of charging external load capacitance. The
state of the entire system can be described by the states of all components.

Our idea is to describe the system state and correlate energy consumption with
the system state transitions. This also gives us the ability to correlate
energy consumption from a lower abstraction level with high-level state
transitions. Figure 3 shows exemplary how the result of gate level simulations
can be utilized to develop an energy model on logic level. This model can be
used together with a cycle accurate system simulator to estimate the energy
consumption.

### Model Constants and Functions

The state of a system can be described at all abstraction levels. An overview
of abstraction levels shows Table 1. On circuit level the system state is
described by differential equations of the transistors. Similarly,
instructions define the system state at the higher cycle abstraction level.

Table 1. Abstraction levels Abstraction Level | System State | Granularity
---|---|---
Circuit | Transistors | Differential Equations
Logic | Gates | Boolean Functions
RTL | Registers | RTL Spec.
Cycle | Instructions | Instruction
Task | Functions | System Spec.
Application | Program | Application Spec.

Without additional information we could now perform energy measurements and
correlate them with state transitions. A comprehensive model containing all
possible state transformations could be created by performing
$n_{meas}=n_{states}\times n_{states}$ measurements. But unfortunately an
incalculable amount of ways to describe the system state exist.

This correlation is visualized in Figure 4 on the basis of the CoreVA-MPSoC.
Figure 4 (a) show a sparse model which only distinguishes whether a component
was used (coloured) or unused (grey) during the period of measurement. In
contrast the model in (c) breaks down the state of each component to a very
fine-grained state description, visualized by multicolored bar charts on the
CPU components and pie charts on the NoC routers. For instance this model
considers the variation of power consumption for different CPU instructions. A
trade-off between both attempts depicts model (b), which considers the time a
used components was in active or idle state. In this case the state
description of a component contains no redundant states.

Figure 4. Abstraction levels of system state description per period of
measurement

On cycle level system state transitions might be caused by a sequence of
instructions as a result of an execution of a function. The energy consumption
per instruction forms the model constants on cycle level (cf. Figure 3).
Setting up the model includes the creation of system behavior, the capturing
of the system’s energy consumption during simulation and the final model
determination. Accordingly, the setup time for a simple CPU core model would
include the generation of micro benchmarks, gate level simulations and fitting
the model constants.

On each abstraction level the system state is modeled with a different
granularity. Therefore, we added a function to the model to transform the
state description from simulation into the state description used by the
model. It is able to add, change and remove state information and provides the
necessary flexibility while it may also include additional system assumptions
to simplify e.g. reduce the system description. This reduction results in
faster model evaluation and setup time. The same model constants can be reused
on different abstraction levels by transforming the system state with
different model functions.

The model is evaluated by executing the model function and combining the
resulting state with the model parameters. Since no simulation is required the
model evaluation is only a fast computation.

### Exemplary MPSoC Component Models

As an example the following paragraph describes the creation of an model for
an ideal NoC communication. The NoC of the CoreVA-MPSoC is organized as a
2D-mesh, described in Section 2.1. A transmission can be characterized by the
coordinates of the involved CPUs and the size of the transmitted buffer. So a
naive model would contain all possible permutations of the coordinates and all
possible packet sizes resulting in a huge model.

The path on which a packet passes through the NoC can be characterized by a
sequence of independent hops (routers). That way, the foregoing model can be
drastically reduced by removing the redundancy of individual links. The
reduced model would then contain all number of hops and possible packet sizes
together with a function to define the new distance metric number of hops, for
example by $n_{\text{hops}}=\text{manhattan\\_dist}(x,y)$.

More knowledge about the hardware could be incorporated into the model to
reduce it even further. For example, one could assume that the energy
consumption will increase linearly with the packet size (under the assumption
of no congestion). A model containing all available number of hops together
with two constants for the linear functions would be reduced again but
possibly at the expense of an increasing estimation error. More detailed
results and further analysis follow in Section 6.

##  5\. Energy Modeling Framework

The contribution of our work is a framework to estimate the energy consumption
at an arbitrary abstraction level without the need to provide further
information about the system. However, the framework allows for adding
additional knowledge to the system to improve estimation accuracy or to
accelerate the model setup time.

Back annotated simulations on gate level could be employed to create a model,
as long as one is in possession of the RTL (VHDL, Verilog) sources of the
system and a process design kit or corresponding IP-Cores. Afterwards the
simulation results can be analyzed by commercial tools like Cadence© Voltus™
to obtain an energy estimation. Unfortunately, comprehensive gate level
simulations for a complex system like the configurable CoreVA-MPSoC
architecture would imply a large number of simulations. Accordingly, these
require a lot of simulation time and resources as well as tool licenses.
Therefore, our framework enables the developer to obtain an energy estimation
faster and with an acceptable loss of accuracy with the help of reducing the
number of necessary simulations by including system knowledge into the model.

We require a set of system states and energy information to find the
correlation parameters for the specific energy model. The system states (the
system behavior) are created by a previously selected application or by
microbenchmarks and a gate level simulation. This simulation provides the
necessary state information together with the occurred switching activity data
that is fed into the tools for energy measurement. We further provide a trace
infrastructure, implemented up to the C level, in order to provide a method to
create additional, high level system state information that could be used in
the model or by the model function. We developed a tool to automatically
create a broad range of microbenchmarks based on user preferences and on the
available information in the toolchain. The tool is able to variate the use
pattern of all SoC components that have a specified API description. This
includes for example all different communication and memory components. It is
further able to create microbenchmarks for the CoreVA CPU cores. This is done
automatically by parsing the instruction set description already available in
the software toolchain. The tool always generates setup code in front of every
test case in order to bring the CoreVA into a defined state to have a
reproductive measuring environment.

##  6\. Validation

There are many effects in a System on Chip that have an impact on the energy
consumption. This work focuses on a many-core system with a 2D-mesh NoC which
connects multiple clusters each containing four fully connected CoreVA CPUs.
The CoreVA CPU cores include two VLIW slots with one load/store and
multiplication/division unit. Each CPU is configured with 16 kB of instruction
memory organized in two banks and 16 kB of local data memory equally
organized. Each cluster additionally contains 64 kByte of shared data memory,
tightly coupled to each CPU and organized in eight banks. Each memory bank in
the system contains 2 k words with 32 bit word-size. Each network interface
(NI) supports 128 synchronous channels with an effective bandwidth of 64 bit
per clock cycle (cf. Section 2.1). The cluster interconnects are carried out
as crossbars with a width of 64 bit.

Figure 5. Physical layout of a CPU cluster node. Left: Placed design, Right:
Routed design

One CPU cluster occupies 0.817 mm2 of chip area using a 28 nm FD-SOI standard
cell technology: STMicroelectronics, 10 metal layer, worst case corner 1.0 V,
125 ${}^{\circ}\text{C}$. Allows for maximum clock frequency of 700 MHz. One
CoreVA CPU dissipates 4.5 mW during idle state and up to 30 mW during active
computation. A full CPU cluster has a total power consumption of up to 150 mW.
The layout corresponds with the design shown in Figure 5. In addition to the
CPUs, memory and the crossbar, each cluster layout includes NoC components,
like the NI and a router with connection ports to all four directions.
Multiple of these clusters can be directly strung together on a top level P&R
step to build a scalable MPSoC with a 2D-Mesh NoC, as we have presented in (Ax
et al., 2017). As we pointed out earlier, many scenarios exist where accurate
energy metrics are important for the final result. This already begins in an
early design stage, e.g. during the specification of the CPU’s feature and
instruction set. With our toolchain it is possible to get a fast indicator for
the energy consumption of the instruction set. We automatically create
microbenchmarks, which variate just the instructions while setting all
registers and the memory content to a specific value. For instance, the
description of CoreVA CPU instruction set (cf. Section 5) contains 155
different definitions of instructions, which results in 20,093 different
instruction groups for the two VLIW slots (not all instructions could be
placed on an arbitrary slot). All instruction groups were analyzed for three
data patterns.

Figure 6. Average power consumption of CoreVA CPU instruction groups for two
VLIW slots divided into CPU components

This part of the model leads to Figure 6 which shows the mean power
consumption of the CoreVA CPU instructions for the given configuration. The
energy consumption is shown in pico joule and broken down into the CPU
components core, instruction and data memory as well as a small constant
unclassified section. Instructions of similar behavior are categorized into
the groups NOP, SIMD, Load/Store and Miscellaneous due to illustrative
purposes. An empty VLIW slot is depicted as Empty. Non SIMD arithmetic
instructions are included in the Miscellaneous instruction group.

As one might expect, the NOP instruction leads to the lowest energy
consumption of the CPU core. Whereas the most energy is consumed when a SIMD
instruction is executed due to the highest core activity. The data memory
shown in blue has obviously an impact only on the power consumption of
Load/Store instructions. The worst case difference between the three data
patterns mentioned above is 2.6 mW. The energy consumption of the instruction
memory depends on the VLIW instruction format. A VLIW instruction for a CoreVA
CPU with two slots may contain instructions for both or only one slot. The
latter case results in instruction compression to save space in memory and
energy on access as seen in Figure 6. The provided hardware macro of a single
CoreVA CPU has a minimum power usage of 15.3 mW and a maximum of 38.0 mW. The
model creation needs 55 h of gate level simulation and 70 days of power
analysis time to run 60,279 different tests on a single Intel® Xeon® E5-1650
v4 CPU core. Actually, the tests are distributed onto 26 CPU cores in
parallel.

Many detailed effects on the power consumption were left aside, like the
position of the instruction or the data dependency. Therefore, we investigated
how the specific instruction memory address effects the energy consumption of
the instruction memory. We create microbenchmarks where the position of the
instruction is variated while the rest of the system state remains constant.
An 16 kb instruction memory is able to store 2048 2-Slot and up to 4096 1-Slot
CoreVA instructions because of instruction compression. The energy consumption
for the first 800 instruction memory addresses is shown in Figure 7. All
addresses were filled with the same instruction hence the energy consumption
only depends on the accessed address and is not a function on the memory
content.

Figure 7. Spatial correlation between energy consumption of nop instruction
and address in instruction memory

The position creates an energy consumption range of 1.7 pJ for two VLIW slot
instructions and 3.2 pJ for a single slot, which can be explained by the
internal interconnect of the memory blocks. We assume that the arrangement of
the data points is caused by the tree like structure of the memory address
decoder.

A communication transfer between two CPU cores across the NoC involves the two
CPU cores itself, routers, the interfaces between CPU cluster and NoC (NI) and
other communication infrastructure like the bus interconnect within a CPU
cluster. Our toolchain is able to generate benchmarks with the tool presented
in Section 5 to analyze communication characteristics in an automated way
based on an API description. Therefore, it is possible to variate the packet
size for a NoC communication between two cores on different clusters. Figure 8
shows the energy consumption of specific transfer sizes. The packet size
variates from 4 to 1024 byte in 4 byte increments. The zoomed view includes 36
data points from 374 to 446 in 2 byte increments. Due to the 64 bit flit size,
into which all NoC packets are segmented (cf. Section 2.1), the plot shows a
staircase function. One bar depicts the energy consumption in nJ and is
separated into the shares of the involved components. CPU synchronization,
data memories, bus, router and the NoC interface are shown in red, blue,
green, bright red and gray. The unclassified part shown in yellow contains,
among others, especially the parts of the CPU cluster interconnect, which are
not assigned to the bus part. The costs for CPU synchronization are software
costs to synchronize the communication channel. More details about these
software costs for communication are presented in (Ax et al., ress).

Figure 8. Energy consumption in nJ of NoC transmissions with varying packet
size separated into components

Generating the results of 256 data points for the plot of Figure 8 took 40 h
of simulation time and 170 h for power analysis. The data was further
approximated with a linear and a staircase function as shown in Section 4. We
used 16 data points around the center for the fit in order to reduce the
simulation time to 3 h and the power analysis time to 11 h. The linear
function shown as a blue line had a maximum error of 62 pJ while the staircase
function shown as a green line had a maximum error of 52 pJ.

We validated the applicability of our approach by estimating the energy
consumption for a set of five many-core benchmarks from the StreamIt benchmark
suite (Thies et al., 2002). The analyzed benchmarks are chosen from the signal
processing, matrix multiplication and sorting algorithm domain. Compared to a
gate level simulation the average estimation error is about 4% with a maximum
absolute error of 10% for a simplified model. The model considers all
instruction types and three data pattern but not the instruction position in
the memory. However, the accuracy depends on how many effects are taken into
account when establishing the model.

##  7\. Conclusion

This paper presents an approach to develop and use energy models for a design
space exploration of embedded many-core systems. The contribution of our work
is a framework to estimate the energy consumption at an arbitrary abstraction
level without the need to provide further information about the system.
Therefore, we developed an automatic tool flow to fill the energy model with
energy information gathered from gate level simulations of basic system states
of different components of a many-core system. Once the energy model is
filled, it can be used on higher abstraction levels to get fast and accurate
energy estimations for specific applications and hardware configurations of
the many-core system.

To validate our framework, we used our, at design time, highly configurable
many-core system CoreVA-MPSoC as a target architecture. Compared to a gate
level simulation of the CoreVA-MPSoC in a 28 nm FD-SOI standard cell
technology, the average estimation error of our energy model framework is
about 4 %. Compared to other approaches, like Ortiz et al. (Ortiz et al.,
2017), we are able to analyze the energy consumption more fine-grained without
the restriction to measure an existing ASIC. Hence, our model contains
internal MPSoC components like instruction memory and VLIW slots as well.
However, Ortiz et al. achieve a more accurate total energy estimation at
instruction level. The design process of highly configurable embedded systems
like the CoreVA-MPSoC benefits from such a sophisticated model. Applications
may be specifically optimized for tailored systems by modeling the different
variants and characteristics of system components. A trade-off can be found
based on comprehensive knowledge about application scenario and hardware
behavior. In future work, our MPSoC compiler will use the presented energy
model to automatically map streaming applications to our CoreVA-MPSoC,
optimized for energy.

## Acknowledgment

The reported research was supported through project grants KogniHome (German
Federal Ministry of Education and Research (BMBF) Grant No. 16SV7054K), the
Cluster of Excellence Cognitive Interaction Technology ”CITEC” (EXC 277) at
Bielefeld University and the BMBF Leading-Edge Cluster ”Intelligent Technical
Systems OstWestfalenLippe” (it’s OWL).

## References

  * (1)
  * Ax et al. (ress) Johannes Ax et al. 2015\. In press.  System-Level Analysis of Network Interfaces for Hierarchical MPSoCs. In NoCArc. ACM.  https://doi.org/10.1145/2835512.2835513
  * Ax et al. (2017) Johannes Ax et al. 2017\.  Comparing synchronous, mesochronous and asynchronous NoCs for GALS based MPSoC. In MCSoC ’17.
  * Brooks et al. (2000) D. Brooks, V. Tiwari, and M. Martonosi. 2000.  Wattch: a framework for architectural-level power analysis and optimizations. In Proceedings of 27th International Symposium on Computer Architecture (IEEE Cat. No.RS00201).
  * de Dinechin et al. (2013) Benoit Dupont de Dinechin et al. 2013\.  A clustered manycore processor architecture for embedded and accelerated applications. In HPEC ’13. IEEE.  https://doi.org/10.1109/HPEC.2013.6670342
  * Flasskamp et al. (2016) Martin Flasskamp et al. 2016\.  Performance estimation of streaming applications for hierarchical MPSoCs. In RAPIDO ’16, Vol. 18-January. ACM Press.  https://doi.org/10.1145/2852339.2852342
  * Hesse et al. (2016) Marc Hesse, Michael Adams, Timm Hormann, and Ulrich Rückert. 2016.  Towards a comprehensive power consumption model for wireless sensor nodes. In BSN ’16.  https://doi.org/10.1109/BSN.2016.7516293
  * Jordans et al. (2013) R. Jordans, R. Corvino, L. Józwiak, and H. Corporaal. 2013\.  An Efficient Method for Energy Estimation of Application Specific Instruction-Set Processors. In Euromicro DSD ’13.  https://doi.org/10.1109/DSD.2013.120
  * Kahng et al. (2009) A. B. Kahng, Bin Li, L. S. Peh, and K. Samadi. 2009\.  ORION 2.0: A fast and accurate NoC power and area model for early-stage design space exploration. In DATE ’09.  https://doi.org/10.1109/DATE.2009.5090700
  * Kelly et al. (2014) Wayne Kelly et al. 2014\.  A Communication Model and Partitioning Algorithm for Streaming Applications for an Embedded MPSoC. In SoC ’14. IEEE.  https://doi.org/10.1109/ISSOC.2014.6972436
  * Leng et al. (2013) Jingwen Leng et al. 2013\.  GPUWattch: Enabling Energy Optimizations in GPGPUs. In ISCA ’13.  https://doi.org/10.1145/2485922.2485964
  * Li et al. (2009) S. Li, J. H. Ahn, R. D. Strong, J. B. Brockman, D. M. Tullsen, and N. P. Jouppi. 2009\.  McPAT: An integrated power, area, and timing modeling framework for multicore and manycore architectures.  MICRO ’09.  https://doi.org/10.1145/1669112.1669172
  * Lütkemeier et al. (2013) Sven Lütkemeier et al. 2013\.  A 65 nm 32 b Subthreshold Processor With 9T Multi-Vt SRAM and Adaptive Supply Voltage Control.  IEEE Journal of Solid-State Circuits 48, 1 (2013).
  * Ortiz et al. (2017) Gabriel Ortiz, Lars Svensson, Erik Alveflo, and Per Larsson-Edefors. 2017.  Instruction level energy model for the Adapteva Epiphany multi-core processor. In CF ’17. ACM Press.  https://doi.org/10.1145/3075564.3078892
  * Pallister et al. (2013) James Pallister et al. 2013\.  Identifying Compiler Options to Minimise Energy Consumption for Embedded Platforms.  CoRR abs/1303.6485 (2013).
  * Rosenfeld et al. (2011) P. Rosenfeld et al. 2011\.  DRAMSim2: A Cycle Accurate Memory System Simulator.  IEEE L-CA ’11 10, 1 (2011).  https://doi.org/10.1109/L-CA.2011.4
  * Rossello et al. (2005) J.L. Rossello, C. de Benito, and J. Segura. 2005.  A compact gate-level energy and delay model of dynamic CMOS gates.  IEEE Transactions on Circuits and Systems II: Express Briefs 52\.  https://doi.org/10.1109/TCSII.2005.851992
  * Sievers et al. (2015a) Gregor Sievers et al. 2015a.  Comparison of Shared and Private L1 Data Memories for an Embedded MPSoC in 28nm FD-SOI. In MCSoC ’15.  https://doi.org/10.1109/MCSoC.2015.25
  * Sievers et al. (2015b) Gregor Sievers et al. 2015b.  Evaluation of Interconnect Fabrics for an Embedded MPSoC in 28nm FD-SOI. In ISCAS ’15. IEEE.
  * Sievers et al. (2017) Gregor Sievers et al. 2017\.  The CoreVA-MPSoC: A Multiprocessor Platform for Software-Defined Radio.  Springer International Publishing, Cham.  https://doi.org/10.1007/978-3-319-49679-5_3
  * Thies et al. (2002) William Thies et al. 2002\.  StreamIt: A Language for Streaming Applications. In CC ’02. Springer.  https://doi.org/10.1007/3-540-45937-5_14
  * Vasilakis (2015) Evangelos Vasilakis. 2015\.  An Instruction Level Energy Characterization of ARM Processors.  Technical Report March.
  * Wilton and Jouppi (1996) S. J. E. Wilton and N. P. Jouppi. 1996.  CACTI: an enhanced cache access and cycle time model.  IEEE Journal of Solid-State Circuits 31, 5 (May 1996).  https://doi.org/10.1109/4.509850
  * Wolkotte et al. (2005) P.T. Wolkotte, G.J.M. Smit, N. Kavaldjiev, J.E. Becker, and Jürgen Becker. 2005.  Energy Model of Networks-on-Chip and a Bus. In SoC ’05. IEEE.  https://doi.org/10.1109/ISSOC.2005.1595650


# Woodbury Transformations for
Deep Generative Flows

You Lu
Department of Computer Science
Virginia Tech
Blacksburg, VA
<EMAIL_ADDRESS>
&Bert Huang
Department of Computer Science
Tufts University
Medford, MA
<EMAIL_ADDRESS>

###### Abstract

Normalizing flows are deep generative models that allow efficient likelihood
calculation and sampling. The core requirement for this advantage is that they
are constructed using functions that can be efficiently inverted and for which
the determinant of the function’s Jacobian can be efficiently computed.
Researchers have introduced various such flow operations, but few of these
allow rich interactions among variables without incurring significant
computational costs. In this paper, we introduce _Woodbury transformations_ ,
which achieve efficient invertibility via the Woodbury matrix identity and
efficient determinant calculation via Sylvester’s determinant identity. In
contrast with other operations used in state-of-the-art normalizing flows,
Woodbury transformations enable (1) high-dimensional interactions, (2)
efficient sampling, and (3) efficient likelihood evaluation. Other similar
operations, such as 1x1 convolutions, emerging convolutions, or periodic
convolutions allow at most two of these three advantages. In our experiments
on multiple image datasets, we find that Woodbury transformations allow
learning of higher-likelihood models than other flow architectures while still
enjoying their efficiency advantages.

##  1 Introduction

Deep generative models are powerful tools for modeling complex distributions
and have been applied to many tasks such as synthetic data generation [27,
39], domain adaption [40], and structured prediction [33]. Examples of these
models include autoregressive models [13, 28], variational autoencoders [21,
31], generative adversarial networks [11], and normalizing flows [6, 30, 7,
22]. Normalizing flows are special because of two advantages: They allow
efficient and exact computation of log-likelihood and sampling.

Flow-based models are composed of a series of invertible functions, which are
specifically designed so that their inverse and determinant of the Jacobian
are easy to compute. However, to preserve this computational efficiency, these
functions usually cannot sufficiently encode dependencies among dimensions of
a variable. For example, affine coupling layers [6] split a variable to two
parts and require the second part to only depend on the first. But they ignore
the dependencies among dimensions in the second part.

To address this problem, Dinh et al. [6, 7] introduced a fixed permutation
operation that reverses the ordering of the channels of pixel variables.
Kingma and Dhariwal [22] introduced a 1$\times$1 convolution, which are a
generalized permutation layer, that uses a weight matrix to model the
interactions among dimensions along the channel axis. Their experiments
demonstrate the importance of capturing dependencies among dimensions.
Relatedly, Hoogeboom et al. [15] proposed emerging convolution operations, and
Hoogeboom et al. [15] and Finz et al. [9] proposed periodic convolution. These
two convolution layers have $d\times d$ kernels that can model dependencies
along the spatial axes in addition to the channel axis. However, the increase
in representational power comes at a cost: These convolution operations do not
scale well to high-dimensional variables. The emerging convolution is a
combination of two autoregressive convolutions [10, 23], whose inverse is not
parallelizable. To compute the inverse or determinant of the Jacobian, the
periodic convolution requires transforming the input and the convolution
kernel to Fourier space. This transformation is computationally costly.

In this paper, we develop _Woodbury transformations_ for generative flows. Our
method is also a generalized permutation layer and uses spatial and channel
transformations to model dependencies among dimensions along spatial and
channel axes. We use the Woodbury matrix identity [37] and Sylvester’s
determinant identity [35] to compute the inverse and Jacobian determinant,
respectively, so that both the training and sampling time complexities are
linear to the input variable’s size. We also develop a memory-efficient
variant of the Woodbury transformation, which has the same advantage as the
full transformation but uses significantly reduced memory when the variable is
high-dimensional. In our experiments, we found that Woodbury transformations
enable model quality comparable to many state-of-the-art flow architectures
while maintaining significant efficiency advantages.

##  2 Deep Generative Flows

In this section, we briefly introduce the deep generative flows. More
background knowledge can be found in the appendix.

A normalizing flow [30] is composed of a series of invertible functions
$\mathbf{f}=\mathbf{f}_{1}\circ\mathbf{f}_{2}\circ...\circ\mathbf{f}_{K}$,
which transform $\mathbf{x}$ to a latent code $\mathbf{z}$ drawn from a simple
distribution. Therefore, with the _change of variables_ formula, we can
rewrite the log-likelihood $\log p_{\theta}(\mathbf{x})$ to be

$\log p_{\theta}(\mathbf{x})=\log
p_{Z}(\mathbf{z})+\sum_{i=1}^{K}\log\left|\det\left(\frac{\partial\mathbf{f}_{i}}{\partial\mathbf{r}_{i-1}}\right)\right|,$
(1)

where $\mathbf{r}_{i}=\mathbf{f}_{i}(\mathbf{r}_{i-1})$,
$\mathbf{r}_{0}=\mathbf{x}$, and $\mathbf{r}_{K}=\mathbf{z}$.

Flow-based generative models [6, 7, 22] are developed on the theory of
normalizing flows. Each transformation function used in the models is a
specifically designed neural network that has a tractable Jacobian determinant
and inverse. We can sample from a trained flow $\mathbf{f}$ by computing
$\mathbf{z}\sim
p_{Z}(\mathbf{z}),\quad\mathbf{x}=\mathbf{f}^{-1}(\mathbf{z})$.

There have been many operations, i.e., layers, proposed in recent years for
generative flows. In this section, we discuss some commonly used ones, and
more related works will be discussed in Section 4.

Actnorm layers [22] perform per-channel affine transformations of the
activations using scale and bias parameters to improve training stability and
performance. The actnorm is formally expressed as
$\mathbf{y}_{:,i,j}=\mathbf{s}\odot\mathbf{x}_{:,i,j}+\mathbf{b}$, where both
the input $\mathbf{x}$ and the output $\mathbf{y}$ are $c\times h\times w$
tensors, $c$ is the channel dimension, and $h\times w$ are spatial dimensions.
The parameters $\mathbf{s}$ and $\mathbf{b}$ are $c\times 1$ vectors.

Affine coupling layers [6, 7] split the input $\mathbf{x}$ into two parts,
$\mathbf{x}_{a},\mathbf{x}_{b}$. And then fix $\mathbf{x}_{a}$ and force
$\mathbf{x}_{b}$ to only relate to $\mathbf{x}_{a}$, so that the Jacobian is a
triangular matrix. Formally, we compute

$\displaystyle\mathbf{x}_{a},\mathbf{x}_{b}=\text{split}(\mathbf{x}),\quad\quad\quad\quad~{}~{}\mathbf{y}_{a}=\mathbf{x}_{a},$
$\displaystyle\mathbf{y}_{b}=\mathbf{s}(\mathbf{x}_{a})\odot\mathbf{x}_{b}+\mathbf{b}(\mathbf{x}_{a}),\quad\mathbf{y}=\text{concat}(\mathbf{y}_{a},\mathbf{y}_{b}),$

where $\mathbf{s}$ and $\mathbf{b}$ are two neural networks with
$\mathbf{x}_{a}$ as input. The split and the concat split and concatenate the
variables along the channel axis. Usually, $s$ is restricted to be positive.
An additive coupling layer is a special case when $\mathbf{s}=\mathbf{1}$.

Actnorm layers only rescale the dimensions of $\mathbf{x}$, and affine
coupling layers only relate $\mathbf{x}_{b}$ to $\mathbf{x}_{a}$ but omit
dependencies among different dimensions of $\mathbf{x}_{b}$. Thus, we need
other layers to capture local dependencies among dimensions.

Invertible convolutional layers [22, 15, 9] are generalized permutation layers
that can capture correlations among dimensions. The 1$\times$1 convolution
[22] is $\mathbf{y}_{:,i,j}=\mathbf{Mx}_{:,i,j}$, where $\mathbf{M}$ is a
$c\times c$ matrix. The Jacobian of a 1$\times$1 convolution is a block
diagonal matrix, so that its log-determinant is $hw\log|\det(\mathbf{M})|$.
Note that the 1$\times$1 convolution only operates along the channel axis and
ignores the dependencies along the spatial axes.

(a) Flow step

(b) Multi-scale arch.

Figure 1: Overview of architecture of generative flows. We can design the flow
step by selecting a suitable convolutional layer and a coupling layer based on
the task. Glow [22] uses 1$\times$1 convolutions and affine coupling.

Emerging convolutions [15] combine two autoregressive convolutions [10, 23].
Each autoregressive convolution masks out some weights to force an
autoregressive structure, so that the Jacobian is a triangular matrix and
computing its determinant is efficient. One problem of emerging convolution is
the computation of inverse is non-parallelizable, so that is inefficent for
high-dimensional variables.

Periodic convolutions [15, 9] transform the input and kernel to the Fourier
domain using discrete Fourier transformations, so the convolution function is
an element-wise matrix product with a block-diagonal Jacobian. The
computational cost of periodic convolutions is
$\mathcal{O}(chw\log(hw)+c^{3}hw)$. Thus, when the input is high-dimensional,
both training and sampling are expensive.

Multi-scale architectures [7] compose flow layers to generate rich models,
using _split layers_ to factor out variables and _squeeze layers_ to shuffle
dimensions, resulting in an architecture with $K$ flow steps and $L$ levels.
See Fig. 1.

##  3 Woodbury Transformations

In this section, we introduce Woodbury transformations as an efficient means
to model high-dimensional correlations.

###  3.1 Channel and Spatial Transformations

Suppose we reshape the input $\mathbf{x}$ to be a $c\times n$ matrix, where
$n=hw$. Then the 1$\times$1 convolution can be reinterpreted as a matrix
transformation

$\mathbf{y}=\mathbf{W}^{(c)}\mathbf{x},$ (2)

where $\mathbf{y}$ is also a $c\times n$ matrix, and $\mathbf{W}^{(c)}$ is a
$c\times c$ matrix. For consistency, we will call this a channel
transformation. For each column $\mathbf{x}_{:,i}$, the correlations among
channels are modeled by $\mathbf{W}^{(c)}$. However, the correlation between
any two rows $\mathbf{x}_{:,i}$ and $\mathbf{x}_{:,j}$ is not captured.
Inspired by Eq. 2, we use a spatial transformation to model interactions among
dimensions along the spatial axis

$\mathbf{y}=\mathbf{xW}^{(s)},$ (3)

where $\mathbf{W}^{(s)}$ is an $n\times n$ matrix that models the correlations
of each row $\mathbf{x}_{i,:}$. Combining Equation 2 and Equation 3, we have

$\displaystyle\mathbf{x}_{c}=\mathbf{W}^{(c)}\mathbf{x},~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{y}=\mathbf{x}_{c}\mathbf{W}^{(s)}.$
(4)

For each dimension of output $\mathbf{y}_{i,j}$, we have
$\mathbf{y}_{i,j}=\sum_{v=1}^{c}\left(\sum_{u=1}^{n}\mathbf{W}^{(c)}_{i,u}\cdot\mathbf{x}_{u,v}\right)\cdot\mathbf{W}^{(s)}_{v,j}$.

Therefore, the spatial and channel transformations together can model the
correlation between any pair of dimensions. However, in this preliminary form,
directly using Eq. 4 is inefficient for large $c$ or $n$. First, we would have
to store two large matrices $\mathbf{W}^{c}$ and $\mathbf{W}^{s}$, so the
space cost is $\mathcal{O}(c^{2}+n^{2})$. Second, the computational cost of
Eq. 4 is $\mathcal{O}(c^{2}n+n^{2}c)$—quadratic in the input size. Third, the
computational cost of the Jacobian determinant is $\mathcal{O}(c^{3}+n^{3})$,
which is far too expensive in practice.

###  3.2 Woodbury Transformations

We solve the three scalability problems by using a low-rank factorization.
Specifically, we define

$\displaystyle\mathbf{W}^{(c)}=\mathbf{I}^{(c)}+\mathbf{U}^{(c)}\mathbf{V}^{(c)},~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{W}^{(s)}=\mathbf{I}^{(s)}+\mathbf{U}^{(s)}\mathbf{V}^{(s)},$

where $\mathbf{I}^{(c)}$ and $\mathbf{I}^{(s)}$ are $c$\- and $n$-dimensional
identity matrices, respectively. The matrices $\mathbf{U}^{c}$,
$\mathbf{V}^{c}$, $\mathbf{U}^{s}$, and $\mathbf{V}^{s}$ are of size $c\times
d_{c}$, $d_{c}\times c$, $n\times d_{s}$, and $d_{c}\times n$, respectively,
where $d_{c}$ and $d_{s}$ are constant latent dimensions of these four
matrices. Therefore, we can rewrite Equation 4 as

$\displaystyle\mathbf{x}_{c}=(\mathbf{I}^{(c)}+\mathbf{U}^{(c)}\mathbf{V}^{(c)})\mathbf{x},~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{y}=\mathbf{x}_{c}(\mathbf{I}^{(s)}+\mathbf{U}^{(s)}\mathbf{V}^{(s)}).$
(5)

We call Eq. 5 the Woodbury transformation because the Woodbury matrix identity
[37] and Sylvester’s determinant identity [35] allow efficient computation of
its inverse and Jacobian determinant.

Woodbury matrix identity.111A more general version replaces $\mathbf{I}^{(n)}$
and $\mathbf{I}^{(k)}$ with arbitrary invertible $n\times n$ and $k\times k$
matrices. But this simplified version is sufficient for our tasks. Let
$\mathbf{I}^{(n)}$ and $\mathbf{I}^{(k)}$ be $n$\- and $k$-dimensional
identity matrices, respectively. Let $\mathbf{U}$ and $\mathbf{V}$ be $n\times
k$ and $k\times n$ matrices, respectively. If $\mathbf{I}^{(k)}+\mathbf{VU}$
is invertible, then
$(\mathbf{I}^{(n)}+\mathbf{UV})^{-1}=\mathbf{I}^{(n)}-\mathbf{U}(\mathbf{I}^{k}+\mathbf{VU})^{-1}\mathbf{V}$.

Sylvester’s determinant identity. Let $\mathbf{I}^{(n)}$ and
$\mathbf{I}^{(k)}$ be $n$\- and $k$-dimensional identity matrices,
respectively. Let $\mathbf{U}$ and $\mathbf{V}$ be $n\times k$ and $k\times n$
matrices, respectively. Then,
$\det(\mathbf{I}^{(n)}+\mathbf{UV})=\det(\mathbf{I}^{(k)}+\mathbf{VU})$.

Based on these two identities, we can efficiently compute the inverse and
Jacobian determinant

$\displaystyle\mathbf{x}_{c}=\mathbf{y}(\mathbf{I}^{(s)}-\mathbf{U}^{(s)}(\mathbf{I}^{(d_{s})}+\mathbf{V}^{(s)}\mathbf{U}^{(s)})^{-1}\mathbf{V}^{(s)}),$
$\displaystyle\mathbf{x}=(\mathbf{I}^{(c)}-\mathbf{U}^{(c)}(\mathbf{I}^{(d_{c})}+\mathbf{V}^{(c)}\mathbf{U}^{(c)})^{-1}\mathbf{V}^{(c)})\mathbf{x}_{c},$
(6)

and

$\displaystyle\log\left|\det\left(\frac{\partial\mathbf{y}}{\partial\mathbf{x}}\right)\right|$
$\displaystyle=n\log\left|\det(\mathbf{I}^{(d_{c})}+\mathbf{V}^{(c)}\mathbf{U}^{(c)})\right|+c\log\left|\det(\mathbf{I}^{(d_{s})}+\mathbf{V}^{(s)}\mathbf{U}^{(s)})\right|,$
(7)

where $\mathbf{I}^{(d_{c})}$ and $\mathbf{I}^{(d_{s})}$ are $d_{c}$\- and
$d_{s}$-dimensional identity matrices, respectively.

A Woodbury transformation is also a generalized permutation layer. We can
directly replace an invertible convolution in Figure 1 with a Woodbury
transformation. In contrast with 1$\times$1 convolutions, Woodbury
transformations are able to model correlations along both channel and spatial
axes. We illustrate this in Figure 2. To implement Woodbury transformations,
we need to store four weight matrices, i.e.,
$\mathbf{U}^{(c)},\mathbf{U}^{(s)},\mathbf{V}^{(c)}$, and $\mathbf{V}^{(s)}$.
To simplify our analysis, let $d_{c}\leq d$ and $d_{s}\leq d$, where $d$ is a
constant. This setting is also consistent with our experiments. The size of
$\mathbf{U}^{(c)}$ and $\mathbf{V}^{(c)}$ is $\mathcal{O}(dc)$, and the size
of $\mathbf{U}^{(c)}$ and $\mathbf{V}^{(c)}$ is $\mathcal{O}(dn)$. The space
complexity is $\mathcal{O}(d(c+n))$.

For training and likelihood computation, the main computational bottleneck is
computing $\mathbf{y}$ and the Jacobian determinant. To compute $\mathbf{y}$
with Equation 4, we need to first compute the channel transformation and then
compute the spatial transformation. The computational complexity is
$\mathcal{O}(dcn)$. To compute the determinant with Equation 7, we need to
first compute the matrix product of $\mathbf{V}$ and $\mathbf{U}$, and then
compute the determinant. The computational complexity is
$\mathcal{O}(d^{2}(c+n)+d^{3})$.

For sampling, we need to compute the inverse transformations, i.e., Equation
6. With the Woodbury identity, we actually only need to compute the inverses
of $\mathbf{I}^{(d_{s})}+\mathbf{V}^{(s)}\mathbf{U}^{(s)}$ and
$\mathbf{I}^{(d_{c})}+\mathbf{V}^{(c)}\mathbf{U}^{(c)}$, which are computed
with time complexity $\mathcal{O}(d^{3})$. To implement the inverse
transformations, we can compute the matrix chain multiplication, so we can
avoid computing the product of two large matrices twice, yielding cost
$\mathcal{O}(c^{2}+n^{2})$. For example, for the inverse spatial
transformation, we can compute it as
$\mathbf{x}_{c}=\mathbf{y}-((\mathbf{yU}^{(s)})(\mathbf{I}^{(d_{s})}+\mathbf{V}^{(s)}\mathbf{U}^{(s)})^{-1})\mathbf{V}^{(s)}$,
so that its complexity is $\mathcal{O}(d^{3}+cd^{2}+cnd)$. The total
computational complexity of Equation 6 is $\mathcal{O}(dcn+d^{2}(n+c)+d^{3})$.

In practice, we found that for a high-dimensional input, a relatively small
$d$ is enough to obtain good performance, e.g., the input is $256\times
256\times 3$ images, and $d=16$. In this situation, $nc\geq d^{3}$. Therefore,
we can omit $d$ and approximately see the spatial complexity as
$\mathcal{O}(c+n)$, and the forward or inverse transformation as
$\mathcal{O}(nc)$. They are all linear to the input size.

We do not restrict $\mathbf{U}$ and $\mathbf{V}$ to force $\mathbf{W}$ to be
invertible. Based on analysis by Hoogeboom et al. [15], the training maximizes
the log-likelihood, which implicitly pushes $\det(\mathbf{I}+\mathbf{VU})$
away from $0$. Therefore, it is not necessary to explicitly force
invertibility. In our experiments, the Woodbury transformations are as robust
as other invertible convolution layers.

(a) 1$\times$1 convolution

(b) Woodbury

(c) ME-Woodbury

Figure 2: Visualization of three transformations. The 1$\times$1 convolution
only operates along the channel axis. The Woodbury transformation operates
along both the channel and spatial axes, modeling the dependencies of one
channel directly via one transformation. The ME-Woodbury transformation
operates along three axes. It uses two transformations to model spatial
dependencies.

###  3.3 Memory-Efficient Variant

In Eq. 4, one potential challenge arises from the sizes of $\mathbf{U}^{(s)}$
and $\mathbf{V}^{((s))}$, which are linear in $n$. The challenge is that $n$
may be large in some practical problems, e.g., high-resolution images. We
develop a memory-efficient variant of Woodbury transformations, i.e., ME-
Woodbury, to solve this problem. The ME version can effectively reduce space
complexity from $\mathcal{O}(d(c+hw))$ to $\mathcal{O}(d(c+h+w))$.

The difference between ME-Woodbury transformations and Woodbury
transformations is that the ME form cannot directly model spatial
correlations. As shown in Figure 2(c), it uses two transformations, for height
and width, together to model the spatial correlations. Therefore, for a
specific channel $k$, when two dimensions $\mathbf{x}_{k,i,j}$ and
$\mathbf{x}_{k,u,v}$ are in two different heights, and widths, their
interaction will be modeled indirectly. In our experiments, we found that this
limitation only slightly impacts ME-Woodbury’s performance. More details on
ME-Woodbury transformations are in the appendix.

##  4 Related Work

Rezende and Mohamed [30] developed planar flows for variational inference
$\mathbf{z}_{t+1}=\mathbf{z}_{t}+\mathbf{u}\delta(\mathbf{w}^{T}\mathbf{z}_{t}+b)$,
where $\mathbf{z}$, $\mathbf{w}$, and $\mathbf{u}$ are $d$-dimensional
vectors, $\delta()$ is an activation function, and $b$ is a scalar.

Berg et al. [3] generalized these to Sylvester flows
$\mathbf{z}_{t+1}=\mathbf{z}_{t}+\mathbf{QR}\delta(\tilde{\mathbf{R}}\mathbf{Q}^{T}\mathbf{z}_{t}+\mathbf{r})$,
where $\mathbf{R}$ and $\tilde{\mathbf{R}}$ are upper triangular matrices,
$\mathbf{Q}$ is composed of a set of orthonormal vectors, and $\mathbf{r}$ is
a $d$-dimensional vector. The resulting Jacobian determinant can be
efficiently computed via Sylvester’s identity, just as our methods do.
However, Woodbury transformations have key differences from Sylvester flows.
First, Berg et al. only analyze their models on vectors. The inputs to our
layers are matrices, so our method operates on high-dimensional input, e.g.,
images. Second, though Sylvester flows are inverse functions, computing their
inverse is difficult. One possible way is to apply iterative methods [2, 34,
5] to compute the inverse. But this research direction is unexplored. Our
layers can be inverted efficiently with the Woodbury identity. Third, our
layers do not restrict the transformation matrices to be triangular or
orthogonal. In fact, Woodbury transformations can be seen as another
generalized variant of planar flows on matrices, with
$\delta(\mathbf{x})=\mathbf{x}$, and whose inverse is tractable. Roughly
speaking, Woodbury transformations can also be viewed as applying the planar
flows sequentially to each row of the input matrix. After this work was
completed and submitted, we learned that the TensorFlow software [1] also uses
the Woodbury identity in their affine bijector.

Normalizing flows have also been used for variational inference, density
estimation, and generative modeling. Autoregressive flows [23, 29, 17, 25]
restrict each variable to depend on those that precede it in a sequence,
forcing a triangular Jacobian. Non-linear coupling layers replace the affine
transformation function. Specifically, spline flows [26, 8] use spline
interpolation, and Flow++ [14] uses a mixture cumulative distribution function
to define these functions. Flow++ also uses variational dequantization to
prevent model collapse. Many works [22, 15, 9, 18] develop invertible
convolutional flows to model interactions among dimensions. MintNet [34] is a
flexible architecture composed of multiple masked invertible layers. I-ResNet
[2, 5] uses discriminative deep network architecture as the flow. These two
models require iterative methods to compute the inverse. Discrete flows [36,
16] and latent flows [41] can be applied to discrete data such as text.
Continuous-time flows [4, 12] have been developed based on the theory of
ordinary differential equations.

##  5 Experiments

In this section, we compare the performance of Woodbury transformations
against other modern flow architectures, measuring running time, bit per-
dimension ($\log_{2}$-likelihood), and sample quality.

Figure 3: Running time comparison. Sampling with emerging convolutions is
slow, since their inverses are not parallelizable. Periodic convolutions are
costly for larger inputs. Both 1$\times$1 convolutions and Woodbury
transformations are efficient in training and sampling.

Running Time We follow Finz et al. [9] and compare the per-sample running time
of Woodbury transformations to other generalized permutations: 1$\times$1
[22], emerging [15], and periodic convolutions [15, 9]. We test the training
time and sampling time. In training, we compute (1) forward propagation, i.e.,
$\mathbf{y}=\mathbf{f}(\mathbf{x})$, of a given function $\mathbf{f}()$, (2)
the Jacobian determinant, i.e.,
$\det\left(\left|\frac{\partial\mathbf{y}}{\partial\mathbf{x}}\right|\right)$,
and (3) the gradient of parameters. For sampling, we compute the inverse of
transformation $\mathbf{x}=\mathbf{f}^{-1}(\mathbf{y})$. For emerging and
periodic convolutions, we use $3\times 3$ kernels. For Woodbury
transformations, we fix the latent dimension $d=16$. For fair comparison, we
implement all methods in Pytorch and run them on an Nvidia Titan V GPU. We
follow Hoogeboom et al. [15] and implement the emerging convolution inverse in
Cython, and we compute it on a 4 Ghz CPU (the GPU version is slower than the
Cython version). We first fix the spatial size to be $64\times 64$ and vary
the channel number. We then fix the channel number to be $96$ and vary the
spatial size.

The results are shown in Figure 3. For training, the emerging convolution is
the fastest. This is because its Jacobian is a triangular matrix, so computing
its determinant is much more efficient than other methods. The Woodbury
transformation and ME-Woodbury are slightly slower than the 1x1 convolution,
since they contain more transformations. Emerging convolutions, Woodbury
transformations, and 1x1 convolutions only slightly increase with input size,
rather than increasing with $\mathcal{O}(c^{3})$. This invariance to input
size is likely because of how the GPU parallelizes computation. The periodic
convolution is efficient only when the input size is small. When the size is
large, it becomes slow, e.g., when the input size is $96\times 64\times 64$,
it is around $30$ times slower than Woodbury transformations. In our
experiments, we found that the Fourier transformation requires a large amount
of memory. According to Finz et al. [9], the Fourier step may be the
bottleneck that impacts periodic convolution’s scalability. A more efficient
implementation of Fourier transformation, e.g., [18], may improve its running
time.

For sampling, both 1$\times$1 convolutions and Woodbury transformations are
efficient. The 1$\times$1 convolution is the fastest, and the Woodbury
transformations are only slightly slower. Neither is sensitive to the change
of input size. Emerging convolutions and periodic convolutions are much slower
than Woodbury transformations, and their running time increases with the input
size. When the input size is $96\times 128\times 128$, they are around $100$
to $200$ times slower than Woodbury transformations. This difference is
because emerging convolutions cannot make use of parallelization, and periodic
transformations require conversion to Fourier form. Based on these results, we
can conclude that both emerging convolution and periodic convolution do not
scale well to high-dimensional inputs. In contrast, Woodbury transformations
are efficient in both training and sampling.

Quantitative Evaluation We compare Woodbury transformations with state-of-the-
art flow models, measuring bit per-dimension (bpd). We train with the CIFAR-10
[24] and ImageNet [32] datasets. We compare with three generalized permutation
methods—1$\times$1 convolution, emerging convolution, and periodic
convolution—and two coupling layers—neural spline coupling [8] and MaCow [25].
We use Glow (Fig. 1, [22]) as the basic flow architecture. For each method, we
replace the corresponding layer. For example, to construct a flow with
Woodbury transformations, we replace the 1$\times$1 convolution with a
Woodbury transformation, i.e., Eq. 4. For all generalized permutation methods,
we use affine coupling. For each of the coupling layer baselines, we
substitute it for the affine coupling. We tune the parameters of neural spline
coupling and MaCow so that their sizes are close to affine coupling. We follow
Hoogeboom et al. [15] and test the performance of small models. For $32\times
32$ images, we set the number of levels to $L=3$ and the number of steps per-
level to $K=8$. For $64\times 64$ images, we use $L=4$ and $K=16$. More
details are in the appendix.

Table 1: Quantitative evaluation results. | Quantitative measure (bpd) | Model sizes (# parameters)
---|---|---
| CIFAR-10 | ImageNet | ImageNet | 32x32 images | 64x64 images
| 32x32 | 32x32 | 64x64 |  |
1$\times$1 convolution | 3.51 | 4.32 | 3.94 | 11.02M | 37.04M
Emerging | 3.48 | 4.26 | 3.91 | 11.43M | 40.37M
Periodic | 3.49 | 4.28 | 3.92 | 11.21M | 38.61M
Neural spline | 3.50 | 4.24 | 3.95 | 10.91M | 38.31M
MaCow | 3.48 | 4.34 | 4.15 | 11.43M | 37.83M
ME-Woodbury | 3.48 | 4.22 | 3.91 | 11.02M | 36.98M
Woodbury | 3.47 | 4.20 | 3.87 | 11.10M | 37.60M

Woodbury-Glow
Iteration 50,000

Woodbury-Glow
Iteration 100,000

Woodbury-Glow
Iteration 600,000

Glow
Iteration 50,000

Glow
Iteration 100,000

Glow
Iteration 600,000

Figure 4: Random samples $64\times 64$ drawn from models trained on CelebA
with temperature $0.7$.

The test-set likelihoods are listed in Table 1 left. Our scores are worse than
those reported by Kingma and Dhariwal [22], Hoogeboom et al. [15] because we
use smaller models and train each model on a single GPU. Based on the scores,
1$\times$1 convolutions perform the worst. Emerging convolutions and periodic
convolutions score better than the 1$\times$1 convolutions, since they are
more flexible and can model the dependencies along the spatial axes. Neural
spline coupling works well on $32\times 32$ images, but do slightly worse than
1$\times$1 convolution on $64\times 64$ images. MaCow does not work well on
ImageNet. This trend demonstrates the importance of permutation layers. They
can model the interactions among dimensions and shuffle them, which coupling
layers cannot do. Without a good permutation layer, a better coupling layer
still cannot always improve the performance. The Woodbury transformation
models perform the best, likely because they can model the interactions
between the target dimension and all other dimensions, while the invertible
convolutions only model the interactions between target dimension its
neighbors. ME-Woodbury performs only slightly worse than the full version,
showing that its restrictions provide a useful tradeoff between model quality
and efficiency.

We list model sizes in Table 1 (right). Despite modeling rich interactions,
Woodbury transformations are not the largest. With $32\times 32$ images, ME-
Woodbury and 1$\times$1 convolution are the same size. When the image size is
$64\times 64$, ME-Woodbury is the smallest. This is because we use the multi-
scale architecture, i.e., Fig. 1, to combine layers. The squeeze layer doubles
the input variable’s channels at each level, so larger $L$ suggests larger
$c$. The space complexities of invertible convolutions are
$\mathcal{O}(c^{2})$, while the space complexity of ME-Woodbury is linear to
$c$. When $c$ is large, the weight matrices of invertible convolutions are
larger than the weight matrices of ME-Woodbury.

Table 2: Evaluation of different $d$ (bpd). | Woodbury | ME-Woodbury
---|---|---
$d=2$ | 3.54 | 3.53
$d=4$ | 3.51 | 3.51
$d=8$ | 3.48 | 3.48
$d=16$ | 3.47 | 3.48
$d=32$ | 3.47 | 3.48

Latent Dimension Evaluation We test the impact of latent dimension $d$ on the
performance of Woodbury-Glow. We train our models on CIFAR-10, and use bpd as
metric. We vary $d$ within $\\{2,4,8,16,32\\}$. The results are in Table 2.
When $d<8$, the model performance will be impacted. When $d>16$, increasing
$d$ will not improve the bpd. This is probably because when $d$ is too small,
the latent features cannot represent the input variables well, and when $d$ is
too big, the models become hard to train. When $8\leq d\leq 16$, the Woodbury
transformations are powerful enough to model the interactions among
dimensions. We also test two values of $d$, i.e., $16,32$, of Woodbury-Glow on
ImageNet $64\times 64$. The bpds of both $d$ are $3.87$, which are consistent
with our conclusion.

Figure 5: Learning curves on CelebA-HQ 64x64. The NLL of Woodbury Glow
decreases faster than Glow.

Sample Quality Comparisons We train Glow and Woodbury-Glow on the CelebA-HQ
dataset [19]. We use $5$-bit images and set the size of images to be $64\times
64$, $128\times 128$, and $256\times 256$. Due to our limited computing
resources, we use relatively small models in our experiments. We follow Kingma
and Dhariwal [22] and choose a temperature parameter to encourage higher
quality samples. Detailed parameter settings are in the appendix. We compare
samples from Glow and Woodbury-Glow during three phases of training, displayed
in Fig. 4. The samples show a clear trend where Woodbury-Glow more quickly
learns to generate reasonable face shapes. After 100,000 iterations, it can
already generate reasonable samples, while Glow’s samples are heavily
distorted. Woodbury-Glow samples are consistently smoother and more realistic
than samples from Glow in all phases of training. The samples demonstrate
Woodbury transformations’ advantages. The learning curves in Figure 5 also
show that the NLL of Woodbury Glow decreases faster, which is consistent to
the sample comparisons. In the appendix, we show analogous comparisons using
higher resolution versions of CelebA data, which also exhibit the trend of
Woodbury-Glow generating more realistic images than Glow at the same training
iterations.

##  6 Conclusion

In this paper, we develop Woodbury transformations, which use the Woodbury
matrix identity to compute the inverse transformations and Sylvester’s
determinant identity to compute Jacobian determinants. Our method has the same
advantages as invertible $d\times d$ convolutions that can capture
correlations among all dimensions. In contrast to the invertible $d\times d$
convolutions, our method is parallelizable and the computational complexity of
our methods are linear to the input size, so that it is still efficient in
computation when the input is high-dimensional. One potential limitation is
that Woodbury transformations do not have parameter sharing scheme as in
convolutional layers, so one potential future research is to develop partially
Woodbury transformations that can share parameters. We test our models on
multiple image datasets and they outperform state-of-the-art methods.

## Broader Impact

This paper presents fundamental research on increasing the expressiveness of
deep probabilistic models. Its impact is therefore linked to the various
applications of such models. By enriching the class of complex deep models for
which we can train with exact likelihood, we may enable a wide variety of
applications that can benefit from modeling of uncertainty. However, a
potential danger of this research is that deep generative models have been
recently applied to synthesize realistic images and text, which can be used
for misinformation campaigns.

## Acknowledgments

We thank NVIDIA’s GPU Grant Program and Amazon’s AWS Cloud Credits for
Research program for their support. The work was completed while both authors
were affiliated with the Virginia Tech Department of Computer Science. Bert
Huang was partially supported by an Amazon Research Award and a grant from the
U.S. Department of Transportation Safe-D Program for work on separate projects
not directly related to this paper.

## 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 _12th USENIX Symposium on Operating Systems Design and Implementation_ , pages 265–283, 2016.
  * Behrmann et al. [2018] Jens Behrmann, Will Grathwohl, Ricky TQ Chen, David Duvenaud, and Jörn-Henrik Jacobsen.  Invertible residual networks.  _arXiv preprint arXiv:1811.00995_ , 2018.
  * Berg et al. [2018] Rianne van den Berg, Leonard Hasenclever, Jakub M Tomczak, and Max Welling.  Sylvester normalizing flows for variational inference.  _arXiv preprint arXiv:1803.05649_ , 2018.
  * Chen et al. [2018] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.  Neural ordinary differential equations.  In _Advances in Neural Information Processing Systems_ , pages 6571–6583, 2018.
  * Chen et al. [2019] Tian Qi Chen, Jens Behrmann, David K Duvenaud, and Jörn-Henrik Jacobsen.  Residual flows for invertible generative modeling.  In _Advances in Neural Information Processing Systems_ , pages 9913–9923, 2019.
  * Dinh et al. [2014] Laurent Dinh, David Krueger, and Yoshua Bengio.  NICE: Non-linear independent components estimation.  _arXiv preprint arXiv:1410.8516_ , 2014.
  * Dinh et al. [2016] Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio.  Density estimation using real NVP.  _arXiv preprint arXiv:1605.08803_ , 2016.
  * Durkan et al. [2019] Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios.  Neural spline flows.  _arXiv preprint arXiv:1906.04032_ , 2019.
  * Finz et al. [2019] Marc Finz, Pavel Izmailov, Wesley Maddox, Polina Kirichenko, and Andrew Gordon Wilson.  Invertible convolutional networks.  In _ICML Workshop on Invertible Neural Networks and Normalizing Flows_ , 2019.
  * Germain et al. [2015] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle.  Made: Masked autoencoder for distribution estimation.  In _International Conference on Machine Learning_ , pages 881–889, 2015.
  * Goodfellow et al. [2014] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  In _Advances in Neural Information Processing Systems_ , pages 2672–2680, 2014.
  * Grathwohl et al. [2018] Will Grathwohl, Ricky TQ Chen, Jesse Betterncourt, Ilya Sutskever, and David Duvenaud.  Ffjord: Free-form continuous dynamics for scalable reversible generative models.  _arXiv preprint arXiv:1810.01367_ , 2018.
  * Graves [2013] Alex Graves.  Generating sequences with recurrent neural networks.  _arXiv preprint arXiv:1308.0850_ , 2013.
  * Ho et al. [2019] Jonathan Ho, Xi Chen, Aravind Srinivas, Yan Duan, and Pieter Abbeel.  Flow++: Improving flow-based generative models with variational dequantization and architecture design.  _arXiv preprint arXiv:1902.00275_ , 2019.
  * Hoogeboom et al. [2019a] Emiel Hoogeboom, Rianne van den Berg, and Max Welling.  Emerging convolutions for generative normalizing flows.  _arXiv preprint arXiv:1901.11137_ , 2019a.
  * Hoogeboom et al. [2019b] Emiel Hoogeboom, Jorn WT Peters, Rianne van den Berg, and Max Welling.  Integer discrete flows and lossless compression.  _arXiv preprint arXiv:1905.07376_ , 2019b.
  * Huang et al. [2018] Chin-Wei Huang, David Krueger, Alexandre Lacoste, and Aaron Courville.  Neural autoregressive flows.  _arXiv preprint arXiv:1804.00779_ , 2018.
  * Karami et al. [2019] Mahdi Karami, Dale Schuurmans, Jascha Sohl-Dickstein, Laurent Dinh, and Daniel Duckworth.  Invertible convolutional flow.  In _Advances in Neural Information Processing Systems_ , pages 5636–5646, 2019.
  * Karras et al. [2017] Tero Karras, Timo Aila, Samuli Laine, and Jaakko Lehtinen.  Progressive growing of gans for improved quality, stability, and variation.  _arXiv preprint arXiv:1710.10196_ , 2017.
  * Kingma and Ba [2014] Diederik Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational Bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Kingma and Dhariwal [2018] Durk P Kingma and Prafulla Dhariwal.  Glow: Generative flow with invertible 1x1 convolutions.  In _Advances in Neural Information Processing Systems_ , pages 10215–10224, 2018.
  * Kingma et al. [2016] Durk P Kingma, Tim Salimans, Rafal Jozefowicz, Xi Chen, Ilya Sutskever, and Max Welling.  Improved variational inference with inverse autoregressive flow.  In _Advances in Neural Information Processing Systems_ , pages 4743–4751, 2016.
  * Krizhevsky et al. [2009] Alex Krizhevsky, Geoffrey Hinton, et al.  Learning multiple layers of features from tiny images.  Technical report, Citeseer, 2009.
  * Ma et al. [2019] Xuezhe Ma, Xiang Kong, Shanghang Zhang, and Eduard Hovy.  Macow: Masked convolutional generative flow.  In _Advances in Neural Information Processing Systems_ , pages 5891–5900. 2019.
  * Müller et al. [2019] Thomas Müller, Brian McWilliams, Fabrice Rousselle, Markus Gross, and Jan Novák.  Neural importance sampling.  _ACM Transactions on Graphics (TOG)_ , 38(5):1–19, 2019.
  * Oord et al. [2016a] Aaron van den Oord, Sander Dieleman, Heiga Zen, Karen Simonyan, Oriol Vinyals, Alex Graves, Nal Kalchbrenner, Andrew Senior, and Koray Kavukcuoglu.  Wavenet: A generative model for raw audio.  _arXiv preprint arXiv:1609.03499_ , 2016a.
  * Oord et al. [2016b] Aaron van den Oord, Nal Kalchbrenner, and Koray Kavukcuoglu.  Pixel recurrent neural networks.  _arXiv preprint arXiv:1601.06759_ , 2016b.
  * Papamakarios et al. [2017] George Papamakarios, Theo Pavlakou, and Iain Murray.  Masked autoregressive flow for density estimation.  In _Advances in Neural Information Processing Systems_ , pages 2338–2347, 2017.
  * Rezende and Mohamed [2015] Danilo Jimenez Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  _arXiv preprint arXiv:1505.05770_ , 2015.
  * Rezende et al. [2014] Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic backpropagation and approximate inference in deep generative models.  _arXiv preprint arXiv:1401.4082_ , 2014.
  * Russakovsky et al. [2015] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  _International Journal of Computer Vision_ , 115(3):211–252, 2015.
  * Sohn et al. [2015] Kihyuk Sohn, Honglak Lee, and Xinchen Yan.  Learning structured output representation using deep conditional generative models.  In _Advances in Neural Information Processing Systems_ , pages 3483–3491, 2015.
  * Song et al. [2019] Yang Song, Chenlin Meng, and Stefano Ermon.  Mintnet: Building invertible neural networks with masked convolutions.  In _Advances in Neural Information Processing Systems_ , pages 11002–11012, 2019.
  * Sylvester [1851] James Joseph Sylvester.  On the relation between the minor determinants of linearly equivalent quadratic functions.  _The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science_ , 1(4):295–305, 1851.
  * Tran et al. [2019] Dustin Tran, Keyon Vafa, Kumar Krishna Agrawal, Laurent Dinh, and Ben Poole.  Discrete flows: Invertible generative models of discrete data.  _arXiv preprint arXiv:1905.10347_ , 2019.
  * Woodbury [1950] Max A Woodbury.  Inverting modified matrices.  1950\.
  * Yu et al. [2015] Fisher Yu, Yinda Zhang, Shuran Song, Ari Seff, and Jianxiong Xiao.  Lsun: Construction of a large-scale image dataset using deep learning with humans in the loop.  _arXiv preprint arXiv:1506.03365_ , 2015.
  * Yu et al. [2017] Lantao Yu, Weinan Zhang, Jun Wang, and Yong Yu.  Seqgan: Sequence generative adversarial nets with policy gradient.  In _Thirty-First AAAI Conference on Artificial Intelligence_ , 2017\.
  * Zhu et al. [2017] Jun-Yan Zhu, Taesung Park, Phillip Isola, and Alexei A Efros.  Unpaired image-to-image translation using cycle-consistent adversarial networks.  In _Proceedings of International Conference on Computer Vision_ , pages 2223–2232, 2017.
  * Ziegler and Rush [2019] Zachary M Ziegler and Alexander M Rush.  Latent normalizing flows for discrete sequences.  _arXiv preprint arXiv:1901.10548_ , 2019.

##  Appendix A More Background

In this section, we introduce more detailed background knowledge.

###  A.1 Normalizing Flows

Let $\mathbf{x}$ be a high-dimensional continuous variable. We suppose that
$\mathbf{x}$ is drawn from $p^{*}(\mathbf{x})$, which is the true data
distribution. Given a collected dataset
$\mathcal{D}=\\{\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{D}\\}$, we are
interested in approximating $p^{*}(\mathbf{x})$ with a model
$p_{\theta}(\mathbf{x})$. We optimize $\theta$ by minimizing the negative log-
likelihood

$\mathcal{L}(\mathcal{D})=\sum_{i=1}^{D}-\log p_{\theta}(\mathbf{x}_{i}).$ (8)

For some settings, variable $\tilde{\mathbf{x}}$ is discrete, e.g., image
pixel values are often integers. In these cases, we dequantize
$\tilde{\mathbf{x}}$ by adding continuous noise $\bm{\mu}$ to it, resulting in
a continuous variable $\mathbf{x}=\tilde{\mathbf{x}}+\bm{\mu}$. As shown by Ho
et al. [14], the log-likelihood of $\tilde{\textbf{x}}$ is lower-bounded by
the log-likelihood of $\mathbf{x}$.

Normalizing flows enable computation of $p_{\theta}(\mathbf{x})$, even though
it is usually intractable for many other model families. A normalizing flow
[30] is composed of a series of invertible functions
$\mathbf{f}=\mathbf{f}_{1}\circ\mathbf{f}_{2}\circ...\circ\mathbf{f}_{K}$,
which transform $\mathbf{x}$ to a latent code $\mathbf{z}$ drawn from a simple
distribution. Therefore, with the _change of variables_ formula, we can
rewrite $\log p_{\theta}(\mathbf{x})$ to be

$\log p_{\theta}(\mathbf{x})=\log
p_{Z}(\mathbf{z})+\sum_{i=1}^{K}\log\left|\det\left(\frac{\partial\mathbf{f}_{i}}{\partial\mathbf{r}_{i-1}}\right)\right|,$
(9)

where $\mathbf{r}_{i}=\mathbf{f}_{i}(\mathbf{r}_{i-1})$,
$\mathbf{r}_{0}=\mathbf{x}$, and $\mathbf{r}_{K}=\mathbf{z}$.

###  A.2 Invertible $d\times d$ Convolutions

Emerging convolutions [15] combine two autoregressive convolutions [10, 23].
Formally,

$\displaystyle\mathbf{M}^{\prime}_{1}=\mathbf{M}_{1}\odot\mathbf{A}_{1},~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{M}^{\prime}_{2}=\mathbf{M}_{2}\odot\mathbf{A}_{2},~{}~{}~{}~{}~{}~{}~{}~{}\mathbf{y}=\mathbf{M}^{\prime}_{2}\star(\mathbf{M}^{\prime}_{1}\star\mathbf{x}),$

where $\mathbf{M}_{1},\mathbf{M}_{2}$ are convolutional kernels whose size is
$c\times c\times d\times d$, and $\mathbf{A}_{1},\mathbf{A}_{2}$ are binary
masks. The symbol $\star$ represents the convolution operator.222In practice,
a convolutional layer is usually implemented as an aggregation of cross-
correlations. We follow Hoogeboom et al. [15] and omit this detail. An
emerging convolutional layer has the same receptive fields as standard
convolutional layers, which can capture correlations between a target pixel
and its neighbor pixels. However, like other autoregressive convolutions,
computing the inverse of an emerging convolution requires sequentially
traversing each dimension of input, so its computation is not parallelizable
and is a computational bottleneck when the input is high-dimensional.

Periodic convolutions [15, 9] use discrete Fourier transformations to
transform both the input and the kernel to Fourier domain. A periodic
convolution is computed as

$\mathbf{y}_{u,:,:}=\sum_{v}\mathcal{F}^{-1}(\mathcal{F}(\mathbf{M}^{(p)}_{u,v,:,:})\odot\mathcal{F}(\mathbf{x}_{v,:,:})),$

where $\mathcal{F}$ is a discrete Fourier transformation, and
$\mathbf{M}^{(p)}$ is the convolution kernel whose size is $c\times c\times
d\times d$. The computational complexity of periodic convolutions is
$\mathcal{O}(c^{2}hw\log(hw)+c^{3}hw)$. In our experiments, we found that the
Fourier transformation requires a large amount of memory. These two problems
impact the efficiency of both training and sampling when the input is high-
dimensional.

##  Appendix B Memory-Efficient Woodbury Transformations

Memory-Efficient Woodbury transformations can effectively reduce the space
complexity. The main idea is to perform spatial transformations along the
height and width axes separately, i.e., a height transformation and a width
transformation. The transformations are:

$\displaystyle\mathbf{x}_{c}$ $\displaystyle=$
$\displaystyle(\mathbf{I}^{(c)}+\mathbf{U}^{(c)}\mathbf{V}^{(c)})\mathbf{x},$
$\displaystyle\mathbf{x}_{w}$ $\displaystyle=$
$\displaystyle\text{reshape}(\mathbf{x}_{c},(ch,w)),$
$\displaystyle\mathbf{x}_{w}$ $\displaystyle=$
$\displaystyle\mathbf{x}_{c}(\mathbf{I}^{(w)}+\mathbf{U}^{(w)}\mathbf{V}^{(w)}),$
$\displaystyle\mathbf{x}_{h}$ $\displaystyle=$
$\displaystyle\text{reshape}(\mathbf{x}_{w},(cw,h)),$
$\displaystyle\mathbf{y}$ $\displaystyle=$
$\displaystyle\mathbf{x}_{h}(\mathbf{I}^{(h)}+\mathbf{U}^{(h)}\mathbf{V}^{(h)}),$
$\displaystyle\mathbf{y}$ $\displaystyle=$
$\displaystyle\text{reshape}(\mathbf{y},(c,hw)),$ (10)

where $\text{reshape}(\mathbf{x},(n,m))$ reshapes $\mathbf{x}$ to be an
$n\times m$ matrix. Matrices $\mathbf{I}^{(w)}$ and $\mathbf{I}^{(h)}$ are
$w$\- and $h$-dimensional identity matrices, respectively. Matrices
$\mathbf{U}^{(w)},\mathbf{V}^{(w)},\mathbf{U}^{(h)}$, and $\mathbf{V}^{(h)}$
are $w\times d_{w}$, $d_{w}\times w$, $w\times d_{w}$, and $d_{w}\times w$
matrices, respectively, where $d_{w}$ and $d_{h}$ are constant latent
dimensions.

Using the Woodbury matrix identity and the Sylvester’s determinant identity,
we can compute the inverse and Jacobian determinant:

$\displaystyle\mathbf{y}$ $\displaystyle=$
$\displaystyle\text{reshape}(\mathbf{y},(cw,h)),$
$\displaystyle\mathbf{x}_{h}$ $\displaystyle=$
$\displaystyle\mathbf{y}(\mathbf{I}^{(h)}-\mathbf{U}^{(h)}(\mathbf{I}^{(d_{h})}+\mathbf{V}^{(h)}\mathbf{U}^{(h)})^{-1}\mathbf{V}^{(h)}),$
$\displaystyle\mathbf{x}_{w}$ $\displaystyle=$
$\displaystyle\text{reshape}(\mathbf{x}_{h},(ch,w)),$
$\displaystyle\mathbf{x}_{w}$ $\displaystyle=$
$\displaystyle\mathbf{x}_{w}(\mathbf{I}^{(w)}-\mathbf{U}^{(w)}(\mathbf{I}^{(d_{w})}+\mathbf{V}^{(w)}\mathbf{U}^{(w)})^{-1}\mathbf{V}^{(w)}),$
$\displaystyle\mathbf{x}_{c}$ $\displaystyle=$
$\displaystyle\text{reshape}(\mathbf{x}_{w},(c,hw)),$
$\displaystyle\mathbf{x}$ $\displaystyle=$
$\displaystyle(\mathbf{I}^{(c)}-\mathbf{U}^{(c)}(\mathbf{I}^{(d_{c})}+\mathbf{V}^{(c)}\mathbf{U}^{(c)})^{-1}\mathbf{V}^{(c)})\mathbf{x}_{c},$
(11)
$\displaystyle\log\left|\det(\frac{\partial\mathbf{y}}{\partial\mathbf{x}})\right|$
$\displaystyle=hw\log\left|\det(\mathbf{I}^{(d_{c})}+\mathbf{V}^{(c)}\mathbf{U}^{(c)})\right|+ch\log\left|\det(\mathbf{I}^{(d_{w})}+\mathbf{V}^{(w)}\mathbf{U}^{(w)})\right|$
(12)
$\displaystyle+cw\log\left|\det\left(\mathbf{I}^{(d_{h})}+\mathbf{V}^{(h)}\mathbf{U}^{(h)}\right)\right|,$

where $\mathbf{I}^{(d_{w})}$ and $\mathbf{I}^{(d_{h})}$ are $d_{w}$\- and
$d_{h}$-dimensional identity matrices, respectively. The Jacobian of the
$\text{reshape}()$ is an identity matrix, so its log-determinant is $0$.

We call Equation 10 the memory-efficient Woodbury transformation because it
reduces space complexity from $\mathcal{O}(c+hw)$ to $\mathcal{O}(c+h+w)$.
This method is effective when $h$ and $w$ are large. To analyze its
complexity, we let all latent dimensions be less than $d$ as before. The
complexity of forward transformation is $\mathcal{O}(dchw)$; the complexity of
computing the determinant is $\mathcal{O}(d(c+h+w)+d^{3})$; and the complexity
of computing the inverse is $\mathcal{O}(dchw+d^{2}(c+ch+cw)+d^{3})$. The same
as Woodbury transformations, when the input is high dimensional, we can omit
$d$. Therefore, the computational complexities of the memory-efficient
Woodbury transformation are also linear with the input size.

We list the complexities of different methods in Table 3. We can see that the
computational complexities of Woodbury transformations are comparable to other
methods, and maybe smaller when the input is high-dimensional, i.e., the
$c,h,w$ are big.

Table 3: Comparisons of computational complexities.  Method | Forward | Backward
---|---|---
1x1 convolution | $\mathcal{O}(c^{2}hw+c^{3})$ | $\mathcal{O}(c^{2}hw)$
Periodic conolution | $\mathcal{O}(chw\log(hw)+c^{3}hw)$ | $\mathcal{O}(chw\log(hw)+c^{2}hw)$
Emerging convolution | $\mathcal{O}(c^{2}hw)$ | $\mathcal{O}(c^{2}hw)$
ME-Woodbury transformation | $\mathcal{O}(dchw)$ | $\mathcal{O}(dchw)$
Woodbury transformation | $\mathcal{O}(dchw)$ | $\mathcal{O}(dchw)$

##  Appendix C Parameter Settings

In this section, we present additional details about our experiments to aid
reproducibility.

###  C.1 Experiments of Quantitative Evaluation

In the experiments of qualitative evalution, we compare Woodbury
transformations with $3$ permutation layer baselines, i.e., 1x1 convolution,
emerging convolution, and periodic coupling, and $2$ coupling layer baselines,
i.e., neural spline coupling, and MaCow. For all generalized permutation
methods, we use affine coupling, which is composed of $3$ convolutional
layers, and the $2$ latent layers have $512$ channels. For the neural spline
coupling, we set the number of spline bins to $4$. The spline parameters are
generated by a neural network, which is also composed of convolutional layers.
For $32\times 32$ images, we set the number of channels to $256$, and for
$64\times 64$ images, we set it to $224$. Ma et al. [25] used steps containing
a MaCow unit, i.e., $4$ autoregressive convolution coupling layers, and a full
Glow step. For fair comparison, we directly use the MaCow unit to replace the
affine coupling. For $32\times 32$ images, we set the convolution channel to
$384$, and for $64\times 64$ images, we set it to $296$.

We run each method to fixed number of iterations and test it every $10,000$
iterations. The bpds reported in our main paper are the best bpds obtained by
each method. The bpds are single-run results. This is because each run of the
experiment requires 3 to 5 days, and running each model multiple times is a
major cost. We found in our experiments that for the same model and parameter
settings, the bpds’ standard deviation of multiple runs are very small, i.e.,
around $0.003$, so single run results are sufficient for comparing bpd.

###  C.2 Hyper-parameter Settings

We use Adam [20] to tune the learning rates, with $\alpha=0.001$,
$\beta_{1}=0.9$, and $\beta_{2}=0.999$. We use uniform dequantization. The
sizes of models we use, and mini-batch sizes for training in our experiments
are listed in Table 4.

Table 4: Model sizes and mini-batch sizes.  Dataset | Mini-batch size | Levels(L) | Steps(K) | Coupling channels
---|---|---|---|---
CIFAR-10 32x32 | 64 | 3 | 8 | 512
ImageNet 32x32 | 64 | 3 | 8 | 512
ImageNet 64x64 | 32 | 4 | 16 | 512
LSUN Church 96x96 | 16 | 5 | 16 | 256
CelebA-HQ 64x64 | 8 | 4 | 16 | 512
CelebA-HQ 128x128 | 4 | 5 | 24 | 256
CelebA-HQ 256x256 | 4 | 6 | 16 | 256

###  C.3 Latent Dimension Settings

In all our experiments, we set the latent dimensions of Woodbury
transformations, and ME-Woodbury transformations as in Table 5.

Table 5: Latent dimensions of Woodbury transformations and ME-Woodbury transformations. The numbers in the brackets represent the latent dimension used in that level. For example, the $d_{c}:\\{8,8,16\\}$, represents that the settings of $d_{c}$ at the three levels are $8$, $8$, and $16$.  Dataset | Woodbury | ME-Woodbury
---|---|---
CIFAR-10 32x32 | $d_{c}:\\{8,8,16\\}$ | $d_{c}:\\{8,8,16\\}$
| $d_{s}:\\{16,16,8\\}$ | $d_{h}:\\{16,16,8\\}$
|  | $d_{w}:\\{16,16,8\\}$
ImageNet 32x32 | $d_{c}:\\{8,8,16\\}$ | $d_{c}:\\{8,8,16\\}$
| $d_{s}:\\{16,16,8\\}$ | $d_{h}:\\{16,16,8\\}$
|  | $d_{w}:\\{16,16,8\\}$
ImageNet 64x64 | $d_{c}:\\{8,8,16,16\\}$ | $d_{c}:\\{8,8,16,16\\}$
| $d_{s}:\\{16,16,8,8\\}$ | $d_{h}:\\{16,16,8,8\\}$
|  | $d_{w}:\\{16,16,8,8\\}$
LSUN Church 96x96 | $d_{c}:\\{8,8,16,16,16\\}$ | —
| $d_{s}:\\{16,16,16,8,8\\}$ |
CelebA-HQ 64x64 | $d_{c}:\\{8,8,16,16\\}$ | —
| $d_{s}:\\{16,16,8,8\\}$ |
CelebA-HQ 128x128 | $d_{c}:\\{8,8,16,16,16\\}$ | —
| $d_{s}:\\{16,16,16,8,8\\}$ |
CelebA-HQ 256x256 | $d_{c}:\\{8,8,16,16,16,16\\}$ | —
| $d_{s}:\\{16,16,16,16,8,8\\}$ |

##  Appendix D Sample Quality Comparisons

We compare the samples generated by Woodbury-Glow and Glow models trained on
the CelebA-HQ dataset. We follow Kingma and Dhariwal [22] and randomly hold
out 3,000 images as a test set. We use $5$-bits images. We use $64\times 64$,
$128\times 128$, $256\times 256$ images. Due to the our limited computing
resources, we use relatively small models. The model sizes and other settings
are listed in Table 4 and Table 5. We generate samples from the models during
different phases of training and display them in Figure 6, and Figure 7 (The
results of $64\times 64$ images are shown in the main paper). For the
$128\times 128$ images, both Glow and Woodbury-Glow generate distorted images
at iteration 100,000, but Woodbury-Glow seems to improve in later stages,
stabilizing the shapes of faces and structure of facial features. Glow,
continues generating faces with distorted overall shapes as training
continues. For the $256\times 256$ images, neither model ever trains
sufficiently to generate highly realistic faces, but Woodbury-Glow makes
significantly more progress in these 300,000 iterations than Glow. Glow’s
samples at 300,000 are still mostly random swirls with an occasional
recognizable face, while almost all of Woodbury-Glow’s samples look like
faces, though distorted. Due to limits on our computational resources, we
stopped the higher resolution experiments at 300,000 iterations (rather than
running to 600,000 iterations as we did for the $64\times 64$ experiments in
the main paper). With a larger model and longer training time, it seems
Woodbury-Glow would reach higher sample quality much faster than Glow.

The likelihoods of test set under the trained model are listed in Table 3. For
the $64\times 64$ and $128\times 128$ images, Woodbury-Glow scores higher
likelihood than Glow. For the $256\times 256$ images, their likelihoods are
almost identical, and are better than the score reported in [22]. This may be
due to three possible reasons: (1) We use affine coupling rather than additive
coupling, which is a non-volume preserving layer and may improve the
likelihoods; (2) Since the test set is randomly collected, it is different
from the one used in [22]; And (3) The model used in [22] is very large, so it
may be somewhat over-fitting. Surprisingly, the clear difference in sample
quality is not reflected by the likelihoods. This discrepancy may be because
we use $5$-bit images, and the images are all faces, so the dataset is less
complicated than other datasets such as ImageNet. Moreover, even though Glow
cannot generate reasonable $256\times 256$ samples, the colors of these
samples already match the colors of real images well, so these strange samples
may non-intuitively be equivalently likely as the face-like samples from
Woodbury-Glow.

Table 6: Bit per-dimension results on CelebA-HQ Size of images | Glow | Woodbury-Glow
---|---|---
$64\times 64$ | 1.27 | 1.23
$128\times 128$ | 1.09 | 1.04
$256\times 256$ | 0.93 | 0.93

Woodbury-Glow
Iteration 100,000

Woodbury-Glow
Iteration 200,000

Woodbury-Glow
Iteration 300,000

Glow
Iteration 100,000

Glow
Iteration 200,000

Glow
Iteration 300,000

Figure 6: Random samples of $128\times 128$ images drawn with temperature
$0.7$ from a model trained on CelebA data.

Woodbury-Glow
Iteration 150,000

Woodbury-Glow
Iteration 220,000

Woodbury-Glow
Iteration 300,000

Glow
Iteration 150,000

Glow
Iteration 220,000

Glow
Iteration 300,000

Figure 7: Random samples of $256\times 256$ images drawn with temperature
$0.7$ from a model trained on CelebA data.

##  Appendix E Additional Samples

In this section, we include additional samples from Woodbury-Glow models
trained on our various datasets. These samples complement our quantitative
analysis. We train our models on CIFAR-10 [24], ImageNet [32], the LSUN church
dataset [38], and the CelebA-HQ dataset [19]. Specifically, for ImageNet, we
use $32\times 32$ and $64\times 64$ images. For the LSUN dataset, we use the
same approach as Kingma and Dhariwal [22] to resize the images to be $96\times
96$. For the CelebA-HQ dataset, we use $64\times 64$, $128\times 128$, and
$256\times 256$ images. For LSUN and CelebA-HQ datasets, we use $5$-bit
images. The parameter settings of our models are in Table 4 and Table 5. The
samples are in Figures 8, 9, 10, 11, 12, 13, and 14.

Figure 8: CIFAR-10 $32\times 32$ Woodbury-Glow samples.

Figure 9: ImageNet $32\times 32$ Woodbury-Glow samples.

Figure 10: ImageNet $64\times 64$ Woodbury-Glow samples.

Figure 11: LSUN church $96\times 96$ Woodbury-Glow samples (temperature
$0.875$).

Figure 12: CelebA-HQ $64\times 64$ Woodbury-Glow samples (temperature $0.7$).

Figure 13: CelebA-HQ $128\times 128$ Woodbury-Glow samples (temperature
$0.5$).

Figure 14: Selected CelebA-HQ $256\times 256$ Woodbury-Glow samples
(temperature $0.5$).