# Scalable Coherent Optical Crossbar Architecture using PCM for AI
Acceleration
††thanks: Acknowledgement Place Holder

Dan Sturm  Electrical and Computer Engineering
University of Washington
Seattle, USA
<EMAIL_ADDRESS>Sajjad Moazeni  Electrical and Computer Engineering
University of Washington
Seattle, USA
<EMAIL_ADDRESS>
###### Abstract

Optical computing has been recently proposed as a new compute paradigm to meet
the demands of future AI/ML workloads in datacenters and supercomputers.
However, proposed implementations so far suffer from lack of scalability,
large footprints and high power consumption, and incomplete system-level
architectures to become integrated within existing datacenter architecture for
real-world applications. In this work, we present a truly scalable optical AI
accelerator based on a crossbar architecture. We have considered all major
roadblocks and address them in this design. Weights will be stored on chip
using phase change material (PCM) that can be monolithically integrated in
silicon photonic processes. All electro-optical components and circuit blocks
are modeled based on measured performance metrics in a 45nm monolithic silicon
photonic process, which can be co-packaged with advanced CPU/GPUs and HBM
memories. We also present a system-level modeling and analysis of our chip’s
performance for the Resnet-50V1.5, considering all critical parameters,
including memory size, array size, photonic losses, and energy consumption of
peripheral electronics. Both on-chip SRAM and off-chip DRAM energy overheads
have been considered in this modeling. We additionally address how using a
dual-core crossbar design can eliminate programming time overhead at practical
SRAM block sizes and batch sizes. Our results show that a $128\times 128$
proposed architecture can achieve inference per second (IPS) similar to Nvidia
A100 GPU at $15.4\times$ lower power and $7.24\times$ lower area.

###### Index Terms:

Optical Neural Networks, AI Accelerator, Crossbar, Phase Change Material,
System-level Optimization

##  I Introduction

Recent advancements in artificial intelligence (AI) and machine learning (ML)
have been challenging our conventional computing paradigms by demanding
enormous computing power at a dramatically faster pace than Moore’s law [1].
We can compare the performance of today’s AI/ML processors from two key
aspects of compute power in terms of Tera operations per second (TOPS) and
energy-efficiency (TOPS/W) as illustrated in Fig. 1. Despite the promising
success of neuromorphic and analog-based computing in electrical domains for
low TOPS applications (edge-computing), these approaches can not satisfy the
requirements of datacenters and super computers. Due fundamental bandwidth
limitations, they cannot achieve high throughput. Optical neural networks
(ONNs) can potentially overcome this barrier by providing tens of GHz
bandwidths and ultra-low losses of photonic integrated circuits [2, 3].
However, realizing a practical ONN-based AI accelerator requires a holistic
system design that considers devices, circuits, chip architectures, and
algorithms.

Figure 1: Comparison of state-of-the-art AI/ML processors.

In this paper, we present a novel architecture for an ONN accelerator in which
the multiply-and-accumulate (MAC) operation is performed on a coherent
photonic crossbar with programmable phase-change materials (PCM). PCM enables
low power photonic computing by storing the weights on-chip in a nonvolatile
fashion. This design provides a compact and scalable solution for the first
time, which minimally relies on thermo-optic phase-shifters effect. We have
considered all critical circuit blocks and parameters, including memory size,
array size, photonic losses, and energy consumption of peripheral electronics
including analog-to-digital converters (ADCs), digital-to-analog converters
(DACs), serializers, and clocking. We develop a custom simulation framework
based on existing cycle-accurate simulation tools to model the compute cycles,
programming cycles, and DRAM accesses for a given neural network running on a
specific set of parameters (including size of SRAM, array, and batch).

Presented work focuses on inference of convolutional neural networks (CNN)
such as ResNet50 v1.5, which is used as a benchmark to compare our proposed
accelerator performance with state-of-the-art. Additionally, while precision
and process variation are major factors in all analog-based computers, here we
assume a INT6 precision for all the components as it has been shown to be
sufficient for neural networks with high accuracy [4, 5].

This paper is organized as follows: We briefly describe related work in
Section II. In Section III, we explain the principles of performing the MAC
operation in this work. We present an overview of overall chip architecture
and CNN operation in Section IV, and section V explains our custom simulation
methodology. Finally, we present the results of our fully optimized design in
Section VI, and compare those with state-of-the-art in Section VIII.

##  II Related Work

Researchers have recently proposed a variety of methods to realize an ONN.
Most of these works focus only on the physics and devices rather than
providing a system-level solution and analysis. In this presented discussion,
we only consider integrated solutions, as free-space solutions [6] lack the
reconfigurability that is an essential part of any “computer” and
compatibility with mainstream CMOS technology. Furthermore, we note that a
suitable application space for ONNs can be datacenters as opposed to edge
computing according to Fig. 1. Below we discuss the most promising solutions
so far from the perspective of three critical factors that we believe have
been addressed in this work:

(1) Scalability: While on-chip photonics provide high bandwidths, their
footprints are fundamentally orders of magnitude larger than advanced nm-scale
CMOS. With only one or two routing layers, building an ONN processor with big
dimensions has remained elusive. This will become even more challenging
considering the need for compact ADCs and DACs. Mach-Zehnder Interferometer
(MZI)-based coherent architectures such as [2] have large chip areas and
large-scale realizations end up exceeding a few $cm^{2}$. Non-coherent PCM-
based crossbars have been also proposed [7], however they require many
wavelengths for large matrix operations and that is impractical. Time-
multiplexed coherent arrays [8, 9] also require 2D arrays of free-space
detectors and yet only perform vector-vector multiplication in one clock
cycle.

(2) Monolithic Integration: Since any practical computing system will
eventually require high-density CMOS electronics for I/O and memory, we argue
that electronics and photonic should be monolithically integrated on a single
chip using processes such as GF 45CLO [10]. While 3D integration is typically
proposed as an alternative, existing advanced 3D integration technologies
(micro-bumps at $55\mu m$ pitch ) cannot provide density for optical computing
applications.

(3) Full Architecture-level Modeling and Optimization: Practical AI
accelerators, including optical processors, should be modeled at the system
level. This has been previously discussed in [11], however, accessing DRAM
through a PCIe switch will have large energy and latency overheads. We
elaborate on this aspect here, and model the system with co-packaged high-
bandwith memory (HBM), similar to state-of-the-art AI accelerators. In
addition, we discuss impacts and trade-offs between the programming time,
batch size, and multiple cores in this work.

##  III Proposed Crossbar Design for Optical MAC

The overall proposed crossbar ONN is shown in Fig. 2. Below, we describe two
key components of this design:

Figure 2: Photonic crossbar array with peripheral electronics for transmitter
and receiver

###  III-A Analog-based Optical MAC Core

The cross-bar design is an $N\times M$ array of PCM-based unit cells. Details
of these unit cells and how the array perform MAC is briefly described below.

####  III-A1 PCM-Based Unit Cell

Each unit cell multiplies an input electric field (E-field) by the weight
programmed into the PCM section, and adds this product to an externally-
inputted electric field through each column. To do so, a portion of light from
the row waveguide (E-field into each row is denoted by $|E_{in,i}|$ in Fig. 3)
is partially coupled into a bended waveguide via a directional coupler (DC)
with a cross-coupling ratio of $k_{in,j}$ (input coupling strength is column-
dependent). The portion of the E-field that does not couple into the unit cell
passes through a multi-mode interference (MMI) waveguide crossing junction and
enters the next column of unit cells to the right. Each bended waveguide has a
$\mu m$-long section covered with PCM. Individual PCM cells can be programmed
electrically to be either in the amorphous or crystalline state, or somewhere
in between, in a non-volatile fashion [7, 8]. Programming energy is estimated
to be around $100pJ$ [7, 8]. The state of PCM changes the absorption
coefficient, and hence it can change the amplitude of E-field. Consequently,
if the PCM’s programmed transmission in E-field domain is $w_{i,j}$, the
E-field at the end of each bended waveguide will be $|E_{in,i}|\times
k_{in,j}\times w_{i,j}$, which we refer to as $E_{p,(i,j)}$. This outcome will
be coupled into a column waveguide via another DC (with $k_{out,i}$ coupling
ratio), and it will travel into the next row in the bottom, where the coherent
summation with another row of products will occur through the DC region across
each column. The output coupling strength is row-dependent.

Figure 3: Principle of optical MAC operation in a coherent crossbar array (an
example of $2\times 2$ unit cell array).

####  III-A2 Crossbar Array of Unit Cells

In order to perform a MAC operation on our $N\times M$ crossbar array (Fig.
3), we need to encode the input data vector ($v_{in}$) on each row’s E-field
amplitude. To do so, we use a splitter tree structure and a set of optical
DACs (ODACs) to generate $v_{in,i}\times E_{Laser}/\sqrt{N}$ at each row.
Correctly designed input coupling coefficients in each unit cell ($k_{in,j}$)
distribute this equally throughout the row, so each unit cell receives
$v_{in,i}\times E_{Laser}/\sqrt{NM}$ in its bended waveguide.

Similarly, output coupling coefficients ($k_{out,i}$) are designed correctly
so that outputted light equally represents the product computed by every unit
cell in a given column. Although this adds an electrical field loss of
$1/\sqrt{N}$, this enables the entire array to operate on a single frequency
of light in a compact footprint (unlike previous works like [2]), which is
critical for scaling the array for high-performance compute. We can calculate
the resultant E-field at the end of each column as:

$E_{j}^{c}=\frac{E_{Laser}}{N\sqrt{M}}\sum_{i=0}^{N-1}|v_{in,i}|\times
w_{i,j}$ (1)

Thus, the overall $E^{c}$ vector will be equal to the matrix-vector
multiplication of input data ($v_{in}$) by a weights matrix of $w$. We can
convert this back into the electrical domain using coherent detection by
coupling each signal with a certain portion of the input laser light in a DC.
Coupler outputs enter balanced photodiodes (PD), with
$I_{out,j}\propto|E_{Laser}||E_{j}^{c}|$. We note that our scheme relies on
maintaining coherency across the entire array, which requires precise optical
path lengths and phase shift matchings. In order to adjust for potential phase
errors due to process variations or random phase variations, etc., we propose
adding a small thermal-phase shifter in each unit cell across the column
waveguides.

The photonic elements in our crossbar array present numerous sources of loss,
which are critical for system modeling:

  * •

Grating coupler: 2 dB [12, 10]

  * •

Splitting tree: 0.8 dB [13]

  * •

MMI Crossing Junction - 1.8 dB/junction [14]

  * •

Waveguide loss: 3 dB/cm [10]

  * •

Effective loss due to ODAC optical modulation amplitude (OMA) - 4 dB [15]

  * •

Laser wall-plug efficiency: 15%

###  III-B Peripheral Electronic Circuitry

While proposed crossbar array can perform optical MAC operation efficiently in
a compact footprint with on-chip PCM-based unit cells, the premise of ultra-
fast operations relies on energy-efficient and high-speed electro-optical
conversions. Here we elaborate on suitable design choices for this aim. We
have added area and power estimates of critical peripheral electronics for a
10GHz MAC operation. Note that some assumptions may involve extrapolation from
papers reporting work in other processes or operation frequencies.

####  III-B1 Optical On-chip Transmitter

As mentioned, the amplitude of input E-field at each row should be modulated
with the input vector data. Notice that coherent operation required that phase
remain constant with data. This has been previously proposed to be done via
MZIs, however, the mm-long footprint of MZIs impose unacceptable area energy
overheads for scalable ONN applications. In this design, we propose using
micro-ring assisted MZIs (RAMZIs). RAMZI devices has been previously proposed
for high linearity [16], but here we propose using them to perform constant-
phase PAM modulation. This can be done by placing one ring-resonator based
ODAC in each arm. Such ODACs have been demonstrated in the 45nm silicon
photonics with +20GS/s at low power with up to 6-bit accuracy for data
modulation [15]. Each ODAC driver consumes $168fJ$ power and occupies a
$0.0012mm^{2}$ area at 10GS/s [15]. We also add a thermal tuning over head of
$0.72mW$ per ring-resonator [15].

####  III-B2 Optical On-chip Coherent Receiver

Output photocurrents ($I_{out,j}$) should be amplified via a trans-impedance
amplifier (TIA) and digitized using an ADC (Fig. 3). Similar coherent receiver
designs have been demonstrated in 45nm silicon photonics with $2.25mW$ per TIA
[17]. ADC power and area at a 10 GHz sample rate is estimated to be $25mW$ and
$0.0475mm^{2}$, respectively, in 45nm CMOS [18].

####  III-B3 SerDes and Clocking

