11institutetext: Ho Chi Minh City University of Technology (HCMUT)
22institutetext: Vietnam National University - Ho Chi Minh City (VNU-HCM)
22email: {trung.pham.ktmt,bao.bachbbace12,lam.luu1602,minh.nguyen207bk,
<EMAIL_ADDRESS>
# An FPGA-based Solution for Convolution Operation Acceleration

Trung Pham-Dinh  Bao Bach-Gia  Lam Luu-Trinh  Minh Nguyen-Dinh  Hai Pham-Duc
Khoa Bui-Anh  Xuan-Quang Nguyen  Cuong Pham-Quoc

###### Abstract

Hardware-based acceleration is an extensive attempt to facilitate many
computationally-intensive mathematics operations. This paper proposes an FPGA-
based architecture to accelerate the convolution operation - a complex and
expensive computing step that appears in many Convolutional Neural Network
models. We target the design to the standard convolution operation, intending
to launch the product as an edge-AI solution. The project’s purpose is to
produce an FPGA IP core that can process a convolutional layer at a time.
System developers can deploy the IP core with various FPGA families by using
Verilog HDL as the primary design language for the architecture. The
experimental results show that our single computing core synthesized on a
simple edge computing FPGA board can offer 0.224 GOPS. When the board is fully
utilized, 4.48 GOPS can be achieved.

###### Keywords:

Convolution Operation, FPGA, Hardware Acceleration, IP core, Edge Computing

##  1 Introduction

The prominence of Convolutional Neural Networks (CNN) has sparked computer
vision and AI technology. As the demand for these applications grows,
scientists have proposed many CNN models and fitted them to many aspects of
the digital world. Soon enough, embedded devices find themselves the urge to
incorporate the AI capability into their functionality to ease the exhausted
computation on the server. However, micro-controllers have limitations
regarding the resources needed to perform AI computation. Therefore, FPGAs
(Field Programmable Gate Arrays) dominate the battlefield, along with ASICs
(Application Specific Integrated Circuits) [1]. The FPGA technology has been
rapidly growing, with each new generation integrating more Intellectual
Properties (IP) to expand the chip’s capability. These IP cores are the
engineers’ resources to integrate any CNN functionality and inferences to the
embedded devices.

This paper proposes our initial attempt to bring AI functionality to edge
devices. Our work focuses on convolution operation acceleration using FPGA
technology. Convolution is long known for its intensive and time-consuming
computing process, and it occupies around 90% in nearly every CNN model.
Therefore, the solution to this acceleration problem will undoubtedly lay the
ground for further AI-at-the-edge developments.

The rest of the paper is organized as follows. The convolution operation and
related work are presented in Section 2. We describe the overview of the
proposed architecture in Section 3. Section 4 and 5 will continue to discuss
the implementation of the IP core and its behavior under simulation
environment, along with its synthesis report on different Xilinx’s FPGA
families. Finally, Section 6 concludes the paper.

##  2 Background and related work

In this section, we introduce an overview of the convolution operation that we
are targeting as our hardware-based IP core. We then present related work for
the FPGA-based CNN computing cores.

###  2.1 Convolution Operation

Convolution is an operation performing on two arguments referred as input
(images) and filters. The output of the operation is known as the feature map.
Consider our input image as a two-dimension matrix. The convolution operation
is performed as shown in Equation 1. We apply the element-wise multiplication
of the Kernel matrix and a partial matrix of the input Image and take the sum
of all the products. We continue to apply the same computation by ”sliding”
and ”weighted-summing” the Kernel on the Image until we yield all the outputs
of the Feature Map.

$F(i,j)=(I\circledast K)(i,j)=\sum\limits_{m}\sum\limits_{n}I(i+m,j+n)\times
K(m,n)$ (1)

In practical usage, most images come in RGB format, encoded to three matrices
corresponding to each dimension of the picture - Red, Green, and Blue. In this
case, the Kernel expands to a 3D tensor. Equation 2 describes the resulting
Feature Map retrieved by performing the Convolution Operation subjecting to
this change.

$F(i,j)=(I\circledast
K)(i,j,d)=\sum\limits_{d}\sum\limits_{m}\sum\limits_{n}I(i+m,j+n,d)\times
K(m,n,d)$ (2)

###  2.2 Related work

In recent years, many studies in the literature have been proposed to
accelerate CNN in FPGA. Some of them try to improve parallelism by exploiting
the huge amount of hardware resources in FPGA such as research in [3], [7],
[8], [9]. Another approach is proposals that reduce the complexity of CNN by
using lightweight models like BNN such as [14], [2], [5], [10], [11], [12],
[13]. However, all these systems focus on high-performance FPGA computing
platforms. Meanwhile, in this work, we target edge computing platforms with
less amount of hardware resources and low energy consumption.

##  3 IP Core Architecture

Figure 1 conveys the primary functionality of our IP core. Because the core
takes on one layer at a time, it always expects a set of feature maps
containing C channels and the associated sets of K kernels (each kernel
includes C channels) as inputs. After the accelerated computation, the IP core
produces another set of feature maps, including K channels. These operations
sum up the overall functionality of the architecture.

Figure 1: IP Core Functionality

The proposed solution in this paper exploits the profound advantages of SoC
architectural design. The integrated on-chip Processing System (PS) will
efficiently transfer the input data down the stream. The IP initially loads
the two sets of input data into its internal groups of BRAMs (Section 4.1 will
further elaborate on the organizing of data among the RAM pool). Since the
amount of data is typically large, we use a direct memory access controller,
or DMA, to handle the transfer; hence cutting down the workload on the PS.

Figure 2: Dataflow

The IP has two groups of BRAMs to store the input image and the input kernel
separately. The IP core’s internal logic enables the Computing Core module
when the DMA has finished transferring the input data. In the most abstract
sense, the Computing Core is a convolutional computation unit of the
architecture. When fed with an image channel and its associated kernels, it
spits out the resulting feature map after applying the convolutional
operation. The detailed implementation to achieve this functionality is far
more complicated and is discussed in Section 4.2.

After computing the feature map, the Computing Core will store it in another
set of BRAMs. The DMA is continually responsible for transferring the
calculated results back to the PS from this set of BRAM. All the
communications between the DMA and the BRAMs mentioned so far are through AXI4
interfaces - a communication interface protocol that Xilinx has adopted as the
standard for their IPs communication.

Figure 2 illustrates the dataflow discussed so far. The final sub-module to
mention here is the Controller unit. To perform a correct convolution
operation, it will receive the information needed from the PS (for example,
the dimension of the input image and the input kernel). Section 4.2 will
discuss the Computing core in depth.

##  4 Implementation

This section presents in detail the architecture of the IP core. We first
start describing our memory management. Then, the introduction of the IP’s
central computing core, including sub-components, will be followed.

###  4.1 Memory Management

Block RAMs (BRAMs) in FPGA are common resources for storing data as they
provide random access to memory elements. Xilinx provides IP Block Memory
Generator (BMG) to help designers access BRAM efficiently under various
configurations. However, BMG has only two ports for concurrently reading and
writing; therefore, the architecture will distribute and organize data into
multiple BMGs to exploit this concurrency behavior.

####  4.1.1 Input BRAMs

The architecture has a set of four BMGs for each input image (Figure 3). Each
BMG in the set will store one fourth of the image channels. Because the depth
(the number of channels) per image varies, the size of each BMG must be large
enough to hold the largest possible image. This also means that for small
images, there will be redundant, or meaningless, slots in the BMGs. Here we
use 4 BMGs but not some other arbitrary number because the majority of CNN
inferences like the AlexNet [6] or the MobileNet [4] architecture has one
common characteristic. For these inferences, all the the produced feature maps
are divisible by 4, except for the first input image. By separating the
channels into four sets, we can process them in parallel, hence speeding up
the computation of the operation.

Figure 3: The set of Image BRAMs. Here, $B$ is the largest possible feature
map size divided by 4, and $C$ is the total number of channels of this layer’s
input image (or feature map).

The same divisible-by-4 idea applies to the kernel inputs. However, one
standard convolutional step typically requires more than one multiple-channel
kernel. Therefore, we organize the kernel data according to the number of
kernels. Furthermore, because the number of channels is a multiple of 4, the
number of kernels associated with each convolution layer must also be
divisible by 4. Therefore, for each BMG that contains one-fourth of the
channels of every kernel, we further divide it into four smaller BMGs, each
holding one-fourth of the whole kernels. This distribution scheme guarantees
that each BMG portion of the image has an associated set of four BMGs,
breaking the computing process down to four separate tasks independent of each
other.

####  4.1.2 Output BRAMs

Since an output feature map of one convolutional layer is the required input
for the next one in a typical CNN model, the distribution of the computed
feature map among the output BRAMs is identical to that of the input image
BRAMs. It means that each one-fourth of the output channels will be stored in
a separate BGM. These output BGMs will then be used to either read the result
out (through the DMA) or compute the next convolution layer in the network.

###  4.2 Computing Core

Computing Core is the essential module of the proposed IP core architecture.
As stated before, this module is responsible for producing the weighted sum
(or partial sum) when feeding the input image and weight data from the BRAMs.
Subsequent sections discuss the design and behavior of this module from the
multi-channel view to the multi-kernel execution.

####  4.2.1 Multi-Channel Architecture

Figure 4: Computing Cores Diagram

The proposed architecture comprises four computing cores; each performs the
convolution operation on one-fourth of the depth of the input feature maps.
With the organization of the input data among the BRAMs discussed previously,
the four computing cores can execute the operation totally in parallel, called
multi-channel architecture. Each computing core has one associated image BRAM
and four corresponding weight BRAMs - complying exactly with the number of
BRAMs in the architecture. Figure 4 illustrates the proposed design. In this
architecture, computed PSUM values of each core get accumulated continually
into the output BRAMs until the processing depth of images is finished. The
computing cores will continue to repeat the process but with another set of
kernels. It keeps processing until all the kernels have been filtered.

####  4.2.2 Multi-Kernel Computing Core

Figure 5: A Single Computing Core. In the diagram, $W[k][c]$ is the channel
$c^{th}$ of kernel $k^{th}$, $i[.][c]$ is the partial input image at the
channel $c^{th}$, $\circledast$ is the convolution operator, $comp.core[i]$ is
the $i^{th}$ computing core. In general, $k$ iterates from $0$ to $K/4-1$,
where $K$ is the total number of input kernels; and $c$ iterates from $i*C/4$
to $(i+1)*C/4-1$, for each iteration of $k$, where $C$ is total number of
input kernels.

Each computing core in the architecture produces four PSUM (Partial Sums)
values due to four kernels received, as shown in Figure 5. The produced PSUM
values belong to the one-fourth output feature map channel, speeding up the
computing process four times more. A single computing core contains four more
sub-module called PCOREs. The transfer of data from the BRAMs to the PCOREs is
through the intermediate loaders. Each PCORE computes a PSUM value according
to the weight input it receives from the Weight Loader, while the Image Loader
holds a set of nine pieces of input values for all the four PCOREs. This
computing model is weight stationery since the weight data persists in the
weight loader.

In contrast, the image loader continually fetches different input images after
each computed set of PSUMs. The internal logic of a PCORE is simple. Each
contains a set of MAC (Multiply and Accumulate) units and adder modules to
perform a weighted-sum operation.

####  4.2.3 Pipeline

Essentially, the computing process of the computing cores can be divided into
two stages. The first stage transfers sufficient data from the BRAMs to the
loaders while the second stage performs the convolution operation and
accumulates the computed PSUMs into the output BRAMs. The architecture
pipelines these two stages to accelerate the computation further, effectively
cutting down the wasted cycles.

####  4.2.4 Bias Handling

We utilize a handful feature from Xilinx’s BMG IP to handle the bias input.
The input bias is first to get initialized into the output BRAMs through the
PS. Since the computed PSUMs from the computing core get accumulated into
these BRAMs, it effectively causes the same effect as if the biases were added
to the weighted sum. Therefore there is no logic needed to handle the bias.

##  5 Experimental Results

In this section, we present our synthesis results with three different Xilinx
FPGA devices, including xc7z020clg400, xc7z020clg484, and xzcu3eg-sbva484-1.
We then anlyze simulation results to claim the throughput of our computing
core.

###  5.1 Synthesis Results

Table 1 show the resources used on different FPGAs including the number of
Look-up Table (LUTs), Flip-flop (FF), and percentage of these resources.
Maximum working frequencies for operation are calculated based on ”Data Path
Delay” values of the synthesis. According to the table, our computing core
consumes less than 5% hardware resources of the Pynq Z2 board, one of the most
suitable FPGA boards for edge computing. In other words, we can deploy up to
20 cores concurrently to further improve computing performance.

Table 1: Synthesis result on different FPGAs FPGA | #LUTs | #FF | Max frequency
---|---|---|---
xc7z020clg400-1 | 5027 (9.45%) | 4959 (4.66%) | 112 MHz
xc7z020clg484-1 | 5243 (9.86%) | 5054 (4.75%) | 93 MHz
xzcu3eg-sbva484-1-i | 11917 (16.89%) | 14522 (10.29%) | 161 MHz

###  5.2 Simulation Results

Figure 6 describes a part of the simulation waveform of one Computing core on
Vivado. Each feature signal (feature 0, feature1, feature2) concatenates three
8-bit feature entries to create a 3x3 tile of the feature. Each weight signal
(weight0, weight1, weight2, weight3) links nine 8-bit weight entries, so each
weight signal represents a channel of input weight. Because each Computing
core can compute convolution on four different kernels simultaneously, four
weight signals correspond to four different kernels. As a result, there are
four partial sum signals (psum0, psum1, psum2, psum3), and each signal results
from a single weight channel overlapping over a single 3x3 tile of the
feature. The computing core needs eight clock cycles to compute four psum
values and accumulate them to BRAM outputs. The system has four Computing
cores to compute sixteen psum values for each eight clock cycles.

The system consists of four Computing cores; a single core can compute four
psum values for each eight clock cycles. With the input feature of [224x224x8]
and input weight of [8x3x3x8], the system needs to compute 3,154,176 psum
values. Based on Synthesis Result, the IP core can operate on 112MHz clock
frequency with the Pynq Z2 FPGA board. Therefore, we can deduce the theory
time needed for computing this sample, which is 0.01408 seconds. In other
words, the throughput of a single core is 0.224 GOPS. With the Pynq Z2 FPGA
board for edge computing devices mentioned above, when 20 cores are deployed,
our computing system can offer up to 4.48 GOPS.

Figure 6: One part of the waveform from the simulation of a single Computing
core

##  6 Conclusion

With the rapid growth in AI, the demand of real time processing is
significantly high. Convolution operation is the core of convolution neural
network and take up most of the computation time. Therefore, this paper
proposes a FPGA-based architecture for accelerating convolution operation. The
architecture is implemented by Verilog and be able to operate on different
FPGA devices. This architecture can compute on multiple channels of feature
and multiple kernels of filter simultaneously, which increases the parallelism
on the computation. In addition, load stage and computation stage are
pipelined, which significantly reduces the computation time and reachs 0.224
GOPS on a single core. This core is synthesized on different FPGAs and reachs
maximum frequency of 161MHz with reasonable used resources. Future
optimization on this architecture is try to make it more flexible and able to
operate on a higher frequency to prove the great prospect of FPGA in AI
processing.

## Acknowledgement

This research is funded by Vietnam National University - Ho Chi Minh City
(VNU-HCM) under grant number B2021-20-02. We acknowledge the support of time
and facilities from Ho Chi Minh City University of Technology (HCMUT), VNU-HCM
for this study.

## References

  * [1] DeMicheli, G., Sami, M.: Hardware/Software Co-Design. Nato Science Series E:, Springer Netherlands (1996), https://www.springer.com/gp/book/9780792338833
  * [2] Ghasemzadeh, M., Samragh, M., Koushanfar, F.: Rebnet: Residual binarized neural network. In: 2018 IEEE 26th Annual International Symposium on Field-Programmable Custom Computing Machines (FCCM). pp. 57–64. IEEE Computer Society, Los Alamitos, CA, USA (may 2018)
  * [3] Guo, K., Han, S., Yao, S., Wang, Y., Xie, Y., Yang, H.: Software-hardware codesign for efficient neural network acceleration. IEEE Micro 37(2), 18–25 (2017)
  * [4] Howard, A.G., Zhu, M., Chen, B., Kalenichenko, D., Wang, W., Weyand, T., Andreetto, M., Adam, H.: Mobilenets: Efficient convolutional neural networks for mobile vision applications (2017)
  * [5] Jiao, L., Luo, C., Cao, W., Zhou, X., Wang, L.: Accelerating low bit-width convolutional neural networks with embedded fpga. In: 2017 27th International Conference on Field Programmable Logic and Applications (FPL). pp. 1–4 (2017)
  * [6] Krizhevsky, A., Sutskever, I., Hinton, G.E.: Imagenet classification with deep convolutional neural networks. In: Pereira, F., Burges, C.J.C., Bottou, L., Weinberger, K.Q. (eds.) Advances in Neural Information Processing Systems. vol. 25. Curran Associates, Inc. (2012), https://proceedings.neurips.cc/paper/2012/file/c399862d3b9d6b76c8436e924a68c45b-Paper.pdf
  * [7] Li, H., Fan, X., Jiao, L., Cao, W., Zhou, X., Wang, L.: A high performance fpga-based accelerator for large-scale convolutional neural networks. In: 2016 26th International Conference on Field Programmable Logic and Applications (FPL). pp. 1–9 (2016)
  * [8] Lin, X., Yin, S., Tu, F., Liu, L., Li, X., Wei, S.: Lcp: a layer clusters paralleling mapping method for accelerating inception and residual networks on fpga. In: 2018 55th ACM/ESDA/IEEE Design Automation Conference (DAC). pp. 1–6 (2018)
  * [9] Liu, Z., Dou, Y., Jiang, J., Xu, J.: Automatic code generation of convolutional neural networks in fpga implementation. In: 2016 International Conference on Field-Programmable Technology (FPT). pp. 61–68 (2016)
  * [10] Moss, D.J.M., Nurvitadhi, E., Sim, J., Mishra, A., Marr, D., Subhaschandra, S., Leong, P.H.W.: High performance binary neural networks on the xeon+fpga™ platform. In: 2017 27th International Conference on Field Programmable Logic and Applications (FPL). pp. 1–4 (2017)
  * [11] Nakahara, H., Fujii, T., Sato, S.: A fully connected layer elimination for a binarizec convolutional neural network on an fpga. In: 2017 27th International Conference on Field Programmable Logic and Applications (FPL). pp. 1–4 (2017)
  * [12] Nurvitadhi, E., Sheffield, D., Sim, J., Mishra, A., Venkatesh, G., Marr, D.: Accelerating binarized neural networks: Comparison of fpga, cpu, gpu, and asic. In: 2016 International Conference on Field-Programmable Technology (FPT). pp. 77–84 (2016)
  * [13] Prost-Boucle, A., Bourge, A., Pétrot, F., Alemdar, H., Caldwell, N., Leroy, V.: Scalable high-performance architecture for convolutional ternary neural networks on fpga. In: 2017 27th International Conference on Field Programmable Logic and Applications (FPL). pp. 1–7 (2017)
  * [14] Yang, Y., Huang, Q., Wu, B., Zhang, T., Ma, L., Gambardella, G., Blott, M., Lavagno, L., Vissers, K., Wawrzynek, J., Keutzer, K.: Synetgy: Algorithm-hardware co-design for convnet accelerators on embedded fpgas. In: Proceedings of the 2019 ACM/SIGDA International Symposium on Field-Programmable Gate Arrays. p. 23–32. FPGA ’19, ACM, New York, NY, USA (2019)


# Expressing linear equality constraints in feedforward neural networks

Anand Rangarajan, Pan He, Jaemoon Lee, Tania Banerjee and Sanjay Ranka

###### Abstract

We seek to impose linear, equality constraints in feedforward neural networks.
As top layer predictors are usually nonlinear, this is a difficult task if we
seek to deploy standard convex optimization methods and strong duality. To
overcome this, we introduce a new saddle-point Lagrangian with auxiliary
predictor variables on which constraints are imposed. Elimination of the
auxiliary variables leads to a dual minimization problem on the Lagrange
multipliers introduced to satisfy the linear constraints. This minimization
problem is combined with the standard learning problem on the weight matrices.
From this theoretical line of development, we obtain the surprising
interpretation of Lagrange parameters as additional, penultimate layer hidden
units with fixed weights stemming from the constraints. Consequently, standard
minimization approaches can be used despite the inclusion of Lagrange
parameters—a very satisfying, albeit unexpected, discovery. Examples ranging
from multi-label classification to constrained autoencoders are envisaged in
the future.

Dept. of CISE

University of Florida, Gainesville, FL

##  1 Introduction

The current resurgence of interest and excitement in neural networks, mainly
driven by high performance computing, has tended to shift the research focus
away from applied math development. While this has enabled a dizzying array of
practical applications and has firmly brought artificial intelligence and
machine learning into the mainstream, fundamental issues remain. One such
issue is the enforcement of constraints in neural network information
processing. By far the most common architecture is the feedforward neural
network with its implicit hidden units, nonlinearities and divergence measures
linking predictors to targets. However, if we seek to impose mathematically
specified constraints on the predictors since the targets satisfy them, there
is a paucity in the literature to address this issue. A survey of the
literature brings out methods that softly impose constraints using penalty
functions and there is a smattering of work focused on constraint
satisfaction, but in the main, there is no work that makes the satisfaction of
hard equality constraints central to its mission. Addressing this lack is the
motivation behind the paper.

The principal difficulty in imposing linear (or other) equality constraints
lies in the nonlinear nature of the predictor of the output. In a standard
feedforward network, the predictor is often a pointwise nonlinearity applied
to a linear transformation of the output of the penultimate layer. Enforcing
linear constraints on the predictor using standard Lagrangian methods is not
straightforward as they effectively result in a set of nonlinear constraints
on the top layer weights—difficult therefore (at first glance) to cast using
standard convex optimization methods. Not taking advantage of convex
optimization and duality would be unfortunate given the richness (and
simplicity) of the methods at hand. On the other hand, the softmax
nonlinearity is a clear cut example of a hard (sum to one) constraint being
enforced, usually at the top layer to obtain one-of-$K$ multi-class soft
memberships. This seems to run counter to our argument that enforcing
constraints on nonlinear predictors is difficult. We therefore ask (and answer
in the affirmative) the following question: can this example be generalized to
other constraints?

We show that linear equality constraints can be expressed in feedforward
neural networks using standard convex optimization methods. To achieve this,
we introduce auxiliary variables (one set for each training set instance) and
express the constraints on them. It turns out that we obtain a _maximization_
problem on the auxiliary variables (which take the form of predictors) and
even more crucially, we obtain a convex dual _minimization_ problem on the
Lagrange parameter vectors upon elimination of the auxiliary set. Furthermore,
if the form of the output nonlinearity is carefully matched with the choice of
the output divergence measure, the Lagrange parameters can be interpreted as
additional penultimate layer hidden units which connect to the output via
fixed weights (depending on the constraints). This is a crucial discovery
since these additional variables can be coupled by standard means to lower
layers (finding application in compression autoencoders for example). No
saddle-point methods are needed and while it may be necessary to explicitly
solve for the Lagrange parameters in certain applications, we can deploy
standard stochastic gradient descent methods here—a huge computational
advantage.

The development of the basic idea proceeds as follows. We begin with the
multi-class logistic regression (MC-LR) negative log-likelihood (which can
serve as the top layer) and note (as mentioned above) that the softmax
nonlinearity of the predictor satisfies a sum to one constraint. Then, we use
an Ansatz (an educated guess) for the form of the saddle-point Lagrangian with
auxiliary variables, which when eliminated results in a dual minimization on
the Lagrange parameter used to set up the softmax one-of-$K$ membership
constraint. Subsequently, we show that a difference between squared $\ell_{2}$
norms can serve as a principle for setting up a saddle-point Lagrangian with
linear constraints—albeit one with a linear output layer. Moving forward to
discover a principle, the next step in the technical development is crucial:
we identify the difference of squared $\ell_{2}$ norms as a special case of
the difference of _Bregman_ _divergences_ with a direct relationship between
the choice of the output nonlinearity and the form of the Bregman divergence
[Bregman, 1967]. Once this identification has been made, the remainder of the
development is entirely straightforward. Auxiliary variables (arranged in the
form of predictors) are eliminated from the difference of Bregman divergences
saddle-point Lagrangian leaving behind a dual _minimization_ w.r.t. the
Lagrange parameter vectors while continuing to obtain a standard minimization
problem w.r.t. the network weights. Matching the nonlinearity to the choice of
the Bregman divergence is important and we further require the output
nonlinearity to be continuous, differentiable and monotonic. Provided these
conditions are met, we obtain penultimate layer hidden unit semantics for the
Lagrange parameters—a surprising (but obvious in hindsight) result.

##  2 A summary of the constraint satisfaction approach

We present a very brief summary of the basic result. Assume a feedforward
neural network wherein we have the situation that the targets
$\left\\{\boldsymbol{y}_{n}\right\\}_{n=1}^{N}$ with
$\boldsymbol{y}_{n}\in\mathbb{R}^{K}$ satisfy the constraints
$A\boldsymbol{y}_{n}=\boldsymbol{b}_{n}$where
$\boldsymbol{b}_{n}\in\mathbb{R}^{I}$. Here the matrix of constraint weights
$A\equiv\left[\boldsymbol{a}_{1},\ldots,\boldsymbol{a}_{K}\right]$. We would
therefore like the nonlinear predictor from the top layer of the network also
satisfy (or approximately satisfy) these constraints. Assume top layer weights
$W\equiv\left[\boldsymbol{w}_{1},\ldots,\boldsymbol{w}_{K}\right]$ and
_penultimate_ layer hidden units $\boldsymbol{x}\in\mathbb{R}^{J}$. (It should
be clear that $K>I$ and that usually $J>K$.) Then, consider the following loss
function of the neural network (written in terms of the top layer weights and
for penultimate layer hidden units and for a single instance):

$\displaystyle\ell(W,\boldsymbol{\lambda})$
$\displaystyle=-\boldsymbol{y}^{T}\left(W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}\right)+\boldsymbol{e}^{T}\boldsymbol{\phi}(W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda})$
(1)
$\displaystyle=\sum_{k=1}^{K}\left[-y_{k}\left(\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\boldsymbol{a}_{k}^{T}\boldsymbol{\lambda}\right)+\phi\left(\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\boldsymbol{a}_{k}^{T}\boldsymbol{\lambda}\right)\right]$

