# Network Traffic Analysis of Medical Devices

Nowfel Mashnoor  Computer Science and Engineering
University of Nevada, Reno
<EMAIL_ADDRESS>Batyr Charyyev  Computer Science and Engineering
University of Nevada Reno
<EMAIL_ADDRESS>
###### Abstract

The availability of medical devices such as glucose level and blood pressure
measuring devices is continuously increasing. These devices have gained user
interest as they are easy to use. However, medical devices introduce extra
complexity to the network by being numerous, heterogeneous, and more
vulnerable to cyber-attacks. For better network management and overall network
security, it is important to understand the network traffic characteristics of
the devices. Thus, in this paper, we analyze in detail the traffic
characteristics of 8 medical devices both at the device level and at the level
of individual functionality of each device. We collect and share both network
and bluetooth traffic from a total of 51 functionalities of the devices. Our
analysis includes different metrics such as protocols, amount of
incoming/outgoing traffic, DNS queries, and analysis of traffic destinations.
We observed that devices have unique network and bluetooth traffic
characteristics, that might be useful in developing networking tools for
medical devices.

###### Index Terms:

Internet of Things, Networking, Medical Devices, Network Traffic Analysis.

Table I: Medical devices used in the study.

Device name | Connectivity | # of functions | Description
---|---|---|---
Wifi | Bluetooth
BabyMonitor | y | y | 7 | Real-time baby monitoring device through a video stream, manufactured by Sense-U.
EarWax | y | n | 5 | Device for ear wax removal with high-definition camera and LED lighting.
FitBitFitness | y | y | 5 | Fitness tracker that tracks daily activity, heart rate, and sleep patterns.
|  |  |  | Provides 24-hour emergency support for seniors, including a wearable alert
GuardianAlert | y | n | 7 | button and two-way communication for immediate assistance in a crisis.
KardiaMobile | y | n | 6 | Device that can record a medical electrocardiogram (EKG) on the smartphone.
KetoScanMeter | y | y | 5 | Breath-acetone analyzer that measures the amount of ketones in a person’s breath.
WithingsBPM | y | y | 8 | Blood pressure monitor that can be synced with a smartphone app for tracking.
WyzeScale | y | y | 8 | Scale that measures weight and body mass indexes.

##  I Introduction

The Internet of Things (IoT) devices are continuously being integrated into
our daily lives. These devices have various utilities that can improve quality
of life through automation and monitoring. For instance, an IoT device ring
doorbell enables you to monitor the front door of your house or any facility,
notifies you when there is a visitor through motion/sound detection, and
enables you to watch the environment remotely by providing live stream from
its camera. These devices are also different from regular devices such as
laptops and mobile phones as IoT devices have more simpler nature with a
limited/dedicated set of functionalities [1]. For instance, smart blood
pressure measuring device should only have functionalities related to blood
pressure measurement and provide some statistical reports about the results.
On the other hand, IoT devices have a heterogeneous nature as there are
different types of devices and numerous manufacturers following different
design principles [2].

In general, IoT devices come in the form of surveillance devices (ring
doorbells, security cameras), home appliances (smart plugs, lights), and
medical devices such as blood pressure measuring devices, and fitness
monitoring devices. In particular medical devices are gaining popularity among
all the people as they enable the user to monitor basic health status such as
blood pressure, and body temperature without visiting a physician or doctor,
and in this paper, we will mainly focus on medical devices. While medical
devices have many use cases, cybersecurity, and privacy are the main factors
that hinder its wide adoption [2]. Poor design principles and errors by
manufacturers may lead to the production of malfunctioning devices with
security vulnerabilities. This might be due to global competition in the IoT
industry and manufacturers trying to push their devices to market as soon as
possible [2].

Network management becomes important with the integration of medical devices
into the network as these devices are more vulnerable to cyber-attacks.
Moreover, due to the limited computational resources of IoT devices, it is
challenging to install protective software on a device [3]. Thus developing
networking tools tailored for medical devices enables better network
management and security. Since network traffic-based monitoring tools such as
device identification [4], and intrusion detection [5] mainly use machine
learning with features extracted from network data, it becomes crucial to
identify protocols and characteristics of network traffic features from the
devices. Towards this direction, it is necessary to understand the network
characteristics of these devices prior to designing the dedicated networking
tools and software. Thus in this paper, we first collect network traffic data
and bluetooth data from 8 medical devices. Then we analyze these data on both
the device level and functionality level of each device. Specifically, we
explore how traffic data changes from device to device, and additionally how
data changes for each functionality of the device. Our analysis includes
statistical attributes such as activity cycles, port numbers, and amount of
incoming/outgoing traffic data. We also explore application/transport layer
protocols used by the devices, DNS characteristics, and contacted IP
destinations. We believe that findings from this study will provide a better
understanding of the traffic characteristics of medical devices. The detailed
results of our analysis and data collected in our laboratory will be publicly
available in our GitHub repository (github.com/iotsec-lab/medical).

##  II Related Works

In this section, we provide related works on medical devices, network traffic,
and bluetooth traffic analysis.

Previous studies on medical devices focus on the adoption of medical devices
in healthcare [2], healthcare systems that utilize wearable electrocardiogram
(ECG) monitoring devices based on bluetooth technology [6], and the potential
of integrating bluetooth devices and big data analysis in healthcare networks
[7]. There exist studies that focus on the ZigBee mesh protocol and enhancing
the efficacy of routine monitoring procedures with medical devices [8].

There exist studies on network traffic analysis of traditional devices such as
personal computers [9] and mobile applications [10]. There also exist studies
on traffic analysis of IoT devices [11, 3], and drones [12]. Additional to
generic traffic characterization [13, 11, 3], most of these studies focus on
fingerprinting the network traffic for OS identification [9], device
identification [4], user interaction classification [1], and intrusion
detection [5]. Studies mostly use machine learning [11, 9], natural language
processing [14, 15], and hashing-based [4] methods to make traffic analysis.

Existing studies on bluetooth traffic analysis focus on analyzing
electrocardiogram data from wearable devices [16], identification of wearable
devices through bluetooth data [17], and estimating travel time and speed on a
freeway with bluetooth traffic [18]. The methodologies to intercept bluetooth
traffic are explored in [19]. There also exist studies focusing on privacy
implications of bluetooth traffic sniffing and designing the countermeasures
to sniffing [20]. Researchers also explored bluetooth based contact tracing
[21].

Compared to existing studies on medical devices [2, 6] our work focuses on
analyzing both network traffic and bluetooth traffic of medical devices.
Compared to studies in network traffic analysis [9, 12, 13] and bluetooth
traffic analysis [16, 17, 22], our work explores how traffic characteristics
vary for different functionalities of the device, and provide in-depth
analysis on features and traffic destinations.

##  III Methodology

In this section, we will discuss the data collection and cleaning process. We
collected data from 8 medical devices, and the summary of devices used in our
study is provided in Table I. For each device, we first identified a set of
functionalities that the device has. For instance, the set of functionalities
for WithingsBPM device includes turning on/off, taking blood pressure,
displaying results, syncing with the server, etc. Then we identified if the
device uses WiFi or Bluetooth connectivity. After that, we identified which
functionalities of the device generate traffic data. The number of
functionalities for each device that produces traffic data is provided in
Table I. Then for each functionality of the device and for each possible type
of connectivity (WiFi or Bluetooth), we collected data for 10 measurements. It
is important to mention that not all device functionalities generate traffic,
thus we provided only a number of functionalities that generate traffic data.
Also, note that not all devices use Bluetooth connectivity.

###  III-A WiFi Data

To collect the WiFi data we connected the devices to a Linksys router flushed
with OpenWrt firmware. Then we used PCAPDroid with automated Python scripts to
collect the data. We repeated each functionality of the device 10 times and
collected the data during the process. Some of the IoT devices generate the
traffic data while they are in idle state [13]. However, we did not observe
any traffic being generated when devices are not being interacted. Thus our
network traffic is limited to the active state of the devices. Overall we
collected data for 510 measurements that include data from 8 devices and 51
different functionalities in total.

###  III-B Bluetooth Data

We collected Bluetooth data from two phases of the device namely advertisement
and interaction. In the advertisement phase devices send advertisement packets
to connect with the base station. Once connection is established devices start
to transmit data packets while the user is interacting with the medical device
through a controller (mobile phone). In order to collect the data from medical
devices, we used a nRF Sniffer [23], a tool designed for intercepting
Bluetooth Low Energy (BLE) packets. When the nRF Sniffer is used with
Bluetooth Dongle connected to the data collection server, Wireshark can
identify the sniffer as a separate network interface enabling us to collect
the bluetooth traffic. First, for each device, we collected 10 advertisement
data, each measurement lasting 20 seconds. Then, for each device, we collected
10 interaction data where we performed a set of interactions for each device
in a particular order. Overall we collected bluetooth traffic from 100
measurements that include advertisement and interaction data from 5 medical
devices.

##  IV Evaluation Results

In this section, we will first present an analysis of WiFi data and then an
analysis of Bluetooth data.

###  IV-A Flow duration and volume

We first evaluated the duration, volume, and packet count of traffic flows
generated by each device using the methodology from [11, 3] for each metric.
Here traffic flow is network traffic data generated by a single functionality
of the device. The flow duration is basically the time between the first and
the last packet in a flow, the volume is the total download and upload bytes,
and packet count is the number of sent and received packets by the device when
it is performing the functionality [11, 3]. Fig. LABEL:average_dur_vol_pc
presents the average duration, volume, and packet count for all devices. We
can see that the average duration of flows is mostly less than 10 seconds, but
for some devices such as BabyMonitor and EarWax, the average duration can
reach up to 59 and 21 seconds respectively. When we look at the values for
individual functionalities of the devices Fig. LABEL:individual_dur_vol_pc_1,
we see that BabyMonitor functionalities 2-5 have an average duration of less
than 10s and functionality 1 (turn on) and 7 (turn off) can reach up to 60s.
While for all devices most of the functionalities have a flow duration of less
than 10s, for each device there exist certain functionalities that have a flow
duration greater than 10s (Fig. LABEL:individual_dur_vol_pc_1). In terms of
traffic volume and packet count, we see that they follow similar patterns as
expected. Overall devices exchange small amounts of data, we observed that 95%
of flows have a volume of less than 300 KB and 61% of flows have a volume of
less than 100 KB. Interestingly, the duration and volume of the data don’t
show any correlation. For instance, in terms of duration, BabyMonitor has the
highest average duration of the flows whereas its volume and packet count are
less than WyzeScale which suggests that WyzeScale has a higher traffic rate
compared to BabyMonitor. In terms of individual functionalities, BabyMonitor
function 7 has the highest duration of the traffic flow and lowest packet
count (Fig. LABEL:individual_dur_vol_pc). Overall it suggests that devices do
not generate a large amount of data but the traffic rate can change for each
device and its functionality, indicating that we can observe a high burst of
network traffic for the devices.

###  IV-B Protocols and port numbers

We analyzed the transport layer and application layer protocols of the devices
and port numbers that they use. Fig. 9 shows for each device the percentage of
transport and application layer protocols. In the figure, each device has two
columns associated with it, where the first column shows the percentage of
transport layer protocols and the second column shows the percentage for the
application layer. In the transport layer, most of the devices use TCP except
for EarWax and WyzeScale. The percentage of UDP packets in WyzeScale is around
11% and this may correspond to DNS queries or application layer QUIC protocol
that uses UDP as underlying transport layer protocol. For Earwax the
percentage of UDP packets is higher. However, the operation of EarWax is
slightly different compared to other devices. It opens its local hotspot to
connect the device and mobile app then it streams the data through that
channel. Thus, EarWax doesn’t have much communication with the rest of the
world, as we will see in the next subsection on incoming/outgoing traffic
analysis. As a result, data generated by the EarWax mostly doesn’t leave the
local network. However, it generates some HTTP and DNS packets as an
application layer protocols.