Operating the ONN crossbar at high clock rates of +10GHz requires a set of
serializer/deserializer (SerDes) to interface the optical transmitter and
receiver data with the digital backend (SRAM in this case). The serialization
ratio depends on the speed of MAC operation and readout speed from SRAM
memory. For a 10GHz ONN operation and a $\sim$1GHz backend clock-rate, we
assume a 10:1 ratio. Estimated power is roughly $100fJ$ per bit [15]. The
transmit and receive blocks require high speed clock generation and
distribution, which we assume will take $200fJ$ power with $0.005mm^{2}$ per
row/column [15].

##  IV Overall Accelerator Architecture

Figure 4: Architecture of proposed photonic crossbar array AI accelerator.

In this work we focus on CNNs for our modeling and benchmarking, as they are
one of the most popular deep learning algorithms due to their strong
performance in image classification and object detection [19]. Performing
inference using a CNN includes multiple convolutional layers using various
sets of filters. Our crossbar array capitalizes on the data reuse inherent to
the algorithm by minimizing on- and off-chip data access and executes a CNN as
follows: The weights in a given 2D filter are flattened and embedded into the
crossbar array in the form of PCM programming weights. Since the PCM can only
absorb light, all the weights are mapped to a value from between 0 and 1 over
64 levels (6-bits). Once the array has been programmed, the crossbar and
peripheral electronics iterate through all subsets of input data (across all
batches) and complete proper MAC operations. For any layer, the array will
likely have to be programmed multiple times to fully fit all the filters. An
accumulator is used at the output of the ADC and deserializer to hold partial
sums and add them to new partial sums from the crossbar array. Once complete
MAC operations have been computed, they are passed to an activation unit which
performs a non-linear transformation to achieve the final output (Fig. 4).

Data can be stored on chip in SRAM blocks (one each for the input, filters,
and output) for rapid and low-power access. SRAM area is estimated to be
$0.45mm^{2}$ per $1MB$ in 45nm CMOS [20]. On-chip SRAM cannot store the entire
datasets and all network parameters for practical CNNs. This necessitates an
off-chip DRAM memory, which typically consumes up to $15pJ/b$ [21] and can
dominate the total accelerator power. Here, similar to state-of-the-art AI
accelerators, we assume that our chip-scale system can be co-packaged with HBM
DRAM stack using advanced packaging [22] (fig. 4), which only consumes
$3.9pJ/bit$ [21]. In order to minimize DRAM access, data can be sent directly
from output SRAM to input SRAM at the end of a full layer computation. For
modeling, we assumed an SRAM access energy of $50fJ/bit$ [20].

Like many non-volatile memories, the PCM in this work has from low programming
speed, which can impact crossbar array performance. For instance, PCM
programming time can be around $100ns$ [7], which is $1000\times$ slower than
our target 10GHz MAC operation speed. Since the array cannot perform any
computations during PCM programming, this has the potential to significantly
slow down the effective inference/MAC speed. We propose a ”dual core” design
featuring two copies of all photonic elements (crossbar array, optical
transmitter, and optical receiver) to mitigate this issue. While one core
performs computation, the other core can be programmed with the next set of
weights. Assuming that programming the array takes less time than running all
input data through it, this effectively hides the programming latency, so that
computation runs at all times. We note that this solution can be practical
only if total area of photonic components will be a small portion of total
area. A single laser source can be shared among two crossbar cores. We discuss
this scheme further in Section VI.

##  V Simulation Framework

Figure 5: Simulation framework for estimating IPS and IPS/W.

We modeled our proposed system architecture comprehensively to measure key
performance metrics including inferences per second (IPS), inferences per
second per watt (IPS/W), and total chip area based on the 45nm monolithic
silicon photonics technology. For this aim, we developed a custom simulation
framework based on SCALE-sim [23] as shown in Fig. 5. There are two main steps
in estimating the key performance metrics values:

(1) Calculating the runtime specs: This includes the number of SRAM/DRAM
accesses, crossbar programming cycles, and number of MAC compute cycles.

(2) Calculating high-level metrics: This is estimated using the runtime specs
as well as the specs of all chip components - programming time, photonic
losses, ADC power, DRAM access power, etc. (See Sections III and IV)

Calculating runtime specs must be done for a specific neural network and
dataset. We achieve this using SCALE-Sim, a cycle-accurate tool for simulating
CNN accelerators [23]. We made slight modifications to SCALE-Sim to account
for non-unity batch sizes, the presence of an accumulator, and output SRAM
data reuse.

In step two of our framework, the number of array programming cycles and MAC
compute operation cycles are combined with MAC operation speed and programming
time to calculate inference time and corresponding IPS. Total chip power is
calculated based on clock speed, number of SRAM and DRAM accesses (from step
one), plus the energy of optical and electrical components. Chip area is
estimated similarly from SRAM sizes, digital accumulator and activation block
areas, photonic devices footprints, and area of peripheral electronics,
including ADCs, DACs, and clocking.

##  VI System-Level Performance Trends and Optimization

Our holistic simulation framework enables us to search the many-dimensional
chip design space to find the optimal configuration for ResNet50 inference. We
do so by studying the trends in terms of IPS and IPS/W for various design
parameters. Here we present observed critical trends and our approach to
choosing the optimal design. While we investigate the effect of array size,
batch size, dual core scheme, and SRAM size, we hold a constant MAC operation
speed of $10GHz$, since the power of many of the peripheral electronics
(specifically ADCs) may rise steeply at higher speeds.

Figure 6: IPS/W as a function of crossbar rows and columns.

###  VI-A Trends in Array Performance

Major trends are described below. All figures and data referenced in this
section assume the following default chip parameters except when clearly being
swept: Array size: 32 x 32, SRAM sizing: 26.3 MB (input), 0.75 MB (output),
0.75 MB (filter), and 0.75 MB (accumulator), dual-core architecture, and batch
size: 32.

####  VI-A1 Dual vs. Single Core

Dual core design increases the IPS, but power consumption is also consistently
higher since computing and programming happen simultaneously. As a result,
IPS/W is the same regardless of the core count.

####  VI-A2 Array Size

IPS always increases approximately linearly with the array size ($N\times M$).
This is no surprise because a larger array can process larger
matrices/vectors, so less programming and compute cycles are required to reach
the final sum. However, the effect of array size on IPS/W is more complicated.
As the array size increases, the power of peripheral electronics (such as ODAC
drivers or ADCs) grows. Also, higher laser power is required for larger arrays
because the power input to each row must spread out to more columns and rows.
However, both of these effects increase the power less than linearly since
many components are unaffected (such as DRAM accesses). However, photonic
losses in the crossbar array, such as those due to MMI waveguide crossings or
waveguides themselves, grow exponentially. Ultimately this causes power to
grow more than linearly with array size, so IPS/W drops. We observe a peak
IPS/W at array size of 128-256 rows and 64-128 columns (Fig. 6).

####  VI-A3 Batch Size and SRAM Size

Running CNNs with large batch sizes offers the possibility of amortizing the
cost of array programming time over a larger set of input data, increasing
IPS. However, if not all of the input data across all batches can fit on the
input SRAM, data will routinely have to be loaded on and off the chip every
time the array is reprogrammed. This effect can be seen in 7a, where DRAM
access energy rises steeply between batch sizes of 32 and 64. For any batch
size there is a critical input SRAM size that must be met for optimal IPS/W,
and any larger input SRAM size does not increase the performance further (Fig.
7b).

It is important to note that the area overhead of increasing SRAM size is
significant due to the chip area limitations. While smaller batch size is more
favorable from this perspective becuase it reduces the need for large SRAM, it
can limit IPS. This can be alleviated using our dual-core scheme which can
hide the array’s programming latency, so IPS becomes only a function of
computation time. Fig. 7c shows the large increase in IPS a dual-core
architecture can have for small-batch systems.

Figure 7: (a, b) Effect of the SRAM size and batch size on chip power and
IPS/W, (c) Effect of dual-core architecture on relationship between batch size
and IPS.

###  VI-B Proposed Approach for Optimization

Building on these observed trends, we develop the following flow for
optimizing system performance: Firstly, we find the smallest batch size that
is large enough to avoid a low IPS caused by a large programming time that
cannot be hidden by our dual core architecture. We do not need to be overly
concerned with increasing the batch size beyond this minimal amount because it
will have minimal effect on IPS/W, even for infinitely large SRAM (Fig. 7b).
Next, we maximize the SRAM size without exceeding a practical chip size ($\sim
1cm^{2}$ in here), assuming that SRAM dominates chip area, validated in Fig.
8. Finally, we find the the array size that maximizes IPS/W. When there are
multiple array sizes with similar IPS/W, we pick that largest array as it
yields a higher IPS, even though exact optimal size may depend on the CNN.

##  VII Results and Comparison to State-of-the-Art

Using our optimization framework, we arrive at the following system parameters
for optimal performance of our ONN, with dual-core architecture and a $10GHz$
MAC operation speed.

  * •

Array size: 128 rows $\times$ 128 columns

  * •

SRAM sizes: 26.3 MB (input), 0.75 MB (output), 0.75 MB (filter), and 0.75 MB
(accumulator)

  * •

Batch Size: 32

The performance of this system is compared with the state-of-the-art NVIDIA
A100 GPU (INT8 mode with a batch of 128) for ResNet50 in the table below [24].
Our proposed system achieves similar IPS with $15.4\times$ lower power and
$7.24\times$ lower area. Overall power and area breakdowns are illustrated in
Fig. 8. While the power consumption is dominated by DRAM accesses, the area is
mainly dominated by the SRAM blocks.

System | IPS | IPS/W | Power | Area
---|---|---|---|---
This work | 36382 | 1196 | 30W |  121 $mm^{2}$
Nvidia A100 | 29,733 | 75 | 396W |  826 $mm^{2}$
Figure 8: Power and area breakdown of our proposed accelerator.

##  VIII Conclusion

We have presented an optical AI accelerator using the coherent crossbar design
with on-chip PCM weight storage. The key performance metrics of this design
are holistically modeled and estimated using practical and realistic reported
results based on 45nm monolithic silicon photonic technology. We have
discussed our modeling and system optimization framework. Our approach can be
used for optimizing other types of emerging AI accelerators as well.
Furthermore, we proposed a dual core system design that can hide programming
latency. Our results show that overall system performances greater than the
state-of-the-art can be achieved. This work can be a first step towards
designing and optimizing large-scale ONN accelerators for real-world
applications by considering all the crucial aspects including SRAM/DRAM
memory, electronic-photonic integration, packaging, and all peripheral
electronics.