where $\boldsymbol{e}$ is the vector of all ones,
$\boldsymbol{\lambda}\in\mathbb{R}^{I}$ and $\boldsymbol{\phi}$ is a vector of
convex, differentiable functions
$\boldsymbol{\phi}(\boldsymbol{u})\equiv\left[\phi(u_{1}),\ldots,\phi(u_{K})\right]^{T}$.
We assume the output nonlinearity corresponding to the nonlinear predictor is
pointwise and chosen to be $z=\sigma(u)=\phi^{\prime}(u)$.

We now minimize the loss function w.r.t. $\boldsymbol{\lambda}$ in (1) by
differentiating and setting the result to the origin to get

$\nabla_{\boldsymbol{\lambda}}\ell(W,\boldsymbol{\lambda})=-A\boldsymbol{y}+A\boldsymbol{z}=\boldsymbol{0}\Rightarrow-\boldsymbol{b}+A\boldsymbol{z}=\boldsymbol{0}$
(2)

where
$\boldsymbol{z}=\left[\phi^{\prime}(\boldsymbol{w}_{1}^{T}\boldsymbol{x}+\boldsymbol{a}_{1}^{T}\boldsymbol{\lambda}),\ldots,\phi^{\prime}(\boldsymbol{w}_{K}^{T}\boldsymbol{x}+\boldsymbol{a}_{K}^{T}\boldsymbol{\lambda})\right]^{T}$.
This implies that the loss function in (1) is an objective function w.r.t.
$\boldsymbol{\lambda}$ as well as the weights of the feedforward neural
network. Furthermore, there exists one $\boldsymbol{\lambda}$ vector per input
instance and when we minimize the loss function w.r.t. $\boldsymbol{\lambda}$,
we obtain a nonlinear predictor $\boldsymbol{z}$ which satisfies the output
constraints $A\boldsymbol{z}=\boldsymbol{b}$. The price of admission therefore
is an additional set of $\boldsymbol{\lambda}$ units (with fixed weights $A$)
which need to be updated in each epoch. In the rest of the paper, we develop
the theory underlying the above loss function.

