# Multi-modal imaging using a cascaded microscope design

Xi Yang  Duke University, Department of Biomedical Engineering  Mark Harfouche
Ramona Optics, Inc.  Kevin C. Zhou  Duke University, Department of Biomedical
Engineering  Lucas Kreiss  Duke University, Department of Biomedical
Engineering  Shiqi Xu  Duke University, Department of Biomedical Engineering
Kanghyun Kim  Duke University, Department of Biomedical Engineering  Roarke
Horstmeyer  Duke University, Department of Biomedical Engineering

###### Abstract

We present a new _Multimodal Fiber Array Snapshot Technique (M-FAST)_ , based
on an array of 96 compact cameras placed behind a primary objective lens and a
fiber bundle array. which is capable of large-area, high-resolution, multi-
channel video acquisition. The proposed design provides two key improvements
to prior cascaded imaging system approaches: a novel optical arrangement that
accommodates the use of planar camera arrays, and the new ability to acquire
multi-modal image data acquisition. M-FAST is a multi-modal, scalable imaging
system that can acquire snapshot dual-channel fluorescence images as well as d
phase contrast measurements over a large 8 x 10 mm2 FOV at 2.2 $\mu$m full-
pitch resolution.

_K_ eywords Computational microscopy, Quantitative phase imaging, Fluorescence
microscopy

Large area, high-resolution microscope imaging plays an increasingly important
role in basic research, clinical workflows, and industrial inspection. Large
field-of-view (FOV) microscopes, capable of resolving more than hundreds of
millions pixels per snapshot, are needed to efficiently image whole-mount
tissue specimens [1], e.g., to monitor different neurons over an entire mouse
cortex [2, 3], and to observe freely moving organisms at cellular level detail
[4, 5, 6, 7]. Standard microscopes that were historically designed for the
human visual system and based on single objective lenses often cannot meet the
needs of such novel applications. Moreover, there is an increasing demand for
acquisition of _multi-modal_ image data, including multi-channel fluorescence,
phase contrast information, and polarimetric data, which is challenging to
integrate into traditional systems on a large scale. It is well-known[8] that
the number of resolved pixels per image, otherwise referred to as the system
space bandwidth product (SBP), is practically limited by two components within
a standard digital microscope: the microscope objective lens and the digital
image sensor. Over the past few years, there have been several computational
imaging-based approaches to increase an imaging system’s SBP to simultaneously
provide high resolution over a large FOV, including Fourier ptychography [9]
and its vectorial variation [10, 11], large-scale scanning light-sheet
microscopes [12] and multi-view microscopes[13]. While with high SBP, such
approaches typically must acquire multiple snapshots over time and cannot
resolve rapidly moving specimens or transient dynamics.

For single-snapshot high-SBP acquisition, explored strategies include direct
optimization of the objective lens itself [14] and the use of multiple imaging
systems, or "multi-scale" systems [15], that are arranged in a cascaded
fashion to acquire data in parallel. A recent work applied such parallelized
"multiscale" concepts [8] to create the RUSH system [3], a video-rate large-
SBP microscope with a custom-designed objective lens, which forms a magnified
spherical intermediate image. A number of individual compact sCMOS cameras
then re-image unique segments of this intermediate surface, and the resulting
data is stitched into a final composite. While the use of a curved
intermediate image plane in such prior designs avoids the introduction of
multiple sources of aberration [15], it also introduces several critical
engineering challenges, including tilted cameras overlapping in a common focal
plane design and tight camera packing density. As we will detail, these two
aspects present significant limitations to the goal of multi-modal acquisition
(e.g., the ability to capture multiple fluorescence emission wavelengths per
snapshot).

Figure 1: (a) Schematic of the M-FAST microscope. (b) Image of micro camera
array (top) and of fiber bundle array plate (bottom). (c) Captured 96
($8\times 12$) raw images of a lung histology sample with
$>$$50\text{\,}\mathrm{\char 37\relax}$ FOV overlap in one dimension (black
box), to be stitched in (d). (e) Large-FOV stitched composite, high resolution
composite of thin-tissue specimen.

In this work, we build upon such prior designs to produce a scalable, cascaded
microscope design that facilitates multi-modal image acquisition, including
both phase contrast as well as dual-channel fluorescence imaging. Our new
"Multimodal Fiber Array Snapshot Technique" (M-FAST) is composed of three main
components: a primary objective lens, a flat integrated micro-camera array
board [4] and a fiber bundle array plate. The primary objective lens (finite
conjugate lens: RMA Electronics, FL-CC0814-5M, 0.36 numerical aperture (NA) )
offers an optical SBP of several hundred megapixels [16] and magnifies the
sample onto an intermediate plane ($8\times 10\text{mm}^{2}$ FOV, $16\times
24\text{cm}^{2}$ magnified image size). The intermediate plane contains the
fiber faceplate array (FFA) (details below). Light emerging from the fiber
array is then recorded by a micro-camera array. Each micro-camera is
configured to focus on the intermediate image plane and capture data from a
unique FOV that overlaps with those of neighboring cameras, as sketched in
Fig. 1(a). The cameras’ magnification was designed to have an approximately
$10\text{\,}\mathrm{\char 37\relax}$ FOV overlap along the short side of the
rectangular sensor (Fig. 1(d) red box) and $52\text{\,}\mathrm{\char
37\relax}$ overlap along the longer dimension (Fig. 1(d) black box), ensuring
each point within the intermediate image is observed at least by two micro-
cameras for multi-channel acquisition.

Figure 2: (a) Photo of $8\times 12$ micro-camera array. (a1-a6) Bright-field
images from both center (first row) and corner (second row) of three unique
micro-camera images, each marked with colored box. (b1-b3) Bright-field images
of resolution target with fiber faceplate raw, time-averaging, and frame
averaging. (c1) and (c2) shows impressive vignetting alleviation with the FFA.
Red circle in (c1) shows the vignetting effect while being eliminated in (c2).
Translating the FFA $1\text{\,}\mathrm{mm}$ within a $0.2\text{\,}\mathrm{s}$
exposure time leads to the example time-averaged result in Fig. 2(c3). Figure
3: (a) Dual channel fluorescence imaging set up schematic showing
$55\text{\,}\mathrm{\char 37\relax}$ overlap along $x$ and alternating
red/green emission filters inserted over adjacent micro-cameras to
concurrently detect red/green fluorescence over entire FOV. (b) Snapshot dual-
channel fluorescence image acquired by the M-FAST ($8\times 10\text{mm}^{2}$
FOV) of red (b1)/green (b2) fluorescent microspheres with zoomed version in
(d)/(e). (d1)(e1) shows the raw results with speckle noise while (d2)(e2)
shows the time-averaged results. (c) Plot of required number of micro-cameras
as a function of desired number of multi-mode channels for standard M-FAST
parameters assuming a primary lens with 11cm by 11cm FOV (blue) and 22cm by
22cm FOV megapixel (orange) SBP.

A key component of our new design is the FFA, which is a thin, large plate of
fused glass single-mode fibers. The FFA accounts for the primary objective
lens’ non-telecentricity. Ideally, light cones at the intermediate image plane
formed by an image-side telecentric primary lens would align with the optical
axis of each micro-array camera, thus producing an image without vignetting.
Unfortunately, it is significantly challenging and expensive to design
microscope objective lenses with a high SBP and image-side telecentricity. To
decouple multiscale microscopes from such a challenging requirement, we
instead adopted a non-telecentric primary objective lens whose chief rays on
the image side are not parallel (see sketch in Fig. 1(a)). Recent works show
sub-micrometer resolution imaging[17], wide angle imaging[18] with a fiber
bundle based system. In addition, it has been established that speckle sensing
methods can be used to image dynamics [19, 20, 21]. In a configuration without
FFA, off-axis light from the primary lens is traveling at an angle that can
exceed the NA of micro-cameras when the micro-cameras are located in the
peripheral array. This leads to image signal reduced or vanished. By
introducing a FFA at the intermediate image plane, light is re-directed to
enter the lens of each single micro-camera of the array to form an image with
minimal vignetting. Fig. 2(c1)(c2) highlights the benefit of adopting the FFA.
Without the array, vignetting severely limits the observable micro-camera FOV
of this exemplary 10 $\mu$m fluorescent microsphere calibration specimen. In
the current version of our M-FAST microscope, we use a
$6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ size fiber faceplate (SCHOTT, glass
type: 47A) with 5 mm thickness and per-fiber NA=1. We placed four such arrays,
each $8\leavevmode\nobreak\ \times\leavevmode\nobreak\ 12$ cm in size,
directly against one another to cover the M-FAST intermediate image plane.
Fig. 1(b2) shows the $6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ fiber faceplate
under a $40\times$ microscope.

The micro-camera array contains 96 individual imagers ($8\times 12$ grid,
$19\text{\,}\mathrm{mm}$ pitch), each consisting of a rectangular CMOS sensor
(Omnivision 108238, 4320 x 2432 pixels with a pixel size of 1.4 um) that are
jointly controlled by a custom electronics arrangement to simultaneously
acquire image data from all sensors
($0.96\text{\,}\mathrm{g}\mathrm{i}\mathrm{g}\mathrm{a}\mathrm{p}\mathrm{i}\mathrm{x}\mathrm{e}\mathrm{l}\mathrm{s}$
acquired per snapshot). Each micro-camera utilizes a $25\text{\,}\mathrm{mm}$
focal length lens ($f/2.5$, Edmund Optics) configured to image at 0.14
magnification. Each micro-camera’s NA is 0.032, which was selected to
approximately match the image-side NA of the primary objective lens. Stitching
the images captured by all 96 micro-cameras (via PTGui software [22]) produces
a final full-FOV composite. Unlike prior multiscale microscope designs, the
M-FAST micro-cameras are arranged to capture FOVs that exhibit a large degree
of overlap with their immediate neighbors to facilitate snapshot multimodal
acquisition. In our current design, we aimed to capture two modalities per
snapshot: each point on the intermediate image plane is observed by at least
two micro-cameras. Within this configuration, the half-pitch resolution of the
micro camera array at the intermediate image plane was 9 $\mu$m, when directly
imaging a target at the intermediate image plane without a primary lens. The
complete M-FAST system has 3.81 magnification. The entire FOV of the complete
system is 8x10 mm. To assess the resolution of our M-FAST microscope. We
translated a USAF resolution target (Ready Optics) via a manual 3-axis stage
to different FOV areas. Fig. 2 (a) shows bright-field imaging results of the
resolution target within the FOV of three different micro-cameras with
uniquely colored outlines. At the center of the FOV of the central micro-
camera, we measured a $2.2\text{\,}\mathrm{\SIUnitSymbolMicro m}$ two-point
resolution, which drops to $2.76\text{\,}\mathrm{\SIUnitSymbolMicro m}$ at the
edge of the central micro-camera’s FOV. We observed a
$3.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and
$3.48\text{\,}\mathrm{\SIUnitSymbolMicro m}$ two-point resolution at the FOV
center and edge, respectively, for the micro-camera at the array edge.

There are two key issues that impact image fidelity within our unique optical
arrangement that we address using software post-processing. First, vignetting
effects introduced by the primary objective lens and camera array lenses lead
to a fixed intensity fall-off towards the FOV edges of each camera as well as
of the composite system. We calibrate the images based on a Gaussian-blurred
blank reference image. The second issue is the introduction of fixed, high-
spatial frequency modulation by the fiber faceplate, whose fibers (6 $\mu$m
average diameter) are smaller than the M-FAST system’s resolution at the
intermediate image plane (9 $\mu$m half-pitch resolution), yet still lead to a
speckle-like modulation effect, as observable in Fig. 2(b1). To remove these
modulation effects, We utilized a motorized stage to vibrate the fiber
faceplate [23] during the finite image exposure. We adopted two methods to
achieve fiber faceplate vibration: time averaging Fig. 2(b2) and frame
averaging Fig. 2(b3). For the time averaging method, we vibrated the fiber
faceplate within a fixed exposure time, while for frame averaging we captured
and computed the average of $10$ frames, each with the plate at a random
displacement ($1\text{\,}\mathrm{mm}3\text{\,}\mathrm{mm}$)(see Fig. 2(b3)).
In the future, calculating the transmission matrix of the fiber faceplate
could enable effective removal without moving parts in the system [24].