## References

  * [1] D. Amodei and D. Hernandez, “AI and Compute. Blog post, OpenAI,” 2018.
  * [2] Y. Shen, N. C. Harris _et al._ , “Deep learning with coherent nanophotonic circuits,” _Nature Photonics_ , vol. 11, no. 7, pp. 441–446, 2017.
  * [3] X. Xu, M. Tan _et al._ , “11 TOPS photonic convolutional accelerator for optical neural networks,” _Nature_ , vol. 589, no. 7840, pp. 44–51, 2021\.
  * [4] N. Wang, J. Choi _et al._ , “Training Deep Neural Networks with 8-bit Floating Point Numbers,” in _Advances in Neural Information Processing Systems_ , S. Bengio, H. Wallach _et al._ , Eds., vol. 31. Curran Associates, Inc., 2018.
  * [5] B. Darvish Rouhani, D. Lo _et al._ , “Pushing the Limits of Narrow Precision Inferencing at Cloud Scale with Microsoft Floating Point,” in _Advances in Neural Information Processing Systems_ , H. Larochelle, M. Ranzato _et al._ , Eds., vol. 33. Curran Associates, Inc., 2020, pp. 10 271–10 281.
  * [6] X. Lin, Y. Rivenson _et al._ , “All-optical machine learning using diffractive deep neural networks,” _Science_ , vol. 361, no. 6406, pp. 1004–1008, 2018.
  * [7] J. Feldmann, N. Youngblood _et al._ , “Parallel convolutional processing using an integrated photonic tensor core,” _Nature_ , vol. 589, no. 7840, pp. 52–58, 2021.
  * [8] N. Youngblood, “Coherent Photonic Crossbar Arrays for Large-Scale Matrix-Matrix Multiplication,” _IEEE Journal of Selected Topics in Quantum Electronics_ , pp. 1–1, 2022.
  * [9] R. Hamerly, A. Sludds _et al._ , “Large-Scale Optical Neural Networks based on Photoelectric Multiplication,” _ArXiv_ , vol. abs/1812.07614, 2019\.
  * [10] M. Rakowski, C. Meagher _et al._ , “45nm CMOS — Silicon Photonics Monolithic Technology (45CLO) for Next-Generation, Low Power and High Speed Optical Interconnects,” in _2020 Optical Fiber Communications Conference and Exhibition (OFC)_ , 2020, pp. 1–3.
  * [11] C. Demirkiran, F. Eris _et al._ , “An electro-photonic system for accelerating deep neural networks,” _CoRR_ , vol. abs/2109.01126, 2021.
  * [12] Y. Luo, Z. Nong _et al._ , “Low-loss two-dimensional silicon photonic grating coupler with a backside metal mirror,” _Opt. Lett._ , vol. 43, no. 3, pp. 474–477, Feb 2018.
  * [13] A. Zanzi, A. Brimont _et al._ , “Compact and low-loss asymmetrical multimode interference splitter for power monitoring applications,” _Opt. Lett._ , vol. 41, no. 2, pp. 227–229, Jan 2016.
  * [14] Y. Ma, Y. Zhang _et al._ , “Ultralow loss single layer submicron silicon waveguide crossing for soi optical interconnect,” _Opt. Express_ , vol. 21, no. 24, pp. 29 374–29 382, Dec 2013.
  * [15] S. Moazeni, S. Lin _et al._ , “A 40-Gb/s PAM-4 Transmitter Based on a Ring-Resonator Optical DAC in 45-nm SOI CMOS,” _IEEE Journal of Solid-State Circuits_ , vol. 52, no. 12, pp. 3503–3516, 2017.
  * [16] J. Cardenas, h. A. Morton _et al._ , “Linearized silicon modulator based on a ring assisted Mach Zehnder inteferometer,” _Optics Express_ , vol. 21, no. 19, p. 22549, 2013.
  * [17] N. Mehta, S. Lin _et al._ , “A laser-forwarded coherent 10gb/s bpsk transceiver using monolithic microring resonators in 45nm soi cmos,” in _2019 Symposium on VLSI Circuits_ , 2019, pp. C192–C193.
  * [18] M. Zhang, Y. Zhu _et al._ , “16.2 a 4$\times$ interleaved 10GS/s 8b time-domain ADC with 16$\times$ interpolation-based inter-stage gain achieving 37.5dB SNDR at 18GHz input,” in _2020 IEEE International Solid-State Circuits Conference-(ISSCC)_. IEEE, 2020, pp. 252–254.
  * [19] J. Gu, Z. Wang _et al._ , “Recent advances in convolutional neural networks,” _Pattern recognition_ , vol. 77, pp. 354–377, 2018.
  * [20] G. K. Chen, “Power management and sram for energy-autonomous and low-power systems,” Ph.D. dissertation, University of Michigan, 2011.
  * [21] M. O’Connor, N. Chatterjee _et al._ , “Fine-grained DRAM: Energy-efficient DRAM for extreme bandwidth systems,” in _Proceedings of the 50th Annual IEEE/ACM International Symposium on Microarchitecture_ , 2017, pp. 41–54.
  * [22] K. Hosseini, E. Kok _et al._ , “5.12 Tbps Co-Packaged FPGA and Silicon Photonics Interconnect I/O,” in _2022 IEEE Symposium on VLSI Technology and Circuits (VLSI Technology and Circuits)_ , 2022, pp. 260–261.
  * [23] A. Samajdar, Y. Zhu _et al._ , “Scale-sim: Systolic cnn accelerator simulator,” _arXiv preprint arXiv:1811.02883_ , 2018.
  * [24] “Nvidia data center deep learning product performance,” Sep 2022.


# (H)airborne:
two spinning black holes balanced by their synchronised scalar hair

Carlos A. R. Herdeiro, Eugen Radu  Departamento de Matemática da Universidade
de Aveiro and CIDMA, Campus de Santiago, 3810-183 Aveiro, Portugal

(May 2023)

###### Abstract

General Relativity minimally coupled to a massive, free, complex scalar field,
is shown to allow asymptotically flat solutions, non-singular on and outside
the event horizon, describing two spinning black holes (2sBHs) in equilibrium,
with co-axial, aligned angular momenta. The 2sBHs configurations bifurcate
from solutions describing dipolar spinning boson stars. The BHs emerge at
equilibrium points diagnosed by a test particle analysis and illustrated by a
Newtonian analogue. The individual BH “charges” are mass and angular momentum
only. Equilibrium is due to the scalar environment, acting as a (compact)
dipolar field, providing a lift against their mutual attraction, making the
2sBHs (h)airborne. We explore the 2sBHs domain of solutions and its main
features.

Introduction. The equilibrium 2-body problem provides a theoretical window
into the non-linear interactions between compact objects, in particular black
holes (BHs), in General Relativity (GR).

Vacuum BHs – the gravitational atoms – have only mass and angular momentum, as
macroscopic degrees of freedom Carter (1971); Robinson (1975); Chrusciel  _et
al._ (2012), measured at their horizon, $(M_{H},J_{H})$. It turns out that
such “atoms” cannot be balanced in vacuum GR, without introducing naked
singularities Costa  _et al._ (2009); Neugebauer and Hennig (2012); Hennig and
Neugebauer (2011). Placing them in an external global gravitational field, on
the other hand, can indeed balance gravitational atoms locally Astorino and
Vigano (2021), but at the cost of losing asymptotic flatness, in fact creating
asymptotic singularities 111The latter pathology is circumvented by replacing
the external field by a positive cosmological constant Dias  _et al._ (2023);
these solutions are asymptotically de Sitter, rather than asymptotically
flat.. This suggests, however, that the right smooth, local environment could
balance gravitational atoms 222Balance can also be achieved by endowing the
individual BHs with extra macroscopic repulsive degrees of freedom - gauge
charges - $Q_{H}$. The simplest such balanced configuration is the Majumdar-
Papapetrou Majumdar (1947); Papapetrou (1948) solution of electrovacuum. Here
we focus on neutral BHs, with $(M_{H},J_{H})$ as the only charges measured at
the horizon.. This is corroborated by a recent example Herdeiro and Radu
(2023), wherein, however, the environment was made up of unphysical “matter”.

In this letter we show that equilibrium gravitational molecules \- two BHs
with only $(M_{H},J_{H})$ individual degrees of freedom and asymptotically
flat -, exist in a smooth local environment, as exact solutions of GR
minimally coupled to simple and physical “matter”. More generically, we unveil
a mechanism, and lay out a methodology, to construct other families of such
non-linear molecules, prone to different applications.

The model and the environment. Consider the Einstein-(massive, complex) scalar
model described by the action $\mathcal{S}=(4\pi)^{-1}\int
d^{4}x\sqrt{-g}\mathcal{L}$, with 333We use geometrized units where $G=1=c$.

$\displaystyle\mathcal{L}=\frac{R}{4}-g^{\alpha\beta}\Phi_{,\,\alpha}^{*}\Phi_{,\,\beta}-\mu^{2}\Phi^{*}\Phi\
,$ (1)

where $R$ is the Ricci scalar, $g_{\alpha\beta}$ the spacetime metric, with
determinant $g$, $\Phi$ is a complex scalar field with mass $\mu$, and the
asterisk denotes complex conjugation.

Model (1), for appropriate ranges of $\mu$, is a popular fuzzy dark matter
model, to which different scalar self-interactions can be added Suárez  _et
al._ (2014); Hui  _et al._ (2017). One class of solutions are scalar
(mini-)boson stars (BSs) Schunck and Mielke (2003), a family of self-
gravitating solitons, often described as macroscopic Bose-Einstein condensates
Liebling and Palenzuela (2012), that can have a multipolar structure, akin to
hydrogenic orbitals Herdeiro  _et al._ (2021). Here we focus on dipolar
spinning boson stars (DsBSs) Kleihaus  _et al._ (2008); Kunz  _et al._ (2019)
\- see also Yoshida and Eriguchi (1997a); Cunha  _et al._ (2022) for the
static case - with the metric

$ds^{2}=-e^{2F_{0}}dt^{2}+e^{2F_{1}}\left(dr^{2}+dz^{2}\right)+e^{2F_{2}}\rho^{2}(d\varphi-
Wdt)^{2},$ (2)

where $\\{\rho,z,\varphi\\}$ are cylindrical coordinates and $F_{i},W$ are
functions of $(\rho,z)$, $i=0,1,2$. The geometry of DsBSs is
$\mathbb{Z}_{2}$-even, $F_{i}(\rho,-z)=F_{i}(\rho,z)$, $W(\rho,-z)=W(\rho,z)$,
with an equatorial plane at $z=0$. On the other hand, the scalar field,
$\Phi=\phi(\rho,z)e^{i(m\varphi-\omega t)}$, is $\mathbb{Z}_{2}$-odd,
$\phi(\rho,z)=-\phi(\rho,-z)$, where $\omega>0$ and $m\in\mathbb{Z}$ are the
frequency and azimuthal harmonic index. The $\pi$ phase difference between the
north and south hemispheres yields a repulsive scalar interaction between the
two constituents of the DsBSs Palenzuela  _et al._ (2007); Cunha  _et al._
(2022). DsBSs can be interpreted as two individual spinning BSs balanced by
their short (long) range scalar repulsion (gravitational attraction). The
domain of existence of the DsBSs is shown in the inset of Fig. 2 in an ADM
mass $vs.$ frequency diagram 444In the following, dimensionful quantities,
such as frequency $\omega$ and mass $M$ are presented in units of $\mu$.. The
stars exist for $0.6835\lesssim\omega\leqslant 1$ and the diagram shows the
characteristic spiral shape of other bosonic star solutions.

Bifurcation points. Individual spinning mini-BSs are mass tori in GR Schunck
and Mielke (1998); Yoshida and Eriguchi (1997b). To gain intuition, consider a
mass torus in Newtonian gravity. The simplest case is an infinitely thin
torus: a ring of radius $\mathcal{R}$ and constant mass density $\chi$. In a
cylindrical chart $(\rho,z,\phi)$ on $\mathbb{R}^{3}$ 555I.e. the spatial
metric (2) with $F_{i}=W=0$., take the ring on the $z=z_{r}\geqslant 0$ plane
and centred at $\rho=0$; its gravitational potential at $(\rho,z,\phi)$ is

$\Psi_{z_{r}}^{\rm
ring}=-\int_{0}^{2\pi}\frac{\chi\mathcal{R}d\varphi}{\sqrt{\rho^{2}-2\rho\mathcal{R}\cos(\phi-\varphi)+\mathcal{R}^{2}+(z-z_{r})^{2}}}\,.$
(3)

The gravitational force ${\bf F}=-\nabla\Psi$ has a $z$-component always
directed towards the plane of the ring (for particles outside that plane). Its
magnitude is not monotonic. Along $\rho=0$, it attains a maximum at a critical
distance $|\tilde{z}-\tilde{z}_{r}|=\tilde{z}^{\rm crit}\equiv{1}/{\sqrt{2}}$,
where $\tilde{z}\equiv z/\mathcal{R}$. This turns out to determine the
equilibrium points of a 2-rings system, a Newtonian analogue to the DsBSs
solutions - Fig. 1.

Figure 1: Equilibrium points in: (left) a Newtonian 2-rings system (diamonds)
and the corresponding gravitational potential along the $z$-axis; (right) the
GR DsBSs system (disks) and a corresponding (illustrative) metric $F_{0}$
function. $L$ ($\Delta_{\rm BS}$) will denote the proper distance between the
horizons (BSs).

Consider two parallel thin rings, both with radius $\mathcal{R}$,
symmetrically placed with centres at $\rho=0$, on the planes $z=\pm z_{r}$ \-
Fig. 1 (top left). The corresponding Newtonian potential is $\Psi^{\rm
2\,rings}_{-z_{r},z_{r}}=\Psi_{-z_{r}}^{\rm ring}+\Psi_{z_{r}}^{\rm ring}$.
The number of equilibrium points on the $z$-axis now depends on the
dimensionless ratio $\tilde{z}_{r}$: for $\tilde{z}_{r}<\tilde{z}_{r}^{\rm
crit}$, there is a single equilibrium point, at $z_{0}^{\rm eq}=0$. This is
stable against vertical displacements, since a test particle displaced from
the origin is attracted to the furthest ring; for
$\tilde{z}_{r}>\tilde{z}_{r}^{\rm crit}$, 2 new equilibrium points emerge,
$0<|z_{\pm}^{\rm eq}|<|z_{r}|$. They bifurcate from $z_{0}^{\rm eq}=0$ and are
symmetric with respect to $z=0$ 666These are determined by $(z_{\pm}^{\rm
eq}-z_{r})[\mathcal{R}^{2}+(z_{\pm}^{\rm eq}+z_{r})^{2}]^{3/2}+(z_{\pm}^{\rm
eq}+z_{r})[\mathcal{R}^{2}+(z_{\pm}^{\rm eq}-z_{r})^{2}]^{3/2}=0$. \- see an
illustrative $\Psi$ in Fig. 1 (bottom left). They are stable against vertical
displacements, whereas $z_{0}^{\rm eq}=0$ now becomes unstable. In the limit
$\tilde{z}_{r}\gg 1$, the new equilibrium points tend to the centre of each of
the rings, $z_{\pm}^{\rm eq}\rightarrow\pm z_{r}$. The upshot is that the
emergence of the new equilibrium points relies on the existence of a maximum
of the vertical force of each ring.