##  3 Formal development of saddle-point Lagrangians and their dual

###  3.1 The Multi-class Logistic Regression (MC-LR) Lagrangian

The negative log-likelihood for multi-class logistic regression is

$-\log\Pr(\boldsymbol{Y}=\boldsymbol{y}|W,\boldsymbol{x})=-\sum_{k=1}^{K}y_{k}\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\log\sum_{k=1}^{K}\exp\left\\{\boldsymbol{w}_{k}^{T}\boldsymbol{x}\right\\}$
(3)

where the weight matrix $W$ comprises a set of weight vectors
$W=\left[\boldsymbol{w}_{1},\ldots,\boldsymbol{w}_{K}\right]$, the class label
column vector $\boldsymbol{y}=\left[y_{1},\ldots,y_{K}\right]^{T}$ and the
log-likelihood is expressed for a generic input-output training set pair
$\left\\{\boldsymbol{x},\boldsymbol{y}\right\\}$.

Rather than immediately embark on a tedious formal exposition of the basic
idea, we try and provide intuition. We begin with an Ansatz (an informed guess
based on a result shown in [Yuille and Kosowsky, 1994]) that leads to (3) and
then provide a justification for that educated guess. Consider the following
Lagrangian Ansatz for multi-class logistic regression:

$\displaystyle\mathcal{L}_{\mathrm{MC-
LR}}^{(\mathrm{total})}(W,\left\\{\boldsymbol{z}_{n}\right\\},\left\\{\lambda_{n}\right\\})$
$\displaystyle=\sum_{n=1}^{N}\left[\left(\boldsymbol{z}_{n}-\boldsymbol{y}_{n}\right)^{T}W^{T}\boldsymbol{x}_{n}-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z}_{n})+\lambda_{n}\left(\boldsymbol{e}^{T}\boldsymbol{z}_{n}-1\right)\right]$
$\displaystyle=\sum_{n=1}^{N}\Big{[}\left(z_{nk}-y_{nk}\right)\sum_{j}w_{jk}x_{nj}-\sum_{k}\psi(z_{nk})+\lambda_{n}\Big{(}\sum_{k}z_{nk}-1\Big{)}\Big{]}$
(4)

where $\boldsymbol{z}=\left[z_{1},\ldots,z_{K}\right]^{T}$ is the auxiliary
variable _predictor_ of the output,
$\boldsymbol{\psi}(\boldsymbol{z})\equiv\left[\psi(z_{1}),\ldots,\psi(z_{K})\right]^{T}$
is the _negative entropy_ and $\boldsymbol{e}$ is the column vector of all
ones with

$\psi(z)=z\log z-z.$ (5)

(We use the term negative entropy for $\psi$ since it returns the negative of
the Shannon entropy for the standard sigmoid.) The one-of-$K$ encoding
constraint is expressed by a Lagrange parameter $\lambda$ (one per training
set instance). In the context of the rest of the paper, the specific
nonlinearity (in place of the sigmoid nonlinearity) used for MC-LR is
$\sigma(u)\equiv\exp\left\\{u\right\\}$, its indefinite integral—the sigmoid
integral—$\phi(u)\equiv\exp\left\\{u\right\\}$, the inverse sigmoid
$u(z)\equiv\sigma^{-1}(z)=\left(\phi^{\prime}\right)^{-1}(z)=\log z$ and the
negative entropy $\psi(z)\equiv zu(z)-\phi(u(z))=z\log z-z$ [effectively a
Legendre transform of $\phi(u)$ [Legendre, 1787, Mjolsness and Garrett,
1990]]. The notation used in (4) and in the rest of the paper is now
described.

#### Notation:

$\sigma(u),\,$$\phi(u),\,$$u(z)$ and $\psi(z)$ are the generalized sigmoid,
sigmoid integral, inverse sigmoid and negative entropy respectively and
correspond to individual variables ($u$ and $z$ both in $\mathbb{R}).$ Please
see Table 1 for examples. Boldfaced quantities always correspond to vectors
and/or vector functions as the case may be. For example,
$\boldsymbol{\sigma}(\boldsymbol{u})$ is not a tensor and is instead unpacked
as a vector of functions
$\left[\sigma(u_{1}),\ldots,\sigma(u_{K})\right]^{T}$. The same applies for
$\boldsymbol{u}(\boldsymbol{z}),$
$\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))$ and
$\boldsymbol{\psi}(\boldsymbol{z})$. Likewise,
$\nabla\boldsymbol{\psi}(\boldsymbol{v})$ is not a gradient tensor but refers
to the vector of derivatives
$\left[\psi^{\prime}(v_{1}),\ldots,\psi^{\prime}(v_{K})\right]^{T}$.
$\boldsymbol{e}$ is the vector of all ones with implied dimensionality.
Similarly, $\nabla\boldsymbol{u}(\boldsymbol{z})$ is a vector of derivatives
$\left[u^{\prime}(z_{1}),\ldots,u^{\prime}(z_{K})\right]^{T}$ and
$\boldsymbol{z}\odot\nabla\boldsymbol{u}(\boldsymbol{z})$ is a vector with
$\odot$ denoting component-wise multiplication, i.e.
$\left[z_{1}u^{\prime}(z_{1}),\ldots,z_{K}u^{\prime}(z_{K})\right]^{T}$. In an
identical fashion,
$\nabla_{\boldsymbol{u}}\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))$ is
a vector of derivatives
$\left[\phi^{\prime}(u(z_{1})),\ldots,\phi^{\prime}(u(z_{K}))\right]^{T}$ and
$\nabla\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))\odot\nabla\boldsymbol{u}(\boldsymbol{z})$
denotes the vector formed by component-wise multiplication
$\left[\phi^{\prime}(u(z_{1}))u^{\prime}(z_{1}),\ldots,\phi^{\prime}(u(z_{K}))u^{\prime}(z_{K})\right]^{T}$.
The Bregman divergence
$B(\boldsymbol{z}||\boldsymbol{v})\equiv\sum_{k=1}^{K}B(z_{k}||v_{k})$—the sum
of component-wise Bregman divergences. Uppercase variables (such as $W,\,A$)
denote matrices with entries $\left\\{w_{kj}\right\\}$ and $\\{a_{ki}\\}$
respectively and $O$ denotes the origin. Lagrangian notation $\mathcal{L}$ is
strictly reserved for (optimization or saddle-point) problems _with
constraints_. Individual training set instance pairs
$(\boldsymbol{x},\boldsymbol{y})$ are often referred to with the instance
subscript index $n$ suppressed.

Generalized sigmoid $\sigma(u)$ |  Sigmoid integral $\phi(u)$ |  Inverse sigmoid $u(z)$ |  Negative entropy $\psi(z)$ |  Bregman divergence $B(z||v)$
---|---|---|---|---
Name | Expression
exp | $\exp\left\\{u\right\\}$ | $\exp\left\\{u\right\\}$ | $\log z$ | $z\log z-z$ | $z\log\frac{z}{v}-z+v$
sig | $\frac{1}{1+\exp\left\\{-u\right\\}}$ | $\log\left(1+\exp\left\\{u\right\\}\right)$ | $\log\frac{z}{1-z}$ |  $z\log z+$ $(1-z)\log(1-z)$ |  $z\log\frac{z}{v}+$ $(1-z)\log\frac{1-z}{1-v}$
tanh | $\tanh u$ | $\log\cosh u$ | $\frac{1}{2}\log\frac{1+z}{1-z}$ |  $\frac{1+z}{2}\log(1+z)+$ $\frac{1-z}{2}\log(1-z)$ |  $\frac{1+z}{2}\log\frac{1+z}{1+v}+$ $\frac{1-z}{2}\log\frac{1-z}{1-v}$
SRLU | $\frac{au\exp\left\\{\frac{au^{2}}{2}\right\\}+bu\exp\left\\{\frac{bu^{2}}{2}\right\\}}{\exp\left\\{\frac{au^{2}}{2}\right\\}+\exp\left\\{\frac{bu^{2}}{2}\right\\}}$ |  $\log\Big{(}\exp\left\\{\frac{au^{2}}{2}\right\\}$ $+\exp\left\\{\frac{bu^{2}}{2}\right\\}\Big{)}$ | - | - | -
$\ell_{2}$ | $u$ | $\frac{u^{2}}{2}$ | $z$ | $\frac{z^{2}}{2}$ | $\frac{1}{2}\left(z-v\right)^{2}$
Table 1: Generalized sigmoid to Bregman divergence table. SRLU stands for soft
rectified linear unit. (1$u(z)$ for the choice
$\sigma(u)=\frac{au\exp\left\\{\frac{au^{2}}{2}\right\\}+bu\exp\left\\{\frac{bu^{2}}{2}\right\\}}{\exp\left\\{\frac{au^{2}}{2}\right\\}+\exp\left\\{\frac{bu^{2}}{2}\right\\}}$
is left implicit since the inverse is not available in closed form. The same
applies for the negative entropy and the Bregman divergence.)

The Lagrangian in (4) expresses a _saddle-point_ problem—a minimization
problem on $W$ and a maximization problem on the predictors
$\left\\{\boldsymbol{z}_{n}\right\\}$. We will show that the Lagrangian is
concave w.r.t. $\left\\{\boldsymbol{z}_{n}\right\\}$ leading to a convex dual
minimization on $\left\\{\lambda_{n}\right\\}$. It turns out to be more
convenient to work with a single instance Lagrangian rather than (4) which is
on the entire training set. This is written as

$\mathcal{L}_{\mathrm{MC-
LR}}(W,\boldsymbol{z},\lambda)=\left(\boldsymbol{z}-\boldsymbol{y}\right)^{T}W^{T}\boldsymbol{x}-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z})+\lambda\left(\boldsymbol{e}^{T}\boldsymbol{z}-1\right).$
(6)

######  Lemma 1.

Eliminating $\boldsymbol{z}$ via maximization in (6) and solving for the
Lagrange parameter $\lambda$ via minimization yield the multi-class logistic
regression negative log-likelihood in (3).

###### Proof.

Since $\boldsymbol{\psi}(\boldsymbol{z})$ is continuous and differentiable, we
differentiate and solve for $\boldsymbol{z}$ in (6). Note that (6) contains a
sum of linear and concave (since $\psi(z)=z\log z-z$ is convex) terms on
$\boldsymbol{z}$ and is therefore concave. The negative of the Lagrangian is
convex and differentiable w.r.t. $\boldsymbol{z}$ and with a linear constraint
on $\boldsymbol{z}$, we are guaranteed that the dual is convex (and without
any duality gap) [Boyd and Vandenberghe, 2004, Slater, 2014]. We obtain