Figure 9: Protocol percentages for application and transport layer.

In terms of application layer protocol, we observed that HTTPS is the dominant
protocol in most of the devices. Every device has a certain amount of DNS
communication and for WyzeScale the QUIC protocol makes around 14% of its
application layer protocol. Only in BabyMonitor, we observe that application
layer protocol is split between HTTPS and HTTP being 51.9% and 42.0% of
traffic flows respectively. In terms of destination port number, we observed
that devices reach traditional port numbers such as port number 443 (used for
HTTPS services over TCP) and port number 80 (used for HTTP over TCP). Other
than traditional port we observed that BabyMonitor reaches port number 8883 in
10% of its traffic flow and WyzeScale has a destination port number 9870 in
0.4% of its traffic flow. While these port numbers might be used for updates
and patches they also might be used for other purposes which are not critical
for device functionality.

###  IV-C Traffic direction

We also analyzed the direction of the traffic generated by the devices. The
outbound traffic is the traffic that leaves the local network and is sent to
the internet. The inbound traffic is the traffic that is sent from the
internet to the device. We observed that in some devices such as EarWax most
of the communications are between the device and its mobile app thus data does
not leave the local network. Thus for consistency, we considered those traffic
as part of the inbound traffic. Fig. LABEL:user_detection_for_each_method
shows average inbound/outbound traffic for all devices and traffic for each
functionality for the sample of devices. Overall we observed that for most of
the devices, the traffic is inbound, meaning that communication is happening
in the local network (Fig. LABEL:alldevice_direction). This is expected as
medical devices are mostly used with their mobile app and most of the time
data doesn’t leave the local network except for the cases with DNS queries,
firmware update checks, or fetching the data from remote servers. For each
device functionality, overall we observed that GuardianAlert, KardiaMobile,
and KetoScanMeter have a similar pattern as BabyMonitor (shown in Fig.
LABEL:babymonitor_direction) where inbound and outbound traffic are split
equally, and FitbitFitness shows a similar pattern to WithingsBPM (shown in
Fig. LABEL:WithingsBPM_direction) where inbound traffic is slightly higher for
each device functionality.

###  IV-D DNS queries and traffic destinations

Then we focused on the destination of the traffic sent to the internet and the
source of the traffic coming from the internet. We started our analysis with
DNS queries. Fig. 10 shows the average DNS query count observed across all
device functionalities and for each functionality of a device. Figure shows
EarWax is a device that has the lowest average DNS query count (0.46) and
WyzeScale has the highest average DNS query count (9.45). We also observed
that not all device functionalities make DNS queries. For instance, while
BabyMonitor and WithingsBPM make DNS queries for each of their
functionalities, FitBitFitness, GuardianAlert, KardiaMobile, KetoScanMeter,
and WyzeScale do not make any DNS queries in their functionality 1, and only
functionalities 1 and 2 of EarWax make DNS queries. Then we analyzed the
unique DNS queries for each device, results are shown as a secondary axis in
Fig. 10. We observed that most of the devices have 4-6 unique DNS queries,
EarWax has only 1 unique DNS query which is android.bugly.qq.com and WyzeScale
has the highest number of unique DNS queries with the count of 26.

Table II: Destinations contacted by the devices, top domain name based on the
frequency of the contact.

Country | Organization | Contacting devices | Top Domain
---|---|---|---
United States | Amazon | BabyMonitor | a1.tuyaus.com
| GuardianAlert | ws-us2.pusher.com
| KardiaMobile | us-telekardia-production.alivecor.com
| KetoScanMeter | api.channel.io
| WyzeScale | wyze-platform-service.wyzecam.com
| Google | BabyMonitor | firebaselogging-pa.googleapis.com
|  | FitbitFitness | scone-pa.googleapis.com
|  | GuardianAlert | clients4.google.com
|  | KardiaMobile | api.mixpanel.com
|  | WithingsBPM | prod.rudderstack.withings.net
|  | WyzeScale | googletagmanager.com
| Alibaba | BabyMonitor | v2.sense-u.com
| Cloudflare | WyzeScale | use.fontawesome.com
| Fastly | WithingsBPM | static.withings.com
|  | WyzeScale | sdk.iad-03.braze.com
| Meta Platforms IRL | KetoScanMeter | graph.facebook.com
| Microsoft Azure | GuardianAlert | mymedicalguardian.com
| NBP NCP | KetoScanMeter | kminiapp.sentechkorea.com
Canada | Shopify | WyzeScale | wyzecom.myshopify.com
China | CNNIC | BabyMonitor | amdcopen.m.taobao.com
China Mobile | BabyMonitor | utoken.umeng.com
Hangzhou Alibaba Adv | BabyMonitor | offmsg.umeng.com
France | BSO Network Solutions | WithingsBPM | scalews.withings.net
Hong Kong | Tencent Cloud Computing | EarWax | android.bugly.qq.com
BabyMonitor
Singapore | Alibaba | BabyMonitor | ulogs.umeng.com
South Korea | Amazon | KetoScanMeter | ws.channel.io

Then to find the country and organization of servers that the device is
contacting, we analyzed the IP addresses returned by the DNS servers using the
WHOIS tool.

Figure 10: Average DNS query count and unique DNS count for each device and
functionality.

Table-II shows the country and names of the organizations that devices are
contacting and the top domain name for the organization. For the organization
name, we provided parent organizations such as Google for Google Cloud and
Amazon for AWS. Overall we observe that devices reach cloud providers most of
the time. This is expected as medical devices have limited computational power
thus, they need to offload computations to remote servers. Cloud providers
include Google, Amazon, Microsoft, and Alibaba, and 7 devices out of 8
contacts at least Amazon or Google. The only device that doesn’t contact these
providers is EarWax which uses Tencent Cloud a cloud provider based in China.
Even though devices are being used in the US we can see that 5 devices out of
8 reach 6 overseas countries. Some devices send traffic to countries where
they are manufactured such as WithingsBPM to France and KetoScanMeter to South
Korea. However, BabyMonitor is made by Sense-U, a company based in California
US but we can see that it sends traffic overseas to countries such as China
and Hong Kong. Interestingly all of the organizations in China are third
parties such as CNNIC (China Internet Network Information Center) a government
agency in China, a telecommunication company China Mobile, and Hangzhou
Alibaba Advertising. From a privacy perspective, this is concerning as the
device being used in the US and made by a company based in the US, sends data
to a government agency and telecommunication company of another country. Even
sending data to cloud (in particular Amazon, Google, and Microsoft) is privacy
concerning as it allows them to potentially profile consumers. For example,
these companies not only can learn the types of devices in a household, but
also how/when they are used, simply by analyzing the traffic from devices to
their cloud services [13].

###  IV-E Bluetooth communication

In addition to the WiFi connection, some medical devices (5 out of 8) also
have Bluetooth connections. In general bluetooth devices first initiate
connection through advertisement then after establishing the connection they
send data packets to enable interaction. We collected data for both the
advertisement and interaction phases and the packets are in the form of
Bluetooth Low Energy Link Layer protocol [24].

Table III: Average packet size (bytes), traffic rate (bps), and inter-arrival
timing (ms) for the Bluetooth traffic of medical devices.

Metric (Avg) | BabyMonitor | FitBitFitness | KetoScanMeter | WithingsBPM | WyzeScale
---|---|---|---|---|---
Advertisement
packet size | 51.70 | 46.17 | 48.16 | 53.95 | 48.66
inter-arrival time | 0.04 | 0.21 | 0.04 | 0.02 | 0.11
traffic rate | 2345.79 | 718.90 | 1379.10 | 2358.56 | 547.21
Interaction
packet size | 30.60 | 42.08 | 30.72 | 38.56 | 35.16
inter-arrival time | 0.06 | 0.08 | 0.04 | 0.03 | 0.04
traffic rate | 1177.66 | 680.53 | 816.67 | 1575.61 | 908.67

Figure 11: Percentage of observed PDU values for each device in advertisement
flows.

Our goal in this subsection is to provide general characteristics of bluetooth
connectivity of medical devices and explore traffic features that might be
useful to better fingerprint the devices. Table III presents the
characteristics of the bluetooth traffic in terms of the average packet size,
inter-arrival timing, and traffic rate, as these features are widely used by
previous studies to fingerprint the bluetooth traffic [17]. Overall we
observed that bluetooth traffic is different from WiFi traffic having a lower
traffic volume. We observed that both in the advertisement and interaction
phase, BabyMonitor and WithingsBPM devices have the highest traffic rate,
whereas FitBitFitness and WyzeScale have the lowest. Devices have unique
characteristics in terms of basic traffic features such as packet size and
traffic rate. Also, advertisement and interaction traffic flows also differ in
terms of basic characteristics even if they belong to the same device. For
instance, we can see that the traffic rate of BabyMonitor in the advertisement
phase (2345 bps) is almost double the rate in the interaction phase (1177
bps). Such differences in basic traffic characteristics can enable
fingerprinting the device for device identification and device state
identification and thus used by previous studies [17, 16] for fingerprinting
purposes.

Next, we explored if there exist additional features that can be used to
fingerprint the medical devices through bluetooth connection. To that we
extracted all 262 features listed in [24] and explored the ones that have high
variance and can be used to fingerprint the devices. Overall we observed that
most of the features have low variance or attributes corresponding to those
features are left empty in the packet. However, there exist certain features
that can be extracted from the header field of the bluetooth packet and can be
used for traffic analysis. One of those features is PDU type which can be
extracted from the header field of the packet. PDU type is a Protocol Data
Unit that can be found in both advertisement and interaction packets. In
general, PDU type is a 2-257 octet value and defines the type of the packet.
For instance, in advertisement packets, PDU value 0000b corresponds to the
Advertising Indications (ADV_IND) packet and 0001b corresponds to the Directed
Advertising (ADV_DIRECT_IND) packet. Fig. 11 shows for each device the
percentage of packets with specific PDU values. Overall we can see that PDU
type 0x00, 0x03, and 0x04 are the most commonly observed values in all
devices. We observed the same behavior in interaction flows as well. However,
when we consider other PDU values we can see that there is a unique
characteristic for each device in the amount of different PDU values. For
instance, in advertisement flows of BabyMonitor, we can observe all PDU types
except 0x0d whereas for WyzeScale we can observe only 7 out of 14 PDU types.
Also, we observed that PDU values for BabyMonitor are different for
advertisement and interaction flows (not shown in the paper), as in
advertisement flows we observed 13 values out of 14, and in interaction flows
we observed only 7 values out of 14. This suggests that in addition to classic
packet size and inter-arrival timing features, the PDU types extracted from
bluetooth packet header can be used to identify a device and its state. In
addition to PDU values, we also observed certain header features such as
header.flags, header.tx_power, and header.llid have high discriminatory power
and can be used for various traffic analysis such as fingerprinting for
network management and intrusion detection to prevent malicious activities.

##  V Conclusion

In this paper, we analyzed the network traffic characteristics of medical
devices. We explored both wifi and bluetooth data of these devices at both
device and functionality levels. Overall we observed that devices have unique
behavior in terms of basic traffic characteristics such as volume of data, and
signaling patterns. In terms of traffic destinations, we observed that some
devices pose privacy concerns, as there are some devices contacting third
parties such as advertising and overseas telecommunication companies. In
bluetooth data, we observed that there are certain traffic features such as
PDU type, header tx power, and header flags that have high discriminatory
power. Since these features can be easily extracted from packet headers they
can be good alternatives and supplements to basic traffic features such as
packet size and inter-arrival timing used in traffic fingerprinting and
analysis.