Let us now turn to the (timelike) equilibrium points on the DsBSs geometry,
$cf.$ eq. (2), of the form $(\rho,z,\varphi)=(0,z^{\rm
eq},\varphi_{0});(\dot{\rho},\dot{z},\dot{\varphi})=0=(\ddot{\rho},\ddot{z},\ddot{\varphi})$,
where derivatives are with respect to proper time. The geodesic equations,
$\ddot{x}^{\alpha}+\Gamma^{\alpha}_{\beta\gamma}\dot{x}^{\beta}\dot{x}^{\gamma}=0$,
yield the Christoffel symbols conditions
$\Gamma^{\rho}_{tt}=\Gamma^{z}_{tt}=\Gamma^{\varphi}_{tt}=0$, at the
equilibrium point. Considering the geometry (2), the last condition is obeyed
and the first follows from smoothness at the axis of the metric functions. The
remaining condition ($\Gamma^{z}_{tt}=0$) yields $\partial_{z}F_{0}=0$ at
equilibrium points. We scanned the parameter space of DsBSs Kunz  _et al._
(2019) \- Fig. 2 \- and observed always three equilibrium points, $z_{0}^{\rm
eq}=0$, $z_{\pm}^{\rm eq}\neq 0$, as in the Newtonian 2-rings system for
$\tilde{z}_{r}>\tilde{z}_{r}^{\rm crit}$ \- Fig. 1 (right). The fact that
$\Delta_{\rm BS}$ (Fig. 2) never nears zero may explain the absence of DsBSs
with a single equilibrium point.

Figure 2: Domain of existence of DsBSs (inset - red line) and the distance
between the two constituents $\Delta_{\rm BS}$ (main panel) $vs.$ $\omega$.
$\Delta_{\rm BS}$ is the proper distance along the $z$-axis between the
extrema of $\Phi$, attained, however, for $\rho\neq 0$, due to the toroidal
morphology. $L$ is the proper distance between $z_{\pm}^{\rm eq}$.

Static equilibrium points for a timelike test particle in a BH/soliton
spacetime suggest (but do not guarantee) a bifurcation towards a new family of
solutions, wherein the test particle is promoted to a back-reacting object. A
familiar example is the extremal Reissner-Nordström BH which bifurcates into
the Majumdar-Papapetrou multi-BH spacetimes, which can be diagnosed by the
equilibrium points of (charged) geodesics around the former. Another example,
closer to our discussion, are spinning BSs, that have an equilibrium point at
their centre. A spinning BH with horizon angular velocity $\Omega_{H}$ can be
added therein under the synchronisation condition Herdeiro and Radu (2014)

$w=m\Omega_{H}\ .$ (4)

The resulting hairy BHs interpolate between spinning BSs and the Kerr solution
Kerr (1963). On the other hand, the bifurcation from the latter end is
diagnosed by a different equilibrium: the existence of a zero mode of a test
bosonic field, at the threshold of superradiance Hod (2012).

We shall now show that two BHs can be grown from the equilibrium points of the
DsBSs. They can either both emerge from $z^{\rm eq}_{0}$ or one from each of
the two $z^{\rm eq}_{\pm}$ points. Equilibrium between each horizon and the
scalar environment requires condition (4). Equilibrium between the two
horizons relies on the extra opposing force of the two DsBS constituents,
themselves readjusting their distance $\Delta_{\rm BS}$ as the BHs are grown.
One anticipates that equilibrium requires BHs carrying a small enough fraction
of the total mass. This is confirmed below. A complementary perspective is
that each horizon-BS pair is a BH with synchronised scalar hair Herdeiro and
Radu (2014), with the repulsive force due to the $\pi$-phase difference
between the two bosonic clouds balancing the system.

Generalized Weyl framework. Axially symmetric, stationary, multi-BH solutions
in electrovacuum can be analytically constructed as Weyl solutions Weyl
(1917); Emparan and Reall (2002). At the heart of this construction is the
notion of rod structure Harmark (2004). In Herdeiro and Radu (2023), a
generalized Weyl construction for numerical implementation in models wherein
integrability is lost, such as $\eqref{action}$, was introduced, still based
on a rod-structure. We shall follow this construction, generalizing it for
stationary (rather than static) solutions.

In our case, the rod structure introduces 2 positive parameters: $\Delta z$,
(roughly) measuring the distance between horizons and $z_{H}$, determining the
horizon(s) size; it reads: $(i)$ a $1^{st}$ semi-infinite spacelike rod,
$[-\infty,-z_{H}]$, in the $\varphi$-direction; $(ii)$ a $1^{st}$ (finite)
timelike rod, $[-z_{H},-\Delta z/2]$; $(iii)$ a finite spacelike rod,
$[-\Delta z/2,\Delta z/2]$; $(iv)$ a $2^{nd}$ (finite) timelike rod $[\Delta
z/2,z_{H}]$; $(v)$ a $2^{nd}$ semi-infinite spacelike rod, $[z_{H},\infty]$.

This rod structure results in a set of boundary conditions on the $z$-axis,
$\rho=0$, conveniently expressed in terms of $f_{0}\equiv
e^{2F_{0}},f_{1}\equiv e^{2F_{1}},f_{2}\equiv\rho^{2}e^{2F_{2}}$. These are as
follows. For $m\neq 0$, $\phi=0$. For the timelike rods, corresponding to the
BH horizons, $f_{0}=\partial_{\rho}f_{1}=\partial_{\rho}f_{2}=0,W=\Omega_{H},$
together with $\lim_{\rho\to 0}\rho^{2}f_{1}/f_{0}=$const. and
$\partial_{\rho}W=0$. The horizon angular velocity of the two BHs is the same,
$\Omega_{H}=W(-z_{H}\leqslant z\leqslant-\Delta z/2)=W(\Delta z/2\leqslant
z\leqslant z_{H})$; the two BHs are corotating. As with the single BH case
Herdeiro and Radu (2014), assuming the existence of a power series expansion
of the solutions near a horizon results in (4). Each horizon has a spherical
topology, but deviating from a round sphere. They have the same area
$A_{H}=2\pi\int_{\Delta z/2}^{z_{H}}dz\sqrt{f_{1}(0,z)f_{2}(0,z)}$ and Hawking
temperature, $T_{H}=\lim_{\rho\to
0}\sqrt{{f_{0}(\rho,z)}/[\rho^{2}f_{1}(\rho,z)}]/2\pi$. Thus the horizons are
in thermodynamical equilibrium. The total entropy is twice that of a single
BH, $S=A_{H}/2$.

On a spacelike rod we impose (again at $\rho=0)$
$\partial_{\rho}f_{0}=\partial_{\rho}f_{1}=f_{2}=\partial_{\rho}W=0$;
$\lim_{\rho\to 0}\rho^{2}f_{1}/f_{2}=1$, additionally, imposes the absence of
conical singularities. The proper length of the finite $\varphi$-rod measures
the inter-horizons distance, $L=\int_{-\Delta z/2}^{\Delta
z/2}dz\sqrt{f_{1}(0,z)}$ \- Fig 1 (right).

Flat spacetime is approached for large $(\rho,|z|)$, wherein $\phi\rightarrow
0$. The ADM mass $M$ and angular momentum $J$ can be read off from the
asymptotic metric components: $-g_{tt}\simeq 1-{2M}/{\sqrt{\rho^{2}+z^{2}}}$,
$g_{\varphi t}\simeq-{2J\rho^{2}}/{\left(\rho^{2}+z^{2}\right)^{3/2}}$. As
usual in (asymptotically flat) BH mechanics, the temperature, entropy and the
global charges are related by a Smarr mass formula Townsend (1997),
$M=2T_{H}S+2\Omega_{H}(J-mQ)+M_{\Phi}$, with $M_{\Phi}$ the energy $outside$
the BHs stored in the scalar field. $Q$ is the conserved Noether charge
associated to the global $U(1)$ symmetry of (1) 777Physical quantities are
also connected via the first law, which, for balanced configurations, takes
the usual form $dM=T_{H}dS+\Omega_{H}dJ$.. The total mass is then
$M=M_{H}+M_{\Phi}$, with $M_{H}$ the (total) mass of the horizons, that can be
computed via Komar integrals Herdeiro and Radu (2015).

Finally, the vacuum, co-rotating double Kerr solution Kramer and Neugebauer
(1980) can be described by the above formalism, but with rather involved
expressions - see $e.g.$ Costa  _et al._ (2009); Cabrera-Munguia  _et al._
(2017); Manko and Ruiz (2017).

Numerics. We have solved numerically the set of five coupled non-linear
elliptic partial differential equations resulting from (1), with the ansatz
(2), subject to the above boundary conditions. For $\Phi$ we take $m=1$ 888We
have found solutions with $m=2$ as well, not reported here.. The input
parameters are $\\{\Omega_{H}\ {\rm or}\ w;\Delta z,z_{H}\\}$. The balanced
two spinning BHs (2sBHs) are constructed by scanning the domain of existence
of general unbalanced solutions at constant $\Omega_{H}$, and varying $z_{H}$
and $\Delta z$. Additionally fixing (say) $z_{H}$, regular configurations may
exist for some discrete values of $\Delta z$ only 999The accurate construction
of solutions with large frequencies $\omega\gtrsim 0.9$ remains a numerical
challenge.. Details on the approach, in particular on the coordinates better
suited for numerics, can be found in Herdeiro and Radu (2023). All solutions
reported here are free of conical singularities, and no pathologies were
observed on and outside the event horizon. Typical numerical errors are $\sim
10^{-3}$.

Results. DsBSs have no timelike rods: $z_{H}=\Delta z/2$. $\omega$ is the only
input parameter ($\Delta z$ is arbitrary). Taking $\Delta z=0,z_{H}\neq 0$
adds a single BH horizon at $z_{0}^{\rm eq}=0$ for any DsBS along the spiral
in Fig 2, yielding a family of BHs with parity-odd synchronised hair, studied
in Kunz  _et al._ (2019).

Any DsBS appears to also bifurcate into a 2sBH configuration, obeying (4). As
for the single BH case Kunz  _et al._ (2019), increasing the horizon size -
via $z_{H}$ -, 2BHs emerge from the seed DsBS at $z_{0}^{\rm eq}=0$. Along the
branch with the corresponding fixed frequency $\omega=\Omega_{H}$, the two BHs
initially increase in size ($A_{H}$). They depart from the equatorial plane -
$L$ grows - becoming (h)airborne. The sequence of solutions can terminate in
two possible ways. Away from the minimal frequency, $\omega\gtrsim 0.735$, $L$
never vanishes again and the sequence ends on the $z_{\pm}^{\rm eq}\neq 0$
equilibrium points of the same DsBS - this is illustrated by two sequences
labelled $(i)$ in Fig. 3. Close to the minimal frequency,
$0.6835\lesssim\omega\lesssim 0.735$,

Figure 3: Horizon area $A_{H}$ $vs.$ inter-horizon distance $L$ for
qualitatively different sequences of (balanced) 2sBHs. The inset shows the
mass fraction stored in the scalar field.

$L$ does vanish again and the sequence ends on the $z_{0}^{\rm eq}=0$
equilibrium point of another DsBS (with the same $\omega$ \- Fig. 2) - this is
labelled sequence $(iia)$ in Fig. 3. Additionally, these very same two DsBSs
can also be connected by a sequence departing from the $z_{\pm}^{\rm eq}\neq
0$ equilibrium points of one and arriving at the analogous points of the
other; here $L$ never vanishes along the whole sequence, labelled $(iib)$ in
Fig. 3.

Fig. 3 (inset) also reveals that the horizon mass is never larger than a few
percent - in our scanning we found a maximum of $\sim 8\%$ total horizon mass
(twice that of an individual BH) along a type $(i)$ sequence, seeded on the
upper DsBSs branch. The largest BHs are also the ones with the lowest Hawking
temperature - Fig. 4 (inset). The main panel of this figure also makes clear
that the two DsBSs connecting sequences $(iia)/(iib)$ are different, since
they have different masses.

Figure 4: Mass $M$ $vs.$ inter-horizon distance $L$ for the same sequences as
in Fig. 3. The inset shows the BH temperature.

Fig. 5 (top panel) provides a different illustration of the sequences
$(i),(iia)$ and $(iib)$. It zooms on the lower frequency end of the
$M$-$\omega$ domain of existence of DsBSs - inset of Fig. 2 -, since the
transition between sequences $(i)$ and $(ii)$ occurs in this frequency region.
Adding $L$ as a third dimension, makes transparent the sequences of 2sBHs
(black solid lines), branching off from the DsBSs at $z_{0}^{\rm eq}=0$ \-
from the red solid line - or at $z_{\pm}^{\rm eq}\neq 0$ \- from the blue
dotted line, which yields $L$. From the drawn sequences one can extrapolate
the surface describing the domain of existence of the 2sBHs: two tendentially
vertical (curved) surfaces for the larger $\omega$ that join together and
become two tendentially horizontal (curved) surfaces for the smaller $\omega$.
Additionally, one can conceive another surface, with $L=0$, bounded by the
DsBS line, describing the single BH solutions found in Kunz  _et al._ (2019).