$\log z_{k}=\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda\Rightarrow
z_{k}=\exp\left\\{\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda\right\\}.$ (7)

(Note that the predictor $z_{k}$ has a nonlinear dependence on both $W$ and
$\lambda$ through the choice of sigmoid: $\exp$.) We eliminate
$\boldsymbol{z}$ to obtain the dual w.r.t. $\lambda$ written as

$\displaystyle R_{\mathrm{MC-LR}}(W,\lambda)$
$\displaystyle=-\sum_{k=1}^{K}y_{k}\boldsymbol{w}_{k}^{T}\boldsymbol{x}-\lambda+\sum_{k=1}^{K}\exp\left\\{\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda\right\\}$
$\displaystyle=-\sum_{k=1}^{K}y_{k}\boldsymbol{w}_{k}^{T}\boldsymbol{x}-\lambda+\sum_{k=1}^{K}\phi(\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda)$
$\displaystyle=-\sum_{k=1}^{K}y_{k}\left(\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda\right)+\sum_{k=1}^{K}\phi(\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda)$
(8)

where we have used the fact that $\sum_{k=1}^{K}y_{k}=1$. Also
$\phi(u)=\exp\left\\{u\right\\}$ is the sigmoid integral corresponding to the
choice of sigmoid $z(u)=\exp\left\\{u\right\\}$. Please see section 3.1 and
Table 1 for more information on the notation used. It should be understood
that we continue to have a minimization problem on $W$ but now augmented with
a convex minimization problem on $\lambda$. Differentiating and solving for
$\lambda$ in (8), we get

$-1+\sum_{k=1}^{K}\exp\left\\{\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\lambda\right\\}=0\Rightarrow\lambda=-\log\sum_{k=1}^{K}\exp\left\\{\boldsymbol{w}_{k}^{T}\boldsymbol{x}\right\\}.$
(9)

Substituting this solution for $\lambda$ in (8), we get

$R_{\mathrm{MC-
LR}}(W,\lambda(W))=-\sum_{k=1}^{K}y_{k}\boldsymbol{w}_{k}^{T}\boldsymbol{x}+\log\sum_{k=1}^{K}\exp\left\\{\boldsymbol{w}_{k}^{T}\boldsymbol{x}\right\\}+\mathrm{constant}$
(10)

which is identical (up to additive constant terms) to (3). ∎

We can write the dual (with the optimization on $W$) as a new MC-LR objective
function:

$E_{\mathrm{MC-
LR}}^{(\mathrm{total})}(W,\left\\{\lambda_{n}\right\\})=\sum_{n=1}^{N}\left[-\boldsymbol{y}_{n}^{T}W^{T}\boldsymbol{x}_{n}-\lambda_{n}+\boldsymbol{e}^{T}\boldsymbol{\phi}\left(W^{T}\boldsymbol{x}_{n}+\lambda_{n}\boldsymbol{e}\right)\right].$
(11)

It should be clear from Lemma 1 (that this is a _minimization_ problem on
$(W,\left\\{\lambda_{n}\right\\})$ despite the presence of the set of Lagrange
parameters $\left\\{\lambda_{n}\right\\}$. The structure of the objective
function strongly indicates that the Lagrange parameter $\lambda_{n}$ is
almost on the same footing as the set of penultimate layer hidden units
$\boldsymbol{x}_{n}$: it appears inside the nonlinearity in the form of a new
unit with a fixed weight of one and also connects to the loss in the same
manner. Noting that $\boldsymbol{e}^{T}\boldsymbol{y}_{n}=1,$ we rewrite (11)
as

$E_{\mathrm{MC-
LR}}^{(\mathrm{total})}(W,\left\\{\lambda_{n}\right\\})=\sum_{n=1}^{N}\left[-\boldsymbol{y}_{n}^{T}\left(W^{T}\boldsymbol{x}_{n}+\lambda_{n}\boldsymbol{e}\right)+\boldsymbol{e}^{T}\boldsymbol{\phi}\left(W^{T}\boldsymbol{x}_{n}+\lambda_{n}\boldsymbol{e}\right)\right]$
(12)

which is similar to (1) but with the specialization to multi-class
classification. This simple observation serves as an _underpinning_ for the
entire paper.

###  3.2 Difference of squared $\ell_{2}$ norms Lagrangian

We now embark on clarifying the intuition that led to the Ansatz in (4). Since
divergence measures are closely coupled to nonlinearities in neural networks,
it can be quite difficult to guess at the right Lagrangian which allows us to
move from a saddle-point problem to a dual minimization (on weight matrices
and Lagrange parameters). To (initially) remove nonlinearities completely from
the picture, consider the following saddle-point Lagrangian with a
minimization problem on $W$, a maximization problem on $\boldsymbol{z}$ and
the constraint $A\boldsymbol{z}=\boldsymbol{b}$ for a single training set pair
$(\boldsymbol{x},\boldsymbol{y}).$

$\mathcal{L}_{\ell_{2}}(W,\boldsymbol{z},\boldsymbol{\lambda})=\frac{1}{2}\|\boldsymbol{y}-W^{T}\boldsymbol{x}\|_{2}^{2}-\frac{1}{2}\|\boldsymbol{z}-W^{T}\boldsymbol{x}\|_{2}^{2}+\boldsymbol{\lambda}^{T}(A\boldsymbol{z}-\boldsymbol{b}).$
(13)

Without the constraints, the predictor attains its maximum for
$\boldsymbol{z}=W^{T}\boldsymbol{x}$, the optimal solution and we recover the
original unconstrained learning problem on the weights. With the constraints,
the maximization problem is not (usually) allowed to reach its unconstrained
value with the predictor likewise not allowed to reach $W^{T}\boldsymbol{x}$
due to constraint satisfaction. With the constrained maximization problem
nested inside the minimization problem, the intuition here is that the
maximization problem leads to a predictor $\boldsymbol{z}$ which depends only
on the input $\boldsymbol{x}$ (and not directly on $\boldsymbol{y}$) and
satisfies the constraints $(A\boldsymbol{z}=\boldsymbol{b})$. Furthermore,
elimination of $\boldsymbol{z}$ leads to a dual minimization problem on
$\boldsymbol{\lambda}$ which can be combined with the minimization problem on
$W$ in exactly the same way as in MC-LR above. After this elimination step,
$W$ is tuned to best fit to the “labels” $\boldsymbol{y}$ (but with a
regularization term obtained from the maximization problem, one that penalizes
deviation from the constraints). By this mechanism, the weight gradient of the
objective function gets projected into the subspace spanned by the
constraints. Consequently, once initial conditions are set up to satisfy the
constraints, the network cannot subsequently deviate from them. Since each
training set instance $(\boldsymbol{x},\boldsymbol{y})$ is deemed to satisfy
$A\boldsymbol{y}=\boldsymbol{b}$, we can separately evaluate the gradient term
for each training set pair for constraint satisfaction.

######  Lemma 2.

If $A\boldsymbol{y}=\boldsymbol{b}$, then
A$\nabla_{W}\mathcal{\mathcal{L}}_{\ell_{2}}(W,\boldsymbol{z}(W),\boldsymbol{\lambda}(W))=O,$
or in other words the gradient w.r.t. $W$ of the Lagrangian in (13) with
$\boldsymbol{z}$ and $\boldsymbol{\lambda}$ set to their optimal values is in
the nullspace of $A$.

###### Proof.

Since the Lagrangian in (13) is quadratic w.r.t. $\boldsymbol{z}$, we
differentiate and solve to get

$\boldsymbol{z}=W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}.$ (14)

Substituting this solution for $\boldsymbol{z}$ in (13), we obtain the dual
objective function (while continuing to have a minimization problem on $W$):

$\displaystyle R_{\ell_{2}}(W,\boldsymbol{\lambda})$
$\displaystyle=\frac{1}{2}\|\boldsymbol{y}-W^{T}\boldsymbol{x}\|_{2}^{2}+\frac{1}{2}\boldsymbol{\lambda}^{T}AA^{T}\boldsymbol{\lambda}+\boldsymbol{\lambda}^{T}AW^{T}\boldsymbol{x}-\boldsymbol{\lambda}^{T}\boldsymbol{b}$
$\displaystyle\propto-\boldsymbol{y}^{T}W^{T}\boldsymbol{x}-\boldsymbol{b}^{T}\boldsymbol{\lambda}+\frac{1}{2}\boldsymbol{x}^{T}WW^{T}\boldsymbol{x}+\boldsymbol{x}^{T}WA^{T}\boldsymbol{\lambda}+\frac{1}{2}\boldsymbol{\lambda}^{T}AA^{T}\boldsymbol{\lambda}$
$\displaystyle=-\boldsymbol{y}^{T}W^{T}\boldsymbol{x}-\boldsymbol{b}^{T}\boldsymbol{\lambda}+\frac{1}{2}\|W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}\|_{2}^{2}.$
(15)

This objective function is similar to the previous one obtained for MC-LR in
(8) and to help see the intuition, we rewrite it as

$\displaystyle S_{\ell_{2}}(W,\boldsymbol{z},\boldsymbol{\lambda})$
$\displaystyle=-\boldsymbol{y}^{T}W^{T}\boldsymbol{x}-\boldsymbol{y}^{T}A^{T}\boldsymbol{\lambda}+\boldsymbol{y}^{T}A^{T}\boldsymbol{\lambda}-\boldsymbol{b}^{T}\boldsymbol{\lambda}+\frac{1}{2}\|W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}\|_{2}^{2}$
$\displaystyle=-\boldsymbol{y}^{T}\boldsymbol{z}+\frac{1}{2}\boldsymbol{z}^{T}\boldsymbol{z}+\boldsymbol{\lambda}^{T}(A\boldsymbol{y}-\boldsymbol{b})\,\,\,\,\,\,\,\,\,\,\,\,(\mathrm{since}\,\boldsymbol{z}=W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda})$
$\displaystyle=-\boldsymbol{y}^{T}\boldsymbol{z}+\frac{1}{2}\boldsymbol{z}^{T}\boldsymbol{z}$
(16)

since $A\boldsymbol{y}=\boldsymbol{b}$ and with $\boldsymbol{z}$ as in (14).
This clearly shows the constrained predictor
$\boldsymbol{z}=W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}$ attempting to
fit to the labels $\boldsymbol{y}$.

We now enforce the constraint $A\boldsymbol{z}=\boldsymbol{b}$ in (14) to get

$AW^{T}\boldsymbol{x}+AA^{T}\boldsymbol{\lambda}=\boldsymbol{b}\Rightarrow\boldsymbol{\lambda}=\left(AA^{T}\right)^{-1}\left(\boldsymbol{b}-AW^{T}\boldsymbol{x}\right)$
(17)

provided $(AA^{T})^{-1}$ exists. Using the solution in (17) and eliminating
$\boldsymbol{\lambda}$ from the objective function in (13), we obtain

$R_{\ell_{2}}(W,\boldsymbol{\lambda}(W))\propto-\boldsymbol{y}^{T}W^{T}\boldsymbol{x}+\boldsymbol{x}^{T}WQ\boldsymbol{b}+\frac{1}{2}\boldsymbol{x}^{T}WPW^{T}\boldsymbol{x}$
(18)

where $P\equiv(I-A^{T}(AA^{T})^{-1}A)$ is a projection matrix (with $P^{2}=P$)
and $Q\equiv A^{T}(AA^{T})^{-1}$ and with properties $PQ=O,$ $AP=O$ and $AQ=I$
which can be readily verified. The objective function w.r.t. $W$ is clearly
bounded from below by virtue of the last term in $(\ref{eq:RL2_final}).$
Taking the gradient of $R_{\ell_{2}}$ w.r.t. $W$, we obtain

$\nabla_{W}R_{\ell_{2}}(W,\boldsymbol{\lambda}(W))=\left[-\boldsymbol{y}+Q\boldsymbol{b}+PW^{T}\boldsymbol{x}\right]\boldsymbol{x}^{T}.$
(19)

From this, we get

$A\nabla_{W}R_{\ell_{2}}(W,\boldsymbol{\lambda}(W))=\left[-A\boldsymbol{y}+AQ\boldsymbol{b}+APW^{T}\boldsymbol{x}\right]\boldsymbol{x}^{T}=\left(-A\boldsymbol{y}+\boldsymbol{b}\right)\boldsymbol{x}^{T}=O$
(20)

since $AQ=I,\,AP=O$ and provided the target $\boldsymbol{y}$ satisfies
$A\boldsymbol{y}=\boldsymbol{b}$ which was our starting assumption. ∎

We have shown that the gradient of the effective objective function w.r.t. $W$
for each training set pair $(\boldsymbol{x},\boldsymbol{y})$ obtained by
eliminating $\boldsymbol{z}$ and $\boldsymbol{\lambda}$ from the Lagrangian in
(13) lies in the nullspace of $A$. This observation may be useful mainly for
analysis since actual constraint satisfaction may be computationally more
efficient in the dual objective function (written for all training set
instances for the sake of completion):