Figure 4: (a) M-FAST differential phase contrast configuration, with a
illumination LED source(I), condensers(C), sample(S) and primary lens(PL). The
FFA and micro-camera array are not shown here. (b) Bright field image of dense
collection of 10$\mu$m microspheres immersed in n=1.56 oil. DPC results shown
for zoomed-in areas at FOV center (b1)-(b2) and corner (b3)-(b4). (c) 1D
surface profiles of microspheres within DPC image (averaged over 20 beads) at
FOV corner (blue) and FOV center(red), along with the 1D simulated surface
profile from known sphere thickness profile. (d) DPC image of Oral epithelial
cells for zoomed-in areas at FOV center (d1)-(d2) and corner (d3)-(d4).

As a first demonstration of snapshot multi-modal acquisition, we configured
the M-FAST setup to acquire dual-channel fluorescence video over thin samples,
since thick sample causes the ambiguity in the lateral shift. As sketched in
Fig. 3(a), we inserted alternating red fluorescent emission filter
(ET610/75m(Chroma), 12 mm) and green fluorescent emission filters
(ET510/20m(Chroma), 12 mm) over successive micro-camera columns, allowing the
M-FAST to simultaneously capture both spectral channels across the full system
FOV (except at the left and right edges). As a proof-of-concept experiment, we
imaged a mixture of $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$ red
(590nm/605nm) and $6\text{\,}\mathrm{\SIUnitSymbolMicro m}$ yellow-green
microspheres (441nm/486nm). For fluorescence excitation, we used a high power
$470\text{\,}\mathrm{nm}$ blue LED (Chanzon, $100\text{\,}\mathrm{W}$ Blue)
with a short-pass excitation filter (Thorlabs, FES500). We stitched the two
fluorescence channels separately to produce the composites in Fig. 3(b1-b2).
M-FAST can clearly resolve the microsphere mixture, both spectrally and
spatially.

As presented above, two-channel imaging is achieved by inserting different
filters in neighboring cameras with an overlapping FOV. More channels can be
designed by increasing the amount of overlap between neighboring cameras (Fig.
3(a)). For four fluorescent channels, there are two options: (1) use
rectangular image sensors with 75$\%$ FOV overlap along one dimension, or (2)
using square sensors with 50$\%$ FOV overlap along both dimensions across the
object plane. In general, the amount of FOV overlap, and thus the number of
multimodal channels $c$ that an M-FAST system simultaneously captures, is a
function of four key parameters: the sensor size ($s$), magnification of the
micro-camera lens ($M$), desired number of cameras ($n$) and the total system
FOV($F$). Within a simplified 1D analysis, the number of cameras along $x$ is
given by $n_{x}=\frac{F_{x}M}{s_{x}}c_{x}-1$. The required camera array pitch
can be shown to be given by $\frac{p_{i}}{s_{i}}=\frac{1}{Mc_{i}}>1,i=x,y$. We
assume the resolution is limited by the camera pixel size $\delta_{x}$. When
giving a specific SBP image and a desired number of fluorescence channels, we
can get required number of cameras using
$n_{x}=\frac{SBP_{x}}{SBP_{cam,x}}c_{x}-1$ by substituting
$SBP_{x}=\frac{FM}{\delta_{x}},SBP_{cam,x}=\frac{s_{x}}{\delta_{x}}$.
Extension to a 2D camera array is achieved by multiplying the 1D expressions –
the total number of cameras is thus $n=n_{x}n_{y}$. Assuming square image
sensors and reasonable M-FAST parameters ($SBP_{cam}$ = 10 MP), we plot $n$
against $c$ in Fig. 3(c). The blue line and orange line shows the relationship
assuming 200-MP and 400-MP composite FOVs.

As a final demonstration of the flexibility of our novel multiscale
microscope, we modified the illumination unit to include the ability to
perform differential phase contrast (DPC) imaging [25, 26]. Here,a 16-mm focal
length condenser (ACL25416U, Thorlabs) and two subsequent condensers
(ACL7560U, Thorlabs) are used for magnifying the active area of the LED
source, shown in Fig. 4(a). We inserted an absorption mask at the illumination
source conjugate plane to produce DPC modulation[27]. In this design, the
specimen is located at the Fourier plane of the DPC mask [28]. To first
quantitatively evaluate M-FAST DPC performance, we imaged a large collection
of 10-µm diameter polystyrene microspheres (refractive index: $n=1.59$)
immersed in $n=1.58$ oil. To provide DPC modulation, we inserted four unique
$0.8$ NA half-circle masks oriented along four cardinal directions and
captured four snapshots. The resulting phase contrast maps are in Fig. 4(b),
with average 1D traces through ten randomly selected beads shown in Fig. 4(c),
both for beads at the micro-camera FOV center and edge. As reference we
analytically derived the phase contrast profile of the microsphere, captured
by a simplified single-objective DPC imaging with matching specifications. The
direct comparison between our experimental results and this analytical
reference highlights the accurate performance of our system in this DPC mode.
Next, we tested the effectiveness of M-FAST DPC imaging of buccal epithelial
cells samples, which are mostly transparent and thus exhibit minimal
structural contrast under a standard brightfield microscope. For this task, we
employed ring-shaped DPC masks with 0.7 outer NA and 0.5 inner NA. Fig. 4(d)
shows the zoom-ins of the associated DPC maps, which reveal detailed phase-
specific features within cells.

In summary, we have presented M-FAST, a new multi-modal, scalable imaging
technique featuring a novel cascaded microscope design based on a finite-
conjugate objective, a camera array and a fiber faceplate. This system
provides fast acquisition of large-FOV images without scanning parts and
without image vignetting. Moreover, our current configuration is versatile for
different imaging modalities and features two fluorescence channels and a DPC
module. We demonstrated that the system is scalable to increase the number of
simultaneously observable fluorescence channels by tuning the camera array
geometry. Our technique thus opens the door to high-speed, high-resolution,
large-FOV, multi-modal microscopy.

##  1 Funding

We acknowledge funding from NIH ORIP and NIH NIEHS for SBIR R44-OD24879, NIH
NCI for SBIR 1R44CA250877, NIH NIBIB for SBIR 1R43EB030979

##  2 Disclosures

MF: Ramona Optics Inc. (F,I,P,S), RH: Ramona Optics Inc. (F,I,P,S), KCZ:
Ramona Optics Inc. (C).

## References

  * [1] A. Y. Sdobnov, M. Darvin, E. Genina, A. Bashkatov, J. Lademann, and V. Tuchin, “Recent progress in tissue optical clearing for spectroscopic application,” _Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy_ , vol. 197, pp. 216–229, 2018.
  * [2] C. K. Kim, S. J. Yang, N. Pichamoorthy, N. P. Young, I. Kauvar, J. H. Jennings, T. N. Lerner, A. Berndt, S. Y. Lee, C. Ramakrishnan _et al._ , “Simultaneous fast measurement of circuit dynamics at multiple sites across the mammalian brain,” _Nature methods_ , vol. 13, no. 4, pp. 325–328, 2016\.
  * [3] J. Fan, J. Suo, J. Wu, H. Xie, Y. Shen, F. Chen, G. Wang, L. Cao, G. Jin, Q. He _et al._ , “Video-rate imaging of biological dynamics at centimetre scale and micrometre resolution,” _Nature Photonics_ , vol. 13, no. 11, pp. 809–816, 2019.
  * [4] E. Thomson, M. Harfouche, P. Konda, C. W. Seitz, K. Kim, C. Cooke, S. Xu, R. Blazing, Y. Chen, W. S. Jacobs _et al._ , “Gigapixel behavioral and neural activity imaging with a novel multi-camera array microscope,” _bioRxiv_ , 2021.
  * [5] D. Krishnamurthy, H. Li, F. B. du Rey, P. Cambournac, A. Larson, and M. Prakash, “Scale-free vertical tracking microscopy: Towards bridging scales in biological oceanography,” _bioRxiv_ , p. 610246, 2019.
  * [6] “Probabilistic models of larval zebrafish behavior reveal structure on many scales,” _Current Biology_ , vol. 30, no. 1, pp. 70–82.e4, 2020.
  * [7] J. P. Nguyen, F. B. Shipley, A. N. Linder, G. S. Plummer, M. Liu, S. U. Setru, J. W. Shaevitz, and A. M. Leifer, “Whole-brain calcium imaging with cellular resolution in freely behaving caenorhabditis elegans,” _Proceedings of the National Academy of Sciences_ , vol. 113, no. 8, pp. E1074–E1081, 2016.
  * [8] J. Park, D. J. Brady, G. Zheng, L. Tian, and L. Gao, “Review of bio-optical imaging systems with a high space-bandwidth product,” _Advanced Photonics_ , vol. 3, no. 4, p. 044001, 2021.
  * [9] P. C. Konda, L. Loetgering, K. C. Zhou, S. Xu, A. R. Harvey, and R. Horstmeyer, “Fourier ptychography: current applications and future promises,” _Optics express_ , vol. 28, no. 7, pp. 9603–9630, 2020.
  * [10] X. Dai, S. Xu, X. Yang, K. C. Zhou, C. Glass, P. C. Konda, and R. Horstmeyer, “Quantitative jones matrix imaging using vectorial fourier ptychography,” _Biomedical optics express_ , vol. 13, no. 3, pp. 1457–1470, 2022.
  * [11] S. Xu, X. Dai, X. Yang, K. Zhou, P. Konda, and R. Horstmeyer, “Imaging anisotropy with vectorial fourier ptychography,” in _Computational Optics 2021_ , vol. 11875. SPIE, 2021, p. 118750F.
  * [12] F. F. Voigt, D. Kirschenbaum, E. Platonova, R. A. Campbell, R. Kastli, M. Schaettin, L. Egolf, A. Van der Bourg, P. Bethge, K. Haenraets _et al._ , “The mesospim initiative: open-source light-sheet microscopes for imaging cleared tissue,” _Nature methods_ , vol. 16, no. 11, pp. 1105–1108, 2019.
  * [13] K. He, X. Wang, Z. W. Wang, H. Yi, N. F. Scherer, A. K. Katsaggelos, and O. Cossairt, “Snapshot multifocal light field microscopy,” _Optics Express_ , vol. 28, no. 8, pp. 12 108–12 120, 2020.
  * [14] G. McConnell, J. Trägårdh, R. Amor, J. Dempster, E. Reid, and W. B. Amos, “A novel optical microscope for imaging large embryos and tissue volumes with sub-cellular resolution throughout,” _Elife_ , vol. 5, p. e18659, 2016.
  * [15] D. J. Brady, M. E. Gehm, R. A. Stack, D. L. Marks, D. S. Kittle, D. R. Golish, E. Vera, and S. D. Feller, “Multiscale gigapixel photography,” _Nature_ , vol. 486, no. 7403, pp. 386–389, 2012.
  * [16] G. Zheng, X. Ou, and C. Yang, “0.5 gigapixel microscopy using a flatbed scanner,” _Biomedical optics express_ , vol. 5, no. 1, pp. 1–8, 2014.
  * [17] Y. Liu, Q.-Y. Yu, Z.-M. Chen, H.-Y. Qiu, R. Chen, S.-J. Jiang, X.-T. He, F.-L. Zhao, and J.-W. Dong, “Meta-objective with sub-micrometer resolution for microendoscopes,” _Photonics Research_ , vol. 9, no. 2, pp. 106–115, 2021\.
  * [18] I. Stamenov, A. Arianpour, S. J. Olivas, I. P. Agurok, A. R. Johnson, R. A. Stack, R. L. Morrison, and J. E. Ford, “Panoramic monocentric imaging using fiber-coupled focal planes,” _Optics express_ , vol. 22, no. 26, pp. 31 708–31 721, 2014.
  * [19] D. A. Boas and A. K. Dunn, “Laser speckle contrast imaging in biomedical optics,” _Journal of biomedical optics_ , vol. 15, no. 1, p. 011109, 2010\.
  * [20] S. Xu, X. Yang, W. Liu, J. Jonsson, R. Qian, P. C. Konda, K. C. Zhou, Q. Dai, H. Wang, E. Berrocal _et al._ , “Imaging dynamics beneath turbid media via parallelized single-photon detection,” _arXiv preprint arXiv:2107.01422_ , 2021.
  * [21] S. Xu, W. Liu, X. Yang, J. Jönsson, R. Qian, P. McKee, K. Kim, P. C. Konda, K. C. Zhou, L. Kreiß _et al._ , “Transient motion classification through turbid volumes via parallelized single-photon detection and deep contrastive embedding,” _Frontiers in neuroscience_ , vol. 16, 2022.
  * [22] New House Internet Services BV, “Ptgui.” [Online]. Available: https://ptgui.com/
  * [23] Y. Fujimaki and H. Taniguchi, “Reduction of speckle contrast in multimode fibers using piezoelectric vibrator,” in _Laser Resonators, Microresonators, and Beam Control XVI_ , vol. 8960. International Society for Optics and Photonics, 2014, p. 89601S.
  * [24] G. Huang, D. Wu, J. Luo, Y. Huang, and Y. Shen, “Retrieving the optical transmission matrix of a multimode fiber using the extended kalman filter,” _Optics express_ , vol. 28, no. 7, pp. 9487–9500, 2020.
  * [25] L. Tian and L. Waller, “Quantitative differential phase contrast imaging in an led array microscope,” _Optics express_ , vol. 23, no. 9, pp. 11 394–11 403, 2015.
  * [26] S. Xu, X. Dai, X. Yang, K. C. Zhou, K. Kim, V. Pathak, C. Glass, and R. Horstmeyer, “Tensorial tomographic differential phase-contrast microscopy,” _arXiv preprint arXiv:2204.11397_ , 2022.
  * [27] S. B. Mehta and C. J. Sheppard, “Quantitative phase-gradient imaging at high resolution with asymmetric illumination-based differential phase contrast,” _Optics letters_ , vol. 34, no. 13, pp. 1924–1926, 2009.
  * [28] C. Bonati, T. Laforest, M. Kunzi, and C. Moser, “Phase sensitivity in differential phase contrast microscopy: limits and strategies to improve it,” _Optics Express_ , vol. 28, no. 22, pp. 33 767–33 783, 2020.