Fig. 5 (bottom panel) shows a sequence of type $(i)$ in a horizon mass $M_{H}$
$vs.$ distance $L$ diagram, for a different (larger) frequency. It also
superimposes the horizon shape and size, computed as an isometric embedding in
$\mathbb{E}^{3}$ Smarr (1973). One can appreciate the growth of the horizon
size from both the $z_{0}^{\rm eq}=0$ and the $z_{\pm}^{\rm eq}\neq 0$
branching points. The deviation from sphericity of the horizons is small; the
individual horizons also possess a small $\mathbb{Z}_{2}$ asymmetry, better
seen in the heat map inset.

Figure 5: (Top) Illustrative sequences of 2sBHs branching off from the DsBHs.
(Bottom) Horizon embedding superimposed on an $(M_{H},L)$-diagram. The heatmap
shows the horizon geometry Ricci scalar for the maximal area solution (the
south one - in blue), which varies by $\sim$20% and is
$\mathbb{Z}_{2}$-asymmetric, maximized in the northern hemisphere.

Discussion. The asymptotically flat balanced 2sBHs described rely on two
ingredients: the synchronisation condition (4), allowing equilibrium between a
bosonic environment and a horizon; a bosonic environment whence the BHs
bifurcate with multiple (three) equilibrium points. These ingredients are also
present - possibly with a slightly different rationale - in other models
leading to natural generalization of these balanced 2sBHs solutions. Adding
scalar self-interactions to (1) allows a variety of models Schunck and Mielke
(2003). Judging by the healing properties seen in other contexts Di Giovanni
_et al._ (2020); Siemonsen and East (2021); Sanchis-Gual  _et al._ (2022),
self-interactions may improve the dynamical properties of these 2sBHs systems,
an interesting (but challenging) problem. Similar 2sBH configurations may also
exist for Proca fields, for which BHs with synchronised hair are also known
Santos  _et al._ (2020), but (currently) not dipolar Proca stars, static or
spinning.

A detailed study of the physical and mathematical properties of these balanced
2sBHs systems is in order and will be presented somewhere else.
Phenomenologically, since Eq. (2) describes a dark matter model, one may
entertain the question if there could be realistic dark matter distributions
that would bringing the mutual attraction of a dynamical BH binary to a halt.

Acknowledgements. We thank Pedro Cunha and Denjoe O’Connor for valuable
discussions. This work is supported by the Center for Research and Development
in Mathematics and Applications (CIDMA) through the Portuguese Foundation for
Science and Technology (FCT – Fundação para a Ciência e a Tecnologia),
references UIDB/04106/2020 and UIDP/04106/2020. The authors acknowledge
support from the projects CERN/FIS-PAR/0027/2019, PTDC/FIS-AST/3041/2020, as
well as CERN/FIS-PAR/0024/2021 and 2022.04560.PTDC. This work has further been
supported by the European Union’s Horizon 2020 research and innovation (RISE)
programme H2020-MSCA-RISE-2017 Grant No. FunFiCO-777740 and by the European
Horizon Europe staff exchange (SE) programme HORIZON-MSCA-2021-SE-01 Grant No.
NewFunFiCO-101086251. Computations were performed at the ARGUS cluster at
Aveiro University.

## References

  * Carter (1971) B. Carter, Phys. Rev. Lett. 26, 331 (1971).
  * Robinson (1975) D. C. Robinson, Phys. Rev. Lett. 34, 905 (1975).
  * Chrusciel  _et al._ (2012) P. T. Chrusciel, J. Lopes Costa, and M. Heusler, Living Rev. Rel. 15, 7 (2012), arXiv:1205.6112 [gr-qc] .
  * Costa  _et al._ (2009) M. S. Costa, C. A. R. Herdeiro, and C. Rebelo, Phys. Rev. D 79, 123508 (2009), arXiv:0903.0264 [gr-qc] .
  * Neugebauer and Hennig (2012) G. Neugebauer and J. Hennig, J. Geom. Phys. 62, 613 (2012), arXiv:1105.5830 [gr-qc] .
  * Hennig and Neugebauer (2011) J. Hennig and G. Neugebauer, Gen. Rel. Grav. 43, 3139 (2011), arXiv:1103.5248 [gr-qc] .
  * Astorino and Vigano (2021) M. Astorino and A. Vigano, Phys. Lett. B 820, 136506 (2021), arXiv:2104.07686 [gr-qc] .
  * Note (1) The latter pathology is circumvented by replacing the external field by a positive cosmological constant Dias  _et al._ (2023); these solutions are asymptotically de Sitter, rather than asymptotically flat.
  * Note (2) Balance can also be achieved by endowing the individual BHs with extra macroscopic repulsive degrees of freedom - gauge charges - $Q_{H}$. The simplest such balanced configuration is the Majumdar-Papapetrou Majumdar (1947); Papapetrou (1948) solution of electrovacuum. Here we focus on neutral BHs, with $(M_{H},J_{H})$ as the only charges measured at the horizon.
  * Herdeiro and Radu (2023) C. A. R. Herdeiro and E. Radu, Phys. Rev. D 107, 064044 (2023), arXiv:2302.00016 [gr-qc] .
  * Note (3) We use geometrized units where $G=1=c$.
  * Suárez  _et al._ (2014) A. Suárez, V. H. Robles, and T. Matos, Astrophys. Space Sci. Proc. 38, 107 (2014), arXiv:1302.0903 [astro-ph.CO] .
  * Hui  _et al._ (2017) L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Phys. Rev. D 95, 043541 (2017), arXiv:1610.08297 [astro-ph.CO] .
  * Schunck and Mielke (2003) F. E. Schunck and E. W. Mielke, Class. Quant. Grav. 20, R301 (2003), arXiv:0801.0307 [astro-ph] .
  * Liebling and Palenzuela (2012) S. L. Liebling and C. Palenzuela, Living Rev. Rel. 15, 6 (2012), arXiv:1202.5809 [gr-qc] .
  * Herdeiro  _et al._ (2021) C. A. R. Herdeiro, J. Kunz, I. Perapechka, E. Radu, and Y. Shnir, Phys. Lett. B 812, 136027 (2021), arXiv:2008.10608 [gr-qc] .
  * Kleihaus  _et al._ (2008) B. Kleihaus, J. Kunz, M. List, and I. Schaffer, Phys. Rev. D77, 064025 (2008), arXiv:0712.3742 [gr-qc] .
  * Kunz  _et al._ (2019) J. Kunz, I. Perapechka, and Y. Shnir, Phys. Rev. D 100, 064032 (2019), arXiv:1904.07630 [gr-qc] .
  * Yoshida and Eriguchi (1997a) S. Yoshida and Y. Eriguchi, Phys. Rev. D 55, 1994 (1997a).
  * Cunha  _et al._ (2022) P. Cunha, C. Herdeiro, E. Radu, and Y. Shnir, Phys. Rev. D 106, 124039 (2022), arXiv:2210.01833 [gr-qc] .
  * Palenzuela  _et al._ (2007) C. Palenzuela, I. Olabarrieta, L. Lehner, and S. L. Liebling, Phys. Rev. D 75, 064005 (2007), arXiv:gr-qc/0612067 .
  * Note (4) In the following, dimensionful quantities, such as frequency $\omega$ and mass $M$ are presented in units of $\mu$.
  * Schunck and Mielke (1998) F. E. Schunck and E. W. Mielke, Phys. Lett. A249, 389 (1998).
  * Yoshida and Eriguchi (1997b) S. Yoshida and Y. Eriguchi, Phys. Rev. D56, 762 (1997b).
  * Note (5) I.e. the spatial metric (2) with $F_{i}=W=0$.
  * Note (6) These are determined by $(z_{\pm}^{\rm eq}-z_{r})[\mathcal{R}^{2}+(z_{\pm}^{\rm eq}+z_{r})^{2}]^{3/2}+(z_{\pm}^{\rm eq}+z_{r})[\mathcal{R}^{2}+(z_{\pm}^{\rm eq}-z_{r})^{2}]^{3/2}=0$.
  * Herdeiro and Radu (2014) C. A. R. Herdeiro and E. Radu, Phys.Rev.Lett. 112, 221101 (2014), arXiv:1403.2757 [gr-qc] .
  * Kerr (1963) R. P. Kerr, Phys. Rev. Lett. 11, 237 (1963).
  * Hod (2012) S. Hod, Phys.Rev. D86, 104026 (2012), arXiv:1211.3202 [gr-qc] .
  * Weyl (1917) H. Weyl, Annalen Phys. 54, 117 (1917).
  * Emparan and Reall (2002) R. Emparan and H. S. Reall, Phys. Rev. D 65, 084025 (2002), arXiv:hep-th/0110258 .
  * Harmark (2004) T. Harmark, Phys. Rev. D 70, 124002 (2004), arXiv:hep-th/0408141 .
  * Townsend (1997) P. K. Townsend, (1997), arXiv:gr-qc/9707012 .
  * Note (7) Physical quantities are also connected via the first law, which, for balanced configurations, takes the usual form $dM=T_{H}dS+\Omega_{H}dJ$.
  * Herdeiro and Radu (2015) C. Herdeiro and E. Radu, Class. Quant. Grav. 32, 144001 (2015), arXiv:1501.04319 [gr-qc] .
  * Kramer and Neugebauer (1980) D. Kramer and G. Neugebauer, Physics Letters A 75, 259 (1980).
  * Cabrera-Munguia  _et al._ (2017) I. Cabrera-Munguia, V. E. Ceron, L. A. López, and O. Pedraza, Phys. Lett. B 772, 10 (2017), arXiv:1702.02209 [gr-qc] .
  * Manko and Ruiz (2017) V. S. Manko and E. Ruiz, Phys. Rev. D 96, 104016 (2017), arXiv:1702.05802 [gr-qc] .
  * Note (8) We have found solutions with $m=2$ as well, not reported here.
  * Note (9) The accurate construction of solutions with large frequencies $\omega\gtrsim 0.9$ remains a numerical challenge.
  * Smarr (1973) L. Smarr, Phys. Rev. D 7, 289 (1973).
  * Di Giovanni  _et al._ (2020) F. Di Giovanni, N. Sanchis-Gual, P. Cerdá-Durán, M. Zilhão, C. Herdeiro, J. A. Font, and E. Radu, Phys. Rev. D 102, 124009 (2020), arXiv:2010.05845 [gr-qc] .
  * Siemonsen and East (2021) N. Siemonsen and W. E. East, Phys. Rev. D 103, 044022 (2021), arXiv:2011.08247 [gr-qc] .
  * Sanchis-Gual  _et al._ (2022) N. Sanchis-Gual, C. Herdeiro, and E. Radu, Class. Quant. Grav. 39, 064001 (2022), arXiv:2110.03000 [gr-qc] .
  * Santos  _et al._ (2020) N. M. Santos, C. L. Benone, L. C. B. Crispino, C. A. R. Herdeiro, and E. Radu, JHEP 07, 010 (2020), arXiv:2004.09536 [gr-qc] .
  * Dias  _et al._ (2023) O. J. C. Dias, G. W. Gibbons, J. E. Santos, and B. Way, (2023), arXiv:2303.07361 [gr-qc] .
  * Majumdar (1947) S. Majumdar, Phys.Rev. 72, 390 (1947).
  * Papapetrou (1948) A. Papapetrou, Proc. Roy. Irish Acad. A 52, 11 (1948).