$\displaystyle
E_{\mathrm{\ell_{2}}}^{(\mathrm{total})}(W,\left\\{\boldsymbol{\lambda}_{n}\right\\})$
$\displaystyle=\sum_{n=1}^{N}\left[-\boldsymbol{y}_{n}^{T}W^{T}\boldsymbol{x}_{n}-\boldsymbol{b}_{n}^{T}\boldsymbol{\lambda}_{n}+\frac{1}{2}\|W^{T}\boldsymbol{x}_{n}+A^{T}\boldsymbol{\lambda}_{n}\|_{2}^{2}\right]$
$\displaystyle=\sum_{n=1}^{N}\left[-\boldsymbol{y}_{n}^{T}\left(W^{T}\boldsymbol{x}_{n}+A^{T}\boldsymbol{\lambda}_{n}\right)+\frac{1}{2}\|W^{T}\boldsymbol{x}_{n}+A^{T}\boldsymbol{\lambda}_{n}\|_{2}^{2}\right]\,\,\,\,\,\,(\mathrm{since}\,A\boldsymbol{y}_{n}=\boldsymbol{b}_{n})$
$\displaystyle=\sum_{n=1}^{N}\left[-\boldsymbol{y}_{n}^{T}\boldsymbol{z}_{n}(W,\lambda_{n})+\frac{1}{2}\|\boldsymbol{z}_{n}(W,\lambda_{n})\|_{2}^{2}\right]$
(21)

where the predictor

$\boldsymbol{z}_{n}(W,\lambda_{n})\equiv
W^{T}\boldsymbol{x}_{n}+A^{T}\boldsymbol{\lambda}_{n}$ (22)

as before in (14) but modified to express the dependence on the parameters $W$
and $\boldsymbol{\lambda}_{n}$. The surprising fact that we can augment the
penultimate layer units and the final layer weights in this manner has its
origin in the difference between the squared $\ell_{2}$ norms used back in
(13). Therefore, we have shown that the components of the Lagrange parameter
vector $\boldsymbol{\lambda}$ can be interpreted as additional, penultimate
layer hidden units with fixed weights obtained from $A$. In the remainder of
the paper, we generalize this idea to Bregman divergences allowing us to prove
the same result for any suitable generalized sigmoid (under the presence of
linear equality constraints).

###  3.3 The difference of Bregman divergences Lagrangian

As we move from the specific (MC-LR and $\ell_{2}$) to the general (Bregman
divergence), we define some terms which simplify the notation and development.

######  Definition 3.

The generalized sigmoid vector $\boldsymbol{\sigma}(\boldsymbol{u})$ is
defined as
$\boldsymbol{\sigma}(\boldsymbol{u})=\left[\sigma(u_{1}),\ldots,\sigma(u_{K})\right]^{T}$
where $\sigma(u)$ is continuous, differentiable and monotonic
($\sigma^{\prime}(u)>0,\,\forall u$).

An example (different from the standard sigmoid) is $\sigma(u)=\tanh(u).$ Its
derivative $\sigma^{\prime}(u)=\mathrm{sech}^{2}(u)>0,\,\forall
u\in\left(-\infty,\infty\right)$.

######  Definition 4.

The generalized sigmoid integral is defined as the vector
$\boldsymbol{\phi}(\boldsymbol{u})=\left[\phi(u_{1}),\ldots,\phi(u_{K})\right]^{T}$
where $\phi(u)=\int_{-\infty}^{u}\sigma(u)du$.

Continuing the $\sigma(u)=\tanh(u)$ example, the sigmoid integral for
$\sigma(u)=\tanh(u)$ is $\phi(u)=\log\cosh(u)$.

######  Definition 5.

The generalized inverse sigmoid is defined as a vector
$\boldsymbol{u}(\boldsymbol{z})=\left[u(z_{1}),\ldots,u(z_{K})\right]^{T}$
where $u(z)=\sigma^{-1}(z)=\left(\phi^{\prime}\right)^{-1}(z)$.

The inverse sigmoid for $\sigma(u)=\tanh(u)$ is
$u(z)=\frac{1}{2}\log\frac{1+z}{1-z}$.

######  Definition 6.

The generalized negative entropy
$\boldsymbol{\psi}(\boldsymbol{z})=\left[\psi(z_{1}),\ldots,\psi(z_{K})\right]^{T}$
is defined as
$\boldsymbol{z}^{T}\boldsymbol{u}(\boldsymbol{z})-\boldsymbol{e}^{T}\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))$.

For $\sigma(u)=\tanh(u)$, the negative entropy
$\psi(z)=\frac{1}{2}\left[(1+z)\log(1+z)+(1-z)\log(1-z)\right]$ and
$\phi(u(z))=\frac{1}{2}\log(1-z^{2})$. The generalized negative entropy
$\psi(z)$ is related to the Legendre transform of $\phi(u)$ (and can be
derived from Legendre-Fenchel duality theory [Fenchel, 2014, Mjolsness and
Garrett, 1990]) but this relationship is not emphasized in the paper.

######  Corollary 7.

The generalized negative entropy $\boldsymbol{\psi}(\boldsymbol{z})$ is
convex.

###### Proof.

Since $\boldsymbol{\psi}(\boldsymbol{z})$ is composed of a set of individual
functions $\psi(z)$, we just need to show that $\psi(z)$ is convex. The
derivative of $\psi(z)$ is

$\psi^{\prime}(z)=u(z)+zu^{\prime}(z)-\phi^{\prime}(u(z))u^{\prime}(z)=u(z)+zu^{\prime}(z)-zu^{\prime}(z)=u(z)$
(23)

since $u(z)=\left(\phi^{\prime}\right)^{-1}(z)$. The second derivative
$\psi^{\prime\prime}(z)$ is

$\psi^{\prime\prime}(z)=u^{\prime}(z)=\frac{1}{\sigma^{\prime}(u(z))}>0$ (24)

since $u(z)=\sigma^{-1}(z)$ and $\sigma^{\prime}(u)>0$. Therefore the
generalized negative entropy is convex. ∎

We can directly show the convexity of the negative entropy $\psi(z)$ for
$\sigma(u)=\tanh(u)$ by evaluating
$\psi^{\prime\prime}(z)=u^{\prime}(z)=\frac{d}{dz}\left(\frac{1}{2}\log\frac{1+z}{1-z}\right)=\frac{1}{1-z^{2}}>0,\forall
z\in\left[-1,1\right]$. We now set up the Bregman divergence [Bregman, 1967]
for the convex negative entropy function $\psi(z)$.

######  Definition 8.

The Bregman divergence for a convex, differentiable function $\psi(z)$ is
defined as

$B(z||v)=\psi(z)-\psi(v)-\psi^{\prime}(v)(z-v).$ (25)

This can be extended to the generalized negative entropy
$\boldsymbol{\psi}(\boldsymbol{z})$ but the notation becomes a bit more
complex. For the sake of completeness, the Bregman divergence for
$\boldsymbol{\psi}(\boldsymbol{z})$ is written as

$B(\boldsymbol{z}||\boldsymbol{v})=\boldsymbol{e}^{T}\left(\boldsymbol{\psi}(\boldsymbol{z})-\boldsymbol{\psi}(\boldsymbol{v})\right)-\left(\nabla\boldsymbol{\psi}(\boldsymbol{v})\right)^{T}(\boldsymbol{z}-\boldsymbol{v})$
(26)

where $\nabla\boldsymbol{\psi}(\boldsymbol{v})$ is the gradient vector of the
generalized negative entropy defined as
$\nabla\boldsymbol{\psi}(\boldsymbol{v})\equiv\left[\psi^{\prime}(v_{1}),\ldots,\psi^{\prime}(v_{K})\right]^{T}$.

This definition can be used to build the Bregman divergence for the
generalized negative entropy $\psi(z)$ written as

$\displaystyle B(z||v)$ $\displaystyle=\psi(z)-\psi(v)-\psi^{\prime}(v)(z-v)$
$\displaystyle=zu(z)-\phi(u(z))-vu(v)+\phi(u(v))-u(v)(z-v)$
$\displaystyle=z\left[u(z)-u(v)\right]-\phi(u(z))+\phi(u(v)).$ (27)

Therefore, the Bregman divergence for the choice of $\sigma(u)=\tanh(u)$ is

$\displaystyle B_{\mathrm{tanh}}(z||v)$
$\displaystyle=z\left[\frac{1}{2}\log\frac{1+z}{1-z}-\frac{1}{2}\log\frac{1+v}{1-v}\right]-\frac{1}{2}\log\left(1-z^{2}\right)+\frac{1}{2}\log\left(1-v^{2}\right)$
$\displaystyle=\frac{1}{2}\left[(1+z)\log\frac{1+z}{1+v}+(1-z)\log\frac{1-z}{1-v}\right]$
(28)

which is the Kullback-Leibler (KL) divergence for a binary random variable $Z$
taking values in $\left\\{-1,1\right\\}$ which is what we would expect once we
make the identification that $z=2\Pr(Z=1)-1$ with a similar identification for
$v$.

######  Corollary 9.

The Bregman divergence $B_{\mathrm{gensig}}(z||v)$ is convex in $z$ for the
choice of the generalized negative entropy function $\psi(z)$.

###### Proof.

Since $B_{\mathrm{gensig}}(z||v)$ is a summation of $\psi(z)$ and a term
linear in $z$, it is convex in $z$ due to the convexity of $\psi(z)$ (shown
earlier in Lemma 7). ∎

We now describe the crux of one of the basic ideas in this paper: that the
difference between squared $\ell_{2}$ norms in (13) when generalized to the
difference of Bregman divergences leads to a dual minimization problem wherein
the Lagrange parameters are additional penultimate layer hidden units. In (13)
we showed that nesting a maximization problem on the constrained predictor
$\boldsymbol{z}$ inside the overall minimization problem w.r.t. $W$ had two
consequences: (i) eliminating $\boldsymbol{z}$ led to a minimization problem
on $\boldsymbol{\lambda}$; (ii) the saddle-point problem led to opposing
forces—from the target $\boldsymbol{y}$ and from the input
$\boldsymbol{x}$—which led to projected gradients which satisfied the
constraints. Furthermore, the Lagrange parameters carried the interpretation
of penultimate layer hidden units (albeit in a final layer with a linear
output unit) via the form of the dual minimization problem in (21). We now
generalize the difference of squared $\ell_{2}$ norms to the difference of
Bregman divergences and show that the same interpretation holds for suitable
nonlinearities in the output layer.

From the generalized negative entropy Bregman divergence, we form the
difference of Bregman divergences to get

$\displaystyle B_{\mathrm{gensig}}(y||v)-B_{\mathrm{gensig}}(z||v)$
$\displaystyle=(z-y)u(v)+yu(y)-\phi(u(y))-zu(z)+\phi(u(z)).$
$\displaystyle=(z-y)u(v)+\psi(y)-\psi(z).$ (29)

######  Lemma 10.

The difference of Bregman divergences is linear in the weight matrix $W$ for
the choice of the convex generalized negative entropy vector
$\boldsymbol{\psi}(\boldsymbol{z})$ and for
$\boldsymbol{v}=\boldsymbol{\sigma}(W^{T}\boldsymbol{x})$.

###### Proof.

This is a straightforward consequence of the two choices.

$\displaystyle B_{\mathrm{gensig}}(y||v)-B_{\mathrm{gensig}}(z||v)$
$\displaystyle=y(u(y)-u(v))-\phi(u(y))+\phi(u(v))-\left[z(u(z)-u(v))-\phi(u(z))+\phi(u(v))\right]$
$\displaystyle=(z-y)u(v)+yu(y)-\phi(u(y))-zu(z)+\phi(u(z))$
$\displaystyle=(z-y)u(v)-zu(z)+\phi(u(z))+\mathrm{terms}\,\mathrm{independent}\,\mathrm{of}\,z$
$\displaystyle=(z-y)u(v)-\psi(z)+\mathrm{terms}\,\mathrm{independent}\,\mathrm{of}\,z$
(30)

We write this for a single instance as

$\displaystyle
B_{\mathrm{gensig}}(\boldsymbol{y}||\boldsymbol{v})-B_{\mathrm{gensig}}(\boldsymbol{z}||\boldsymbol{v})$
$\displaystyle\propto(\boldsymbol{z}-\boldsymbol{y})^{T}\boldsymbol{u}(\boldsymbol{v})-\boldsymbol{z}^{T}\boldsymbol{u}(\boldsymbol{z})+\boldsymbol{e}^{T}\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z})).$
(31)
$\displaystyle=(\boldsymbol{z}-\boldsymbol{y})^{T}\boldsymbol{u}(\boldsymbol{v})-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z})$
(32)

Note that the difference of Bregman divergences depends on $\boldsymbol{v}$
only through $\boldsymbol{u}(\boldsymbol{v})$. We now set
$\boldsymbol{v}=\boldsymbol{\sigma}(W^{T}\boldsymbol{x})$ and rewrite the
difference between Bregman divergences using terms that only depend on
$W,\boldsymbol{x}$ and $\boldsymbol{z}$ to get