11institutetext: Department of Statistics, Stanford University 11email:
<EMAIL_ADDRESS>
# FastAD: Expression Template-Based C++ Library for Fast and Memory-Efficient
Automatic Differentiation

James Yang

###### Abstract

Automatic differentiation is a set of techniques to efficiently and accurately
compute the derivative of a function represented by a computer program.
Existing C++ libraries for automatic differentiation (e.g. Adept, Stan Math
Library), however, exhibit large memory consumptions and runtime performance
issues. This paper introduces FastAD, a new C++ template library for automatic
differentiation, that overcomes all of these challenges in existing libraries
by using vectorization, simpler memory management using a fully expression-
template-based design, and other compile-time optimizations to remove some
run-time overhead. Benchmarks show that FastAD performs 2-10 times faster than
Adept and 2-19 times faster than Stan across various test cases including a
few real-world examples.

###### Keywords:

automatic differentiation forward-mode reverse-mode C++ expression templates
template metaprogramming lazy-evaluation lazy-allocation vectorization.

##  1 Introduction

Gradient computation plays a critical role in modern computing problems
surrounding optimization, statistics, and machine learning. For example, one
may wish to compute sensitivities of an Ordinary-Differential-Equation (ODE)
integrator for optimization or parameter estimation [3]. In Bayesian
statistics, advanced Markov-Chain-Monte-Carlo (MCMC) algorithms such as the
Hamiltonian Monte Carlo (HMC) and the No-U-Turn-Sampler (NUTS) rely heavily on
computing the gradient of a (log) joint probability density function to update
proposal samples in the leapfrog algorithm [7][10]. Neural networks rely on
computing the gradient of a loss function during back-propagation to update
the weights between each layer of the network [4].

Oftentimes, the target function to differentiate is extremely complicated, and
it is very tedious and error-prone for the programmer to manually define the
analytical formula for the gradient [9]. It is rather desirable to have a
generic framework where the programmer only needs to specify the target
function to differentiate it. Moreover, computing the derivative may be one of
the more expensive parts of an algorithm if it requires numerous such
evaluations, as it is the case for the three examples discussed above. Hence,
it is imperative that gradient computation is as efficient as possible. These
desires motivated the development of such a framework: automatic
differentiation (AD).

FastAD is a general-purpose AD library in C++ supporting both forward and
reverse modes of automatic differentiation. It is highly optimized to compute
gradients, but also supports computing the full Jacobian matrix. Similar to
the Stan Math Library, FastAD is primarily intended for differentiating scalar
functions, and it is well-known that reverse-mode is more efficient than
forward-mode for these cases [3]. For these reasons, this paper will focus
only on reverse-mode and computing gradients rather than forward-mode or
computing a full Jacobian matrix.

##  2 Overview

Section 3 will first explain the reverse-mode automatic differentiation
algorithm to give context for how FastAD is implemented. With this background,
Section 4 will discuss some implementation detail and draw design comparisons
with other libraries to see how FastAD overcomes the common challenges seen in
these libraries. In Section 5, we show a series of benchmarks with other
libraries including two real-world examples of a regression model and a
stochastic volatility model.

##  3 Reverse-Mode Automatic Differentiation

We briefly discuss the reverse-mode automatic differentiation algorithm for
context and to motivate the discussion for the next sections. For a more in-
depth treatment, we direct the readers to [3][9][5].

As an example, consider the function

$f(x_{1},x_{2},x_{3})=\sin(x_{1})+\cos(x_{2})\cdot x_{3}-\log(x_{3})$ (1)

This function can be represented as an expression graph where each node of the
graph represents a sub-expression. Fig. 1 shows the corresponding expression
graph drawn to evaluate in the same order as defined by the operator
precedence in the C++ standard.

$x_{1}$$x_{2}$$x_{3}$$\sin$$\cos$$\log$$+$$*$$-$ Figure 1: Expression graph
for Eq. 1

Note that, in general, a variable $x_{i}$ can be referenced by multiple nodes.
For example, $x_{3}$ is referenced by the * and log nodes. It is actually more
helpful to convert this expression graph into an expression tree by replacing
all such nodes with multiple parents as separate nodes that have a reference
back to the actual variable. Mathematically,

$\displaystyle f(x_{1},x_{2},x_{3})$
$\displaystyle=\tilde{f}(g(x_{1},x_{2},x_{3}))$ (2)
$\displaystyle\tilde{f}(w_{1},w_{2},w_{3},w_{4})$
$\displaystyle=\sin(w_{1})+\cos(w_{2})\cdot w_{3}-\log(w_{4})$ $\displaystyle
g(x_{1},x_{2},x_{3})$ $\displaystyle=(x_{1},x_{2},x_{3},x_{3})$

Fig. 2 shows the corresponding converted expression tree. This way all nodes
except possibly the $x_{i}$ nodes have exactly one parent and have a unique
path back up to the root, if we view the tree as a directed graph. While this
may seem more complicated than the original expression graph, implementation
becomes much cleaner with this new approach. With the converted expression
tree, we will see momentarily that we may actually start the algorithm from
$w_{i}$ instead of $x_{i}$. In fact, rather than treating $x_{i}$ as nodes of
the expression graph, it is more helpful to instead treat them as containers
that hold the initial values and their adjoints, $\frac{\partial f}{\partial
x_{i}}$. For this reason, we denote the path between $x_{i}$ and $w_{i}$ with
dotted lines, and consider only up to nodes $w_{i}$ in the expression tree.

$x_{1}$$x_{2}$$x_{3}$$w_{1}$$w_{2}$$w_{3}$$w_{4}$$\sin$$\cos$$\log$$+$$*$$-$
Figure 2: Converted expression tree for Eq. 2. Nodes $x_{1},x_{2},x_{3}$ are
now separated from the rest of the graph by a layer of $w$ variables. Note in
particular that $x_{3}$ is now replaced with $w_{3}$ and $w_{4}$ (red
boundary).

The reverse-mode algorithm consists of two passes of the expression graph:
_forward_ -evaluation (not to be confused with forward-mode AD), and
_backward_ -evaluation. During the _forward_ -evaluation, we compute the
expression in the usual fashion, i.e. start at the root, recursively forward-
evaluate from left to right all of its children, then finally take those
results and evaluate the current node. Evaluation of each node will compute
the operation that it represents. For example, after forward-evaluating the
$w_{1}$ node, which is a trivial operation of retrieving the value from the
container $x_{1}$, the sin node will return $\sin(w_{1})=\sin(x_{1})$.

The _backward_ -evaluation for a given node starts by receiving its adjoint
from its parent. We will also refer to this value as its _seed_ to further
distinguish $x_{i}$ from the expression nodes. Since the seed is exactly the
partial derivative of $f$ with respect to the node, the root of the expression
tree will receive the value $\frac{\partial f}{\partial f}=1$. Using the seed,
the node computes the correct seed for all of its children and backward-
evaluates from _right-to-left_ , the opposite direction of forward-evaluation.

The correct seed for each child is computed by a simple chain-rule. Assuming
the current node is represented by $w\in\mathbb{R}$ and one of its children is
$v\in\mathbb{R}$, the seed for $v$ is the following:

$\displaystyle\frac{\partial f}{\partial v}$
$\displaystyle=\frac{df}{dw}\frac{\partial w}{\partial v}$ (3)

Each node $w$ has node-specific information to compute $\frac{\partial
w}{\partial v}$. For example, for the log node in Fig. 2, with $w$ as log and
$v$ as $w_{4}$,

$\displaystyle\frac{\partial w}{\partial
v}=\frac{d\log(w_{4})}{dw_{4}}=\frac{1}{w_{4}}$

In general, if $v\in\mathbb{R}^{m\times n}$ and $w\in\mathbb{R}^{p\times q}$,
then

$\displaystyle\frac{\partial f}{\partial v_{ij}}$
$\displaystyle=\sum\limits_{k=1}^{p}\sum\limits_{l=1}^{q}\frac{\partial
f}{\partial w_{kl}}\frac{\partial w_{kl}}{\partial v_{ij}}$ (4)

In particular, the adjoint will always have the same shape and size as the
value.

For nodes that have references back to the containers $x_{i}$, i.e.
$w_{1},\ldots,w_{4}$, they must increment the adjoints in the containers with
their seeds. For example, nodes $w_{3}$ and $w_{4}$ will take their seeds and
increment the adjoint for $x_{3}$. This is easily seen by chain-rule again:
let $w_{1},\ldots,w_{k}$ denote all of the variables that references $x$ and
for simplicity assume they are all scalar, although the result can be easily
generalized for multiple dimensions. Then,

$\displaystyle\frac{\partial f}{\partial x}$
$\displaystyle=\sum\limits_{i=1}^{k}\frac{\partial f}{\partial
w_{i}}\frac{\partial w_{i}}{\partial x}=\sum\limits_{i=1}^{k}\frac{\partial
f}{\partial w_{i}}$

where $\frac{\partial w_{i}}{\partial x}=1$ because $w_{i}$ is simply the
identity function with respect to $x$. The fully accumulated adjoints for
$x_{1},x_{2},x_{3}$ form the gradient of $f$.

In general, an expression node can be quite general and is only required to
define how to compute its value and the adjoints of its children using Eq. 4.

In the general case where $f:\mathbb{R}^{n}\to\mathbb{R}^{m}$, we can apply
this algorithm for the scalar functions $f_{j}$ for $j=1,\ldots,m$ and save
each gradient as the $j$th row of a matrix. The final matrix is then the
Jacobian of $f$.

##  4 FastAD Implementation

In this section, we cover a few key ideas of our implementation111github page:
https://github.com/JamesYang007/FastAD. In Section 4.1, we first discuss the
benefits of vectorization and the difficulties of integrating it into an AD
system. We then explain how FastAD fully utilizes vectorization and
demonstrate that other libraries do not fully take advantage of it. In Section
4.2, we discuss some memory management and performance issues stemming from
the use of the “tape” data structure. We then explain how FastAD overcomes
these challenges using expression templates and a lazy allocation strategy.
Finally, Section 4.3 covers other compile-time optimizations that can further
maximize performance.

###  4.1 Vectorization