* Li and Leskovec (2022) Pan Li and Jure Leskovec.  The expressive power of graph neural networks.  _Graph Neural Networks: Foundations, Frontiers, and Applications_ , pages 63–98, 2022.
  * Lin and Cohen (2010) Frank Lin and William W Cohen.  Power iteration clustering.  In _ICML_ , 2010.
  * Lofgren et al. (2016) Peter Lofgren, Siddhartha Banerjee, and Ashish Goel.  Personalized pagerank estimation and search: A bidirectional approach.  In _Proceedings of the Ninth ACM International Conference on Web Search and Data Mining_ , pages 163–172, 2016.
  * Loveland et al. (2022) Donald Loveland, Jiayi Pan, Aaresh Farrokh Bhathena, and Yiyang Lu.  Fairedit: Preserving fairness in graph neural networks through greedy graph editing.  _arXiv preprint arXiv:2201.03681_ , 2022.
  * Lu et al. (2019) Lu Lu, Yeonjong Shin, Yanhui Su, and George Em Karniadakis.  Dying relu and initialization: Theory and numerical examples.  _arXiv preprint arXiv:1903.06733_ , 2019.
  * Ma et al. (2021) Jiaqi Ma, Junwei Deng, and Qiaozhu Mei.  Subgroup generalization and fairness of graph neural networks.  _Advances in Neural Information Processing Systems_ , 34:1048–1061, 2021.
  * Mastrandrea et al. (2015) Rossana Mastrandrea, Julie Fournet, and Alain Barrat.  Contact patterns in a high school: a comparison between data collected using wearable sensors, contact diaries and friendship surveys.  _PloS one_ , 10(9):e0136497, 2015.
  * McPherson et al. (2001) Miller McPherson, Lynn Smith-Lovin, and James M Cook.  Birds of a feather: Homophily in social networks.  _Annual review of sociology_ , 27(1):415–444, 2001.
  * Ntoutsi et al. (2020) Eirini Ntoutsi, Pavlos Fafalios, Ujwal Gadiraju, Vasileios Iosifidis, Wolfgang Nejdl, Maria-Esther Vidal, Salvatore Ruggieri, Franco Turini, Symeon Papadopoulos, Emmanouil Krasanakis, et al.  Bias in data-driven artificial intelligence systems—an introductory survey.  _Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery_ , 10(3):e1356, 2020.
  * Ortega et al. (2018) Antonio Ortega, Pascal Frossard, Jelena Kovačević, José MF Moura, and Pierre Vandergheynst.  Graph signal processing: Overview, challenges, and applications.  _Proceedings of the IEEE_ , 106(5):808–828, 2018\.
  * Papadopoulos et al. (2012) Symeon Papadopoulos, Yiannis Kompatsiaris, Athena Vakali, and Ploutarchos Spyridonos.  Community detection in social media.  _Data Mining and Knowledge Discovery_ , 24(3):515–554, 2012.
  * Peel et al. (2017) Leto Peel, Daniel B Larremore, and Aaron Clauset.  The ground truth about metadata and community detection in networks.  _Science advances_ , 3(5):e1602548, 2017.
  * Rahman et al. (2019) Tahleen A Rahman, Bartlomiej Surma, Michael Backes, and Yang Zhang.  Fairwalk: Towards fair graph embedding.  In _IJCAI_ , pages 3289–3295, 2019.
  * Sandryhaila and Moura (2013) Aliaksei Sandryhaila and José MF Moura.  Discrete signal processing on graphs: Graph filters.  In _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , pages 6163–6166. IEEE, 2013.
  * Schaeffer (2007) Satu Elisa Schaeffer.  Graph clustering.  _Computer science review_ , 1(1):27–64, 2007\.
  * Sen et al. (2008) Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad.  Collective classification in network data.  _AI magazine_ , 29(3):93–93, 2008.
  * Spielman (2012) Daniel Spielman.  Spectral graph theory.  In _Combinatorial scientific computing_ , number 18. Citeseer, 2012\.
  * Stoica et al. (2020) Ana-Andreea Stoica, Jessy Xinyi Han, and Augustin Chaintreau.  Seeding network influence in biased networks and the benefits of diversity.  In _Proceedings of The Web Conference 2020_ , pages 2089–2098, 2020\.
  * Takac and Zabovsky (2012) Lubos Takac and Michal Zabovsky.  Data analysis in public social networks.  In _International scientific conference and international workshop present day trends of innovations_ , volume 1. Present Day Trends of Innovations Lamza Poland, 2012.
  * Tong et al. (2006) Hanghang Tong, Christos Faloutsos, and Jia-Yu Pan.  Fast random walk with restart and its applications.  In _Sixth international conference on data mining (ICDM’06)_ , pages 613–622. IEEE, 2006.
  * Tsagris et al. (2014) Michail Tsagris, Christina Beneki, and Hossein Hassani.  On the folded normal distribution.  _Mathematics_ , 2(1):12–28, 2014.
  * Tsioutsiouliklis et al. (2020) Sotiris Tsioutsiouliklis, Evaggelia Pitoura, Panayiotis Tsaparas, Ilias Kleftakis, and Nikos Mamoulis.  Fairness-aware link analysis.  _arXiv preprint arXiv:2005.14431_ , 2020.
  * Tsioutsiouliklis et al. (2021) Sotiris Tsioutsiouliklis, Evaggelia Pitoura, Panayiotis Tsaparas, Ilias Kleftakis, and Nikos Mamoulis.  Fairness-aware pagerank.  In _Proceedings of the Web Conference 2021_ , pages 3815–3826, 2021\.
  * Tsioutsiouliklis et al. (2022) Sotiris Tsioutsiouliklis, Evaggelia Pitoura, Konstantinos Semertzidis, and Panayiotis Tsaparas.  Link recommendations for pagerank fairness.  In _Proceedings of the ACM Web Conference 2022_ , pages 3541–3551, 2022.
  * Vaccario et al. (2017) Giacomo Vaccario, Matúš Medo, Nicolas Wider, and Manuel Sebastian Mariani.  Quantifying and suppressing ranking bias in a large citation network.  _Journal of informetrics_ , 11(3):766–782, 2017\.
  * Whang et al. (2016) Joyce Jiyoung Whang, David F Gleich, and Inderjit S Dhillon.  Overlapping community detection using neighborhood-inflated seed expansion.  _IEEE Transactions on Knowledge and Data Engineering_ , 28(5):1272–1284, 2016.
  * Wu et al. (2012) Xiao-Ming Wu, Zhenguo Li, Anthony M So, John Wright, and Shih-Fu Chang.  Learning with partially absorbing random walks.  In _Advances in neural information processing systems_ , pages 3077–3085, 2012.
  * Xie et al. (2013) Jierui Xie, Stephen Kelley, and Boleslaw K Szymanski.  Overlapping community detection in networks: The state-of-the-art and comparative study.  _Acm computing surveys (csur)_ , 45(4):43, 2013\.
  * Xu et al. (2018) Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  _arXiv preprint arXiv:1810.00826_ , 2018.
  * Yang and Stoyanovich (2017) Ke Yang and Julia Stoyanovich.  Measuring fairness in ranked outputs.  In _Proceedings of the 29th International Conference on Scientific and Statistical Database Management_ , pages 1–6, 2017.
  * Yoshida (2019) Yuichi Yoshida.  Cheeger inequalities for submodular transformations.  In _Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms_ , pages 2582–2601. SIAM, 2019.
  * Zafar et al. (2019) Muhammad Bilal Zafar, Isabel Valera, Manuel Gomez-Rodriguez, and Krishna P Gummadi.  Fairness constraints: A flexible approach for fair classification.  _J. Mach. Learn. Res._ , 20(75):1–42, 2019.
  * Zehlike et al. (2017) Meike Zehlike, Francesco Bonchi, Carlos Castillo, Sara Hajian, Mohamed Megahed, and Ricardo Baeza-Yates.  Fa* ir: A fair top-k ranking algorithm.  In _Proceedings of the 2017 ACM on Conference on Information and Knowledge Management_ , pages 1569–1578, 2017.
  * Zhang et al. (2022) Tao Zhang, Tianqing Zhu, Mengde Han, Fengwen Chen, Jing Li, Wanlei Zhou, and Philip S Yu.  Fairness in graph-based semi-supervised learning.  _Knowledge and Information Systems_ , pages 1–28, 2022.
  * Zhou et al. (2003) Dengyong Zhou, Olivier Bousquet, Thomas Lal, Jason Weston, and Bernhard Schölkopf.  Learning with local and global consistency.  _Advances in neural information processing systems_ , 16, 2003.

##  A Towards Universal Approximation on the Node Domain

The theoretical analysis of this work reaches beyond fairness-aware graph
filtering and presents the first steps towards universal approximation
theorems in graph neural networks for node attribute prediction (e.g.,
classification). Recent research grapples with the expressive power of node
representations (Xu et al., 2018; Chen et al., 2019; Bécigneul et al., 2020;
Azizian and Lelarge, 2020), but this concept is weaker than the universal
approximation exhibited by traditional neural networks (Li and Leskovec,
2022). In this work, we express universal approximation of node scores in one
specific graph each time and guarantee (local) optimization under the
existence of both sufficiently good filters to be used and a full encoding of
graph signals involved in objective calculations in the priors.

Thanks to the above restrictions, our analysis can universally approximate
graph signal posterior objectives in the traditional sense. Furthermore, it
indicates that practices injecting additional (e.g., random) attributes to
graph neural network features do not necessarily need to capture node
identities, as Li and Leskovec (2022) hypothesize, but only classes of nodes
with different desired posteriors. For example, two nodes with the same
features or representations would be allowed if they also had the same
prediction targets, and despite whether this phenomenon would be a
manifestation of limited expressive power.

At the same time, we provide a first—to our knowledge—formal understanding of
why graph neural networks work so well compared to optimization that does not
acknowledge the graph structure. In theory, for tasks where data are abundant,
graph shift operators could be approximated via neural network components
(e.g., that mimic matrix decomposition) so that posterior objectives could be
optimized via traditional deep architectures, but this is rarely sufficient in
practice (Hoang et al., 2021). The discussion of Subsection 6.1 presents a
more delicate view than brute force optimization; the graph structure plays
the role of a selection mechanism that constraints the optimization to local
minima of only certain types. Thus, if there exist correlations between
structural characteristics captured by graph filters and deep minima, shallow
minima will be avoided by graph neural networks.

Kindred principles have led to the adoption of personalized pagerank under
symmetric adjacency matrix normalization in lieu of regularization mechanisms
that impose similar predictions between neighbors (Zhou et al., 2003; Huang et
al., 2020). Our work explains why this is needed in terms of optimization
theory; given that the imposed regularization corresponds to qualitative
predictions by itself, neural preprocessing can appropriately adjust its
predictions. Furthermore, we generalize to any kind of graph filter at the
role of the regularizer, given that it works as an adequately good predictor
by itself. The need for expressing qualitative predictions despite not
encountering these during training is also reflected in graph neural network
architectures like GCNII (Chen et al., 2020), which artificially impose low-
pass filtering by reducing the importance of many-hop graph diffusion.

Finally, a similar theory to ours is presented by Hoang et al. (2021), who
point out that existing graph neural network research focuses on low-pass
filters and datasets with small eigenvalues and that new approaches should
consider other variations too. In this work, we support generalized filtering,
but also discover pitfalls corresponding to small optimization horizon radii:
some types of high-pass filtering (e.g., combining high-pass filters and
graphs with small eigenvalues) admit too small radii for a small number of
nodes and are therefore only applicable to graphs with many nodes.

##  B Ablation Study

Subsection 5.4 demonstrates the ability of our proposed graph neural network
framework to act as a generic surrogate of fairness-aware graph filtering.
However, the NSGFF system is built by combining several components. In this
appendix, we conduct an empirical ablation study to showcase that each of
these components proactively circumvents issues that arise from misleadingly
similar architectures.

For demonstration purposes, we conduct a small scale comparison between our
approach and variations that give an idea of how things can go wrong when
optimizing the posterior objectives of symmetric filtering. The comparison is
made on a diffusion task defined on the Citeseer graph by assigning either
values $1$ or uniformly random values from the range $[0,1]$ to the first
community’s members and diffusing those. We investigate the following
variations of NSGFF:

  * Mult (baseline). The posterior rebalancing approach described in Subsection 5.2. This effectively removes the learning aspect of our approach.

  * NN (baseline). Directly editing posteriors by removing the graph propagation step of NSGFF to turn it into a multilayer perceptron. This strategy retains the same output activation functions, including those imposing fairness constraints, is trained on the same objective, and runs the same neural depth hyperparameter exploration.

  * APPNP (Klicpera et al., 2018). Adopting a fixed personalized pagerank filter with parameter $a=0.9$ and $10$ (instead of $20$ that we otherwise use in this work) propagation steps as NSGFF’s filter. Selecting a filter different than the one our framework aims to substitute is of unknown efficacy, especially since even the receptive breadth is different. Future works may be interested in investigating different filters of similar structural characteristics but greater optimization horizon radii, but this prospect lies outside the scope of our research.

  * NSGFF (this work). Our proposed system, as summarized in Subsection 5.2.

Results from making the particular diffusion task fairness-aware are presented
in Table 5. The full architecture does not always arrive at smaller utility
losses compared to cruder variations, but this is not unexpected, given that
the parameter $\delta_{0}$ controlling the optimization horizon radius is
selected from only a small range of options and may let original posteriors
lie within the radii of multiple local minima, including shallower ones.
Still, all utility losses of our approach reside in deep minima, i.e., they
are near-best compared to alternatives. Conversely, variations of our scheme
run the risk of ending up with much higher losses.