## References

  * [1] B. Charyyev and M. H. Gunes, “Iot event classification based on network traffic,” in _2020 IEEE INFOCOM WKSHPS_. IEEE, 2020, pp. 854–859.
  * [2] B. Charyyev, M. Mansouri, and M. H. Gunes, “Modeling the adoption of internet of things in healthcare: A systems approach,” in _2021 IEEE ISSE_. IEEE, 2021, pp. 1–8.
  * [3] A. Sivanathan, D. Sherratt, H. H. Gharakheili, A. Radford, C. Wijenayake, A. Vishwanath, and V. Sivaraman, “Characterizing and classifying iot traffic in smart cities and campuses,” in _2017 IEEE INFOCOM WKSHPS_. IEEE, 2017, pp. 559–564.
  * [4] B. Charyyev and M. H. Gunes, “Locality-sensitive iot network traffic fingerprinting for device identification,” _IEEE Internet Things J._ , vol. 8, no. 3, pp. 1272–1281, 2020.
  * [5] A. Hamza, H. H. Gharakheili, T. A. Benson, and V. Sivaraman, “Detecting volumetric attacks on lot devices via sdn-based monitoring of mud activity,” in _Proceedings of the 2019 ACM SOSR_ , 2019, pp. 36–48.
  * [6] J.-H. Kim, S.-C. Lee, B.-G. Lee, and W.-Y. Chung, “Mobile healthcare system based on bluetooth medical device,” _J. Sens. Sci. Technol_ , vol. 21, no. 4, pp. 241–248, 2012.
  * [7] D. V. Dimitrov, “Medical internet of things and big data in healthcare,” _Healthc. Inform. Res._ , vol. 22, no. 3, pp. 156–163, 2016.
  * [8] R. K. Kodali, G. Swamy, and B. Lakshmi, “An implementation of iot for healthcare,” in _2015 IEEE RAICS_. IEEE, 2015, pp. 411–416.
  * [9] A. Aksoy, S. Louis, and M. H. Gunes, “Operating system fingerprinting via automated network traffic analysis,” in _2017 IEEE CEC_.
  * [10] J. Ren, A. Rao, M. Lindorfer, A. Legout, and D. Choffnes, “Recon: Revealing and controlling pii leaks in mobile network traffic,” in _Proceedings of the 14th MobiSys_ , 2016, pp. 361–374.
  * [11] A. Sivanathan, H. H. Gharakheili, F. Loi, A. Radford, C. Wijenayake, A. Vishwanath, and V. Sivaraman, “Classifying iot devices in smart environments using network traffic characteristics,” _IEEE TMC_ , vol. 18, no. 8, pp. 1745–1759, 2018.
  * [12] S. Sciancalepore, O. A. Ibrahim, G. Oligeri, and R. Di Pietro, “Pinch: An effective, efficient, and robust solution to drone detection via network traffic analysis,” _Computer Networks_ , vol. 168, p. 107044, 2020.
  * [13] J. Ren, D. J. Dubois, D. Choffnes, A. M. Mandalari, R. Kolcun, and H. Haddadi, “Information exposure from consumer iot devices: A multidimensional, network-informed measurement approach,” in _Proceedings of the IMC_ , 2019, pp. 267–279.
  * [14] F. Le, J. Ortiz, D. Verma, and D. Kandlur, _Policy-Based Identification of IoT Devices’ Vendor and Type by DNS Traffic Analysis_. Cham: Springer International Publishing, 2019, pp. 180–201.
  * [15] E. Valdez, D. Pendarakis, and H. Jamjoom, “How to discover iot devices when network traffic is encrypted,” in _2019 IEEE ICIOT_ , pp. 17–24.
  * [16] L. Yi, L. Shiyao, and C. Deyuan, “The reliability of bluetooth data transmission in mobile medical information acquisition system,” in _2016 IEEE IMCEC_ , pp. 498–505.
  * [17] H. Aksu, A. S. Uluagac, and E. S. Bentley, “Identification of wearable devices with bluetooth,” _IEEE Trans. Serv. Comput._ , vol. 6, no. 2, pp. 221–230, 2018.
  * [18] J. Barcelö, L. Montero, L. Marqués, and C. Carmona, “Travel time forecasting and dynamic origin-destination estimation for freeways based on bluetooth traffic monitoring,” _Transp. Res. Rec._ , vol. 2175, no. 1, pp. 19–27, 2010.
  * [19] D. Spill and A. Bittau, “Bluesniff: Eve meets alice and bluetooth.” _WooT_ , vol. 7, pp. 1–10, 2007.
  * [20] W. Albazrqaoe, J. Huang, and G. Xing, “Practical bluetooth traffic sniffing: Systems and privacy implications,” in _Proceedings of the 14th MobiSys_ , 2016, pp. 333–345.
  * [21] P. Di Marco, P. Park, M. Pratesi, and F. Santucci, “A bluetooth-based architecture for contact tracing in healthcare facilities,” _J. Sens. Actuator Netw._ , vol. 10, no. 1, p. 2, 2020.
  * [22] M. A. Sofi, “Bluetooth protocol in internet of things (iot), security challenges and a comparison with wi-fi protocol: A review,” _IJETR_ , vol. 5, 2016.
  * [23] “nrf sniffer for bluetooth le,” https://www.nordicsemi.com/Products/ Development-tools/nrf-sniffer-for-bluetooth-le, accessed: 2023-11-10.
  * [24] “Bluetooth low energy link laye,” https://www.wireshark.org/docs/dfref/b/btle.html, accessed: 2023-11-10.


§ INTRODUCTION

A multivariate probability distribution for discrete variables with finite state
spaces can be represented by a multi-dimensional array.  However, for
high-dimensional distributions where each variable may have a large
state space, lack of computer memory can become a problem. For
example, an $80$-dimensional random vector in which each variable
has $10$ levels will lead to a state space with $10^{80}$ cells. Such
a distribution can not be stored in a computer; in fact, $10^{80}$ is
one of the estimates of the number of atoms in the universe. However,
if the array consists of only a few non-zero values, we need only store
these values along with information about their location. That is, a sparse representation of a table. We describe here the R <cit.> package <cit.> for efficient multiplication and marginalization of sparse tables.

The family of graphical models <cit.> is vast, and includes many different models. In this paper, we consider Bayesian networks and decomposable undirected graphical models. Undirected graphical models are also known as Markov random fields (MRFs). Decomposability is a property ensuring a closed form of the maximum likelihood parameters. Graphical models, enjoy the property that conditional independencies can be read off from a
graph consisting of nodes, representing the random variables. In Bayesian networks, edges are directed from a node to another and represents a directed connection. In a MRF the edges are undirected and these should be regarded as associations between pairs of nodes in a broad sense. A decomposable MRF can be turned into a Bayesian network while retaining the correct likelihood function. However, one can no longer interpret conditional independencies in this transformed graph. This transformation is merely a step towards making inference in decomposable MRFs and as such, we shall mainly focus on introducing concepts regarding Bayesian networks.

We illustrate using the package <cit.> which implements the Junction Tree Algorithm (JTA) for discrete variables using the Lauritzen-Spiegelhalter updating scheme <cit.> with table operations relying on . JTA is used for inference in Bayesian networks (BNs). We also show, that multiplication and marginalization of sparse tables is fast compared with standard built-in R functions and is comparable to when the tables are sparse.

In addition to , there are to our knowledge three other packages
for belief propagation in R; gRain <cit.>,
<cit.> and RHugin <cit.> where the
latter is not on the Comprehensive R Archive Network (CRAN). The
only R package on CRAN that has an API for (dense) table
operations is , which gRain depends upon. Some R
packages that rely on gRain and are geneNetBP
<cit.> and bnspatial <cit.>. The
bnclassify package <cit.> has a lower C++
class implementation of conditional probability tables. Our goal is
twofold: Firstly, to provide an efficient back-end for sparse table
operations for other R packages. Secondly, to provide a new implementation of JTA, namely , using as a back-end.

In Section <ref> we introduce basic notation and terminology and in Section <ref> we motivate the usage of sparse tables through JTA. Section <ref> serves as a primer to our novel representation of tables and their algebra given in Section <ref>. In Section <ref> we also demonstrate how to use . Section <ref> outlines how to use and gives specific examples, using two BNs which are well known from the literature. Moreover, we demonstrate the strength of using the <cit.> package to fit a decomposable MRF on real data. In Section <ref> we show, that the trade-off between execution time and memory allocation using is acceptable for small and medium sized tables, and comparable to in high dimensional sparse tables.

Finally, we mention that leverages from the RcppArmadillo package <cit.> by implementing compute-intensive procedures in C++ for better run-time performance.

§ NOTATION AND TERMINOLOGY

Let $p$ be a discrete probability distribution for the random vector $X = X_{V} = (X_{v} \mid v \in V)$ where $V$ is a set of labels of the random variables $X_{v}, v\in V$ and where each $X_{v}$ has finite state space. A realized value $x = (x_{v})_{v\in V}$ is called a cell. If we denote by $I_{v}$ the finite state space of variable $X_{v}$, the domain of $p$ is given by

\[
I = \times_{v\in V} I_{v}.
\]

Given a subset $A$ of $V$, the A-marginal cell of $x_{V}$ is the vector, $x_{A} = (x_v)_{v\in A}$, with state space $I_{A} = \times_{v\in A} I_{v}$. A Bayesian Network can be defined as a directed acyclic graph (DAG), see Figure <ref>, for which each node represents a random variable together with a joint probability of the form

\begin{align}
p(x) = \prod_{v\in V} p(x_{v} \mid x_{pa(v)}),
\end{align}

where $x_{pa(v)}$ denotes the parents of $x_v$; i.e. the set of nodes with an arrow pointing towards $x_v$ in the DAG. Also, $x_v$ is said to be a child of the variables $x_{pa(v)}$. Notice, that $p(x_{v} \mid x_{pa(v)})$ has domain $I_{v} \times I_{pa(v)}$. Hence, we can encode the conditional probabilities in a table, say $\phi(x_{v}, x_{pa(v)})$, of dimension $|I_{v}|\cdot |I_{pa(v)}|$. It is also common in the literature to refer to these tables as potentials and we shall use these terms interchangeably. In general, a potential does not have an interpretation. Sometimes, we also use subscript notation to explicitly show the set of variables for which a potential depends on. That is, $\phi_{A}$ is a potential defined over the variables $X_{A}$. The product, $\phi_{A} \otimes \phi_{B}$, of two generic tables over $A$ and $B$ is defined cell-wise as

\[
(\phi_{A} \otimes \phi_{B})(x_{A\cup B}) := \phi_{A}(x_{A})\phi_{B}(x_{B}).
\]

In other words, the product is defined over the union of the variables of each of the two potentials. Division of two tables, $\phi_{A} \oslash \phi_{B}$, is defined analogously. The marginal table, $\phi_{A}^{\downarrow B}$, over the variables $B\subseteq A$ is defined cell-wise as

\[
\phi_{A}^{\downarrow B}(x_B) := \sum_{x_{A\setminus B} \in I_{A\setminus B}}\phi_{A}(x_{B}, x_{A\setminus{B}}).
\]

Finally, for some $B \subseteq A$, fix $x^{\ast}_{B}$. The $x^{\ast}_{B}$ slice of $\phi_A(x)$ is then given by

\[
\phi_{A}^{x^{\ast}_B}(x_{A\setminus B}) = \phi_{A}(x_{A\setminus B}, x^{\ast}_{B}).
\]

§ MOTIVATION THROUGH MESSAGE PASSING IN BAYESIAN NETWORKS

Consider the simple DAG given in Figure <ref>, from which the joint density can be read off:

\begin{align}\label{eq:pot_rep}
p(x_a, x_b, x_c, x_d, x_e) &= p(x_c)p(x_a \mid x_c)p(x_b \mid x_a, x_d)p(x_d \mid x_c)p(x_e \mid x_c, x_d)% \notag\\
%                             &= \phi(x_c)\phi(x_a, x_c)\phi(x_b, x_a, x_d)\phi(x_d, x_c)\phi(x_e, x_c, x_d).
\end{align}

If, for example, interest is in the joint distribution of $(x_a, x_d)$ we have to sum over $x_b, x_c$ and $x_e$ and exploiting the factorization we could calculate this as

\begin{align*}
p(x_a, x_d) = \sum_{x_c} p(x_c) p(x_a \mid x_c) p(x_d \mid x_c)\sum_{x_e}p(x_e \mid x_c, x_d)\sum_{x_b}p(x_b \mid x_a, x_d).
\end{align*}

= 0.8;
_ = .5cm;
cir/.style=thick, circle, draw = black, minimum size = _, inner sep = 0

[cir, label = above:$a$] (a) at (-1.5*, );
[cir, label = above:$b$] (b) at (1.5*, );
[cir, label = below:$c$] (c) at (-1.5*, -1.5*);
[cir, label = below:$d$] (d) at (1.5*, -1.5*);
[cir, label = below:$e$] (e) at (0*,  -3*);

[-Latex[width=3mm], thick, draw=black] (a) – (b);
[-Latex[width=3mm], thick, draw=black] (c) – (a);
[-Latex[width=3mm], thick, draw=black] (c) – (d);
[-Latex[width=3mm], thick, draw=black] (c) – (e);
[-Latex[width=3mm], thick, draw=black] (d) – (b);
[-Latex[width=3mm], thick, draw=black] (d) – (e);

= 0.8;
_ = .5cm;
cir/.style=thick, circle, draw = black, minimum size = _, inner sep = 0

[cir, label = above:$a$] (a) at (-1.5*, );
[cir, label = above:$b$] (b) at (1.5*, );
[cir, label = below:$c$] (c) at (-1.5*, -1.5*);
[cir, label = below:$d$] (d) at (1.5*, -1.5*);
[cir, label = below:$e$] (e) at (0*,  -3*);

[thick, draw=black] (a) – (b);
[thick, draw=black] (a) – (d);
[thick, draw=black] (c) – (a);
[thick, draw=black] (c) – (d);
[thick, draw=black] (c) – (e);
[thick, draw=black] (d) – (b);
[thick, draw=black] (d) – (e);

= 1;
_ = .7cm;
dot/.style=thick, circle, draw = black, minimum size = _, inner sep = 0, fill = black
cir/.style=thick, circle, draw = black, minimum size = _, inner sep = 2

[cir, label = left:$C_1$] (C1) at (0*, -1.5*)$a,b,d$;
[cir, label = left:$C_2$] (C2) at (0*, -3*)$a,c,d$;
[cir, label = left:$C_3$] (C3) at (0*, -4.5*)$c,d,e$;

[thick, draw=black, ->] (C3) – (C2);
[thick, draw=black, ->] (C2) – (C1);

(a) A DAG. (b) A moralized and triangulated version of (a). (c) A rooted junction tree representation of (b) with root $C_1$.

The Junction Tree Algorithm (JTA) can be seen as an algorithm for automatically factorizing in order to circumvent the direct summation as described in what follows using a minimal example <cit.>: First

* moralize the DAG; i.e. connect nodes that share a common child node,

* remove directions in the DAG to obtain an undirected graph and

* triangulate the resulting graph.

Moralization ensures that the corresponding parent and child nodes are put in the same maximal clique. A clique is a subset of the nodes for which the induced subgraph is complete and it is maximal if it is not contained in any other clique. From here, by clique, we always mean a maximal clique and we refer to these as the cliques of the graph.

A graph is triangulated if it has no cycles of length greater than $3$. If such cycles are present, we must add fill edges to produce a triangulated graph. A triangulated graph is also called decomposable, and hence the connection to decomposable MRFs shows here. Finding an optimal triangulation (in terms of minimizing the number of fill edges) is a NP hard problem, but there exists good heuristic methods, see <cit.>. Finding a good triangulation can have huge impact on the performance of JTA which we detail in Section <ref>. A moralized and triangulated version of the graph in Figure <ref> is shown Figure <ref>, where no fill edges was necessary to make the graph triangulated.

A triangulated graph can always be represented as a junction tree. A junction tree is a tree where the nodes are given by the cliques of the triangulated graph with the property that for two cliques, $C$ and $C'$, the intersection $C \cap C'$ is contained in all clique nodes on the unique path between $C$ and $C'$. The cliques of the graph in Figure <ref> are given as $C_3 = \{c,d,e\}, C_2 = \{a,c,d\}$ and $C_1 = \{a,b,d\}$ where we arbitrarily designate $C_1$ as the root to obtain the rooted junction tree in Figure <ref>. Now, assign each potential in (<ref>) to a clique potential for which the variables conform, e.g.

\begin{align*}
\phi_{C_1}(x_a, x_b, x_d) &\leftarrow p(x_a \mid x_b, x_d), \\
\phi_{C_2}(x_a, x_c, x_d) &\leftarrow p(x_c)p(x_a \mid x_c), \\
\phi_{C_3}(x_c, x_d, x_e) &\leftarrow p(x_d \mid x_c)p(x_e \mid x_c, x_d).
\end{align*}

The clique potentials are now initialized (we also say that the network has been initialized) and note that the clique potentials in general do not have any interpretation at this stage. We have obtained the clique potential representation

\begin{equation}
\label{eq:clique_pot_rep}
p(x_a, x_b, x_c, x_d, x_e) = \phi_{C_1}(x_a, x_b, x_d)\phi_{C_2}(x_a, x_c, x_d)\phi_{C_3}(x_a, x_c, x_d).
\end{equation}

The network is said to be compiled at this stage, i.e. when moralization and triangulization has been performed and the clique potential representation is obtained. In complex networks with large clique potentials it might not be feasible to even initialize the clique potentials due to lack of memory. We give an example of how to overcome this in Section <ref> by entering evidence into the model as described in Section <ref>.

Next, the message passing scheme can now be applied to the junction tree. We describe here, the Lauritzen-Spiegelhalter (LS) scheme which works as follows. Locate a leaf node, here we choose $C_{3}$, and find the intersection, $S_{32} = C_3 \cap C_2 = \{c, d\}$, with its parent clique $C_2$. Then calculate the marginal potential $\phi_{S_{32}}(x_c, x_d) = \sum_{x_e}\phi_{C_3}(x_c, x_d, x_e)$ and perform an inward message by setting

\[
\phi_{C_{2}}(x_a, x_c, x_d) \leftarrow \phi_{C_{2}}(x_a, x_c, x_d)\phi_{S_{32}}(x_c, x_d)
\]

and update the leaf node as

\begin{equation}
\label{eq:update_c2}
\phi_{C_3}(x_c, x_d, x_e) \leftarrow \phi_{C_3}(x_c, x_d, x_e) / \phi_{S_{32}}(x_c, x_d),
\end{equation}

where $0/0:=0$. We say that $C_2$ have collected its messages from all of its children. This procedure must be repeated until the root, $C_1$, has collected all its messages. Hence, we perform another inwards message by setting $\phi_{S_{21}}(x_c) = \sum_{x_a, x_d}\phi_{C_{2}}(x_a, x_c, x_d)$ and update:

\begin{align*}
\phi_{C_{1}}(x_a, x_b, x_d) &\leftarrow \phi_{C_{2}}(x_a, x_c, x_d)\phi_{S_{21}}(x_c, x_d), \\
\phi_{C_{2}}(x_a, x_c, x_d) &\leftarrow \phi_{C_{2}}(x_a, x_c, x_d) / \phi_{S_{21}}(x_c, x_d).
\end{align*}

The inward phase terminates when the root clique potential has been normalized:

\[
\phi_{C_1}(x_a, x_b, x_d) \leftarrow \phi_{C_1}(x_a, x_b, x_d) / \sum_{x_a, x_b, x_d}\phi_{C_{1}}(x_a, x_b, x_d).
\]

To summarize, we have now obtained what is called the set chain representation

\begin{align*}
\label{eq:setchain}
p(x_a, x_b, x_c, x_d, x_e) &= \phi_{C_1}(x_a, x_b, x_d) \phi_{C_2}(x_a, x_c, x_d) \phi_{C_3}(x_c, x_d, x_e)\\
&= p(x_a, x_b, x_d)p(x_a \mid x_c, x_d)p(x_c, x_d \mid x_e),
\end{align*}

where the clique potentials are now conditional probability tables. Notice especially that $\phi_{C_1}(x_a, x_b, x_d) = p(x_a, x_b, x_d)$. In the outward phase we start by sending messages from the root by performing an outward message by letting $\phi_{S_{12}}(x_a, x_d) = \sum_{x_b}\phi_{C_1}(x_a, x_b, x_d)$ and update:

\begin{equation}
\phi_{C_2}(x_a, x_c, x_d) \leftarrow \phi_{C_{2}}(x_a, x_c, x_d)\phi_{S_{12}}(x_c, x_d).
\end{equation}

We say that $C_1$ has distributed evidence to $C_2$. Notice, that $\phi_{C_2}$ is now identical to the probability distribution defined over the variables $x_a, x_c$ and $x_d$. Finally, let $\phi_{S_{23}}(x_c, x_d) = \sum_{x_a}\phi_{C_2}(x_a, x_c, x_d)$ and update $\phi_{C_3}$ as

\[
\phi_{C_3}(x_c, x_d, x_e) \leftarrow \phi_{C_3}(x_c, x_d, x_e) \phi_{S_{23}}(x_c, x_d).
\]

As a consequence we finally obtain the clique marginal representation

\begin{align*}
p(x_a, x_b, x_c, x_d, x_e) &= \frac{\phi_{C_{1}}(x_a, x_b, x_d)\phi_{C_{2}}(x_a, x_c, x_d)\phi_{C_3}(x_c, x_d, x_e)}{\phi_{S_{12}}(x_a, x_d)\phi_{S_{23}}(x_c, x_d)} \\
&= \frac{p(x_a, x_b, x_d)p(x_a, x_c, x_d)p(x_c, x_d, x_e)}{p(x_a, x_d)p(x_c, x_d)},
\end{align*}

where all clique and separator potentials are identical to the marginal probability distribution over the variables involved. Hence we can now find $p(x_a, x_d)$ by locating a clique containing $x_a$ and $x_d$ and sum out all other variables. If we choose $C_2$ we get

\[
p(x_a, x_d) = \sum_{x_c}\phi_{C_{2}}(x_a, x_c, x_d).
\]

Each time we multiply, divide or marginalize potentials, a number of binary operations (addition, multiplication and division) are conducted under the machinery. For a network with $41$ variables and a maximum size of the state space for each variable being $3$, <cit.> recorded a total number of $2,371,178$ binary operations. We do not intend to follow the same analysis here, but for sparse tables the number of necessary binary operations is potentially much smaller.

§.§ Evidence and Slicing

Suppose it is known, before message passing, that $X_{E} = x^{\ast}_{E}$ for some labels $E \subset V$. We refer to $x^{\ast}_{E}$ as evidence. Evidence can be entered into the clique potential representation (<ref>) as follows. For each $v \in E$ choose an arbitrary clique, $C$, where $v\in C$ and set entries in $\phi_{C}$ that are inconsistent with $x^{\ast}_{v}$ equal to zero. The resulting clique potential is then said to be sliced. After message passing, all queries are then conditional on $X_{E} = x^{\ast}_{E}$. Thus entering evidence leads to more zero-cells, and in a sparse setup, the resulting clique potentials will be even more sparse. After message passing, the clique potential $\phi_{C}(x_{C})$ is now equal to the conditional probability $p(x_{C\setminus E} \mid x_E^\ast)$.