Vectorization refers to the parallelization of operations on multiple data at
the hardware level. On a modern Intel 64-bit processor supporting AVX, four
double-precision floating point numbers can be processed simultaneously,
roughly improving performance by a factor of four. While the compiler
optimization is able to vectorize a user’s code sometimes, it is not
guaranteed because vectorization requirements are quite stringent. For
example, vectorization is not guaranteed if memory access is not done in a
contiguous fashion and is impossible if there is any dependency between loop
iterations. This makes it quite challenging to design an AD system that can
always predict compiler optimization to create vectorized code. However,
vectorization can make AD extremely fast, powerful, and practical even in
complex problems. In practice, we come across many examples where operations
can be vectorized during gradient computation. For example, matrix
multiplication, any reduction from a multi-dimensional variable to a scalar
such as summation or product of all elements, and any unary and binary
function that is applied element-wise such as exponential, logarithm, power,
sin, cos, tan, and the usual arithmetic operators.

Since most of the opportunities for vectorization occur in matrix operations,
the goal is to use a well-polished and efficient matrix library such as Eigen,
one of the most popular C++ matrix libraries. However, this is not an easy
task. Adept2.0 notes in their documentation that integrating an expression-
template based matrix library such as Eigen into an AD system can be quite
challenging in design. To circumvent these design issues, Adept2.0 integrates
their own matrix library designed specifically for their AD system. This,
however, only introduces another problem of writing an efficient matrix
library, another daunting task. In fact, the author of Adept notes that matrix
multiplication, one of the most widely used matrix operations, is currently
very slow [8]. Stan mentions that integrating matrix libraries can lead to
unexpected problems that would require extra memory allocations to resolve
[3]. Other libraries such as ADOL-C, CppAD, and Sacado do not integrate any
matrix libraries directly, but they do allow the use of Eigen matrix classes
simply as containers.

The one library among those mentioned previously that ultimately incorporates
Eigen is Stan. Stan provides their own “plug-ins”, which extend Eigen classes
for their AD variables. For example, in Stan, one would use
Eigen::Matrix<stan::math::var, …> as a matrix of AD (univariate) variables,
where the generic matrix object is extended to work differently for
stan::math::var. However, this design still does not take full advantage of
the vectorization that Eigen provides. One reason is that Eigen is only highly
optimized for matrices of primitive types such as double-precision floating
points and integers. For any other class types, Eigen defaults to a less
optimized version with no guarantees of vectorization. Therefore, unless one
extends all Eigen methods with optimizations for their specific class type,
they will not receive much benefit of vectorization. Another reason is that
these matrices now have a heterogeneous structure where each element is an AD
variable which represents a pair of value and adjoint. As a result, it is not
possible to read only the values (and similarly, only the adjoints) of the
matrix in a contiguous fashion. The compiler then cannot guarantee any
automatic vectorization.

Figure 3: The left shows the memory layout for a matrix of univariate AD
variables. This refers to the design pattern “a matrix of pair of pointers to
double” since each element of the matrix contains two pointers pointing to its
value and adjoint. The right shows the memory layout for a single matrix-
shaped AD variable, referring to the reverse pattern “a pair of pointers to
matrix of doubles”. This AD variable has two pointers, each pointing to a
matrix of doubles for the value and adjoint.

FastAD fully utilizes the benefits of Eigen through one simple design
difference, which we refer to as _shape traits_. Stan as well as the other
libraries mentioned previously except Adept2.0 follow the design pattern of “a
matrix of pair of pointers to double” when defining a matrix of AD variables.
Note that a univariate AD variable internally contains two pointers to
doubles, one pointing to the value and one to the adjoint. In contrast, FastAD
follows the reverse pattern: “a pair of pointers to matrix of doubles”. That
is, rather than defining a univariate AD variable, which gets reused as an
element of a matrix, we define a separate matrix AD variable containing a pair
of pointers, each pointing to a matrix of double. Figure 3 shows the memory
layout for each of the design patterns. This subtle difference provides a few
important benefits. Since the value and adjoint are now represented as
matrices of primitive types, matrix operations will be vectorized. The second
benefit is the significant reduction of memory consumption. For other
libraries, a matrix of AD variables contains two pointers _for each element_.
However, with our design, a single matrix AD variable would only store two
pointers regardless of the size of the matrix. Hence, if a matrix is $m\times
n$, the traditional approach has a $O(mn)$ memory consumption for the
pointers, and FastAD approach has a $O(1)$ consumption.

Using templates in C++, it is easy to unify the API for the different AD
variable shapes by providing an additional template parameter:

⬇

ad::Var<double, ad::scl> x; // scalar shape

ad::Var<double, ad::vec> v; // vector shape

ad::Var<double, ad::mat> m; // matrix shape

Shape traits also apply for any arbitrary node because we can deduce the
output shape given the input shapes. The following is a declaration of a
generic node representing an element-wise unary operation, which demonstrates
this idea:

⬇

template <class Unary, class ExprType>

struct UnaryNode:

ValueAdjView<

typename util::expr_traits<ExprType>::value_t,

typename util::shape_traits<ExprType>::shape_t>

{ /*...*/ };

The portion highlighted in red related to shape_t takes the input expression
type ExprType and deduces its shape type. Since an element-wise unary
operation always has the same shape as its input, the unary node takes on the
same shape.

To verify that FastAD vectorizes more than Stan, we performed reverse-AD for a
simple summation function $f(x)=\sum\limits_{i=1}^{n}x_{i}$ and generated the
disassembly for Stan and FastAD 222github page:
https://github.com/JamesYang007/ADBenchmark. We extracted the part that
performs the summation and compared the instructions to see whether
vectorization was taking place. The following is the disassembly for Stan:

⬇

L3178:

movq (%rax), %rdx

addq $8, %rax

vaddsd 8(%rdx), %xmm0, %xmm0

cmpq %rcx, %rax

jne L3178

The instruction used to add is vaddsd, which is an AVX instruction to add
_scalar_ double-precision values. This is not a vectorized instruction, and
hence the addition is not done in parallel. This portion of the disassembly is
related to a specialization of an Eigen class responsible for summation with
_default traversal_ , which is no different from a naive for-loop.

Compare the above disassembly with the one generated for FastAD:

⬇

L3020:

addq $8, %rdx

vaddpd (%rax), %ymm1, %ymm1

vaddpd 32(%rax), %ymm0, %ymm0

addq $64, %rax

cmpq %rdx, %rcx

jg L3020

This portion of the assembly is indeed related to the _linear vectorized
traversal_ specialization of the same Eigen class responsible for the
summation. The instruction used to add is vaddpd, which is an AVX instruction
to add _packed_ double-precision values. This is a vectorized instruction and
the operation is truly done in parallel.

Sometimes, Stan is able to produce vectorized code such as in matrix
multiplication. This is consistent with our benchmark results since Stan came
closest to FastAD for this operation (see Section 5.3). It is also consistent
with how it is implemented, since they allocate extra memory for double values
for each matrix and the multiplication is carried out with these matrices of
primitive types. However, this vectorization does come at a cost of at least 4
times extra memory allocation than what FastAD allocates. Moreover, the
backward-evaluation requires heap-allocating a matrix on-the-fly every time.
FastAD incurs no such cost, only allocates what is needed, and never heap-
allocates during AD evaluation.

###  4.2 Memory Management

Most AD systems manage a data structure in memory often referred to as the
“tape” to store the sequence of operations via function pointers as well as
the node values and adjoints. This tape is modified dynamically and requires
sophisticated memory management to efficiently reuse memory whenever possible.
Stan even writes their own custom memory allocator to alleviate memory
fragmentation, promote data locality, and amortize the cost of memory
allocations [3]. However, the memory is not fully contiguous and may still
over-allocate. For some libraries, on top of memory management of these
operations, a run-time check must be performed at every evaluation to
determine the correct operation [2]. Others like Stan rely on dynamic
polymorphism to look up the vtable to call the correct operation [3].

FastAD is unique in that it uses expression templates to represent the
sequence of operations as a single stack-allocated object. By overloading the
comma operator, we can chain expressions together into a single object. The
following is an example of chaining multiple expressions:

⬇

auto expr = (

x = y * z, // expr 1

w = x * x + 3 * ad::sin(x + y), // expr 2

w + z * x // expr 3

);

Each of the three sub-expressions separated by the commas returns an
expression object containing the necessary information to perform reverse-AD
on their respective structure. Those expression objects are then chained by
the comma operators to build a final expression object that contains the
information to perform reverse-AD on all 3 sub-expressions in the order
presented. This final object is saved into the variable expr.

Expression template makes it possible to build these expression objects
containing the reverse-AD logic. Expression template is a template
metaprogramming technique that builds complex structures representing a
computation at compile-time. For a full treatment of expression templates, we
direct the readers to [12]. As an example, in the following case,

⬇

Var<double, scl> x, y;

auto expr = x + y;

x+y returns a new object of type BinaryNode<Add, Var<double, scl>, Var<double,
scl>>, which represents the addition of two AD variables. In particular, this
object has member functions that define the logic for the forward and backward
evaluations of the reverse-mode AD. This design brings many performance
benefits. Since the compiler now knows the entire sequence of operations for
the expression, it allows for the reverse-AD logic to be inlined with no
virtual function calls, and it removes the need to dynamically manage a
separate vector of function pointers. Additionally, the expression object is
stack-allocated, which is cheaper than being heap-allocated, and its byte size
is proportional to the number of nodes, which is relatively small in practice.

We can optimize the memory management even further with another observation:
an expression can determine the _exact_ number of values and adjoints needed
to represent all of the nodes. If every node can determine this number for
itself, then any expression tree built using these nodes can determine it as
well by induction. It is the case that all nodes defined in FastAD can, in
fact, determine this number. This leads to the idea of _lazy allocation_. Lazy
allocation refers to allocating memory only when the memory is required by the
user. In other words, an expression object does not necessarily need to
allocate memory for the values and adjoints at construction time, since it is
only needed later at differentiation time. Once an expression is fully built
and the user is ready to differentiate it, the user can lazily determine this
number of values and adjoints by interfacing with the expression object,
allocate that amount in a contiguous manner, and “bind” the expression object
to use that region of memory. This solves the issue with the traditional tape
where the values and adjoints are not stored in a fully contiguous manner.
Conveniently, the allocated values and adjoints also do not necessarily need
to be stored and managed by a global tape. Furthermore, the expression objects
can be given additional logic to bind to this region systematically so that
the forward and backward evaluations will read this memory region almost
linearly, which increases cache hits.

While Stan notes that they are more memory-efficient than other popular C++
libraries [3], we noticed a non-negligible difference in memory consumption
between Stan and FastAD. We took the stochastic volatility example in Section
5.6, and compared the memory allocation in number of bytes. For Stan, we took
the member variables from the tape and computed the number of used bytes. We
did not take into account other miscellaneous members for simplicity, and this
estimate serves as a very rough lower bound on the total amount of memory
allocated. For FastAD, we computed the memory allocated for the values,
adjoints, and the expression object. Our result shows that Stan uses 4718696
bytes and FastAD uses 1836216 bytes. This rough estimate shows that Stan uses
at least 2.5 times more memory than FastAD.

With all these optimizations, FastAD removes the need for a global tape by
using expression templates to manage the sequence of operations in an
expression at compile-time, contiguously allocating the exact number of values
and adjoints, and localizing the memory allocations for each expression.

###  4.3 Other Compile-Time Optimizations

Using C++17 features, we can make further compile-time optimizations that
could potentially save tremendous amount of time during run-time.

One example is choosing the correct specialization of an operation depending
on the shapes of the input. As seen in Section 4.1, all nodes are given a
shape trait. Depending on the input shapes, one may need to invoke different
routines for the same node. For example, the normal log-pdf node behaves quite
differently depending on whether the variance parameter is a scalar
$\sigma^{2}$ or a (covariance) matrix $\Sigma$. Namely, if the variance has a
matrix shape, we must perform a matrix inverse to compute the log-pdf, which
requires a different code from the scalar case. Using a C++ design pattern
called Substitution-Failure-Is-Not-An-Error (SFINAE), we can choose the
correct routine at compile-time. The benefit is that there is no time spent
during run-time in choosing the routine anymore, whereas in libraries like
CppAD, they choose the routines at run-time _for every evaluation_ of the node
[2].