| None | Mult | NN | APPNP | NSGFF
---|---|---|---|---|---
| $\mathcal{L}_{util}$ | prule | $\mathcal{L}_{util}$ | prule | $\mathcal{L}_{util}$ | prule | $\mathcal{L}_{util}$ | prule | $\mathcal{L}_{util}$ | prule
Diffusion of random values
PPR0.85Sym | 0.000 | 0.138 | 0.689 | 1.000 | 0.623 | 1.000 | 0.467 | 1.000 | 0.356 | 1.000
PPR0.90Sym | 0.000 | 0.153 | 0.619 | 1.000 | 0.555 | 1.000 | 0.336 | 1.000 | 0.343 | 1.000
HK1Sym | 0.000 | 0.051 | 1.860 | 1.000 | 1.788 | 1.000 | 0.598 | 1.000 | 0.488 | 1.000
HK3Sym | 0.000 | 0.116 | 0.818 | 1.000 | 0.751 | 1.000 | 0.376 | 1.000 | 0.365 | 1.000
Diffusion of binary values
PPR0.85Sym | 0.000 | 0.134 | 0.720 | 1.000 | 0.656 | 1.000 | 0.399 | 1.000 | 0.369 | 1.000
PPR0.90Sym | 0.000 | 0.138 | 0.698 | 1.000 | 0.636 | 1.000 | 0.361 | 1.000 | 0.366 | 1.000
HK1Sym | 0.000 | 0.039 | 2.477 | 1.000 | 2.406 | 1.000 | 0.732 | 1.000 | 0.648 | 1.000
HK3Sym | 0.000 | 0.107 | 0.901 | 1.000 | 0.836 | 1.000 | 0.398 | 1.000 | 0.383 | 1.000
Table 5: Diffusion in the Citeseer graph. Best fairness-aware utility loss is
bolded.

##  C Proofs

In this appendix we present proofs of the theoretical results supporting our
work.

Lemma 2 Let $F(\hat{A})$ be a positive definite graph filter and $p$ a post-
processing vector. If $f(r)$ is a
$\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}$-optimizer of a
differentiable loss $\mathcal{L}(r)$ over graph signal domain
$\mathcal{R}\subseteq\mathbb{R}^{|\mathcal{V}|}$, where
$\lambda_{1},\lambda_{\max}>0$ are the smallest positive and largest
eigenvalues of $F(\hat{A})$ respectively, updating graph signals per the rule:

$\begin{split}&\frac{\partial q(t)}{\partial t}=f(r(t))\\\
&r(t)=\text{diag}(p)F(\hat{A})q(t)\end{split}$ (13)

asymptotically leads to the loss to local optimality if posterior updates are
closed in the domain, i.e. $r(t)\in\mathcal{R}\,\forall
t\in[0,T]\Rightarrow\text{diag}(p)F(\hat{A})\int_{0}^{T}f\big{(}r(t)\big{)}dt\in\mathcal{R}$.

Proof  For non-negative posteriors $r=\text{diag}(p)F(\hat{A})q(t)$, and non-
zero loss gradients, the Cauchy-Shwartz inequality in the bilinear space
$\langle x,y\rangle=x^{T}F(\hat{A})y$ determined by the positive definite
graph filter $F(\hat{A})$ yields:

$\displaystyle\frac{d\mathcal{L}\big{(}r(t)\big{)}}{dt}$
$\displaystyle\quad=\Big{(}\nabla\mathcal{L}\big{(}r(t)\big{)}\Big{)}^{T}\text{diag}(p)F(\hat{A})\frac{dq(t)}{dt}$
$\displaystyle\quad=\Big{(}\nabla\mathcal{L}\big{(}r(t)\big{)}\Big{)}^{T}\text{diag}(p)F(\hat{A})f\big{(}r(t)\big{)}$
$\displaystyle\quad=\Big{(}\nabla\mathcal{L}\big{(}r(t)\big{)}\Big{)}^{T}\text{diag}(p)F(\hat{A})\big{(}-\nabla\mathcal{L}\big{(}r(t)\big{)}+(f\big{(}r(t)\big{)}+\nabla\mathcal{L}\big{(}r(t)\big{)})\big{)}$
$\displaystyle\quad\leq-\min_{v}p[v]\Big{(}\nabla\mathcal{L}\big{(}r(t)\big{)}\Big{)}^{T}F(\hat{A})\nabla\mathcal{L}\big{(}r(t)\big{)}+(\nabla\mathcal{L}(r))^{T}\text{diag}(p)F(\hat{A})\Big{(}f\big{(}r(t)\big{)}+\nabla\mathcal{L}\big{(}r(t)\big{)}\Big{)}$
$\displaystyle\quad\leq-\lambda_{1}\min_{v}p[v]\|\nabla\mathcal{L}\big{(}r(t)\big{)}\|^{2}+\lambda_{\max}\big{\|}f\big{(}r(t)\big{)}+\nabla\mathcal{L}\big{(}r(t)\big{)}\big{\|}\|\text{diag}(p)\nabla\mathcal{L}\big{(}r(t)\big{)}\|$
$\displaystyle\quad<-\min_{v}p[v]\lambda_{1}\big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}\big{\|}^{2}+\lambda_{\max}\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}\big{\|}\max_{v}p[v]\big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}\big{\|}$
$\displaystyle\quad=0$

Therefore, the loss asymptotically converges to a locally optimal point.

Lemma 3 Let us consider graph signal priors $q_{0}$, positive definite graph
filter $F(\hat{A})$ with largest and smallest eigenvalues
$\lambda_{\max},\lambda_{1}$, post-processing vector $p$, differentiable loss
function $\mathcal{L}(r)$ with at least one minimum in domain
$\mathcal{R}\subseteq\mathbb{R}^{|\mathcal{V}|}$, and a differentiable graph
signal generation function
$\mathcal{M}:\mathbb{R}^{K}\to\mathbb{R}^{|\mathcal{M}|}$ for which there
exist parameters $\theta_{0}$ satisfying $\mathcal{M}(\theta_{0})=q_{0}$ and
$\text{diag}(p)F(\hat{A})\mathcal{M}(\theta)\in\mathcal{R}$. If, for any
parameters $\theta$, its Jacobian $\mathbb{J}_{\mathcal{M}}(\theta)$ has
linearly independent rows and satisfies:

$\displaystyle\big{\|}E(\theta)\nabla\mathcal{L}(r)\big{\|}<\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\|\nabla\mathcal{L}(r)\|\quad\text{
for }\nabla\mathcal{L}(r)\neq\textbf{0}$ $\displaystyle
E(\theta)=\mathcal{I}-\mathbb{J}_{\mathcal{M}}(\theta)(\mathbb{J}_{\mathcal{M}}^{T}(\theta)\mathbb{J}_{\mathcal{M}}(\theta))^{-1}\mathbb{J}_{\mathcal{M}}^{T}(\theta)$
$\displaystyle r(\theta)=\text{diag}(p)F(\hat{A})\mathcal{M}(\theta)$

then there exist parameters $\theta_{\infty}$ that make
$\mathcal{L}\big{(}r(\theta_{\infty})\big{)}$ locally optimal.

Proof  Let us consider a differentiable trajectory for parameters $\theta(t)$
for times $t\in[0,\infty)$ that starts from $\theta(0)=\theta_{0}$ and
(asymptotically) arrives at $\theta_{\infty}=\lim_{t\to\infty}\theta(t)$. For
this trajectory, it holds that
$\tfrac{d\mathcal{M}(\theta(t))}{dt}=\mathbb{J}_{\mathcal{M}}\big{(}\theta(t)\big{)}\frac{d\theta(t)}{dt}$.
Then, let us consider the posteriors
$r(t)=\text{diag}(p)F(\hat{A})\mathcal{M}\big{(}\theta(t)\big{)}$ arising from
the graph signal function $\mathcal{M}\big{(}\theta(t)\big{)}$ at times $t$,
as well the least square problem of minimizing the projection of the loss’s
gradient to the row space of
$\mathbb{J}_{\mathcal{M}}\big{(}\theta(t)\big{)}$:

$\text{ minimize
}\big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}-\mathbb{J}_{\mathcal{M}}\big{(}\theta(t)\big{)}x(t)\big{\|}$

The closed form solution to this problem can be found by:

$x(t)=\big{(}{\mathbb{J}_{\mathcal{M}}}^{T}\big{(}\theta(t)\big{)}\mathbb{J}_{\mathcal{M}}\big{(}\theta(t)\big{)}\big{)}^{-1}\mathbb{J}_{\mathcal{M}}^{T}\big{(}\theta(t)\big{)}\nabla\mathcal{L}\big{(}r(t)\big{)}$

Thus, as long as $q(t)=\mathcal{M}\big{(}\theta(t)\big{)}$ (Proposition I),
the theorem’s precondition for posteriors $r(t)=\text{diag}(p)F(\hat{A})q(t)$
and $\nabla\mathcal{L}\big{(}r(t)\big{)}\neq\textbf{0}$ can be written as:

$\big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}-\mathbb{J}_{\mathcal{M}}\big{(}\theta(t)\big{)}x(t)\big{\|}<\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\|\nabla\mathcal{L}\big{(}r(t)\big{)}\|$

Hence, if we set $x(t)$ as parameter derivatives:

$\tfrac{d\theta(t)}{dt}=x(t)\Rightarrow\Big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}-\tfrac{d\mathcal{M}(\theta(t))}{dt}\Big{\|}<\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\big{\|}\nabla\mathcal{L}\big{(}r(t)\big{)}\big{\|}$

which means that, on the selected parameter trajectory,
$\tfrac{\mathcal{M}(\theta(t))}{dt}$ is a
$\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}$-optimizer of the
loss function. As, such, from Lemma 2 it discovers locally optimal posteriors.

As a final step, we now investigate the priors signal editing converges at. To
do this, we can see that the update rule leads to the selection of priors
$q(t)$ at times $t$ for which:

$\displaystyle\tfrac{dq(t)}{dt}=\tfrac{d\mathcal{M}(\theta(t))}{dt}$
$\displaystyle\Rightarrow\lim_{t\to\infty}q(t)=q(0)+\int_{t=0}^{\infty}\tfrac{\mathcal{M}(\theta(t))}{dt}dt=\mathcal{M}\big{(}\theta(t_{\infty})\big{)}-\mathcal{M}\big{(}\theta(0)\big{)}=\mathcal{M}\big{(}\theta(t_{\infty})\big{)}$

We can similarly show that (Proposition I) holds true. Hence, there exists an
optimization path of prior editing parameters (not necessarily the same as the
one followed by the optimization algorithm used in practice) that arrives at
the edited priors $\mathcal{M}(\theta_{\infty})$ for some parameters
$\theta_{\infty}=\theta(t_{\infty})$ that let posteriors exhibit local
optimality.

Theorem 4 Let us consider graph signal priors $q_{0}$, positive definite graph
filter $F(\hat{A})$ with largest and smallest eigenvalues
$\lambda_{\max},\lambda_{1}$ respectively, post-processing vector $p$, and
twice differentiable loss function $\mathcal{L}(r)$, with some unique minimum
at unknown posteriors $r=r_{\infty}$ of a domain
$\mathcal{R}\subseteq\mathbb{R}^{|\mathcal{V}|}$, and whose Hessian matrix
$\mathbb{H}_{\mathcal{L}}(r)$ linearly approximates its gradients
$\nabla\mathcal{L}(r)$ within the domain with error at most
$\epsilon_{\mathbb{H}}$. Let us also construct a node feature matrix $H^{(0)}$
whose columns are graph signals (its rows $H^{(0)}[v]$ correspond to nodes
$v$). If $q_{0}$ and all graph signals involved in the calculation of
$\mathcal{L}(r)$ are columns of $H^{(0)}$, and it holds that:

$\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\|\nabla\mathcal{L}(r)\|>2\epsilon_{\mathbb{H}}\|r-r_{\infty}\|$

for any $r\in\mathcal{R}\setminus\\{r_{\infty}\\}$, then for any
$\epsilon_{\infty}>0$ there exists a deep neural network architecture

$H^{(\ell)}=\phi_{\ell}(H^{(\ell-1)}W^{(\ell)}+b^{(\ell)})$

with activation functions

$\phi_{\ell}(x)=\\{x\text{ for }\ell=L,relu(x)\text{ otherwise}\\}$

learnable weights $W^{(\ell)}$, and biases $b^{(\ell)}$ for layers
$\ell=1,2,\dots,L$ with depth $L\geq 4$, in which using the output of the last
layer $H^{(L)}$ as priors yields posteriors arbitrarily close to the local
minimum $\|diag(p)F(\hat{A})H^{(L)}-r_{\infty}\|<\epsilon_{\infty}$.
Additionally, layers other than the last can have output columns equal to the
number of $\text{columns of }H^{(0)}+2$.

Proof  Let $\mathbb{H}_{\mathcal{L}}(r)$ be the Hessian matrix of
$\mathcal{L}(r)$ with respect to node values of the graph signal
$r\in\mathcal{R}$ of a domain $\mathcal{R}$. Given that its rank is $K(r)$,
there exists a decomposition $\mathbb{H}_{\mathcal{L}}(r)=J(r)D(r)$ with fixed
inner dimensions $K_{\max}=\max_{r\in\mathcal{R}}K(r)$. When $K(r)=K_{\max}$,
the decomposition even becomes unique for the particular $r$.