It suffices to modify a single clique potential such that it is inconsistent with $v \in E$, for all $v$ as described above. However, for sparse tables it is advantageous to enter evidence in all clique potentials containing $v$ since this leads to a higher degree of sparsity. This is how evidence is handled in . In fact, this is the reason why is able to handle very complex networks by exploiting evidence using sparse tables from . It is also possible to enter evidence into the factorization (<ref>). This is the key to handle complex networks that are otherwise infeasible due to lack of memory as we show in Section <ref>. This trick is related to cutset conditioning and if one can enter evidence into variables that breaks cycles the effect can be huge <cit.>. One can exploit the law of total probability to sum out the variables that was conditioned on (since we obtain the probability of the evidence during propagation) and recover whatever probablity is of interest.

§ AN INTUITIVE WAY OF REPRESENTING SPARSE TABLES

Before describing our method for multiplication and marginalization of sparse tables, it is illuminating to describe sparse tables in a standard R language setup. Consider two arrays f and g:

R> dn <- function(x) setNames(lapply(x, paste0, 1:2), toupper(x))
R> d  <- c(2, 2, 2)
R> f  <- array(c(5, 4, 0, 7, 0, 9, 0, 0), d, dn(c("x", "y", "z")))
R> g  <- array(c(7, 6, 0, 6, 0, 0, 9, 0), d, dn(c("y", "z", "w")))

with flat layouts

R> ftable(f, row.vars = "X")

Y y1    y2
Z z1 z2 z1 z2
x1    5  0  0  0
x2    4  9  7  0

R> ftable(g, row.vars = "W")

Y y1    y2
Z z1 z2 z1 z2
w1    7  0  6  6
w2    0  9  0  0

If we convert f and g to data.frame objects and exclude the cases with a value of zero we have:

R> df <- as.data.frame.table(f, stringsAsFactors=FALSE)
R> df <- df[df$Freq != 0,]
R> dg <- as.data.frame.table(g, stringsAsFactors=FALSE)
R> dg <- dg[dg$Freq != 0,]

R> print(df, row.names = FALSE)

X  Y  Z Freq
x1 y1 z1    5
x2 y1 z1    4
x2 y2 z1    7
x2 y1 z2    9

R> print(dg, row.names = FALSE)

Y  Z  W Freq
y1 z1 w1    7
y2 z1 w1    6
y2 z2 w1    6
y1 z2 w2    9

which literally leaves us with two sparse tables, df and dg respectively. In order to multiply df by dg we must, by definition, determine the cases that match on the variables Y,Z that they have in common. For example, row $4$ in df must be multiplied with row $4$ in dg such that (y1, z2, x2, w2) is an element in the product with value $81$. And since the tables are sparse, no multiplication by zero will be performed. The multiplication can be performed with the following small piece of R code (which will be used in Section <ref> in connection with bench-marking):

R> sparse_prod <- function(df, dg)
+    S   <- setdiff(intersect(names(df), names(dg)), "Freq")
+    mrg <- merge(df, dg, by = S)
+    val <- mrg$Freq.x * mrg$Freq.y
+    mrg$val <- val
+    mrg[, setdiff(names(mrg), c("Freq.x", "Freq.y"))]
+  }
\end{Sinput}
\end{Schunk}
The \code{merge} function performs, by default, what is also called an \textit{inner join} or \textit{natural join} in SQL terminology, which is exactly how we defined table multiplication in Section \ref{sec:notation}. Multiplying \code{df} and \code{dg} yields
\begin{Schunk}
\begin{Sinput}
R> sparse_prod(df, dg)
\end{Sinput}
\begin{Soutput}
Y  Z  X  W val
1 y1 z1 x1 w1  35
2 y1 z1 x2 w1  28
3 y1 z2 x2 w2  81
4 y2 z1 x2 w1  42
\end{Soutput}
\end{Schunk}
Marginalization is even more straightforward. Marginalizing out $X$ from \code{df} can for example be done using the built-in \proglang{R} function \texttt{aggregate} (which is also used in Section \ref{sec:benchmark} for benchmarking):
\begin{Schunk}
\begin{Sinput}
R> aggregate(Freq ~ Y + Z, data = df, FUN = sum)
\end{Sinput}
\begin{Soutput}
Y  Z Freq
1 y1 z1    9
2 y2 z1    7
3 y1 z2    9
\end{Soutput}
\end{Schunk}
Thus, we have the necessary tools to implement JTA using sparse tables. So why should we bother redefining sparse tables and algebras on these; because of execution time and memory storage. In Section \ref{sec:benchmark} we show the effect of the effort of going beyond the \code{merge} and \code{aggregate} functions.

\section{Sparse Tables}
\label{sec:sparse_tables}
Let $T$ be a dense table with domain $I = ×_v∈VI_v$ as described in \ref{sec:notation}. Define the \textit{level set} $ℒ :=×_v∈Vℒ_v$ where $ℒ_v = {1, 2, …, |I_v|}$ and let $# : I →ℒ$ be a bijection. We define the sparse table $τ= (Φ, ϕ)$, of $T$ as the pair where $Φ$ is a matrix with columns given by the set of vectors in the \textit{sparse domain} $ℐ := {#(x) |T(x) ≠0, x ∈I}$, consisting of \textit{non-zero cells} and where $ϕ$ is the corresponding vector of values. Thus, a column in $Φ$ represents a cell in $ℐ$ and is written a tuple $i = (i_1, i_2, …, i_|V|; i_v ∈ℒ_v)$ which explicitly determines the ordering of the labels and hence the order of the rows in $Φ$. The order of the columns in $Φ$ is not important as long as it agrees with $ϕ$. We denote by $Φ[j]$ the $j'$th column of $Φ$ and by $ϕ_j$ the corresponding $j'$th value in $ϕ$. The sub-matrix $Φ_S$ defined over the set of labels, $S ⊆V$, is the resulting matrix when rows corresponding to labels in $V∖S$ have been removed. Let $T$ be the table \code{f} from Section \ref{sec:intui}:
\begin{Schunk}
\begin{Soutput}
Y y1    y2
Z z1 z2 z1 z2
x1    5  0  0  0
x2    4  9  7  0
\end{Soutput}
\end{Schunk}
The domain is given by
\[
I = \{x_1, x_2\} \times \{y_1, y_2\} \times \{z_1, z_2\}
\]
and we can choose $#$ as the map $(x_ℓ_1, y_ℓ_2, z_ℓ_3) ↦(ℓ_1, ℓ_2, ℓ_3)$ for $ℓ_1, ℓ_2, ℓ_3 ∈{1,2}$. The non-zero cells can be identified from the table and we have $ℐ = {(1,1,1), (2,1,1), (2,2,1), (2,1,2)}$. Hence
\[
\Phi = \begin{bmatrix}
1&  2& 2& 2\\
1&  1& 2& 1\\
1&  1& 1& 2
\end{bmatrix},
\]
with $ϕ= (5, 4, 7, 9)$, corresponding to \code{df} in Section \ref{sec:intui}. Let $G$ be another dense table with domain $J = ×_u∈UJ_u$ and sparse representation $γ= (Ψ, ψ)$ with sparse domain $𝒥$. We then aim at defining the sparse multiplication $τ⊗γ$ of $T ⊗G$. Let $S = V ∩U$ be the \textit{separator labels} shared between the two sparse tables $τ$ and $γ$. Next, define the map $M_S(Φ)$ which transform $Φ_S$ into a look-up table\footnote{A lookup table is a list arranged as key-value pairs. In \proglang{R} one can think of a look-up table as a named list where the names are the keys and the values are the elements of the list.} as follows: the keys are the unique columns of $Φ_S$ and the value of $M_S(Φ)$ at key $k$ is the set of column indices where column $k$ can be found in $Φ_S$ and hence also in $Φ$ and is given by
\[
M_{S}(\Phi)[k] = \{ j \in \{1, 2, \ldots, |\mathcal{I}| \} : \Phi[j] = k \}.
\]
Let $𝒦$ denote the mutual keys of $M_S(Φ)$ and $M_S(Ψ)$. The number of columns in the matrix of the resulting product $τ⊗γ$ is then given as
\[
N := \sum_{k\in\mathcal{K}} |M_{S}(\Phi)[k]| \cdot |M_{S}(\Psi)[k]|.
\]
This observation is crucial, since the memory storage of the sparse product can then be computed in advance. If $(Π, π)$ is the sparse product of $τ$ and $γ$ we can therefore initialize $Π$ as a matrix with $|V| + |U∖V|$ rows and $N$ columns and $π$ as an $N-$dimensional vector. Finally, $π$ is given by the values $ϕ_j ·ψ_j'$ for $j ∈M_S(Φ)[k]$ and $j' ∈M_S(Ψ)[k]$ for all $k ∈𝒦$. The procedure is formalized in Algorithm \ref{alg:mult}.