Another example is detecting constants in an expression. We can optimize a
node by saving temporary results when certain inputs are constants, which we
can check at compile-time using the C++ language feature if constexpr. These
results can then be reused in subsequent AD evaluations, sometimes changing
orders of magnitude of the performance. As an example, consider an expression
containing a normal log-pdf node for which we differentiate many times.
Suppose also that the input variables to the node are $x$, an $n$-dimensional
vector of constants, and $\mu,\sigma$, which are scalar AD variables. In
general, for every evaluation of the node, the time to forward-evaluate the
log-pdf is $O(n)$, since we must compute

$\displaystyle\log(p(x|\mu,\sigma))$
$\displaystyle=-\frac{\sum\limits_{i=1}^{n}(x_{i}-\mu)^{2}}{2\sigma^{2}}-n\log(\sigma)+C$

However, the normal log-pdf node has an opportunity to make a powerful
optimization in this particular case. Since the normal family defines a two-
dimensional exponential family, it admits a sufficient statistic
$T(x)=\left(\bar{x},S^{2}_{n}(x)\right)$ where

$\displaystyle\bar{x}:=\frac{1}{n}\sum\limits_{i=1}^{n}x_{i},\,S_{n}^{2}(x):=\frac{1}{n}\sum\limits_{i=1}^{n}(x_{i}-\bar{x})^{2}$

Since $x$ is a constant, the sufficient statistic can then be computed _once_
and saved for future use. The normal log-pdf forward-evaluation now only takes
$O(1)$ time given the sufficient statistic, as seen in Eq. 5 below,

$\displaystyle\log(p(x|\mu,\sigma))$
$\displaystyle=-\frac{1}{2\sigma^{2}}\sum\limits_{i=1}^{n}(x_{i}-\mu)^{2}-n\log(\sigma)+C$
$\displaystyle=-\frac{n}{2\sigma^{2}}\left(S_{n}^{2}+(\bar{x}-\mu)^{2}\right)-n\log(\sigma)+C$
(5)

##  5 Benchmarks

In this section, we compare performances of various libraries against FastAD
for a range of examples333github page:
https://github.com/JamesYang007/ADBenchmark. The following is an alphabetical
list of the libraries used for benchmarking:

  * •

Adept 2.0.8 [8]

  * •

ADOL-C 2.7.2 [6]

  * •

CppAD 20200000 [2]

  * •

FastAD 3.1.0

  * •

Sacado 13.0.0 [1]

  * •

Stan Math Library 3.3.0 [3]

All the libraries are at their most recent release at the time of
benchmarking. These libraries have also been used by others [3][9][8].

All benchmarks were run on a Google Cloud Virtual Machine with the following
configuration:

  * •

OS: Ubuntu 18.04.1

  * •

Architecture: x86 64-bit

  * •

Processor: Intel Xeon Processor

  * •

Cores: 8

  * •

Compiler: GCC10

  * •

C++ Standard: 17

  * •

Compiler Optimization Flags: -O3 -march=native

All benchmarks benchmark the case where a user wishes to differentiate a
scalar function $f$ for different values of $x$. This is a very common use-
case. For example, in the Newton-Raphson method, we have to compute
$f^{\prime}(x_{n})$ at every iteration with the updated $x_{n}$ value. In HMC
and NUTS, the leapfrog algorithm frequently updates a quantity called the
“momentum vector” $\rho$ with $\nabla_{\theta}\log(p(\theta,x))$ ($x$ here is
treated as a constant), where $\theta$ is a “position vector” that also gets
frequently updated.

Our benchmark drivers are very similar to the ones used by Stan [3]. As noted
in [9], there is no standard set of benchmark examples for AD, but since Stan
is most similar to FastAD in both design and intended usage, we wanted to keep
the benchmark as similar as possible.

We measure the average time to differentiate a function with an initial input
and save the function evaluation result as well as the gradient. After timing
each library, the gradients are compared with manually-written gradient
computation to check accuracy. All libraries had some numerical issues for
some of the examples, but the maximum proportion of error to the actual
gradient values was on the order of $10^{-15}$, which is negligible. Hence, in
terms of accuracy, all libraries were equally acceptable. Finally, all
benchmarks were written in the most optimal way for every library based on
their documentation.

Every benchmark also times the “baseline”, which is a manually-written forward
evaluation (computing function value). This will serve as a way to measure the
extra overhead of computing the gradient relative to computing the function
value. This approach was also used in [3], however in our benchmark, the
baseline is also optimized to be vectorized. Hence, our results for all
libraries with respect to the baseline may differ from those in the reference.

Sections 5.1-5.4 cover some micro-benchmarks where we benchmark commonly used
functions: summation and product of elements in a vector, log-sum-exponential,
matrix multiplication, and normal log-pdf. Sections 5.5-5.6 cover some macro-
benchmarks where we benchmark practical examples: a regression model and a
stochastic volatility model.

###  5.1 Sum and Product

(a) Sum

(b) Sum (naive, for-loop-based)

(c) Product

(d) Product (naive, for-loop-based)

Figure 4: Sum and product benchmarks of other libraries against FastAD plotted
relative to FastAD average time. Fig. 4(a),4(c) use built-in functions
whenever available. Fig. 4(b),4(d) use the naive iterative-based code for all
libraries.

The summation function is defined as $f_{s}:\mathbb{R}^{n}\to\mathbb{R}$ where
$f_{s}(x)=\sum\limits_{i=1}^{n}x_{i}$, and the product function is defined as
$f_{p}:\mathbb{R}^{n}\to\mathbb{R}$ where
$f_{p}(x)=\prod\limits_{i=1}^{n}x_{i}$. We also tested the case where we
forced all libraries to use a manually-written for-loop-based summation and
product. Fig. 4 shows the benchmark results.

In all four cases, it is clear that FastAD outperforms all libraries for all
values of $N$ by a significant factor. For sum, the next three fastest
libraries, asymptotically, are Adept, Stan, and CppAD, respectively. The trend
stabilizes starting from $N=2^{6}=64$ where Adept is about $10$ times slower
than FastAD and Stan about $15$ times. The naive, for-loop-based benchmark
shown in Fig. 4(b) exhibits a similar pattern, although the performance
difference with FastAD is less extreme: Adept is about $4.5$ times slower and
Stan about $6$ times. Nonetheless, this is still a very notable difference.

For the prod and prod_iter cases, Fig. 4(c),4(d) again show a similar trend as
sum and sum_iter. Overall, the comparison with FastAD is less extreme. For
prod, Adept is about $5$ times slower than FastAD, and Stan about $7$ times.
For prod_iter, Adept is about $3$ times slower and Stan about $4.7$ times.

###  5.2 Log-Sum-Exp

Figure 5: Log-sum-exp benchmark of other libraries against FastAD plotted
relative to FastAD average time.

The log-sum-exp function is defined as $f:\mathbb{R}^{n}\to\mathbb{R}$ where
$f(x)=\log(\sum\limits_{i=1}^{n}e^{x_{i}})$, Fig. 5 shows the benchmark
results.

FastAD outperforms all libraries for all values of $N$. The next three fastest
libraries are Adept, Stan, and CppAD, respectively. The trend stabilizes
starting from $N=2^{6}=64$ where Adept is about $2.5$ times slower than
FastAD, and Stan about $7.7$ times. Note that FastAD is only marginally slower
than the baseline, especially for small to moderate values of $N$, despite
calls to expensive functions like log and exp. In fact, at the largest value
of $N=2^{14}$, FastAD is $1.5$ times slower than the baseline, meaning there
is about $50\%$ overhead from one forward-evaluation to also compute the
gradient. This is because FastAD is optimized such that in this case the
forward-evaluated result for exp expression is reused during the backward-
evaluation.

###  5.3 Matrix Multiplication

Figure 6: Matrix multiplication benchmark of other libraries against FastAD
plotted relative to FastAD average time.

For this benchmark, all matrices are square matrices of the same size. In
order to have a scalar target function, we add another step of adding all of
the entries of the matrix multiplication, i.e.

$f(A,B)=\sum\limits_{i=1}^{K}\sum\limits_{j=1}^{K}{(A\cdot B)}_{ij}$

Fig. 6 shows the benchmark results.

FastAD is still the fastest library for all values of $N$, but Stan performs
much closer to FastAD than in the previous examples. All other libraries
consistently take longer than both FastAD and Stan as $N$ increases. Towards
the end at around $N=2^{14}$, Stan is about $2$ times slower. For moderate
sized $N\in[2^{8},2^{10}]$, Stan is about $3$ times slower.

This example really highlights the benefits of vectorization. As noted in
Section 4.1, this was the one benchmark example where Stan was able to produce
vectorized code, which is consistent with Figure 6 that Stan is the only
library that has the same order of magnitude as FastAD. Other libraries did
not produce vectorized code.

The comparison with the baseline shows that FastAD takes $3.27$ times longer.
Note that forward-evaluation requires one matrix multiplication between two
$K\times K$ matrices, and backward-evaluation additionally requires two matrix
multiplications of the same order, one for each adjoint:

$\displaystyle\frac{\partial f}{\partial A}$ $\displaystyle=\frac{\partial
f}{\partial(A\cdot B)}\cdot B^{T},\,\frac{\partial f}{\partial
B}=A^{T}\cdot\frac{\partial f}{\partial(A\cdot B)}$

Hence, in total, one AD evaluation requires three matrix multiplications
between two $K\times K$ matrices. If we approximate a manually-written
gradient computation to take three times as long as the baseline (one
multiplication), FastAD time relative to this approximated time is
$\frac{3.27}{3}=1.09$. This shows then that FastAD only has about $9\%$
overhead from a manually-written code, which is extremely optimal.

###  5.4 Normal Log-PDF

Figure 7: Normal log-pdf benchmark of other libraries against FastAD plotted
relative to FastAD average time.

The normal log-pdf function is defined up to a constant as:

$f(x)=-\frac{1}{2\sigma^{2}}\sum\limits_{i=1}^{N}\left(x_{i}-\mu\right)^{2}-N\log(\sigma)$

For this benchmark, $\mu=-0.56,\,\sigma=1.37$ and are kept as constants. Fig.
7 shows the benchmark results.

FastAD is the fastest library for all values of $N$. The trend stabilizes at
around $N=2^{7}=128$. Towards the end, Adept is about $9$ times slower, and
Stan about $19$ times slower. Comparing all libraries, the overall difference
we see in this example is the largest we have seen so far, and this is partly
due to how we compute $\log(\sigma)$. Section 4.3 showed that we can check at
compile-time whether a certain variable is a constant, in which case, we can
perform a more optimized routine. In this case, since $\sigma$ is a constant,
it computes the normalizing constant $\log(\sigma)$ once and gets reused over
multiple AD evaluations with no additional cost during runtime, which saves a
lot of time since logarithm is a relatively expensive function.

###  5.5 Bayesian Linear Regression

Figure 8: Bayesian linear regression benchmark of Stan against FastAD plotted
relative to FastAD average time.

This section marks the first macro-benchmark example. We consider the
following Bayesian linear regression model:

$\displaystyle y$ $\displaystyle\sim N\left(X\cdot w+b,\sigma^{2}\right)$
$\displaystyle w$ $\displaystyle\sim N\left(0,1\right)$ $\displaystyle b$
$\displaystyle\sim N\left(0,1\right)$ $\displaystyle\sigma$ $\displaystyle\sim
Unif\left(0.1,10.\right)$

The target function is the log of the joint probability density function (up
to a constant). Fig. 8 shows the benchmark results.

FastAD outperforms Stan by $8.6$ times for the largest $N$ and Adept by $10$
times. The trend stabilizes starting from around $N=70$. It is interesting to
see that around $N=2^{7}$, FastAD is only $2.2$ times slower than the
baseline, despite the model consisting of a large matrix multiplication and
many normal log-pdfs. One of the reasons is that the compiler was able to
optimize-out the backward-evaluation of the matrix constant $X$, since
constants implement a no-op for backward-evaluation.

If we assume that the most expensive operation is the matrix multiplication,
AD evaluation approximately takes two matrix multiplications between a matrix
and a vector. We can then approximate a lower bound for the manually-written
gradient computation time to be two times that of the baseline. The relative
time of FastAD to this approximated time is $1.1$, implying about $10\%$
overhead from a manually-written code.