$B_{\mathrm{gensig}}(\boldsymbol{y}||\boldsymbol{\sigma}(W^{T}\boldsymbol{x}))-B_{\mathrm{gensig}}(\boldsymbol{z}||\boldsymbol{\sigma}(W^{T}\boldsymbol{x}))\propto(\boldsymbol{z}-\boldsymbol{y})^{T}W^{T}\boldsymbol{x}-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z})$
(33)

which is linear in the weight matrix $W$. ∎

######  Corollary 11.

The difference of Bregman divergences for the choice of
$\phi(z)=\frac{z^{2}}{2}$ is identical to the difference of squared $\ell_{2}$
norms in (13).

###### Proof.

We evaluate the difference of Bregman divergences for
$\phi(z)=\frac{z^{2}}{2}$ by substituting this choice of $\phi(z)$ in (29) to
get

$\displaystyle B_{\mathrm{\ell_{2}}}(y||v)-B_{\mathrm{\ell_{2}}}(z||v)$
$\displaystyle=(z-y)v+y^{2}-\frac{y^{2}}{2}-z^{2}+\frac{z^{2}}{2}$
$\displaystyle=\frac{1}{2}(y-v)^{2}-\frac{1}{2}(z-v)^{2}$ (34)

which when extended to the entire set of hidden units is the difference
between the squared $\ell_{2}$ norms as in (13). ∎

The general saddle-point problem can be now be expressed (for a single
training set pair) as

$\min_{W}\max_{\boldsymbol{z}}B_{\mathrm{gensig}}(\boldsymbol{y}||\boldsymbol{\sigma}(W^{T}\boldsymbol{x}))-B_{\mathrm{gensig}}(\boldsymbol{z}||\boldsymbol{\sigma}(W^{T}\boldsymbol{x}))\propto\min_{W}\max_{\boldsymbol{z}}(\boldsymbol{z}-\boldsymbol{y})^{T}W^{T}\boldsymbol{x}-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z})$
(35)

$\mathrm{subject}\,\mathrm{to}\,A\boldsymbol{z}=\boldsymbol{b}.$

The order of the maximization and minimization matters. It is precisely
because we have introduced a maximization problem on $\boldsymbol{z}$ _inside_
a minimization problem on the weights $W$ that we are able to obtain a dual
minimization problem on the Lagrange parameters in exactly the same manner as
in multi-class logistic regression. This will have significant implications
later. The difference of Bregman divergences is _concave_ in $\boldsymbol{z}$.
Note that the weight matrix $W$ appears linearly in (33) due to the choice of
$\boldsymbol{v}=\boldsymbol{\sigma}(W^{T}\boldsymbol{x})$ which produces the
correct generalization of (6). We add the linear constraints
$A\boldsymbol{z}=\boldsymbol{b}$ as before and write the Lagrangian for a
single instance

$\mathcal{L}_{\mathrm{gensig}}(W,\boldsymbol{z},\boldsymbol{\lambda})=\left(\boldsymbol{z}-\boldsymbol{y}\right)^{T}W^{T}\boldsymbol{x}-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z})+\boldsymbol{\lambda}^{T}(A\boldsymbol{z}-\boldsymbol{b})$
(36)

which can be expanded to the entire training set as

$\mathcal{L}_{\mathrm{gensig}}^{(\mathrm{total})}(W,\left\\{\boldsymbol{z}_{n}\right\\},\left\\{\boldsymbol{\lambda}_{n}\right\\})=\sum_{n=1}^{N}\left[\left(\boldsymbol{z}_{n}-\boldsymbol{y}_{n}\right)^{T}W^{T}\boldsymbol{x}_{n}-\boldsymbol{e}^{T}\boldsymbol{\psi}(\boldsymbol{z}_{n})+\boldsymbol{\lambda}_{n}^{T}(A\boldsymbol{z}_{n}-\boldsymbol{b}_{n})\right].$
(37)

where the negative entropy
$\boldsymbol{\psi}(\boldsymbol{z})=\boldsymbol{z}^{T}\boldsymbol{u}(\boldsymbol{z})-\boldsymbol{e}^{T}\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))$.
The final saddle-point Lagrangian contains a minimization problem on $W$ and a
maximization problem on $\left\\{\boldsymbol{z}_{n}\right\\}$ with the unknown
Lagrange parameter vectors $\left\\{\boldsymbol{\lambda}_{n}\right\\}$
handling the linear constraints. Note that
$\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))$ lurking inside
$\boldsymbol{\psi}(\boldsymbol{z})$ is a generalized sigmoid integral (as per
Definition 4) with $\boldsymbol{u}(\boldsymbol{z})$ being the _matched_
generalized inverse sigmoid. The generalized sigmoid integral is not
restricted to the exponential function.

We now show that the generalized dual objective function on
$\left\\{\boldsymbol{\lambda}_{n}\right\\}$ (along with a minimization problem
on $W$) is of the same form as the MC-LR dual objective function (11) and is
written as

$E_{\mathrm{gensig}}^{(\mathrm{total})}(W,\left\\{\boldsymbol{\lambda}_{n}\right\\})=\sum_{n=1}^{N}\left[-\boldsymbol{y}_{n}^{T}\left(W^{T}\boldsymbol{x}_{n}+A^{T}\boldsymbol{\lambda}_{n}\right)+\boldsymbol{e}^{T}\boldsymbol{\phi}(W^{T}\boldsymbol{x}_{n}+A^{T}\boldsymbol{\lambda}_{n})\right].$
(38)

######  Lemma 12.

Eliminating $\boldsymbol{z}$ via maximization in (37) and solving for the
Lagrange parameter vector $\boldsymbol{\lambda}$ via minimization yield the
dual objective function in (38) provided $u(z)$ is set to
$\left(\phi^{\prime}\right)^{-1}(z)$.

###### Proof.

Since $\boldsymbol{\psi}(\boldsymbol{z})$ is continuous and differentiable, we
differentiate and solve for $\boldsymbol{z}$ in (37). This is simpler to write
for the single instance Lagrangian in (36). Differentiating (36) w.r.t.
$\boldsymbol{z}$ and setting the result to $\boldsymbol{0}$, we get

$W^{T}\boldsymbol{x}-\boldsymbol{u}(\boldsymbol{z})-\boldsymbol{z}\odot\nabla\boldsymbol{u}(\boldsymbol{z})+\nabla\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))\odot\nabla\boldsymbol{u}(\boldsymbol{z})+A^{T}\boldsymbol{\lambda}=\boldsymbol{0}.$
(39)

(Please see section 3.1 on the notation used for an explanation of these
terms.) Since $\phi^{\prime}(u(z))=\phi^{\prime}((\phi^{\prime})^{-1}(z))=z$,
we have
$\boldsymbol{z}=\nabla_{\boldsymbol{u}}\boldsymbol{\phi}(\boldsymbol{u}(\boldsymbol{z}))$.
From this we get

$\boldsymbol{u}(\boldsymbol{z})=W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}$
(40)

and therefore

$\boldsymbol{z}^{T}W^{T}\boldsymbol{x}-\boldsymbol{z}^{T}\boldsymbol{u}(\boldsymbol{z})+\boldsymbol{z}^{T}A^{T}\boldsymbol{\lambda}=0.$
(41)

Substituting this in (37), we get the dual objective function (along with a
minimization problem on $W$) written for just a single training set instance
to be

$\displaystyle E_{\mathrm{gensig}}(W,\boldsymbol{\lambda})$
$\displaystyle=\left[-\boldsymbol{y}^{T}W^{T}\boldsymbol{x}-\boldsymbol{b}^{T}\boldsymbol{\lambda}+\boldsymbol{e}^{T}\boldsymbol{\phi}(W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda})\right]$
$\displaystyle=\left[-\boldsymbol{y}^{T}\left(W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}\right)+\boldsymbol{e}^{T}\boldsymbol{\phi}(W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda})\right]\,\,\,\,(\mathrm{since}\,A\boldsymbol{y}=\boldsymbol{b})$
$\displaystyle=\left[-\boldsymbol{y}^{T}\boldsymbol{u}+\boldsymbol{e}^{T}\boldsymbol{\phi}(\boldsymbol{u})\right]$
(42)

where $\boldsymbol{u}\equiv W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}$ is
the linear predictor (with
$\boldsymbol{z}=\boldsymbol{\sigma}(\boldsymbol{u})$ being the nonlinear
predictor). Equation (42) when expanded to the entire training set yields (38)
and can be rewritten in component-wise form as

$\displaystyle
E_{\mathrm{gensig}}^{(\mathrm{total})}\left(W,\left\\{\boldsymbol{\lambda}_{n}\right\\}\right)$
$\displaystyle=\sum_{n=1}^{N}\left[-\sum_{k=1}^{K}y_{nk}\left(\sum_{j=1}^{J}w_{jk}x_{nj}+\sum_{i=1}^{I}a_{ik}\lambda_{ni}\right)+\sum_{k=1}^{K}\phi\left(\sum_{j=1}^{J}w_{jk}x_{nj}+\sum_{i=1}^{I}a_{ik}\lambda_{ni}\right)\right]$
$\displaystyle=\sum_{n=1}^{N}\left[-\sum_{k=1}^{K}y_{nk}u_{nk}+\sum_{k=1}^{K}\phi\left(u_{nk}\right)\right]$
(43)

where

$u_{nk}\equiv\sum_{j=1}^{J}w_{jk}x_{nj}+\sum_{i=1}^{I}a_{ik}\lambda_{ni},$
(44)

the component-wise linear predictor counterpart to
$\boldsymbol{u}=W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}$. We have shown
that we get a minimization problem on both the weights and the Lagrange
parameter vectors (one per training set instance) in a similar manner to the
case of multi-class logistic regression. We stress that $z=\sigma(u)$ requires
$\sigma(u)=\phi^{\prime}(u):$ the chosen nonlinearity and the choice of
$\phi(u)$ in the loss function in (42) are in lockstep. This is a departure
from standard practice in present day deployments of feedforward neural
networks.

The minimization problem in (38) is valid for any generalized sigmoid
$\sigma(u)$ provided it is continuous, differentiable and monotonic. We state
without proof that the dual objective function is convex w.r.t. both the top
layer weights $W$ and the Lagrange parameter vectors
$\left\\{\boldsymbol{\lambda}_{n}\right\\}$. However, since we plan to deploy
standard stochastic gradient descent methods for optimization, this fact is
not leveraged here. ∎

We end by showing that the predictor $\boldsymbol{z}$ satisfies the linear
equality constraints.

######  Corollary 13.

The nonlinear predictor
$\boldsymbol{z}=\boldsymbol{\sigma}(\boldsymbol{u})=\boldsymbol{\sigma}\left(W^{T}\boldsymbol{x}+A^{T}\boldsymbol{\lambda}\right)$
satisfies the constraints $A\boldsymbol{z}=\boldsymbol{b}$ for the optimal
setting of the Lagrange parameter vector $\boldsymbol{\lambda}$.

###### Proof.

We differentiate (38) w.r.t. $\boldsymbol{\lambda}_{n}$ and set the result to
$\boldsymbol{0}$ to obtain

$\sum_{k=1}^{K}a_{ik}\phi^{\prime}\left(\sum_{j=1}^{J}w_{jk}x_{nj}+\sum_{i^{\prime}=1}^{I}a_{i^{\prime}k}\lambda_{ni^{\prime}}\right)=b_{ni},\,\forall
i\in\left\\{1,\ldots,I\right\\}$ (45)

or

$A\boldsymbol{z}_{n}=\boldsymbol{b}_{n},\,\forall n$ (46)

at the optimal value of the Lagrange parameter vector
$\boldsymbol{\lambda}_{n}$ and

$z_{nk}=\phi^{\prime}\left(\sum_{j=1}^{J}w_{jk}x_{nj}+\sum_{i=1}^{I}a_{ik}\lambda_{ni}\right)=\sigma\left(\sum_{j=1}^{J}w_{jk}x_{nj}+\sum_{i=1}^{I}a_{ik}\lambda_{ni}\right),\,\forall
n,k,$ (47)

and we are done. ∎

##  4 Applications

###  4.1 Synthetic Equality Constraints

We now demonstrate the incorporation of equality constraints through a simple
example. To this end, we construct a feedforward network for the XOR problem
using 4 fully connected layers. The last layer (fc4) is a linear classifier
with a single output whereas the previous layer (fc3) converts the fc2 outputs
into new features which are easier to classify than the original coordinates.
Consequently, we visualize the fc3 weights as enacting linear constraints on
the fc2 outputs in the form $A\boldsymbol{x}_{2}=\boldsymbol{b}$ where $A$ is
the set of weights corresponding to the fc3 layer with $\boldsymbol{x}_{2}$
being the output of the fc2 layer and $\boldsymbol{b}$ being the set of final
XOR features. We train this network in PyTorch using a standard approach with
2000 training set samples in $\mathbb{R}^{2}$. The network architecture is as
follows: 10 weights in fc1 (2 units to 5 units), 50 weights in fc2 (5 units to
10 units), 40 weights in fc3 (10 units to 4 units) and 4 weights in fc4 (4
units to 1 unit). Each set of weights in each layer is accompanied by bias
weights as well. After training, we extract the weights corresponding to $A$
(the weights of the fc3 layer) and store them.