\begin{algorithm}[h!]
\caption{Multiplication of Sparse Tables}\label{alg:mult}
\begin{algorithmic}[1]
\Procedure{}{$\tau = (\Phi, \phi)$: sparse table, $\gamma = (\Psi, \psi)$: sparse table}
\State $S := V \cap U$
\State $\mathcal{K}$: Mutual keys of $M_{S}(\Phi)$ and $M_{S}(\Psi)$
\State $N := \sum_{k\in\mathcal{K}} |M_{S}(\Phi)[k]| \cdot |M_{S}(\Psi)[k]|$
\State Initialize the matrix $\Pi$ with $|V| + |U\setminus V|$ rows and $N$ columns
\State Initialize the vector $\pi$ of dimension $N$
\State $l := 1$
\For{$k \in \mathcal{K}$}
\For{$j\in M_{S}(\Phi)[k]$ and $j' \in M_{S}(\Psi)[k]$}
\State $\Pi[l] := (\Phi[j], \Psi_{U\setminus V}[j'])$
\State $\pi_l  := \phi_{j} \cdot \psi_{j'}$ \label{alg:mult_pie}
\State $l = l + 1$
\EndFor
\EndFor
\State\Return $(\Pi, \pi)$
\EndProcedure
\end{algorithmic}
\end{algorithm}

The number of binary operations is smaller than the equivalent dense
table multiplication since every multiplication with zero is
avoided. Moreover, since we only loop over the mutual keys,
$𝒦$, the execution time will depend on the table having the
least unique keys over the separator labels. Trivially, division of
two sparse tables can be obtained by changing line \ref{alg:mult_pie}
of Algorithm \ref{alg:mult} with $π_l = ϕ_j / ϕ_j'$.

Now, let $G$ be the table \code{g} from Section \ref{sec:intui} where the domain is given by
\[
J = \{w_1, w_2\} \times \{y_1, y_2\} \times \{z_1, z_2\}.
\]
Choose the map $(w_ℓ_1, y_ℓ_2, z_ℓ_3) ↦(ℓ_1, ℓ_2, ℓ_3)$ for $ℓ_1, ℓ_2, ℓ_3 ∈{1,2}$. Then, in summary we have the tables
\[
\Phi = \begin{bmatrix}
1&  2& 2& 2\\
1&  1& 2& 1\\
1&  1& 1& 2
\end{bmatrix}, \quad
\Psi = \begin{bmatrix}
1&  2& 2 & 1\\
1&  1& 2 & 2\\
1&  1& 1 & 2
\end{bmatrix},
\]
along with $ϕ= (5,4,7,9)$ and $ψ= (7,6,6,9)$. The separator labels are given by $S = {y, z}$ and the lookup tables of $Φ$ and $Ψ$ are then given by
\begin{align*}
M_{S}(\Phi) &= \{ (1,1):= \{1, 2\} , (2,1): = \{3\}, (1,2):=\{4\} \} \\
M_{S}(\Psi) &= \{ (1,1):= \{1\} , (2,1): = \{2\}, (1,2):=\{4\}, (2,2):=\{3\}\}.
\end{align*}
Above, $y$ corresponds to row two and $z$ corresponds to row three in $Φ$. So for example, $M_S(Φ)[(1,1)] = {1,2}$ means, that the key $(1,1)$ has value ${1,2}$ which means that $(1,1)$ is found in columns $1$ and $2$ in $Φ$.

Therefore, all values $ϕ_j$ for $j∈M_S(Φ)[(1,1)]$ must be multiplied with all values $ψ_j$ for $j∈M_S(Ψ)[(1,1)]$ etc. Hence,
\[
\Pi = \begin{bmatrix}
1&  2& 2& 2 \\
1&  1& 2& 1 \\
1&  1& 1& 2 \\
1&  1& 1& 2
\end{bmatrix},
\]
\[
\pi = (\phi_1 \cdot \psi_1, \phi_2 \cdot \psi_1, \phi_3 \cdot \psi_2, \phi_4 \cdot \psi_4) = (35,28,42,81),
\]
as expected from the result in Section \ref{sec:intui} using \code{sparse_mult}. Notice, that we save any computation with $ψ_3$ since $(2,2)$ is not a key in $M_S(Φ)$.

We mention that addition and subtraction of sparse tables is more demanding since we have to reconstruct zero-cells if one of the tables
have a non-zero cell-value while the other table has a zero-cell in
the corresponding separator cell. Fortunately, these operations are not
needed in JTA.

The marginal sparse table $τ^↓A = (Φ^↓A, ϕ^↓A)$ of $τ$, corresponding to $T^↓A$, can be calculated using the map $M_A(Φ)$ and, for each key $k ∈M_A(Φ)$, sum the corresponding values in $ϕ$. However, for massive tables the memory footprint of $M_A(Φ)$ is unnecessarily large. Instead, we construct the lookup-table $H_A(Φ)$ where the keys are the unique columns of $Φ_A$ as was the case in $M_A(Φ)$. However, the values are themselves pairs where the first element is an index to any of the column indices where the corresponding key can be found in $Φ_A$. The second element is the final cell value in the marginalized table corresponding to the key. The pair corresponding to the key $k$ is therefore on the form
\[
H_{A}(\Phi)[k] = (j, v), \quad v = \sum_{\ell : \Phi_{A}[\ell] = k} \phi_\ell \mbox{ and } \Phi_{A}[j] = k.
\]
The value $v$ can easily be computed iteratively. The point here is, that we never have to store $Φ_A$ since we can deduce all information from $Φ$ on the fly given the row indices corresponding to $A$ in $Φ$. The number of columns in the final matrix $Φ^↓A$, and hence the number of elements in $ϕ^↓A$, is given by $|H_A(Φ)|$. The procedure is formalized in Algorithm \ref{alg:marg}.
\begin{algorithm}[h!]
\caption{Marginalization of Sparse Tables}\label{alg:marg}
\begin{algorithmic}[1]
\Procedure{}{$\tau = (\Phi, \phi)$: sparse table, $A$: Set of labels}
\State Construct $H_{A}(\Phi)$
\State $N = |H_{A}(\Phi)|$
\State Initialize the matrix $\Phi^{\downarrow A}$ with $|A|$ rows and $N$ columns
\State Initialize the vector $\phi^{\downarrow A}$ of dimension $N$
\State Let $\mathcal{K}$ be the keys of $H_{A}(\Phi)$
\State $l:=1$
\For{$k \in \mathcal{K}$}
\State $(j, v) := H_{A}(\Phi)[k]$
\State $\Phi^{\downarrow A}[l] := \Phi_{A}[j]$ \Comment deduced by picking elements from $\Phi[j]$
\State $\phi^{\downarrow A}_l := v$
\State $l = l+1$
\EndFor
\State\Return $(\Phi^{\downarrow A}, \phi^{\downarrow A})$
\EndProcedure
\end{algorithmic}
\end{algorithm}
Consider again the sparse table $τ$ of $T$ and let $A = {y, z}$. Then, the resulting sparse marginal table have two rows corresponding to $y$ and $z$. The construction of $H_A(Φ)$ is as follows. The first column in $Φ$ is extracted and the entry corresponding to $x$ is deleted. Call the resulting vector (key) $k_1$. Now, set $H_A(Φ)[k_1] = (j = 1, v = 5)$ since $ϕ_1 = 5$. Extract now, the second column of $Φ$ and let $k_2$ be the resulting key when removing the entry corresponding to $x$. Since $k_1 = k_2$ and $ϕ_2 = 4$ we update $H_A(Φ)[k_1] = (j = 2, v = 9)$. Proceeding this way, gives
\[
H_{A}(\Phi) = \{(1,1) := (j = 2, v = 9), (2,1):=(j = 3, v = 7), (1,2):= (j = 4, v = 9)\}\}.
\]
\[
\Phi^{\downarrow A} = \begin{bmatrix}
1&  2& 1 \\
1&  1& 2
\end{bmatrix},
\]
and $ψ^↓A = (9, 7, 9)$. For $B ⊂V$, the \textit{$i^{\ast}_{B}-$ slice} of a sparse table $τ= (Φ, ϕ)$ is obtained by removing columns, $k$, in $Φ$ for which $k$ does not agree with $i^∗_B$. We leave out the formal procedure for slicing.

\subsection{Sparse tables in JTA}
Triangulation may lead to cliques for which no potential is assigned during compilation and the clique potential is then regarded as the unity table. In the dense version, this corresponds to an array filled with ones. However, the idea of sparsity will be ruined if we just adapt the dense version. For the general case we define a \textit{sparse unity table} and show how to multiply it with a non-unity table in Appendix \ref{app:sparse_unities}. In the light of message passing, these unity tables are not needed. When a message is sent to a unity, we simply replace it with the message. And as such, we avoid a lot of computational overhead.

\subsection{How to use sparta}
\label{sec:spartatables}

In order to demonstrate the use of \sparta we revisit the example from Section \ref{sec:intui} of the two (dense) tables \code{f} and \code{g} with mutual variables, \code{Y} and \code{Z}
\begin{multicols}{2}
\begin{Schunk}
\begin{Sinput}
R> ftable(f, row.vars = "X")
\end{Sinput}
\begin{Soutput}
Y y1    y2
Z z1 z2 z1 z2
x1    5  0  0  0
x2    4  9  7  0
\end{Soutput}
\begin{Sinput}
R> ftable(g, row.vars = "W")
\end{Sinput}
\begin{Soutput}
Y y1    y2
Z z1 z2 z1 z2
w1    7  0  6  6
w2    0  9  0  0
\end{Soutput}
\end{Schunk}
\end{multicols}
We can convert these to their equivalent \code{sparta} versions as
\begin{Schunk}
\begin{Sinput}
R> sf <- as_sparta(f); sg <- as_sparta(g)
\end{Sinput}
\end{Schunk}
Printing \code{sf}
\begin{Schunk}
\begin{Sinput}
R> print.default(sf)
\end{Sinput}
\begin{Soutput}
[,1] [,2] [,3] [,4]
X    1    2    2    2
Y    1    1    2    1
Z    1    1    1    2
[1] 5 4 7 9
[1] "x1" "x2"

[1] "y1" "y2"

[1] "z1" "z2"

[1] "sparta" "matrix"
\end{Soutput}
\end{Schunk}
where we use the default print method to look under the hood. The columns are the cells in the sparse matrix and the \code{vals} attribute are the corresponding values which can be extracted with the \code{vals} function. Furthermore, domain resides in the \code{dim_names} attribute which can also be extracted using the \code{dim_names} function. From the output, we see that (\code{x2}, \code{y2}, \code{z1}) has a value of $2$. Using the regular print method prettifies things:
% TODO: Make dispatch to generic dimnames instead of dim_names
\begin{Schunk}
\begin{Sinput}
R> print(sf)
\end{Sinput}
\begin{Soutput}
X Y Z val
1 1 1 1   5
2 2 1 1   4
3 2 2 1   7
4 2 1 2   9
\end{Soutput}
\end{Schunk}
where row $i$ corresponds to column $i$ in the sparse matrix. We settled for this print method because printing column wise leads to unwanted formatting when the values are decimal numbers.
Consider the cell (2,1,1). The corresponding named cell is then
\begin{Schunk}
\begin{Sinput}
R> get_cell_name(sf, sf[, 2L])
\end{Sinput}
\begin{Soutput}
X    Y    Z
"x2" "y1" "z1"
\end{Soutput}
\end{Schunk}
where \code{sf[, 2L]} is the second column (row in the output) of \code{sf} which is $(2,1,1)$. The product of \code{sf} and \code{sg} is
\begin{Schunk}
\begin{Sinput}
R> mfg <- mult(sf, sg); mfg
\end{Sinput}
\begin{Soutput}
X Y Z W val
1 2 1 2 2  81
2 2 2 1 1  42
3 1 1 1 1  35
4 2 1 1 1  28
\end{Soutput}
\end{Schunk}
The equivalent dense table has $2^4=16$ entries. However, \code{mfg} stores $20$ values after all, $16$ of which are information about the cells. That is, there is some overhead storing the information about the cells, see Section \ref{sec:when}. Converting \code{sf} into a conditional probability table (CPT) with conditioning variable \code{Z}:
\begin{Schunk}
\begin{Sinput}
R> sf_cpt <- as_cpt(sf, y = "Z"); sf_cpt
\end{Sinput}
\begin{Soutput}
X Y Z   val
1 1 1 1 0.312
2 2 1 1 0.250
3 2 2 1 0.438
4 2 1 2 1.000
\end{Soutput}
\end{Schunk}
We can further slice \code{sf_cpt} to obtain a probability distribution conditioned on \code{Z = z1}
\begin{Schunk}
\begin{Sinput}
R> slice(sf_cpt, s = c(Z = "z1"))
\end{Sinput}
\begin{Soutput}
X Y Z   val
1 1 1 1 0.312
2 2 1 1 0.250
3 2 2 1 0.438
\end{Soutput}
\end{Schunk}
Marginalizing out \code{Y} in \code{sg} yields
\begin{Schunk}
\begin{Sinput}
R> marg(sg, y = c("Y"))
\end{Sinput}
\begin{Soutput}
Z W val
1 2 2   9
2 2 1   6
3 1 1  13
\end{Soutput}
\end{Schunk}
which is in correspondence with the example in Section \ref{sec:sparse_tables}. Finally, we mention that a sparse table can be created using the constructor \code{sparta_struct} which can be necessary to use if the corresponding dense table is too large to have in memory. A table with relevant methods in \sparta is given in Table \ref{tab:sparta_funcs}.

% Magical spacing
% \vspace{2.5cm}

\begin{table}[h!]
\centering
\begin{tabular}{@{}ll@{}}
\midrule
Function name              & Description                                          \\ \midrule
\code{as_<sparta>}         & Convert \code{array}-like object to a \code{sparta}  \\
\code{as_<array/df/cpt>}   & Convert \code{sparta} object to an \code{array/data.frame/CPT} \\
\code{sparta_struct}       & Constructor for \code{sparta} objects                \\
\code{mult, div, marg, slice} & Multiply/divide/marginalize/slice                 \\
\code{normalize}           & Normalize (the values of the result sum to one)      \\
\code{get_val}             & Extract the value for a specific named cell           \\
\code{get_cell_name}       & Extract the named cell                                \\
\code{get_values}          & Extract the values                                   \\
\code{dim_names}           & Extract the domain                                   \\
\code{names}               & Extract the variable names                           \\
\code{max/min}             & The maximum/minimum value                            \\
\code{which_<max/min>_cell}& The column index referring to the max/min value      \\
\code{which_<max/min>_idx} & The configuration corresponding to the max/min value \\
\code{sum}                 & Sum the values                                       \\
\code{equiv}               & Test if two tables are identical up to permutations of the columns \\ \midrule
\end{tabular}
\caption{Description of essential functions from the \sparta package.}
\label{tab:sparta_funcs}
\end{table}

\newpage

\subsection{When to Use Sparta}
\label{sec:when}
As shown in Section~\ref{sec:spartatables}, there is an overhead of
storing the information in a \sparta object. A dense array with $x$
elements takes up $8x$ bytes plus some negligible memory of storing
the variable names and etc. On the contrary, a \sparta object with
$y < x$ elements takes up $y(4k + 8)$ bytes, where $k$ is the number
of variables (these can be stored as integers and hence only requires
$4$ bytes each). In Figure~\ref{fig:when}, we have plotted this
relation for $k = 4, 6$ and $8$, and different levels of
sparsity. That is, a sparsity of $1/2$ implies that $y = x/2$. The
black identity line indicates the number of gigabytes needed to store
the dense table with $x$ elements. The size of the state spaces of the
variables are implicitly reflected by the memory needed to store the
dense table. The more memory needed, the bigger state space of the
variables. However, the more variables and the bigger state space of
these variables will intuitively result in a more sparse table, making
\sparta efficient even for several variables.

The take away message from Figure \ref{fig:when} is that when the
state space of the variables and the sparsity increases the benefit of
storing the tables using \sparta will outweigh the overhead of storing
the additional information.

In connection to JTA, \sparta is favorable when cliques with many variables
imply a high degree of sparsity. In particular, this is often the case for tables in a Bayesian network representing a pedigree. In this case, cliques tend to be small but the state space of the variables can be arbitrarily large due to the large amount of DNA information for each member of the pedigree. In Section \ref{sec:mrf} we fit a decomposable MRF to real data and show that the sparsity of the clique potentials is much favorable towards \sparta.

\begin{figure}[t] \centering
\includegraphics[width=1\textwidth]{img/sparta_when.pdf}
\caption{The black identity line indicates the number of gigabytes
needed to store the dense table with $x$ elements. The colored
lines indicates the number of gigabytes required to store the
equivalent \sparta object with the respective number of variables
and sparsity.}
\label{fig:when}
\end{figure}

\section{Usecases of jti and sparta}
\label{sec:usecases}

In \jti there are two ways of specifying a Bayesian network. Either by a list of CPTs or a dataset together with a DAG. In the latter case, the CPTs are found using maximum likelihood estimates. Here, we describe how to use \jti using the classic Bayesian network \code{asia} \citep{lauritzen1988local} where the corresponding CPTs is part of \jti. The network represents a simplified model to help diagnose the patients arriving at a respiratory clinic. A history of smoking has a direct influence on both whether or not a patient has bronchitis and whether or not a patient has lung cancer. Both lung cancer and bronchitis can result in dyspnoea. An x-ray result depends on the presence of either tuberculosis or lung cancer. Finally, a visit to Asia influence the probability of having tuberculosis. The DAG is depicted in Figure \ref{fig:asia}.
\begin{figure}[!h] \centering
\includegraphics[width=0.3\textwidth]{img/asia_graph.pdf}
\caption{The DAG for the \texttt{asia} network.}
\label{fig:asia}
\end{figure}

We use the version of \texttt{asia} called \texttt{asia2} which is a list of CPTs and which is shipped with \jti. The first step is to call \texttt{cpt\_list} for some initial checks and conversion to \sparta tables:
\begin{Schunk}
\begin{Sinput}
R> cl <- cpt_list(asia2)
R> cl
\end{Sinput}
\begin{Soutput}
List of CPTs
P( asia )
P( tub | asia )
P( smoke )
P( lung | smoke )
P( bronc | smoke )
P( either | lung, tub )
P( xray | either )
P( dysp | bronc, either )

<cpt_list, list>
\end{Soutput}
\end{Schunk}
From the output we see the inferred CPTs corresponds to Figure \ref{fig:asia} giving rise to a factorization in the same way as in (\ref{eq:pot_rep}). The network is now ready for compilation which involves moralization and triangulation. In \jti there are four different choices for triangulation which are all based on the \textit{elimination game} algorithm, see \citep{flores2007review}. One of the most well-known heuristics is \code{min_fill} which tries to minimize the number of fill edges. Evidence can be entered either at compile stage or just before message passing begins. It is always advisable to enter evidence at compile stage since we know from Section \ref{sec:evidence_slicing} that this reduces the number of non-zero elements in the CPTs and hence the memory footprint and execution time. A good strategy might be to locate one or more of the largest cliques and enter evidence on the nodes contained in these. We can investigate the cliques and their state spaces prior to compilation by triangulating the graph as follows:
\begin{Schunk}
\begin{Sinput}
R> tri <- triangulate(cl, tri = "min_fill")
\end{Sinput}
\end{Schunk}
The \code{tri} object is a list containing the triangulated graph, \code{new_graph} as a matrix, a list of \code{fill_edges}, the \code{cliques} and the size of the dense \code{statespace} of each clique. In Section \ref{sec:impact_of_evidence} we show how to exploit information about the cliques with largest state space in order to compile and propagate in a network that is otherwise infeasible on a 'standard' laptop. Now, let for example \code{tub = yes} be the evidence indicating that a given person has tuberculosis. The compiled network is then constructed as

\begin{Schunk}
\begin{Sinput}
R> cp <- compile(cl, evidence = c(tub = "yes"), tri = "min_fill")
R> cp
\end{Sinput}
\begin{Soutput}
Compiled network
Nodes: 8
Cliques: 6
- max: 3
- min: 2
- avg: 2.67
- tub: yes
<charge, list>
\end{Soutput}
\end{Schunk}
The cliques can be extracted from the compiled object with \code{get_cliques(cp)}. The compiled object can now be entered into the message passing procedure as:
\begin{Schunk}
\begin{Sinput}
R> j <- jt(cp)
\end{Sinput}
\end{Schunk}
The junction tree can be visualized by plotting the object as \code{plot(j)}, see Figure \ref{fig:asia_junc}. Finally, we can calculate the probability of a given person having a positive x-ray result, \code{xray = yes}, given that the person has tuberculosis as
\begin{Schunk}
\begin{Sinput}
R> query_belief(j, nodes = "xray")
\end{Sinput}
\begin{Soutput}
yes   no
0.98 0.02
\end{Soutput}
\end{Schunk}

Thus, given that a person has tuberculosis, the probability of observing a positive x-ray result is 0.98. The probability of observing positive x-ray result given that \code{tub = no} can be calculated accordingly and equals 0.1012. Joint queries can be calculated by specifying \code{type = "joint"} in \code{query_belief}.

\begin{figure}[h!] \centering
\includegraphics[width=0.3\textwidth]{img/asia_junction_tree_graph.pdf}
\caption{A junction tree for the asia network.}
\label{fig:asia_junc}
\end{figure}

\newpage

\subsection{Inference in Decomposable Markov Random Fields}
\label{sec:mrf}
In this section we illustrate the gain of using \sparta in connection with the public-domain \texttt{derma} data set (part of the \ess package) originally obtained from the UCI Machine Learning Repository \citep{Dua:2019}. The data set consists of $358$ observations, $35$ variables of which there is one class variable called \texttt{ES}. The class variable has six states; each representing a skin disease. The remaining $34$ clinical variables (all with $4$ states except \texttt{age} which has been discretized into six bins) are used to predict the skin disease. We first fit a decomposable MRF
\begin{Schunk}
\begin{Sinput}
R> g <- ess::fit_graph(derma, q = 1.5, sparse_qic = TRUE)
\end{Sinput}
\end{Schunk}
The fitted graph is plotted in Figure \ref{fig:derma}, where we just notice, that it is a rather complex graph with large cliques. Hence, for such a small data set and a graph with big cliques, the chance of sparse clique potentials is significant.
\begin{figure}[!h] \centering
\includegraphics[width=0.5\textwidth]{img/derma_graph.pdf}
\caption{The fitted decomposable MRF for the derma data set.}
\label{fig:derma}
\end{figure}
We can now convert \code{g} into a Bayesian network simply by entering into the \code{cpt_list} function and from there compile and propagate the network:
\begin{Schunk}
\begin{Sinput}
R> cl <- jti::cpt_list(derma,  ess::as_igraph(g))
R> j  <- jti::jt(jti::compile(cl, root_node = "ES"), propagate = "no")
\end{Sinput}
\end{Schunk}
The argument \code{root_node = "ES"} forces \code{ES} into the root clique. Hence, if we are interested in queries about \code{ES} we only need to conduct the inwards message passing, collect, since the root clique potential is identical to the probability distribution over the variables in the root clique. The argument \code{propagate = "no"} means that messages passing should be postponed. By doing so, initialization is only performed once and the junction tree can be exploited as many times as needed with different evidence.

The largest clique of \code{j} has 6 variables and the mean sparsity of the clique potential is given by
\begin{Schunk}
\begin{Sinput}
R> mean(sapply(j$charge$C, sparta::sparsity))
\end{Sinput}
\begin{Soutput}
[1] 0.9233838
\end{Soutput}
\end{Schunk}
showing that the tables are extremely sparse. The junction tree can now be used for e.g. classification. Consider the first observation in \code{derma}
\begin{Schunk}
\begin{Sinput}
R> z <- unlist(derma[1, -ncol(derma)])
\end{Sinput}
\end{Schunk}
and update the junction tree with evidence corresponding to \code{z} and perform inwards message passing
\begin{Schunk}
\begin{Sinput}
R> jz <- jti::propagate(jti::set_evidence(j, z), prop = "collect")
\end{Sinput}
\end{Schunk}
The probability distribution of the class variable given the evidence \code{z} is then
\begin{Schunk}
\begin{Sinput}
R> q <- jti::query_belief(jz, nodes = "ES", type = "marginal")[[1]]
R> q[which.max(q)]
\end{Sinput}
\begin{Soutput}
seboreic dermatitis
\end{Soutput}
\end{Schunk}
Hence, we classify \code{z} as an observation from class seboreic dermatitis.

\subsection{The Impact of Evidence}
\label{sec:impact_of_evidence}
The Bayesian Network \code{Link} \citep{jensen1999885} is a large Bayesian network with $724$ nodes and $1,125$ arcs. The network has been used for linkage analysis where, previously, only approximate methods has been applied to make inference in the network.

We have verified that \sparta generates a large amount of overhead in the CPTs for the \code{Link} network. In fact, \jti uses approximately $16$ times more memory to store the clique potentials than the dense version. However, in the following we show how \sparta leverages from entering evidence and hence reducing the CPTs.

We picked a handful of the most comprehensive software packages, across different programming languages, for belief propagation and test if they were able to handle \code{Link} on a ``standard'' laptop machine. We used the \proglang{R} package \gRain, the \proglang{R} package \pkg{BayesNetBP}, the \proglang{Python} package \pkg{pyArum} \citep{gonzales2017agrum}, the \proglang{Julia} package \pkg{BayesNet} \citep{bayesnets}. Interestingly, all of these failed due to lack of computer memory on the first authors machine\footnote{See Section \ref{sec:comp_details} for details about this machine.} while \jti succeeded. The \pkg{pyAgrum} package is a high-level interface to the extremely efficient \proglang{C++} library \pkg{aGrum} \citep{gonzales2017agrum}. In \cite{gonzales2017agrum} the authors of \pkg{pyAgrum} themselves failed to propagate in \code{Link} on a computer with $32$Gb memory. An investigation revealed that exploiting the triangulation found by \jti, \pkg{gRain} is now also able to make inference in \code{Link}. We have not tried this with the other packages mentioned but conjecture, that the other packages would also be able to handle the \code{Link} network. Interestingly, this triangulation was the result of one of the most well-known heuristics, namely \code{min_fill}. The reason that triangulation have such a big impact is, that it determines the resulting cliques and hence the size of the clique potentials. If the cliques are too big, and hence the number of elements in the state space is large, it may not be possible to initialize the clique potentials due to lack of memory.

We now compare the cliques with larges state spaces resulting from the two heuristics \code{min_fill} and \code{min_nei} (a variant of \code{min_fill} which also tries to minimize the number of fill-in edges). We have extracted the \code{Link} network from \url{https://www.bnlearn.com/bnrepository/} as a list of CPTs with information about child and parent-node relations. In \jti we just need a list of CPTs which we extract as
\begin{Schunk}
\begin{Sinput}
R> cpts <- bnfit_to_cpts(Link)
\end{Sinput}
\end{Schunk}

Next, we convert these CPTs to their \sparta equivalent using \code{cpt_list} while deducing the underlying network, conducting some sanity checks and triangulate with the two heuristics
\begin{Schunk}
\begin{Sinput}
R> cl <- cpt_list(cpts)
R> tri_min_fill <- triangulate(cl, tri = "min_fill")
R> tri_min_nei  <- triangulate(cl, tri = "min_nei")
\end{Sinput}
\end{Schunk}
The size of the five largest clique state spaces can then be extracted as
\begin{Schunk}
\begin{Sinput}
R> sp_min_fill <- sort(tri_min_fill$statespace, decreasing = TRUE)[1:5]
R> sp_min_nei  <- sort(tri_min_nei$statespace, decreasing = TRUE)[1:5]
\end{Sinput}
\end{Schunk}
and the memory usage, for dense tables, in gigabytes of allocating memory for each of these are
\begin{Schunk}
\begin{Sinput}
R> round(sum(sp_min_fill) / 1e9*8, 2)
\end{Sinput}
\begin{Soutput}
[1] 0.2
\end{Soutput}
\begin{Sinput}
R> round(sum(sp_min_nei) / 1e9*8, 2)
\end{Sinput}
\begin{Soutput}
[1] 26.98
\end{Soutput}
\end{Schunk}

\newpage

This indicates, that any software for message passing should be able to handle the link network on a standard laptop using the \code{min_fill} heuristic from \jti, but also that \code{min_nei} will fail.

We imagine that there is no triangulation of \code{Link}, such that any software packages, including \jti, are able to compile and propagate in \code{Link} on a standard laptop. In what follows we show, that by exploiting evidence and sparsity, we can overcome the initialization step by using \jti. We use the \code{min_nei} heuristic since this was the worst, locate one of the cliques with largest state space and insert evidence on $10$ nodes.

\begin{Schunk}
\begin{Sinput}
R> idx_max  <- which.max(tri_min_nei$statespace)
R> max_vars <- tri_min_nei$cliques[[idx_max]]
R> max_dim  <- jti::dim_names(cl)[max_vars[1:10]]
R> e        <- sapply(max_dim, `[[`, 1L) # get first state of all variables
R> cp       <- jti::compile(cl, e, tri = "min_nei")
R> j        <- jti::jt(cp)
\end{Sinput}
\end{Schunk}
The memory size (Gb), in the dense sense, in the resulting five largest clique potentials is
\begin{Schunk}
\begin{Sinput}
R> dense_sp <- function(x) prod(sapply(sparta::dim_names(x),length))
R> round(sort(sum(sapply(j$charge$C, dense_sp) / 1e9*8), TRUE), 2)
\end{Sinput}
\begin{Soutput}
[1] 0.06
\end{Soutput}
\end{Schunk}
which is negligible compared to 26.98 Gbs. The memory size (Gb) of the sparse tables is
\begin{Schunk}
\begin{Sinput}
R> round(sum(sort(sapply(j$charge$C, ncol), TRUE)[1:5]) / 1e9*8, 2)
\end{Sinput}
\begin{Soutput}
[1] 0.02
\end{Soutput}
\end{Schunk}

\section{Time and Memory Tradeoff in Sparta}
\label{sec:benchmark}
We investigate the tradeoff between memory allocation and execution
time for multiplication and marginalization on \sparta objects. If
\sparta objects do not perform reasonably well for small and medium
sized tables their usage is limited in real usecases.

Thus, we compare three functions for multiplication and three
functions for marginalization: 1) The \code{mult} and \code{marg}
functions from \sparta, 2) the \code{tabMult} and \code{tabMarg}
functions from \gRbase and 3) the \code{sparse_prod} and \code{aggregate} functions from Section \ref{sec:intui}. In the latter case we just refer to \pkg{R} as the package.

We randomly generate pairs of tables such that the number of cells in
the product of the two tables does not exceed $10^6$. We varied the
sparsity of the product of the tables; $0\%$ sparsity (dense tables),
$(1-75]\%$ sparsity, and  $(75, 99]\%$ sparsity. For each pair of tables, we multiply them together and record the memory usage (in megabytes) of the product and the execution time (in seconds). As \gRbase is the standard and most
mature package for graphical models in \proglang{R}, the performance
comparisons are relative to that of \gRbase.  Hence, in
Figure~\ref{fig:relative}, \gRbase performs better for tables with the
relative scores above one (indicated by horizontal dashed lines),
whereas values below one show cases where the alternative approaches
are better. The comparisons are plotted for different ranges of table
sizes (panels) and sparsity of the resulting table (first axis).

\paragraph{Multiplication (size):}

The first row of Figure~\ref{fig:relative} describes the size of the
table resulting from multiplying two tables. It can be seen that
\pkg{R} consistently produces tables of smaller sizes than \sparta for
very small tables with $100$ ($10^2$) cells. For tables with more than
$100$ cells, \sparta consistently produces smaller tables except for a
single case. Increasing the degree of sparsity leads to reduced object
sizes for both \pkg{R} and \sparta, and for tables with more than
$75\%$ sparsity, \sparta outperforms both \pkg{R} and \grbase expect
in the first two panels with small tables.

\paragraph{Multiplication (time):}

The second row of Figure~\ref{fig:relative} describes the computing
time for multiplying two tables. Clearly \sparta outperforms \pkg{R}
by orders of magnitude. For larger tables, (the two rightmost
columns), there is also a clear effect of the degree of sparsity on
the computing time.

\paragraph{Marginalisation (time):}

The third row of Figure~\ref{fig:relative} describes the computing
time for marginalising a table. When the degree of sparsity increases,
the computing time decreases. In the comparison between \gRbase and
\sparta, we see that the marginalization implementation of \sparta is
competitive to that of \gRbase. For tables with $10,000$ ($10^4$) or
fewer cells it is faster irrelevant of the sparsity except in a single
case. For 75\% sparsity \sparta's marginalization is consistently

\begin{figure}[!ht] \centering
\includegraphics[width=\textwidth]{img/relative.pdf}
\caption{\label{fig:relative}Relative comparison of \pkg{R}
(\texttt{sparse\_prod} and \texttt{aggregate}) and \sparta
(\texttt{mult} and \texttt{marg}) to \gRbase (\texttt{tabMult} and
\texttt{tabMarg}) in terms of memory usage and timing. The top row
shows the relative memory usage for multiplication, whereas the
two lower rows shows the relative timing for multiplication and
marginalizaion, respectively.}
\end{figure}

This small, and non-exhaustive, benchmark study indicates, that our
proposed method for table multiplication and marginalization performs
well for small, medium and large tables. However, our real interest is
in the performance on massive tables which is impossible to benchmark
in this paper due to the increased running time and memory usage of
\gRbase and \pkg{R}.

\newpage

\section{Summary}
We have presented a novel method for multiplication and
marginalization of sparse tables. The method is implemented in the
\proglang{R} package \sparta. However the method is generic and we
have provided detailed pseudo algorithms facilitating the extension to
other languages. In addition we presented the companion package \jti
to illustrate some of the advantages of \sparta in connection to the
Junction Tree Algorithm. We hope to explore the benefit from the
\proglang{C} API for working with external pointers in order to reduce
the memory usage for \sparta objects in the future.

The memory footprint of the clique potentials can become prohibitively large when the size of the cliques are large This may not be true in general for sparse tables. As a matter of fact, it may be optimal to have large cliques if they are very sparse and/or if it is common to observe variables in such a clique. Finally, we hope to explore pedigree networks that can potentially benefit immensely from \sparta in the future.

\section*{Computational details}
\label{sec:comp_details}
In addition to the packages already mentioned, the following \proglang{R} packages was used to make the benchmark results

\begin{itemize}
\item \code{dplyr} version $1.0.6$
\item \code{glue} version $1.4.2$
\item \code{tictoc} version $1.0$
\item \code{ggplot2} version $3.3.3$
\end{itemize}

We used \proglang{R} version $4.1.0$. All computations were carried out on a $64$-bit Linux computer with Ubuntu $20.04.2$ and Intel(R) Core(TM) i7-6600U CPU 2.60GHz LTS. The machine has approximately $6$Gb of free memory for use in calculations.

\appendix
\section{Sparse unity tables}
\label{app:sparse_unities}
We define a \textit{sparse unity table}, $\mu$, over the variables $U$ to equal the set of level sets $\{\mathcal{L}_{u}\}_{u\in U}$ induced by the domain $J$. Multiplying a sparse table, $\tau = (\Phi, \phi)$, defined over the set of labels $V$ with sparse domain $\mathcal{I}$ and the sparse unity table, $\mu$, can be performed as follows. Define $R = U\setminus V$ and  $\mathcal{L}_{R} = \times_{r\in R}\mathcal{L}_{r}$ and initialize the product matrix $\Pi$ with $|V| + |R|$ rows and $|\mathcal{I}| \cdot |\mathcal{L}_{R}|$ columns. In other words, $\Pi$ has $|\mathcal{L}_{R}|$ times more columns than $\Phi$. For $j = 1, 2, \ldots, |\mathcal{I}|$ and $i_{R} \in \mathcal{L}_{R}$ append $(\Phi[j], i_{R})$ to $\Pi$ and set the corresponding value in $\pi$ to $\phi_j$. That is, we multiply a sparse table with a sparse unity table on the fly. When $|I| \gg |\mathcal{I}|$ this approach greatly reduces the number of binary operations compared to the dense case; in fact, there are no binary operations. The procedure is summarized in Algorithm \ref{alg:mult_unity}. However, dividing $\mu$ with $\tau$, one must change line \ref{alg:mult_unity_pie} in Algorithm \ref{alg:mult_unity} with $\pi_l := 1 / \phi_j$.

\begin{algorithm}[h!]
\caption{Multiplication of a Sparse Table with a Sparse Unity Table}\label{alg:mult_unity}
\begin{algorithmic}[1]
\Procedure{}{$\tau = (\Phi, \phi)$: sparse table, $\mu$: sparse unity table}
\State $R := U \setminus V$
\State Set $N = |\mathcal{I}| \cdot |\mathcal{L}_{R}|$
\State Initialize the matrix $\Pi$ with $|V| + |R|$ rows and $N$ columns
\State Initialize the vector $\pi$ of dimension $N$
\State $l := 1$
\For{$j = 1, 2, \ldots, |\mathcal{I}|$}
\For{$i_{R} \in \mathcal{L}_{R}$}
\State $\pi[l] := (\Phi[j], i_{R})$
\State $\pi_l := \phi_j$ \label{alg:mult_unity_pie}
\State $l = l + 1$
\EndFor
\EndFor
\State\Return $(\Pi, \pi)$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\newpage
\input{article.bbl}

\end{document}