We now select a prior editing/generation model $\mathcal{M}(\theta)$ with
linearly independent rows and a Jacobian
$\mathbb{J}_{\mathcal{\mathcal{M}}}(\theta)$ that is a deep neural network
approximating well the first of the decomposition elements. That is, for any
chosen small enough (but non-negligible) constant $\epsilon_{\mathbb{J}}>0$
within the domain $\mathcal{R}$ it holds that:

$\big{\|}\mathbb{J}_{\mathcal{\mathcal{M}}}(\theta)-J(r)\big{\|}\leq\epsilon_{\mathbb{J}}$

Throughout this proof, $r$ is a function of $\theta$, but we do not write
$r(\theta)$ for simplicity. The Jacobian does not need to have linearly
independent rows; small perturbations (with cumulative effect that, when added
to the actual approximation error, makes it smaller than
$\epsilon_{\mathbb{J}}$) can be injected in each row to make the rows of its
approximation linearly independent. Additionally, we use the Hessian as a
linear estimator of gradients $\nabla\mathcal{L}(r_{\infty})$ per:

$\|\nabla\mathcal{L}(r_{\infty})-\nabla\mathcal{L}(r)-\mathbb{H}_{\mathcal{L}}(r)(r_{\infty}-r)\|\leq\epsilon_{\mathbb{H}}\|r-r_{\infty}\|$

We now rewrite the loss’s gradient as a linear transformation of a factor
$\mathcal{L}^{\prime}(r)$ with an error term:

$\nabla\mathcal{L}(r)=\mathbb{J}_{\mathcal{\mathcal{M}}}(\theta)\mathcal{L}^{\prime}(r)+\epsilon(r)$

where
$\mathcal{L}^{\prime}(r)=\mathbb{J}_{\mathcal{\mathcal{M}}}^{T}(\theta)\big{(}\mathbb{J}_{\mathcal{\mathcal{M}}}(\theta)\mathbb{J}_{\mathcal{\mathcal{M}}}^{T}(\theta)\big{)}^{-1}\nabla\mathcal{L}(r_{\infty})+D(r)(r-r_{\infty})$,
and the error term admits a small bound:

$\|\epsilon(r)\|\leq\epsilon_{\mathbb{H}}\|r-r_{\infty}\|+\epsilon_{\mathbb{J}}\|D(r)\|\|r-r_{\infty}\|\leq\sup_{r\in\mathcal{R}}\|r-r_{\infty}\|(\epsilon_{\mathbb{H}}+\sup_{r\in\mathcal{R}}\|D(r)\|\epsilon_{\mathbb{J}})$

At this point, we choose
$\epsilon_{\mathbb{J}}=\tfrac{\epsilon_{\mathbb{H}}}{\sup_{r\in\mathcal{R}}\|D(r)\|}$,
which yields $\|\epsilon(r)\|\leq 2\epsilon_{\mathbb{H}}\|r-r_{\infty}\|$

Applying the linear estimation obtained above on the quantities of Lemma 3 we
obtain:

$\displaystyle\big{\|}E\nabla\mathcal{L}(r)\big{\|}=\|\nabla\mathcal{L}(r)-\mathbb{J}_{\mathcal{M}}(\theta)(\mathbb{J}_{\mathcal{M}}^{T}(\theta)\mathbb{J}_{\mathcal{M}}(\theta))^{-1}\mathbb{J}_{\mathcal{M}}^{T}(\theta)\big{(}\mathbb{J}_{\mathcal{M}}(\theta)\mathcal{L}^{\prime}(r)+\epsilon(r)\big{)}\big{\|}$
$\displaystyle=\|\nabla\mathcal{L}(r)-\mathbb{J}_{\mathcal{M}}(\theta)\mathcal{L}^{\prime}(r)+\mathbb{J}_{\mathcal{M}}(\theta)(\mathbb{J}_{\mathcal{M}}^{T}(\theta)\mathbb{J}_{\mathcal{M}}(\theta))^{-1}\mathbb{J}_{\mathcal{M}}^{T}(\theta)\epsilon(r)\|$
$\displaystyle=\|\epsilon(r)-\mathbb{J}_{\mathcal{M}}(\theta)(\mathbb{J}_{\mathcal{M}}^{T}(\theta)\mathbb{J}_{\mathcal{M}}(\theta))^{-1}\mathbb{J}_{\mathcal{M}}^{T}(\theta)\epsilon(r)\|=\big{\|}E\epsilon(r)\big{\|}$

However, the matrix
$\mathcal{I}-E=\mathbb{J}_{\mathcal{M}}(\theta)(\mathbb{J}_{\mathcal{M}}^{T}(\theta)\mathbb{J}_{\mathcal{M}}(\theta))^{-1}\mathbb{J}_{\mathcal{M}}^{T}$
is indempondent (that is, it holds
$(\mathcal{I}-E)(\mathcal{I}-E)=\mathcal{I}-E$) and therefore its eigenvalues
are either $0$ or $1$. Thus:

$\|E\nabla\mathcal{L}(r)\big{\|}\leq\|\epsilon(r)\|\leq
2\epsilon_{\mathbb{H}}\|r-r_{\infty}\|<\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\|\nabla\mathcal{L}(r)\|$

As a final step, and given that we showed the existence of a (deep neural
network) function, $\mathbb{J}_{\mathcal{\mathcal{M}}}(\theta)$, we integrate
it to retrieve $\mathcal{M}(\theta)$ and then approximate the integral with a
new deep neural network. Given the universal approximation theorem’s form
described by Kidger and Lyons (2020) the latter can be the architecture
described by this theorem, used to produce error less than
$\epsilon_{\infty}\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}$.
In this case:

$\displaystyle\|diag(p)F(\hat{A})H^{(L)}-r_{\infty}\|=\|diag(p)F(\hat{A})H^{(L)}-diag(p)F(\hat{A})\mathcal{M}(\theta_{\infty})\|$
$\displaystyle\quad\leq\|diag(p)F(\hat{A})\|\|H^{(L)}-\mathcal{M}(\theta_{\infty})\|<\|\epsilon_{\infty}\|$

We could not immediately assume the existence of that final network, because
we needed to verify that the Jacobian’s properties for the function it
approximates and “hide” methods of tracking equal partial derivatives of
different nodes behind the term $\epsilon_{\mathbb{J}}$.

Theorem 5 For any twice differentiable loss $\mathcal{L}(r)$, there exists
sufficiently large enough parameter
$l_{reg}\in\big{[}0,2\tfrac{\sup_{q}\|q-q_{0}\|}{|\mathcal{V}|}\big{]}$ such
that the loss regularization of Equation 10 satisfies the properties needed by
Theorem 4 within the optimization domain

$\mathcal{R}=\\{r_{\infty}\\}\cup\bigg{\\{}r:\|r-r_{\infty}\|<0.5\max\big{\\{}|\mathcal{V}|,\tfrac{\|\nabla\mathcal{L}(r)\|}{\epsilon_{\mathbb{H}}}\big{\\}}\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}\bigg{\\}}$

where $r_{\infty}$ is the ideal posteriors optimizing the regularized loss. If
the second term of the $\max$ is selected, it suffices to have no
regularization $l_{reg}=0$.

Proof  For ease of notation, during this proof we define
$w=0.5\tfrac{\lambda_{1}\min_{v}p[v]}{\lambda_{\max}\max_{v}p[v]}$. Without
loss of generality, consider $sgn$ to be a closure of the sign function on the
optimization trajectory. If $\tilde{\epsilon}_{\mathbb{H}}$ and
${\epsilon}_{\mathbb{H}}$ are the maximum derivative approximation error of
$\tilde{\mathcal{L}}(r)$ and $\mathcal{L}(r)$ respectively, it suffices to
select $l_{reg}>0$ such that:

$\displaystyle
w\|\nabla\mathcal{L}(r)\|+wl_{reg}|\mathcal{V}|>(\epsilon_{\mathbb{H}}+l_{reg})\sup_{r}\|r-r_{\infty}\|$
$\displaystyle\quad\Leftrightarrow
l_{reg}(w|\mathcal{V}|-\sup_{r}\|r-r_{\infty}\|)>\epsilon_{\mathbb{H}}\sup_{r}\|r-r_{\infty}\|-w\|\nabla\mathcal{L}(r)\|$
$\displaystyle\quad\Leftarrow|\mathcal{V}|w>\sup_{r}\|r-r_{\infty}\|\text{ or
}\epsilon_{\mathbb{H}}\sup_{r}\|r-r_{\infty}\|<w\|\nabla\mathcal{L}(r)\|$

where the first case needs to come alongside the condition

$l_{reg}>\tfrac{\sup_{r}\|r-r_{\infty}\|}{w|\mathcal{V}|-\sup_{r}\|r-r_{\infty}\|}\Leftarrow
l_{reg}\Leftarrow
l_{reg}=\tfrac{\sup_{r}\|r-r_{\infty}\|}{w|\mathcal{V}|}\geq\tfrac{2w\|q-q_{0}\|}{w|\mathcal{V}|}$

but in the second case it suffices to select any $l_{reg}\geq 0$. Given an
appropriate selection, we obtain the necessary criterion:

$w\|\tilde{\mathcal{L}}(r)\|\geq
w\|\nabla\mathcal{L}(r)\|+wl_{reg}|\mathcal{V}|>(\epsilon_{\mathbb{H}}+l_{reg})\sup_{r}\|r-r_{\infty}\|\geq\tilde{\epsilon}_{\mathcal{H}}\|r-r_{\infty}\|$

Theorem 6 Let us consider the case where priors and filter parameters are all
non-negative, and asymmetric adjacency matrix normalization yields positive
edge weights. Equation 12 can express locally optimal asymmetric filtering
posteriors, which can be found by applying the graph neural architecture of
this section on the equivalent symmetric filter to optimize the loss
$\mathcal{L}(f_{\delta}(r))$ from node features $H^{(0)}$ extended with a
fourth dimension $H^{(0)}[v,3]=r_{ns0}[v]$ and $6$ neural layer dimensions.
The optimization tightness of the symmetric filter is the same as if we
minimized
$\mathcal{L}\big{(}diag\big{(}\tfrac{p}{\delta+r_{0}}\big{)}F(\hat{\mathcal{A}})q\big{)}$.

Proof  We first define the receptive breadth of filters $F(\hat{A})$
characterized by real-valued non-negative weights $\\{f_{n}\geq
0|n=0,1,2,\dots\\}$ as the number of maximum possible number of hops away from
non-zero prior values in which non-zero posteriors can be found before post-
processing. We mathematically express this quantity as:

$\text{breadth}(diag(p)\mathcal{F}(\hat{A}))=\max_{n=0,1,2,\dots,dm}\\{n:f_{n}\neq
0\\}$

where $dm$ is the graph’s unweighted diameter. As long as the adjacency matrix
normalization does not assign a zero value to any edge, the breadth does not
depend on the exact type of normalization. Thus, for non-negative priors and
any post-processing vector $p$, it holds that the conditions
$\\{r_{ns0}[v]=0,r_{ns}[v],r_{0}[0]=0,f_{\delta}(r)[v]=0\\}$ can only
simultaneously hold true for node $v$ (this does not necessarily hold true for
$r[v]$, given that prior editing is allowed to procure negative values).

For the rest of nodes, and given that $\delta+r_{0}[v]>0$, there exist
appropriate values
$r[v]=\tfrac{r_{ns}[v]}{r_{ns0}[v]}(\delta+r_{0}[v])-\delta$ corresponding to
the optimal ones $r_{ns}[v]$. Importantly, for $r=r_{0}$ it also holds that
$r_{ns0}=r_{ns}$, which means that the search around the original symmetric
priors is translated to a search around the original asymmetric priors and not
at a random area in $\mathbb{R}^{|\mathcal{V}|}$.

As a final step, we rewrite the desired loss to be minimized as:

$\mathcal{L}_{ns}(r_{ns})=\mathcal{L}_{ns}\big{(}r_{ns0}\tfrac{\delta+r}{\delta+r_{0}}\big{)}=\mathcal{L}_{ns}\big{(}\tfrac{1}{\delta+r_{0}}r+\tfrac{r_{ns0}\delta}{\delta+r_{0}}\big{)}$

This is equivalent to applying the post-processing vector
$p^{\prime}=\tfrac{1}{\delta+r_{0}}p$ (where $p$ is the original post-
processing) and adding a constant additive term to $r$ within the loss
function. Given that the last additive term does not affect the supremum of
loss gradients, our analysis is applicable as-is, with the only difference in
our theoretical results being that attraction radii are equivalent to those
that would have been calculated for post-processing vectors
$\tfrac{p}{\delta+r_{0}}$.