Next, we set up our constraint satisfaction network. The goal of this network
is to begin with fresh inputs in $\mathbb{R}^{2}$ and map them to the outputs
of the previous network’s fc2 layer ($\boldsymbol{x}_{2}$ which are the new
targets $\boldsymbol{y}$) while satisfying the constraints
($A\boldsymbol{x}_{2}=\boldsymbol{b})$. This is explained in greater detail.
Assume we have new inputs $\boldsymbol{x}_{0}^{(\mathrm{new})}$. These are
sent into the XOR network to produce the output of the fc2 layer
$\boldsymbol{x}_{2}$ which we now call $\boldsymbol{y}$, the new targets. But,
we also know that the targets $\boldsymbol{y}$ satisfy additional constraints
in the form of the XOR network’s fc3 layer. Therefore, if we construct a new
network with $\boldsymbol{y}$ being the target of the fc2 layer then the new
network (henceforth XOR Lagrange) has constraints
$A\boldsymbol{x}_{2}^{(\mathrm{new})}=A\boldsymbol{y}=A\boldsymbol{x}_{2}$
which need to be enforced. Since the XOR fc3 layer maps 10 units to 4 units,
there are 4 constraints for each $\boldsymbol{x}_{2}^{(\mathrm{new})}$ in XOR
Lagrange. To enforce these, we set up a custom layer which has additional
fixed weights (mapping 10 units to 4 units) in the form of $A$ and a set of
Lagrange parameters (4 per training or test pattern) which need to be
obtained. We train the new XOR Lagrange network in the same way as the XOR
network but with one crucial difference. In each epoch and for each batch, we
use a constraint satisfaction network to solve for the 4 Lagrange parameters
for each pattern in each batch. A standard gradient descent approach with a
step size solver is used. We use a threshold of $5\times 10^{-3}$ for the loss
function w.r.t. $\boldsymbol{\lambda}$ (the Lagrange parameters)—essentially a
standard BCE loss—and a threshold of $1.0\times 10^{-6}$ for the step-size
parameter (where the search begins with $\alpha$ the step-size parameter set
to 0.1 and then halved each time if too high). If convergence is not achieved,
we stop after 1000 iterations of constraint satisfaction. Note that we do not
need to completely solve for constraint satisfaction in each epoch but have
chosen to do so in this first work.

XOR Lagrange was executed 1000 times on a fresh input (with 2000 patterns)
with the results shown in Figure 1. We discovered a very common pattern in the
optimization. Initially for every fresh input
$\boldsymbol{x}_{0}^{(\mathrm{new})}$, $\boldsymbol{x}_{2}^{(\mathrm{new})}$
is far away from $\boldsymbol{y}$ and the constraint satisfaction solver can
take many iterations to obtain $A\boldsymbol{x}_{2}^{(\mathrm{new})}\approx
A\boldsymbol{y}$ (and note that both $\boldsymbol{x}_{2}^{(\mathrm{new})}$ and
$\boldsymbol{y}$ are the outputs of sigmoids). Then, as optimization proceeds,
the number of constraint satisfaction iterations steadily decreases until
$\boldsymbol{x}_{2}^{(\mathrm{new})}\approx\boldsymbol{y}$ (and we always
execute for 100 epochs). In the figure, we depict this progression for the
maximum, minimum and mean number of iterations. The loss function in each
epoch is also shown (in lockstep) for the sake of comparison. While the mean
number of constraint satisfaction iterations trends lower, we do see instances
of large numbers of iterations in the maximum case. This suggests that we need
to explore more sophisticated optimization schemes in the future.

Figure 1: XOR Lagrange. In each subplot, we show training progress with
constraint satisfaction. In each top subplot, the constraint satisfaction
progress is depicted with its loss counterpart proceeding in lockstep below.
The maximum, minimum and mean iterations (for 1000 runs) are shown in the
left, middle and right subplots respectively.

##  5 Discussion

The motivation behind the work is to endow neural networks with constraint
satisfaction capabilities. We have shown this can be accomplished by appealing
to a saddle-point principle—the difference of Bregman divergences—with
constraints added to a maximization problem which is carefully defined on a
new set of auxiliary predictor variables. When these predictors are eliminated
from the Lagrangian, we obtain a dual minimization problem on the Lagrange
parameters (one per training set instance). The price of admission for
constraint satisfaction therefore is an expansion of the penultimate hidden
layer with additional Lagrange parameter units (and fixed weights
corresponding to the constraints). Solutions for the Lagrange parameters are
required in addition to weight learning with further optimization possibly
needed during testing. The advantage of this formulation is that we have a
single minimization problem (as opposed to a saddle-point problem) with the
unusual interpretation of Lagrange parameters as hidden units. Immediate
future work will contend not only with linear inequality constraints (perhaps
achieved via non-negativity constraints on the Lagrange parameters) but
nonlinear equality constraints which may require modifications of the basic
framework. We envisage a host of applications stemming from this work, driven
by domain scientists who can now attempt to impose more structure during
learning using application-specific constraints.

## References

  * [Boyd and Vandenberghe, 2004] Boyd, S. P. and Vandenberghe, L. (2004).  Convex Optimization.  Cambridge University Press.  7th printing with corrections, 2009, page 226.
  * [Bregman, 1967] Bregman, L. (1967).  The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming.  USSR Computational Mathematics and Mathematical Physics, 7(3):200–217.
  * [Fenchel, 2014] Fenchel, W. (2014).  On conjugate convex functions.  In Giorgi, G. and Kjeldsen, T. H., editors, Traces and Emergence of Nonlinear Programming, pages 125–129. Springer, New York, NY.
  * [Legendre, 1787] Legendre, A. M. (1787).  Mémoire sur l’intégration de quelques équations aux différences partielles.  Imprimerie royale.
  * [Mjolsness and Garrett, 1990] Mjolsness, E. and Garrett, C. (1990).  Algebraic transformations of objective functions.  Neural Networks, 3(6):651–669.
  * [Slater, 2014] Slater, M. (2014).  Lagrange multipliers revisited.  In Giorgi, G. and Kjeldsen, T. H., editors, Traces and Emergence of Nonlinear Programming, pages 293–306. Birkhauser, New York, NY.
  * [Yuille and Kosowsky, 1994] Yuille, A. L. and Kosowsky, J. J. (1994).  Statistical physics algorithms that converge.  Neural Computation, 6(3):341–356.


# Dynamic spreading and infiltration of a molten sand droplet on a porous
surface

Rahul Babu Koneru
Department of Aerospace Engineering
University of Maryland
College Park, MD 20742
<EMAIL_ADDRESS>
&Garrett Foresman
Department of Chemistry and Biochemistry
University of Maryland
College Park, MD 20742
&Alison Flatau
Department of Aerospace Engineering
University of Maryland
College Park, MD 20742
&Zhen Li
Department of Mechanical Engineering
Clemson University
Clemson, SC 29634 &Luis Bravo, Muthuvel Murugan, Anindya Ghoshal
US Army Research Laboratory
Aberdeen Proving Ground, MD 21005 &George Em Karniadakis
Division of Applied Mathematics
Brown University
Providence, RI 02912 Corresponding Author

###### Abstract

Compared to smooth surfaces, droplet spreading on porous surfaces is more
complex and has relevance in many engineering applications. In this work, we
investigate the infiltration dynamics of molten sand droplets on structured
porous surfaces using the multiphase many-body dissipative particle dynamics
(mDPD) method. We carry out three-dimensional simulations with different
equilibrium contact angles and surface porosities. The temporal evolution of
the radius of the wetted area follows a power law, as in the case of a smooth
surface. The infiltration rate on the other hand is dictated by the
competition between spreading and capillary inhibition of the pores.
Additionally, the temporal evolution of the droplet height and the contact
angle on the porous surface is also presented.

_Keywords_ Mesoscopic simulations, mDPD, droplet infiltration, porous surface

##  1 Introduction

Flow through porous media is widespread in nature and in many industrial
applications [1, 2]. On a nonporous surface, the early-time spreading is
dictated by the inertial-capillary forces until viscous dissipation at the
contact line enforces an equilibrium. On a porous surface, the spreading
dynamics are more complex due to the competition between the lateral motion of
the contact line and capillary inhibition of the pores. This depends on the
physical properties of the system and also the topology of the porous
substrate. The spreading of liquid droplets on porous substrates has been
studied widely using theoretical, experimental and numerical approaches [3, 4,
5]. The Lattice Boltzmann (LB) simulations of Frank _et al._ [5] and fully-
resolved volume-of-fluid simulations by Das _et al._ [4] provide a great
insight in to the dynamics of wetting and infiltration on porous surfaces for
low viscosity liquids. Typically in most theoretical studies, the droplet is
assumed to retain the spherical cap shape and the infiltration is treated
using the Darcy’s law. However, there is a substantial body of evidence which
indicate a non-spherical shape of the droplet cap based on the topology of the
surface textures and requires a detailed investigation at a pore level.

Of particular interest is the spreading and infiltration of molten sand
droplets, often referred to as CMAS, on the thermal barrier coating (TBC) of a
gas turbine engine. The molten sand droplets infiltrate the porous TBC coating
and upon solidification alter the thermal properties and cause mechanical
failure. To develop CMAS-resistant TBCs, it is imperative to understand the
underlying interfacial dynamics of the infiltration process. In this work,
three-dimensional (3D) numerical simulations are carried out using the many-
body dissipative particle dynamics (mDPD) method to investigate the spreading
and infiltration of a molten CMAS on porous substrates. The simulations are
carried-out at two equilibrium contact angles: $\theta_{eq}=39^{o}$ and
$70^{o}$. The porous surface is characterized by an isotropic arrangement of
identical cylindrical pillars to represent the top-coat of a TBC. While a
realistic CMAS deposition process is nonisothermal, the droplet is assumed to
be isothermal at $1260\ \rm{{}^{o}C}$ in this work to isolate the effects of
heat transfer.

##  2 Numerical method

The interactions between mDPD particles are considered pairwise and their
motion is governed by the Newton’s second law of motion [6]. The position and
velocity of the $i$-th mDPD particle is tracked using

$\frac{d\mathbf{r}_{i}}{dt}=\mathbf{v}_{i},\
m_{i}\frac{d\mathbf{v}_{i}}{dt}=\mathbf{F}_{i}=\sum_{j\neq
i}\mathbf{F}^{D}_{ij}+\mathbf{F}^{R}_{ij}+\mathbf{F}^{C}_{ij},$ (1)

where $m_{i}$, $\mathbf{r}_{i}$, $\mathbf{v}_{i}$, and $\mathbf{F}_{i}$
denotes the mass, position, velocity, and force vector, respectively, of
particle $i$. The total force $\mathbf{F}_{i}$ is the linear superposition of
the dissipative ($\mathbf{F}^{D}_{ij}$), random ($\mathbf{F}^{R}_{ij}$) and
conservative forces ($\mathbf{F}^{C}_{ij}$) between neighboring particles $i$
and $j$, which are computed by

$\displaystyle\mathbf{F}^{C}_{ij}$
$\displaystyle=A\omega_{c}(r_{ij})\mathbf{e}_{ij}+B\left(\rho_{i}+\rho_{j}\right)\omega_{d}(\mathbf{r}_{ij})\mathbf{e}_{ij},$
(2) $\displaystyle\mathbf{F}^{D}_{ij}$
$\displaystyle=-\gamma\omega_{D}(r_{ij})\left(\mathbf{e}_{ij}\cdot\mathbf{v}_{ij}\right)\mathbf{e}_{ij},$
(3) $\displaystyle\mathbf{F}^{R}_{ij}$
$\displaystyle=\beta\omega_{R}(r_{ij})(dt)^{-1/2}\xi_{ij}\mathbf{e}_{ij}.$ (4)

The pairwise interactions are governed by the weighting function
$\omega_{(*)}$ which depends on the relative distance
$r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ between the particles $i$ and $j$.
These interactions have a compact support and vanish beyond a cut-off distance
$r_{c}$. The constants $\gamma$, $\beta$, $A$ and $B$ determine the strengths
of each individual force. The relative velocity between a pair of particles is
given by $\mathbf{v}_{ij}=\mathbf{v}_{i}-\mathbf{v}_{j}$. The forces between a
pair of particles always lie along the line of centers
$\mathbf{e}_{ij}=\mathbf{r}_{ij}/r_{ij}$. The Gaussian white noise $\xi_{ij}$
is a random variable drawn from a Gaussian distribution with
$\langle\xi_{ij}(t)\rangle=0$ and
$\langle\xi_{ij}(t)\xi_{kl}(t^{\prime})\rangle=\left(\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right)\delta(t-t^{\prime})$
where $\delta_{ij}$ is the Kronecker delta and $\delta(t-t^{\prime})$ is the
Dirac delta function. The system of stochastic ODEs is integrated in time
using a modified velocity-Verlet (MVV) method. The fluctuating-dissipation
theorem (FDT)[6] leads to the following relation between the weights and the
strengths of these two forces