###  5.6 Stochastic Volatility

Figure 9: Stochastic volatility benchmark of Stan against FastAD plotted
relative to FastAD average time.

This section marks the second and last macro-benchmark example. We consider
the following stochastic volatility model taken from the Stan user guide [11]:

$\displaystyle y$ $\displaystyle\sim N(0,e^{h})$ $\displaystyle h_{std}$
$\displaystyle\sim N(0,1)$ $\displaystyle\sigma$ $\displaystyle\sim
Cauchy(0,5)$ $\displaystyle\mu$ $\displaystyle\sim Cauchy(0,10)$
$\displaystyle\phi$ $\displaystyle\sim Unif(-1,1)$ $\displaystyle h$
$\displaystyle=h_{std}\cdot\sigma$ $\displaystyle h[0]$
$\displaystyle=\frac{h[0]}{\sqrt{1-\phi^{2}}}$ $\displaystyle h$
$\displaystyle=h+\mu$ $\displaystyle h[i]$
$\displaystyle=\phi\cdot(h[i-1]-\mu),\,i>0$

The target function is the log of the joint probability density function (up
to a constant) and we wish to differentiate it with respect to
$h_{std},h,\phi,\sigma,\mu$. Fig. 9 shows the benchmark results.

FastAD outperforms Adept by $2$ times and Stan by $2.6$ times for the largest
$N$. The trend seems to stabilize starting from $N=2^{10}$. It is interesting
to see that FastAD actually performs better than the baseline for moderate to
large $N$ values. For a complicated model as such, there are many
opportunities for FastAD to cache certain evaluations for constants as
mentioned in Section 4.3 and 5.4. In particular, the exponential function
$e^{h}$ reuses its forward-evaluated result, and many log-pdfs cache the log
of its parameters such as $\log(\sigma)$ in the Normal log-pdfs and
$\log(\gamma)$ in the Cauchy log-pdfs ($\sigma,\gamma$ are the second
parameters of their respective distributions, which are constant in this
model). Note that this caching is automatically done in FastAD, which would be
tedious to manually code for the baseline. Hence, this shows that due to
automatic caching, FastAD forward and backward-evaluation combined can be
faster than a manually-written forward evaluation only, which puts FastAD at
an optimal performance.

##  6 Conclusion

In this paper, we first introduced the reverse-mode automatic differentiation
algorithm to give context and background on how FastAD is implemented. We then
discussed how FastAD uses vectorization to boost the overall performance,
expression template and lazy allocation strategies to simplify memory
management, and compile-time checks to further reduce run-time. To see how
FastAD performs in practice, we rigorously benchmarked a set of micro and
macro-benchmarks with other popular C++ AD libraries and showed that FastAD
consistently achieved 2 to 19 times faster speed than the next two fastest
libraries across a wide range of examples.

##  7 Further Studies

FastAD is still quite new and does not have full support for all commonly-used
functions, probability density functions, and matrix decompositions. While
FastAD is currently optimized for CPU performance, its design can also be
extended to support GPU. Having support for both processors will allow FastAD
to be well-suited for extremely large-scale problems as well. Although
computing hessian would be easy to implement in FastAD, generalizing to
higher-order derivatives seems to pose a great challenge, and this feature
could be useful in areas such as physics and optimization problems. Lastly,
one application that could potentially benefit greatly from FastAD is
probabilistic programming language such as Stan, which heavily uses automatic
differentiation for differentiating scalar functions.

##  8 Acknowledgements

We would like to give special thanks to the members of the Laplace Lab at
Columbia University, especially Bob Carpenter, Ben Bales, Charles Margossian,
and Steve Bronder, as well as Andrea Walther, Eric Phipps, and Brad Bell, the
maintainers of ADOL-C, Sacado, and CppAD, respectively, for their helpful
comments on the first draft and for taking the time to optimize the benchmark
code for their respective libraries. We also thank Art Owen and our colleagues
Kevin Guo, Dean Deng, and John Cherian for useful feedback and corrections on
the first draft.

## References

  * [1] Bartlett, R.A., Gay, D.M., Phipps, E.T.: Automatic differentiation of c++ codes for large-scale scientific computing. In: Alexandrov, V.N., van Albada, G.D., Sloot, P.M.A., Dongarra, J. (eds.) Computational Science – ICCS 2006. pp. 525–532. Springer Berlin Heidelberg, Berlin, Heidelberg (2006)
  * [2] Bell, B.: Cppad: a package for c++ algorithmic differentiation http://www.coin-or.org/CppAD
  * [3] Carpenter, B., Hoffman, M.D., Brubaker, M., Lee, D., Li, P., Betancourt, M.: The stan math library: Reverse-mode automatic differentiation in c++ http://arxiv.org/abs/1509.07164v1
  * [4] Goodfellow, I., Bengio, Y., Courville, A.: Deep Learning, chap. 6.5 Back-Propagation and Other Differentiation Algorithms, p. 200–220. MIT Press (2016)
  * [5] Griewank, A., Walther, A.: Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation. Society for Industrial and Applied Mathematics (SIAM), 2 edn. (2008)
  * [6] Griewank, A., Juedes, D., Utke, J.: Algorithm 755: Adol-c: A package for the automatic differentiation of algorithms written in c/c++. ACM Trans. Math. Softw. 22(2), 131–167 (Jun 1996). https://doi.org/10.1145/229473.229474, https://doi.org/10.1145/229473.229474
  * [7] Hoffman, M.D., Gelman, A.: The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo http://arxiv.org/abs/1111.4246v1
  * [8] Hogan, R.J.: Fast reverse-mode automatic differentiation using expression templates in c++ 40 (2014)
  * [9] Margossian, C.C.: A review of automatic differentiation and its efficient implementation http://arxiv.org/abs/1811.05031v2
  * [10] Neal, R.M.: Handbook of Markov Chain Monte Carlo, chap. 5: MCMC Using Hamiltonian Dynamics. CRC Press (2011)
  * [11] Team, S.D.: Stan modeling language users guide and reference manual (2018), http://mc-stan.org
  * [12] Vandevoorde, D., Josuttis, N.M.: C++ Templates: The Complete Guide. Addison-Wesley, 2 edn. (2002)


# Self-stabilizing Algorithm for Maximal Distance-2 Independent Set

Badreddine Benreguia‡ Hamouma Moumen‡111Address for correspondence: H. Moumen,
Department of Informatics, University of Batna 2, Fesdis 05078, Batna,
Algeria. Soheila Bouam‡ Chafik Arar‡

‡ University of Batna 2, 05000 Batna, Algeria
<EMAIL_ADDRESS>
###### Abstract

In graph theory, an independent set is a subset of nodes where there are no
two adjacent nodes. The independent set is maximal if no node outside the
independent set can join it. In network applications, maximal independent sets
can be used as cluster heads in ad hoc and wireless sensor networks. In order
to deal with any failure in networks, self-stabilizing algorithms have been
proposed in the literature to calculate the maximal independent set under
different hypothesis.In this paper, we propose a self-stabilizing algorithm to
compute a maximal independent set where nodes of the independent set are far
from each other at least with distance 3. We prove the correctness and the
convergence of the proposed algorithm. Simulation tests show the ability of
our algorithm to find a reduced number of nodes in large scale networks which
allows strong control of networks.

Keywords: Self-stabilizing algorithm, distributed system, network, independent
set.

##  1 Introduction

###  1.1 Context of the study and motivation

Self-stabilization is a fault tolerance approach that allows distributed
systems to achieve a global correct configuration starting from an unknown
initial configuration. Without external intervention, a self-stabilizing
algorithm is able to correct the global configuration of the distributed
system in a finite time. Various self-stabilizing distributed algorithms have
been proposed in the literature using graph theory such as leader election,
nodes coloring, domination problem, identifying the independent set,
constructing the spanning tree. These algorithms have many benefits in the
real-life applications, for example independent sets have been used as cluster
heads in ad hoc and sensor networks [1, 2, 3, 12].

In graphs, independence is commonly defined as follow: let $G=(V,E)$ be a
graph, where $V$ is the set of nodes and $E$ is the set of edges. An
independent set is a subset of nodes $S\subset V$ such that there is no two
adjacent nodes in $S$. The distance between any two nodes in $S$ is greater
than 1. An independent set $S$ is said $maximal$, if there is no superset of
$S$ that could be an independent set. In other words, there is no node outside
the maximal independent set (MIS) that may join MIS. It is well known in graph
literature that MIS is considered also as dominating set because every node
out of MIS has at least a neighbor in MIS (every node outside MIS is dominated
by a node of MIS).

In this paper, we deal with a particular case of independent set. We call $S$
maximal distance-2 independent set (MD2IS), if of $S$ nodes are independent
and the distance between any two nodes among them is strictly greater than 2.
Figure 1 illustrates difference between MIS and MD2IS where green nodes are
independent. Observe that in MIS (a), distance of 2 could be found between
independent nodes. However, the distance between green nodes in MD2IS (b) is
strictly greater than 2.

Figure 1: (a) Maximal independent set. (b) Maximal distance-2 independent set.

Nodes of MIS are used as servers (cluster heads) in ad hoc and wireless sensor
networks to provide important services for other nodes. Each cluster head has
to guarantee services for nodes connected to it, that are called members of
the cluster. Cluster members represent nodes outside of MIS. A cluster head
could serve its members by routing information, providing keys of encryption,
giving names for members,…

Figure 1(b) shows that elements of MD2SI could be used as cluster heads where
members connected to the head could be within distance of 2. However, using
MIS, members could not be located within distance more than 1. Obviously,
MD2IS provides a more reduced number of clusters than MIS. The choice of the
cluster heads is important in order to contribute in extending lifetime of
wireless sensor and ad hoc networks. Using MD2IS rather than MIS as cluster
heads could provide a good alternative in this sense especially that lifetime
is the major problem of these networks. In addition to that and in order to
deal with any possible failure, we use self-stabilizing algorithm that ensures
reconstructing cluster heads, after the failure occurs, allowing the network
still operational.

Finding the maximal independent set (MIS) in graphs using self-stabilization
paradigm was studied in literature for the first time by Shukla et al. in 1995
[16]. Authors have used a straightforward idea based on two rules: (1) a node
$v$ joins the set $S$ (which is under construction) if $v$ has no neighbor in
$S$, and (2) a node $v$ leaves the set $S$ if at least one of its neighbors is
in $S$. Other variants of self-stabilizing algorithms constructing independent
set have been introduced to deal with particular problems which try to
minimize the algorithm complexity [18] or to be suitable for distributed
daemon 222See section 2.1. [5, 8]. Reader can refer to the survey [6] for more
details on MIS self-stabilizing algorithms. Other self-stabilizing algorithms
have been proposed for independent sets imposing additional constraints
besides to the independence condition. For example, [13] has presented an
algorithm to discover the independent set where each node $u$ out of $S$ has
at least a neighbor $v$ in $S$ such that $deg(v)>deg(u)$. In [2] authors
propose a distributed self-stabilizing algorithm to find MIS using two-hop
(distance-2) information in wireless sensor networks.

###  1.2 Related works

[9, 10] has proposed a self-stabilizing algorithm to find the independent
dominating set imposing a distance greater than $k$ between any two nodes of
the independent set. Work [10] is an improvement of the memory management
regarding the first one [9]. Every node outside the independent set is within
distance $k$. [4] presented a self-stabilizing algorithm to compute a
dominating set $S$ (which is not independent) where every node out of $S$ has
to be distant from $S$ at most by $k$. Although the precedent algorithms have
bounded complexity $O(n+k)$ in rounds, authors indicate that these algorithms
might still never converge under the distributed daemon, since the daemon
could ignore an enabled nodes. It is known in literature that: if the round
complexity of a self-stabilizing algorithm is finite, this does not mean it
converges [4]. Therefore, the computation of the convergence time still an
open question [4, 9, 10] for independent (or dominating) set at distance
$k\geq 2$.

###  1.3 Contribution