$\omega_{D}=\omega_{R}^{2},\ \beta^{2}=2\gamma k_{B}T,$ (5)

where $\omega_{D}=\omega_{R}^{2}=\left(1-r_{ij}/r_{D}\right)^{s}$ with
$s=0.5$, $r_{D}=1.45$ and the non-dimensional temperature $k_{B}T=1.0$ where
$k_{B}$ is the Boltzmann constant. The local density $\rho$ of a particle is
computed as weighted sum of its neighbors _i.e._
$\rho_{i}=\sum\omega_{\rho}\left(r_{ij}\right)$. For a more detailed
discussion see [7].

###  2.1 Parameter Mapping

In an mDPD simulation the fluid properties such as $\rho$, $\sigma$ and $\nu$
are computed a posteriori from mDPD simulations as opposed to the continuum-
based methods where they are imposed directly. To accurately represent the
properties of the fluid, the mDPD parameters have to be chosen carefully by
mapping the mDPD parameters to the physical system of interest _i.e._ molten
CMAS. The physical properties of CMAS at $1260\ \rm{{}^{o}C}$ are [8]: density
$\rho_{l}=2690$ $\rm{kg/m^{3}}$, surface tension $\sigma_{l}=0.46$ $\rm{N/m}$
and dynamic viscosity $\mu_{l}=3.6$ $\rm{Pa\cdot s}$. Following Arienti _et
al._ [9], reference length [L], mass [M] and time [T] scales are computed to
non-dimensionalize the system. The reference quantities are defined in terms
of the reference units as

$\rho_{ref}=\frac{[M]}{[L]^{3}};\ \sigma_{ref}=\frac{[M]}{[T]^{2}};\
\nu_{ref}=\frac{[L]^{2}}{[T]}.$ (6)

Following Koneru _et al._ [7], the kinematic viscosity ($\nu$) and the surface
tension ($\sigma$) of the mDPD system are computed using the doubly periodic
Poiseuille flow and the thin-film method. The physical and the reference
quantities obtained are tabulated in Table 1.

Table 1:  Physical and reference properties of molten CMAS. Property | Physical units | mDPD units
---|---|---
Density, $\rho$ | 2690 $\rm{kg/m^{3}}$ | 6.74
Dynamic viscosity, $\mu$ | 3.6 $\rm{Pa\cdot s}$ | 196.08
Kinematic viscosity, $\nu$ |  $1.33\times 10^{-3}$ $\rm{m^{2}/s}$ | 29.093
Surface tension, $\sigma$ | 0.46 $\rm{N/m}$ | 9.287
Length $[L]$ |  $17.017\times 10^{-6}$ $\rm{m}$ | 1.0
Mass $[M]$ | $1.964\times 10^{-8}~{}\rm{kg}$ | 1.0
Time $[T]$ | $6.297\times 10^{-6}~{}\rm{s}$ | 1.0

###  2.2 Simulation setup

The droplet and the porous wall are composed of randomly generated pre-
equilibriated mDPD particles. The initial radius of the droplet is fixed at
$R=10$ mDPD units in the simulations, which corresponds to 0.17 mm in physical
units. The porous surface is comprised of an isotropic arrangement of
cylindrical pillars. The void fraction of the surface $\phi_{v}=1-\pi
R_{c}^{2}/l$ is varied by changing the radius of the pillars while keeping the
number of pillars identical. This way the droplet interaction with the
cylindrical pillars remains consistent between different void fractions. A
schematic of the cross-section of the porous surface is shown in Fig. 1.

Figure 1: (a) A representative lattice of the surface porosity. The shaded
circles denote the cross-section of the cylindrical pillars of radius $R_{c}$
separated by a distance $l$. (b) The contact angle, spreading radius and the
droplet height measured with reference to the top of the porous surface. (c)
Perspective view of the simulation depicting the droplet spreading on top of
the pillars.

In this work, simulations are carried out for $\phi_{v}=0.35,\ 0.47\
\text{and}\ 0.58$. Various experimental studies [10, 8] indicate a shallow
equilibrium contact angle of molten CMAS. The mDPD parameters used in this
work are identical to that in [7] where the wetting behavior of molten CMAS
was investigated on a nonporous surface.

##  3 Results

The Ohnesorge number $Oh=\mu/\sqrt{\sigma\rho R}$ is typically used to
identify the relative strength of the viscous and inertial force which for the
molten CMAS droplet in this paper is 7.85. Given the large $Oh$ indicative of
viscosity-dominated spreading, $\tau_{v}$ is used as the reference timescale.
In Fig. 2, the infiltration of the droplet on a surface with $\phi_{v}=0.58$
is shown at three different nondimensional times $t/\tau_{v}$ for droplets
with an equilibrium contact angle of $\theta_{eq}=70^{o}$ and $39^{o}$ on a
nonporous surface. In the case of $\theta_{eq}=39^{o}$, the droplet
infiltrates completely within as short nondimensional time. Due to combined
lateral and normal inhibition of the droplet, the infiltration is highly
nonhomogeneous and exhibits a departure from the spherical cap shape of the
non-infiltrated portion of the droplet.

Figure 2: Droplet infiltration and spreading for different values of
$\theta_{eq}$ at $\phi_{v}=0.58$.

In Fig. 3, the nondimensional spreading radius, scaled by initial droplet
radius $R$, is plotted as a function of nondimensional time which is scaled by
$\tau_{v}$ for two configurations which produce an equilibrium contact angle
$\theta_{eq}=70^{o}$ and $39^{o}$ on a nonporous surface. The droplet
spreading is clearly characterized by an inertial regime which is soon
followed by a reduction in the spreading radius due to infiltration. For a
smaller $\theta_{eq}$, the rate of infiltration exceeds the rate of spreading
much earlier compared to a larger $\theta_{eq}$. Compared to the case of
$\theta_{eq}=70^{o}$, the transition at the end of the inertial regime is much
sharper for $\theta_{eq}=39^{o}$ due to a larger contact area between the
droplet and the porous surface. In [7], the inertial spreading rate of 0.26
was reported on a nonporous surface. Due to surface porosity, the spreading
rate in the inertial regime is slightly reduced to a value between 0.18 and
0.21.

Figure 3: (a) Nondimensional spreading radius of a molten CMAS droplet with
(a) $\theta_{eq}=39^{o}$ and (b) $\theta_{eq}=70^{o}$.

The nondimensional droplet height measured from the top of the cylindrical
pillars is plotted as a function of nondimensional time in Fig. 4. For
$\theta_{eq}=70^{o}$, the decrease in the droplet height is very similar for
all the surface porosities until about $t/\tau_{v}=2$ after which the case
with $\phi_{v}=0.58$ diverges. The drop height on the porous surface crosses
over the curve on a nonporous surface at about $t/\tau_{v}=10$. This crossover
was also reported by Frank _et al._ [5]. In the case of $\theta_{eq}=39^{o}$,
this cross-over happens much earlier at about $t/\tau_{v}=2$ and a continuous
and complete infiltration is observed. Since the spreading is driven by the
inertia at early time, the differences in the drop height are minimal. In the
later stages when capillary inhibition is dominant, the drop height decreases
faster for the porous surface with the rate of decrease of droplet height
increasing with increasing porosity in agreement with Washburn’s law [11].

Figure 4: (a) Nondimensional droplet height with (a) $\theta_{eq}=39^{o}$ and
(b) $\theta_{eq}=70^{o}$.

The apparent contact angle of the droplet over the porous surface is also
calculated by fitting an sphere through the mDPD particles present on the
surface of the droplet. The temporal evolution of the apparent contact angle
on the porous surface is shown in Fig. 5 along with the data from a nonporous
surface. Interestingly, the porous surface behaves as a Cassie-Baxter surface
for $\theta_{eq}=70^{o}$ wherein the apparent contact angle is larger on the
porous surface than that on a nonporous surface. On the other hand, a Wenzel-
like state is observed for $\theta_{eq}=39^{o}$ wherein the apparent contact
angle is lower than that on a nonporous surface. Due to the large viscosity of
molten CMAS, the droplet takes a very long time to reach an equilibrium on a
nonporous surface. In this work, homogeneous wetting of the pore is observed
wherein the porous surface is completely filled by the liquid.

Figure 5: Temporal evolution of the apparent contact angle for
$\theta_{eq}=39^{o}$ and $\theta_{eq}=70^{o}$.

##  4 Conclusions

In this work, we have investigated the spreading and infiltration of a molten
CMAS droplet on a porous surface using two equilibrium contact angles and
three different porosities. Due to the topology of the porous surface used in
this work, the infiltrated liquid spreads around the cylindrical pillar
resulting in a nonhomogeneous spreading. We also observe a departure from the
spherical cap shape of the droplet. For a given porosity, a lower equilibrium
contact angle induces a much faster drainage. Just like that on a nonporous
surface, we observe an inertial spreading regime but with a reduced slope
between 0.18 and 0.21. The competing effects between lateral spreading and
capillary inhibition of the droplet produce interesting results. For both the
contact angles, the droplet height reduces faster on the porous surface than
that on a nonporous surface after a certain $t/\tau_{v}$. The crossover of the
droplet height between the porous and nonporous occurs at much earlier for the
smaller contact angle since a larger surface area of the droplet is in contact
with the porous surface. We observe that the surface porosity has an opposite
effect on the the evolution of the contact angle. Compared to a nonporous
surface, the porous surface increases the apparent contact angle for a larger
$\theta_{eq}$ whereas the apparent contact angle is reduced for a smaller
$\theta_{eq}$.

### Acknowledgements

Authors A.F and R.B.K acknowledge the support received from the US Army
Research Office Mathematical Sciences Division for this research through grant
number W911NF1910225. L.B, A.G, M.M were supported by the US Army Research
Laboratory 6.1 basic research program in propulsion sciences. Z.L acknowledges
the support from the NSF grant OAC-2103967 and the NASA project 80NSSC21M0153.
The authors gratefully acknowledge the resources and support provided by
Department of Defense Supercomputing Resource Center (DSRC) through use of
"Narwhal" as part of the 2022 Frontier Project, Large-Scale Integrated
Simulations of Transient Aerothermodynamics in Gas Turbine Engines. The views
and conclusions contained in this document are those of the authors and should
not be interpreted as representing the official policies or positions, either
expressed or implied, of the U.S. Army Research Laboratory or the U.S.
Government. The U.S. Government is authorized to reproduce and distribute
reprints for Government purposes notwithstanding any copyright notation
herein.

## References

  * [1] Detlef Lohse.  Fundamental fluid dynamics challenges in inkjet printing.  Annual Review of Fluid Mechanics, 54:349–382, 2022.
  * [2] Tatiana Gambaryan-Roisman.  Liquids on porous layers: wetting, imbibition and transport processes.  Current opinion in colloid & interface science, 19(4):320–335, 2014.
  * [3] Timothy W Munuhe, Ruey-Hung Chen, Liang Zhu, and Ronghui Ma.  Modeling molten droplet spreading and infiltration into non-isothermal thermal barrier coatings.  International Journal of Heat and Mass Transfer, 182:121942, 2022\.
  * [4] Saurish Das, HV Patel, E Milacic, NG Deen, and JAM Kuipers.  Droplet spreading and capillary imbibition in a porous medium: A coupled IB-VOF method based numerical study.  Physics of Fluids, 30(1):012112, 2018.
  * [5] X Frank and Patrick Perré.  Droplet spreading on a porous surface: A lattice Boltzmann study.  Physics of Fluids, 24(4):042101, 2012.
  * [6] Robert D Groot and Patrick B Warren.  Dissipative particle dynamics: Bridging the gap between atomistic and mesoscopic simulation.  The Journal of Chemical Physics, 107(11):4423–4435, 1997.
  * [7] Rahul Babu Koneru, Alison Flatau, Zhen Li, Luis Bravo, Muthuvel Murugan, Anindya Ghoshal, and George Em Karniadakis.  Quantifying the dynamic spreading of a molten sand droplet using multiphase mesoscopic simulations.  Physical Review Fluids, 7(10):103602, 2022.
  * [8] Ravisankar Naraparaju, Juan J Gomez Chavez, Philipp Niemeyer, Kai-Uwe Hess, Wenjia Song, Donald B Dingwell, Siddharth Lokachari, CV Ramana, and Uwe Schulz.  Estimation of CMAS infiltration depth in EB-PVD TBCs: A new constraint model supported with experimental approach.  Journal of the European Ceramic Society, 39(9):2936–2945, 2019\.
  * [9] Marco Arienti, Wenxiao Pan, Xiaoyi Li, and George Karniadakis.  Many-body dissipative particle dynamics simulation of liquid/vapor and liquid/solid interactions.  The Journal of Chemical Physics, 134(20):204114, 2011.
  * [10] YX Kang, Y Bai, GQ Du, FL Yu, CG Bao, YT Wang, and F Ding.  High temperature wettability between CMAS and YSZ coating with tailored surface microstructures.  Materials Letters, 229:40–43, 2018.
  * [11] Edward W Washburn.  The dynamics of capillary flow.  Physical review, 17(3):273, 1921.