In this paper, we propose a self-stabilizing algorithm to find maximal
distance-2 independent set called MD2IS. We prove the correctness and the
convergence of the presented algorithm. Using a serial central daemon, and
starting from an arbitrary configuration, MD2IS reaches the correct
configuration in a limit number of moves. The serial central daemon allows
reaching the correct configuration in the worst case at $2n$ moves. Which
means that MD2IS converges in $O(n)$ moves using a central daemon under
expression distance-2 model. For distance-one model, our algorithm reaches the
correct configuration in $O(nm)$ moves using distributed daemon, where $n$ is
the nodes number and $m$ is the edges number. Proofs and simulation tests
confirm the convergence of the proposed algorithm that provides smaller
independent sets in large scale graphs.

Having independent set more reduced is useful in many applications and allows
more control on large scale networks. For example, the problem of locating an
anonymous source of diffusion [14, 15] needs the placement of few nodes as
observers [17]. A reduced number of nodes that occupy important locations in
graphs, can be used as observers to detect the source of rumors in social
networks.

###  1.4 Organization of the paper

This paper is made up of five sections. Section 2 presents the model and the
terminology used for self-stabilization. In section 3, we introduce the
proposed self-stabilizing algorithm for finding a maximal independent set at
distance-2. Simulation tests are conducted in section 4. Finally, section 5
concludes the paper.

##  2 Model and terminology

Networks and distributed systems are modelled generally as an undirected graph
$G=(V,E)$ where the units of processing represent the set of nodes $V$ and the
links are the set of edges $E$. The neighborhood of a node $v\in V$ is defined
as $N(v)=\\{u\in V:vu\in E\\}$. Usually, we say that two nodes $v,u$ are
$adjacent$ if $u\in N(v)$. We define the neighborhood at distance-2 as
$N(v)_{dist2}=N(v)\cup\\{e\in V:(\exists u\in N(v):e\in N(u))\\}$. It is known
in the graph literature that the distance between two adjacent nodes is 1.
Clearly, the $v$’s neighborhood at distance-2 gathers the set of nodes within
distance of 1 and 2.

A set $S\subset V$ is $independent$ if no two nodes of $S$ are adjacent. In
other word, there is no two nodes in $S$ at distance 1. Generally, a set $S$
of nodes is $distance$-$k$ independent if every node in $S$ is distant at
least $k+1$ to any other node of $S$ [9]. Consequently, a distance-2
independent set is a subset $S$ of $V$ where every two nodes of $S$ are at
distance $>2$. Recall that in the usual case of MIS, each node in the graph is
either independent or dominated $i.e.$ every node of MIS is independent and
every node outside MIS is dominated [7]. In our case, every node $u$ out of
MD2IS is either dominated by a node $v\in$ MD2IS where $dist(u,v)=1$ or is
dominated by $v$ where $dist(u,v)=2$. Note that, a node $u$ out of MD2IS could
be dominated by many nodes of MD2IS, for example it is possible to find a node
$u$ dominated by $v1$ and $v2$ where $dist(u,v1)=1$ and $dist(u,v2)=2$.

Definition : A $distance$-$2$ $independent$ $set$ of a graph $G(V,E)$ is a
subset $S$ of $V$ such that the distance between any two nodes of $S$ is
strictly greater than $2$. $S$ is $maximal$ if no superset of $S$ is also a
distance-2 independent set.

An algorithm is self-stabilizing if it can (1) reach a global correct
configuration called $legitimate$ (2) during a finite time after it starts
from an unknown configuration. When a self-stabilizing algorithm reaches the
correct configuration it must stay inside the correct configuration (known as
$closure$ condition). Hence, to show that an algorithm is self-stabilizing, it
is sufficient to prove its $closure$ for the legitimate configuration and its
$convergence$ to achieve the desired configuration in a finite time. In the
uniform self-stabilizing system, all nodes execute the same code which is a
set of rules having the form: if $guard$ then $statement$ (written as:
$guard\longrightarrow statement$). In this case, nodes use the same local
variables that describe their $state$. The guard is a (or a collection of)
boolean expression. If a guard is evaluated to be true, the corresponding rule
is said $enabled$ (or $priviliged$). We say that a node is $enabled$, if at
least one of its rules is enabled.

Executing the statement of the enabled rule by the node is called a $move$. An
enabled node can make a move only if it is selected by a scheduler called a
$daemon$. A move allows updating the local state (local variables) in order to
make the node more legitimate with its neighborhood.

###  2.1 Daemon notion

The execution of self-stabilizing algorithms is managed by a daemon
(scheduler) that selects the enabled nodes to move from a configuration to
another configuration. Two types of daemons are widely used in literature:
central and distributed daemons. In the central daemons, only one enabled node
is selected to be moved among all the enabled nodes. The central daemon,
called also serial, allows enabled nodes executing a move (one by one) by a
serial order. However, in the distributed daemons, a subset of nodes are
selected among the set of privileged nodes to make a move simultaneously. The
selected subset of nodes to be moved simultaneously forms a $round$. The
distributed daemon is said $synchronous$ when all the enabled nodes are
selected to move simultaneously.

###  2.2 Distance model

Generally, most of the existing self-stabilizing algorithms, use the distance-
one model wherein each node has a partial view on neighbors at distance one
(through an access to the $state$ variable of the neighbors). However, there
is other distance-2 models (like expression model) where every node has to get
information of its neighborhood at distance-2. In the expression model which
is particular case of distance-2 model, access to distance-2 is reached
indirectly through access to $expression$ of neighbors. In the distance-2
model, hypothesis should be constructed carefully. For example, to the best of
our knowledge, there is no algorithm of distance-2 that can operate under
distributed daemon. The existing algorithms of distance-2 have been developed
only under central daemon.

###  2.3 Transformers

A common approach, known in literature [19, 20], allows converting a self-
stabilizing algorithm $A$ operating under a given hypothesis to a new
algorithm $A^{T}$, that operates under other hypothesis different from the
first ones. However, the transformation guarantees that the two algorithms
converge to the same legitimate configuration. Different kinds of transformers
can be found in literature like distance transformers and daemon transformers.
Generally, an overhead of the algorithm complexity is caused by the
transformation which leads generally to a slowdown in the convergence of the
transformed algorithms. In this paper, we use the transformer proposed by [19]
that allows to transform any self-stabilizing algorithm operating under serial
central daemon and expression distance-2 model to a self-stabilizing algorithm
that operates under distributed daemon and distance-one model.

###  2.4 Our execution model

In this paper, we develop a uniform self-stabilizing algorithm. In a first
step, we suppose our algorithm operates under expression distance-2 model
using central daemon. After that, we use transformer proposed by [19] that
allows converting our algorithm to another algorithm that runs under distance-
one model using distributed daemon.

##  3 Self-stabilizing algorithm MD2IS

The proposed self-stabilizing algorithm presented in Algorithm 1, allows
finding the maximal distance-2 independent set. Each node $v$ maintains a
local variable $state$ and an expression $exp$. The $state$ variable could
take one of the values $In$ or $Out$. Once the system reaches the legitimate
configuration, all the nodes are disabled and the set $S=\\{v\in
V:v.state=In\\}$ forms the maximal distance-2 independent set. However, in the
illegitimate configuration, a serial moves of enabled nodes are executed until
the global correct configuration is reached. At every move, an enabled node is
selected randomly by the central daemon.

Every node checks its state regarding its neighborhood using $state$ and
$exp$. Expression $exp$ is used to calculate the number of the neighbors in
$S$. Generally, the expression model allows discovering the neighborhood at
distance-2. Once a node could read expressions of its neighbors, it will have
an information about its neighbors at distance-2.

In our algorithm, R1 shows that if a node $v$ out of $S$ reads the $state$ and
the $exp$ of all its neighbors and finds that all the $state=Out$ and all the
$exp=0$, hence all the neighbors at distance-2 are out, therefore node $v$ has
to join the set $S$. Conversely, R2 illustrates that if a node $v\in S$ knows
that there exists at least one of its neighbors $u$ such that $u.state=In$ or
$u.exp>1$, thus $v$ has to leave $S$. Note that if $v\in S$ and $u.exp>1$
means that $u$ is dominated at least by two nodes: $v$ and another node.
Observe also that if $v.state=IN$ and $u.exp=1$, this means that $u$ is
dominated only by $v$ and it is impossible that $u.exp=0$ because there exists
at least $v\in S$ as a neighbor of $u$.

It is clear that R1 ensures the independence property because node $v$ will be
independent at distance-2 by executing R1 leading to $v.state=In\wedge\forall
u\in N(v):(u.state=Out\wedge u.exp=0)$. However, R2 guarantees that every node
$v$ out of MD2IS is dominated at least by distance-2 because
($v.state=Out\wedge\exists u\in N(v):(u.state=In\vee u.exp>1)$).

Figure 2 shows how our algorithm converges from an initial illegitimate
configuration to a final legitimate configuration using a central serial
daemon. Nodes outlined in red are (privileged) candidates to execute a move.
The central daemon selects in every step only one privileged node to execute a
move. Green nodes represents MD2IS which is reached after 5 moves by the
sequence of rules: R2, R2, R2, R1, R1. Observe that in some cases, a move (for
example transition from configuration c to configuration d) can make other
nodes enabled in the next configuration.

Figure 2: Convergence to the final legitimate MD2IS configuration.

(1) $v.exp::|\\{u\in N(v):u.state=In\\}|$ (2) R1: $v.state=Out\wedge\forall
u\in N(v):(u.state=Out\wedge u.exp=0)\longrightarrow v.state=In$ (3) R2:
$v.state=In\wedge\exists u\in N(v):(u.state=In\vee u.exp>1)\longrightarrow
v.state=Out$

Algorithm 1 Maximal Distance-2 Independent Dominating Set - MD2IS

###  3.1 Closure

###### Lemma 1

When all the nodes are disabled, the set $S=\\{v\in V,v.state=In\\}$ is
maximal distance-2 independent set.

Proof  Suppose that the system is in a legitimate configuration. Since $R2$ is
not enabled for nodes $v$ in $S$, the condition $(\exists u\in
N(v):(u.state=In\vee u.exp>1)$ is false. Thus, $\forall u\in
N(v):(u.state=Out\wedge u.exp\leq 1$). Therefore, at distance 1 from $v$, all
the nodes are out of $S$. And, at distance-2 from $v$ there is no node in $S$
because $u.exp\leq 1$ which is exactly $=1$ ($u$ has at least a neighbor $v$
in $S$, Therefore, all the neighbors of $u$ are out $S$ except $v$. This
implies that all the nodes in the neighborhood of $v$ at distance-2 are out of
$S$).

To show that $S$ is maximal, observe that if we want to add one node $v$ with
$v.state=Out$ then $v$ will have not all the neighbors at distance-2 out of
$S$ ($i.e.$ $v$ has at least a neighbor at distance-2 in $S$). Thus, the
addition of this node will violate distance-2 independence.

$\Box_{Lemma~{}\ref{lem1}}$

###  3.2 Convergence

###### Lemma 2

If a node $v$ executes $R1$ becoming independent, it remains independent, and
every node in its neighborhood at distance-2 still out of $S$ and cannot be
enabled.

Proof  When a node $v$ executes $R1$, that means all the nodes $u$ at
distance-2 from $v$ are out of $S$. It is clear that no one of $u$ could
enabled $R1$ because there exists at least $v$ in $S$ as a node in the
neighborhood of $u$ at distance-2 from $u$. Thus, $v$ will still in $S$ and
all the neighbors at distance-2 still out of $S$. $\Box_{Lemma~{}\ref{lem2}}$

###### Lemma 3

Any node of $V$ will be enabled at most twice by R2 then R1. Thus, MD2IS
terminates in the worst case at $2n$ moves under the expression model using a
central daemon.

Proof  Since Lemma 2 shows that every node executes R1, it cannot move again.
It follows that any node could be enabled (in the worst case) by only R2 and
then R1. Consequently, $2n$ moves is an upper bound to stabilize for $n$
nodes. $\Box_{Lemma~{}\ref{lem3}}$

###### Theorem 1

MD2IS is a self-stabilizing algorithm that constructs Maximal Distance-2
Independent Set in $O(n)$ moves under expression model using a central daemon.

Proof  The proof follows from Lemma 1 and Lemma 3.
$\Box_{Theorem~{}\ref{th1}}$

The last theorem shows that MD2IS stabilizes under the central daemon and
expression model. Now, we use the transformer proposed by [19] that gives a
corresponding self-stabilizing algorithm MD2ISD which operates under
distributed daemon and distance-one model.

###### Theorem 2

MD2ISD converges into legitimate configuration in $O(nm)$ moves in the
distance-one model under a distributed daemon. Where $m$ is the number of
edges.

Proof  Using Theorem 1, the proof follows from Theorem 18 of [19], where $m$
is the number of edges.

$\Box_{Theorem~{}\ref{th2}}$

##  4 Simulation results

In this section, simulation tests are carried out in order to evaluate MD2IS.
We calculate the cardinality of the MD2IS with MIS, and then, we observe how
the cardinality of MD2IS is related when the graph size grows. Although the
comparison with MIS is unfair, there is no other algorithm that we could refer
to evaluate MD2IS cardinality. The algorithms are written in Java using
expression model for MD2IS. For MIS [5], we use the implementation of Lukasz
Kuszner [11] in which we have generated arbitrary graphs with different
density having sizes from 500 nodes to 20000 nodes. For each size of graphs,
we have carried out 5 to 10 executions and then we have taken the average
value.

Figure 3: Cardinality of MD2IS and MIS according density using graphs of 1000
and 3000 nodes.

Figure 3 shows the cardinality of MD2IS and MIS according graphs density. It
is clear that MD2IS gives independent sets smaller than those produced by MIS.
The density of the graphs has a clear impact on the cardinality of the
independent sets. More the density grows, more the cardinality converges to be
smaller. The important observation is that the cardinality of MD2IS will be
close to 1 when the density becomes greater than 0.5. This is a rational
result because for complete graphs (density=1), cardinality of MD2IS is 1.

Figure 4: Cardinality of MD2IS and MIS according graphs size (Density=0.01).

Figure 4 illustrates the cardinality of independent sets according size of
graphs. Using constant density 0.01, curves show that MIS increases
proportionally with the graph sizes. However, cardinality of MD2IS is
inversely proportional to the graph sizes. For graphs of 10000 nodes, the
cardinality of MIS is greater than 100 whereas the cardinality of MD2IS is
less than 10 nodes.

Graph size |  MD2IS Cardinality |  MIS Cardinality |  MD2IS Convergence |  MIS Convergence
---|---|---|---|---
1000 | 601.8 (60.18%) | 687.8 (68.78%) | 400.2 | 415.6
1500 | 703.6 (46.91%) | 926.4 (61.76%) | 652.8 | 636.6
2000 | 758.2 (37.91%) | 1099.4 (54.97%) | 910.8 | 845.4
2500 | 794.4 (31.78%) | 1248.2 (49.93%) | 1229.2 | 1078.4
3000 | 786.6 (26.22%) | 1385.0 (46.17%) | 1570.8 | 1327.4
3500 | 776.4 (22.18%) | 1522.4 (43.50%) | 1877.2 | 1624.6
4000 | 762.2 (19.06%) | 1614.4 (40.36%) | 2202.8 | 1887.4
4500 | 748.8 (16.64%) | 1720.4 (38.23%) | 2483.6 | 2206.0
5000 | 724.2 (14.48%) | 1801.0 (36.02%) | 2785.8 | 2454.8
5500 | 711.0 (12.93%) | 1887.0 (34.31%) | 3048.4 | 2749.2
6000 | 692.2 (11.54%) | 1964.4 (32.74%) | 3335.6 | 3061.6
6500 | 662.8 (10.20%) | 2020.8 (31.09%) | 3551.6 | 3324.6
7000 | 642.6 (9.18%) | 2111.8 (30.17%) | 3825.2 | 3691,4
7500 | 633.8 (8.45%) | 2149.8 (28.66%) | 4098.8 | 3927.2
8000 | 604.6 (7.56%) | 2216.6 (27.71%) | 4333.4 | 4229.2
8500 | 601.8 (7.08%) | 2269.2 (26.70%) | 4586.4 | 4542.8
9000 | 577.6 (6.42%) | 2320.4 (25.78%) | 4836.0 | 4830.8
9500 | 552.8 (5.82%) | 2364.6 (24.89%) | 5097.2 | 5152.4
10000 | 548.4 (5.48%) | 2415.2 (24.15%) | 5350.8 | 5379.2
10500 | 527.4 (5.02%) | 2476.2 (23.58%) | 5645.0 | 5728.4
11000 | 501.6 (4.56%) | 2507.4 (22.79%) | 5817.6 | 6063.2
11500 | 498.2 (4.33%) | 2560.4 (22.26%) | 6091.0 | 6303.6
12000 | 481.2 (4.01%) | 2602.2 (21.69%) | 6300.6 | 6622.0
12500 | 477.8 (3.82%) | 2635.2 (21.08%) | 6555.8 | 6923.8
13000 | 454.6 (3.50%) | 2665.6 (20.50%) | 6790.0 | 7195.0
13500 | 448.0 (3.32%) | 2705.0 (20.04%) | 7057.6 | 7496.6
14000 | 435.2 (3.11%) | 2747.8 (19.63%) | 7251.6 | 7758.6
14500 | 427.4 (2.95%) | 2785.0 (19.21%) | 7512.0 | 8045.0
15000 | 416.2 (2.77%) | 2793.0 (18.62%) | 7798.4 | 8234.0
15500 | 408.2 (2.63%) | 2839.4 (18.32%) | 7988.8 | 8579.8
16000 | 398.8 (2.49%) | 2867.4 (17.92%) | 8293.6 | 8857.8
16500 | 392.0 (2.38%) | 2886.6 (17.49%) | 8496.2 | 9177.6
17000 | 383.4 (2.26%) | 2934.2 (17.26%) | 8798.4 | 9433.6
17500 | 380.2 (2.17%) | 2962.6 (16.93%) | 9043.4 | 9704.0
18000 | 367.6 (2.04%) | 2976.6 (16.54%) | 9249.4 | 9972.0
18500 | 361.0 (1.95%) | 3017.8 (16.31%) | 9481.0 | 10284.0
19000 | 355.2 (1.87%) | 3036.6 (15.98%) | 9748.6 | 10535.8
19500 | 353.0 (1.81%) | 3054.8 (15.67%) | 9998.8 | 10785.0
20000 | 340.4 (1.70%) | 3080.2 (15.40%) | 10207.8 | 11032.4

Table 1: Cardinality and time of convergence of MD2IS (on graphs of
density=0.001).

Further results of simulation are given in table 1 where the time of
convergence is shown to be finite and variates proportionally to the graph
size. Results confirm the lemma 3 that algorithm MD2IS converges at most in
$2n$ moves although the simulations gives smaller number of moves that is
close to $n/2$. In these tests, random graphs have been generated for small
density = 0.001. We use this later value of density in order to be more close
to real networks. For example, for a graph of 10000 nodes, a density of 0.001
will give an average degree = 10 for each node which models a user having 10
friends on social networks. The generated graphs have orders from 1000 nodes
to 20000 nodes. We take the average value after carrying out 5 tests.

##  5 Conclusion

In this paper we proposed a first self-stabilizing MD2IS algorithm that
converges into the correct configuration in $O(n)$ moves using a central
daemon under expression model. In the distance-one model, MD2ISD terminates in
$O(nm)$ moves using a distributed daemon.

The computation of the convergence time for independent sets of $k$ distance
which still an open question is left for the future work. We plan also to
evaluate the use of MD2IS as observers for the problem of locating sources
propagating rumors on real graphs of social networks.

## References

  * [1] O. Arapoglu, V. K. Akram, and O. Dagdeviren. An energy-effcient, self-stabilizing and distributed algorithm for maximal independent set construction in wireless sensor networks. Computer Standards and Interfaces, 62-32; 251 42, 2019.
  * [2] O. Arapoglu and O. Dagdeviren. An Asynchronous Self-Stabilizing Maximal Independent Set Algorithm in Wireless Sensor Networks Using Two-Hop Information, International Symposium on Networks, Computers and Communications (ISNCC), Istanbul, Turkey, 2019, pp. 1-5. doi:10.1109/ISNCC.2019.8909189.
  * [3] D. Bein, A. K. Datta, C. R. Jagganagari, and V. Villain. A self-stabilizing link-cluster algorithm in mobile ad hoc networks. In 8th International Symposium on Parallel Architectures, Algorithms, and Networks, ISPAN, pages 436-441. IEEE Computer Society, 2005.
  * [4] A. K. Datta, S. Devismes, and L. L. Larmore. A silent self-stabilizing algorithm for the generalized minimal k-dominating set problem. Theoretical 259 Computer Science, 753:35-63, 2019.
  * [5] W. Goddard, S. T. Hedetniemi, D. P. Jacobs, and P. K. Srimani. Self-stabilizing protocols for maximal matching and maximal independent sets for ad hoc networks. In 5th IPDPS Workshop on Advances in Parallel and Distributed Computational Models, page 162. IEEE Computer Society, 2003. doi:10.1109/IPDPS.2003.1213302.
  * [6] N. Guellati and H. Kheddouci. A survey on self-stabilizing algorithms for independence, domination, coloring, and matching in graphs. Journal of Parallel and Distributed Computing, 70(4):406-415, 2010.
  * [7] S. Hedetniemi, S. Hedetniemi, D. Jacobs, and P. Srimani. Self-stabilizing algorithms for minimal dominating sets and maximal independent sets. Computer and Mathematics with Applications, 46(56):805-811, 2003.
  * [8] M. Ikeda, S. Kamei, and H. Kakugawa. A space-optimal self-stabilizing algorithm for the maximal independent set problem. In 3rd International Conference on Parallel and Distributed Computing, Applications and Technologies, pages 70-74, 2002.
  * [9] C. Johnen. Fast, silent self-stabilizing distance-k independent dominating set construction. Information Processing Letters, 114:551-555, 2014.
  * [10] C. Johnen. Memory efficient self-stabilizing distance-k independent dominating set construction. In Networked Systems - Third International Conference, pages 354-366. Springer International Publishing, 2015.
  * [11] L. Kuszner. Tools to develop and test self-stabilizing algorithms. http://kaims.eti.pg.gda.pl/ kuszner/self-stab/main.html, 2005.
  * [12] A. Larsson and P. Tsigas. A self-stabilizing (k,r)-clustering algorithm with multiple paths for wireless ad-hoc networks. In 31 International Conference on Distributed Computing Systems, pages 353-362. IEEE, 2011.
  * [13] B. Neggazi, N. Guellati, M. Haddad, and H. Kheddouci. Efficient self-stabilizing algorithm for independent strong dominating sets in arbitrary graphs. International Journal of Foundations of Computer Science, 26(6):751-768, 2015.
  * [14] P. C. Pinto, P. Thiran, and M. Vetterli. Locating the source of diffusion in large-scale networks. Physical Review Letters, 109(6):068702, 2012.
  * [15] E. Seo, P. Mohapatra, and T. Abdelzaher. Identifying rumors and their sources in social networks. In Proceedings of SPIE - The International Society for Optical Engineering, volume 8389, 2012.
  * [16] S. K. Shukla, D. J. Rosenkrantz, and S. S. Ravi. Observations on self-stabilizing graph algorithms for anonymous networks. In 2nd Workshop on Self-Stabilizing Systems, 1995.
  * [17] B. Spinelli, L. E. Celis, and P. Thiran. Observer placement for source localization: The effect of budgets and transmission variance. In 54th Annual Allerton Conference on Communication, Control, and Computing , pages 743-751, 2016.
  * [18] V. Turau. Linear self-stabilizing algorithms for the independent and dominating set problems using an unfair distributed scheduler. Information Processing Letters, 103(3):88-93, 2007.
  * [19] V. Turau. Efficient transformation of distance-2 self-stabilizing algorithms. Journal of Parallel and Distributed Computing, 72(4):603-612, 2012.
  * [20] W. Goddard, P. K. Srimani, Daemon conversions in distributed selfstabilizing algorithms, in: Ghosh S.K., Tokuyama T. (eds) WALCOM: Algorithms and Computation. WALCOM 2013. Lecture Notes in Computer Science, Springer, Berlin, Heidelberg, 2013, pp. 146–157.