# Performance Optimization of Deep Learning Sparse Matrix Kernels on Intel Max
Series GPU

Mohammad Zubair
Old Dominion University
Norfolk, Virginia, USA
<EMAIL_ADDRESS>
& Christoph Bauinger
Intel Corporation
Santa Clara, CA, USA
<EMAIL_ADDRESS>

###### Abstract

In this paper, we focus on three sparse matrix operations that are relevant
for machine learning applications, namely, the sparse-dense matrix
multiplication (SPMM), the sampled dense-dense matrix multiplication (SDDMM),
and the composition of the SDDMM with SPMM, also termed as FusedMM. We develop
optimized implementations for SPMM, SDDMM, and FusedMM operations utilizing
Intel oneAPI’s Explicit SIMD (ESIMD) SYCL extension API. In contrast to CUDA
or SYCL, the ESIMD API enables the writing of explicitly vectorized kernel
code. Sparse matrix algorithms implemented with the ESIMD API achieved
performance close to the peak of the targeted Intel Data Center GPU. We
compare our performance results to Intel’s oneMKL library on Intel GPUs and to
a recent CUDA implementation for the sparse matrix operations on NVIDIA’s V100
GPU and demonstrate that our implementations for sparse matrix operations
outperform either.

_Keywords_ Machine Learning $\cdot$ Optimization $\cdot$ Sparse Matrix
Operations $\cdot$ ESIMD $\cdot$ Intel Data Center GPU

##  1 Introduction

The key to building and utilizing increasingly large machine learning models
is efficient implementations of the training and inferencing on emerging high-
performance computing architectures, mainly graphics processing units (GPUs).
Due to their high throughput in single precision (FP32) and reduced precision
arithmetic (FP8, FP16, BF16, TF32, etc.), which is typically delivered by
specialized hardware built into the devices, GPUs have become the standard
workhorse for machine learning applications [1, 2]. In this context, three
sparse matrix operations are of particular interest. Namely, the sparse-dense
matrix multiplication (SPMM), the sampled dense-dense matrix multiplication
(SDDMM), and the composition of the SDDMM with SPMM, also termed as FusedMM
operation [3, 4, 5, 6, 7, 8]. The SPMM operation involves the multiplication
of a sparse matrix with a dense matrix, resulting in a dense matrix. The
SDDMM, in turn, multiplies a dense matrix with another dense matrix with the
constraint that not all the entries of the resultant matrix are needed. In
other words, the output matrix of the SDDMM operation is sparse. The sparsity
structure is given as input to the SDDMM operation that specifies what entries
of the output matrix need to be computed. The FusedMM operation merges the
SDDMM and SPMM operations into a single operation, where the output of the
SDDMM operation, which is a sparse matrix, is used as input to the following
SPMM operation. Fusing the operations can increase the performance and can
thus be beneficial for applications such as sparse transformer [4], where SPMM
follows the SDDMM operation.

Sparse matrix operations have been extensively studied in scientific computing
[9]. Methods such as finite element (FE), finite volume (FV), or finite
difference (FD) results in sparse matrices with a high degree of sparsity
(i.e., the fraction of zeros in the matrix). More precisely, in these methods,
the density of the matrices, which describes the fraction of non-zero entries
in a matrix, is often less than $1\%$ of the overall number of entries. In
scientific applications, sparse matrices are typically involved in solving
linear sparse systems that require either sparse matrix-vector operations or
sparse matrix factorizations [9, 10, 11]. On the other hand, the sparse
matrices arising in machine learning are relatively dense, where the density
can vary from $10\%$ to $30\%$ [3, 4]. Nevertheless, the matrices are
sufficiently sparse — especially in the cases with a density close to $10\%$ —
to warrant exploration of sparse matrix operations with compressed storage
formats [4, 6].

Several researchers have explored efficient implementations of sparse matrix
operations, including SDDMM, SPMM, and FusedMM, on emerging high-performance
architectures GPUs [4, 6, 5, 12]. The major challenge in optimizing the
performance of sparse matrix operations on GPUs is effectively utilizing the
available memory bandwidth. We found that an implicit SIMD programming
environment such as CUDA [13] or SYCL [14], where the compiler performs the
vectorization implicitly, introduces difficulties when trying to achieve
performance close to the theoretical peak hardware capability. Thus, we
optimize sparse matrix operations utilizing Intel oneAPI’s Explicit SIMD SYCL
extension (ESIMD) API [15]. In contrast to CUDA or SYCL, the ESIMD API enables
the writing of explicitly vectorized kernel code. Due to this, it allows more
precise control over register usage and better handles thread divergence
compared to CUDA or SYCL. In prior contributions [16], it was shown that
kernel code written with the ESIMD API can perform close to the peak of the
targeted Intel Data Center GPU hardware [17]. The main disadvantage of ESIMD,
in contrast to SYCL, is the lack of support for non-Intel GPUs.

In this paper, we develop optimized implementations for the SPMM, SDDMM, and
FusedMM operations for Intel Data Center GPUs utilizing the ESIMD API. In
Section 2, we describe the sparse matrix operations, Intel Data Center GPU,
and Intel ESIMD API. The GPU algorithms for the sparse matrix operations are
covered in Section 3. Section 4 then demonstrates the optimized ESIMD
implementations of the matrix operations. Finally, Section 5 compares our
results with a recent implementation for the two sparse matrix operations on
an NVIDIA V100 GPU presented in [4] and shows that we outperform the reference
implementation [4] on V100 by up to a factor 10. In addition, with our
implementation, we increase the device utilization from up to 27% of the
theoretical single-precision peak of a V100 GPU (cf. [4]) to up to 49% on the
Intel hardware. Finally, we compare our results to the implementations of the
sparse operations available in Intel’s oneMKL library and show that our
performance outperforms the state-of-the-art MKL implementations by up to a
factor of 3.4 on Intel’s Data Center GPU.

##  2 Background

In this section, we give an overview of the relevant sparse matrix operations,
introduce the ESIMD API as well as the targeted Intel Data Center GPU Max
1550.

###  2.1 Sparse Matrix Operations

Table 1 shows the notation used in this paper to describe the sparse matrix
operations. It also lists the matrix sizes investigated in the reference
implementation in [4] and which commonly occur in deep learning networks.

$M\in\mathbb{N}$ |  Value of $M$ ranges from $1024\eqqcolon 1k$ to $32768\eqqcolon 32k$
---|---
$K\in\mathbb{N}$ |  Value of $K$ ranges from $1024\eqqcolon 1k$ to $8192\eqqcolon 8k$
$N\in\mathbb{N}$ |  Value of $N$ is typically $32$ or $128$
$\mathbf{A}\in\mathbb{R}^{M\times K}$ |  A sparse matrix of size $M\times K$
$\mathbf{B}\in\mathbb{R}^{K\times N}$ |  A dense matrix of size $K\times N$
$\mathbf{B}^{T}\in\mathbb{R}^{N\times K}$ |  The transpose of $\mathbf{B}$, with size $N\times K$
$\mathbf{D}\in\mathbb{R}^{K\times N}$ |  A dense matrix of size $K\times N$
$\mathbf{C}\in\mathbb{R}^{M\times N}$ |  A dense matrix of size $M\times N$
$\mathbf{E}\in\mathbb{R}^{M\times N}$ |  A dense matrix of size $M\times N$
$\mathbf{I}_{A}\in\mathbb{R}^{M\times K}$ |  A sparse matrix of size $M\times K$, where all the non-zero entries have a value of 1
$\mathbf{O}_{i}\in\mathbb{R}^{1\times\kappa}$ |  The row $i$, $1\leq i\leq\mu$, of any matrix $\mathbf{O}\in\mathbb{R}^{\mu\times\kappa}$.
$\mathbf{O}_{i,j}\in\mathbb{R}$ |  The element with the row index $i$, $1\leq i\leq\mu$, and the column index $j$, $1\leq j\leq\kappa$, of any matrix $\mathbf{O}\in\mathbb{R}^{\mu\times\kappa}$.
$\alpha\in[0,1)\subset\mathbb{R}$ |  Sparsity. The fraction of matrix elements which are zero.
$nnz\in\mathbb{N}$ |  number of non-zero elements in the matrix $\mathbf{A}$. Computed as $nnz=MK\beta$
Table 1: Summary of the relevant quantities and the notation used in the
description of the algorithms.

#### SDDMM

The SDDMM operation is defined as follows.

$\displaystyle\mathbf{A}$ $\displaystyle=$
$\displaystyle\mathbf{C}\mathbf{B}^{T}\odot\mathbf{I}_{A}$ (1)

For this operation, a dense matrix $\mathbf{C}$ is multiplied with a dense
matrix $\mathbf{B}^{T}$ followed by an element-wise product, which is
indicated by the operation $\odot$, with a matrix $\mathbf{I}_{A}$. In a
typical implementation, a dot product of a row of the matrix $\mathbf{C}$ with
a column of $\mathbf{B}^{T}$ is computed only at a location of a non-zero
entry in $\mathbf{I}_{A}$.

#### SPMM

The SPMM operation, on the other hand, multiplies a sparse matrix $\mathbf{A}$
with a dense matrix $\mathbf{B}$, resulting in a dense matrix $\mathbf{C}$, as
shown below.

$\displaystyle\mathbf{C}$ $\displaystyle=$ $\displaystyle\mathbf{A}\mathbf{B}$
(2)

#### FusedMM

The FusedMM operation is a composition of SDDMM with SPMM. The composite
operation avoids storing the result of the SDDMM operation explicitly and can
have a better performance compared to an implementation with the two
operations applied one after another.

$\displaystyle\mathbf{E}$ $\displaystyle=$
$\displaystyle(\mathbf{C}\mathbf{B}^{T}\odot\mathbf{I}_{A})\mathbf{D}$ (3)

#### Matrix Storage Layout

All dense matrices are stored in row-major order. The sparse matrices are
stored in a compressed sparse row format, which is described in what follows.

The sparsity of the matrix $\mathbf{A}$ is denoted by a value
$1>\alpha\in\mathbb{R}$, which indicates the fraction of zero-entries in the
matrix $\mathbf{A}$. In deep learning networks [3, 4], the sparsity typically
varies from $0.7$ to $0.9$. The compressed sparse row (CSR) format [9] avoids
the explicit storage of zeros in the matrix. In this format, only non-zero
entries in the matrix are stored along with two integer arrays to resolve the
row and column indices of the non-zero values. This format reduces the memory
footprint and operation count for sparse matrix operations at the expense of
indirect addressing. The CSR format for the matrix $\mathbf{A}$ consists of
three one-dimensional arrays: $avalue\in\mathbb{R}^{nnz}$,
$ia\in\mathbb{N}^{M+1}$, and $ja\in\mathbb{N}^{nnz}$. The size of the $avalue$
and $ja$ arrays is $nnz$, the number of non-zeros in the sparse matrix. The
array $avalue$ contains the non-zero entries in $\mathbf{A}$ in row-major
order, and the $ja$ array contains column indices of those values. The array
$ia$ is of size $M+1$ whose $i$-th entry indicates the index in $avalue$ and
$ja$ where the $i$-th row of $\mathbf{A}$ starts. The array $ia$ includes a
fictitious $M+1$-th entry to facilitate easy traversal of the elements through
the last row $M$. Figure 1 shows a sample sparse matrix with the corresponding
CSR arrays.

$\begin{bmatrix}~{}&~{}&1&2\\\ ~{}&~{}&3&~{}\\\ 4&5&~{}&~{}\\\
6&~{}&~{}&~{}\\\ 7&~{}&8&9\end{bmatrix}$

$\displaystyle ia$ $\displaystyle=$ $\displaystyle\left[0,2,3,5,6,9\right]$
$\displaystyle ja$ $\displaystyle=$
$\displaystyle\left[2,3,2,0,1,0,0,2,3\right]$ $\displaystyle avalue$
$\displaystyle=$ $\displaystyle\left[1,2,3,4,5,6,7,8,9\right]$

Figure 1: An example of a sparse matrix of size $M\times K$ with $M=5$ and
$K=4$ on the left-hand side. Note that only non-zero entries are shown. The
three CSR arrays $ia$, $ja$, and $avalue$ required for storing the sparse
matrix are presented on the right-hand side.

###  2.2 Intel Data Center GPU Max 1550

The code in the present contribution targets the nascent Intel Data Center GPU
Max 1550 [17]. In what follows, this GPU is briefly introduced. The Intel Data
Center GPU Max 1550, which is shortened to Intel GPU in what follows, is a HPC
accelerator with 128 gigabyte (GB) high-bandwidth memory (HBM), and a
theoretical peak FP32, and FP64 throughput of approximately 54 tera floating-
point operations per second (Tflops/s) which is delivered by 1024 so-called
vector engines. For machine learning applications the device includes so-
called matrix engines (XMX) delivering a theoretical peak performance of 832
Tflops/s on the bfloat16 data type.

From a performance optimization point of view, it is important to note that
the device is not one monolithic piece of hardware but consists of two Xe
stacks. This is relevant considering that, i) the interconnect between the
stacks is slower than the access to HBM memory, and ii) cross-stack accesses
are not cached. It is therefore recommended for most applications to scale
explicitly to the two stacks using, e.g., MPI or multiple queues. Thus, in
what follows, we focus on the performance of a single stack of the Intel GPU.
For more information on the Intel hardware, we refer to [17, 18]

###  2.3 Intel oneAPI/ESIMD

The ESIMD API [15, 19] is an extension to the SYCL standard developed
specifically for Intel GPUs. It is based on Intel’s GPU Instruction Set
Architecture (ISA). Whereas SYCL relies on the compiler for the vectorization
along the work-items within a sub-group [20] (cf. “warp” in CUDA), in ESIMD
the code uses simd objects (see Fig. 2) for explicit vectorization. These simd
objects enable the vectorization over SIMD-sizes different from the sub-group
size. Thus, ESIMD offers finer control over the vectorization compared to
standard SYCL. Since the simd objects are mapped to the registers, ESIMD
permits fine control over register usage. Additionally, ESIMD provides APIs
for explicit memory load, store, and prefetch operations with parameters to
control the caching behavior, and it simplifies the management of divergent
branches in kernel code.

In contrast to SYCL, ESIMD does not have the concept of sub-groups. Since each
work-item in ESIMD utilizes explicit vectorization, the notion of sub-groups
is not required and each work-item in ESIMD represents, in some sense, a SYCL
sub-group with variable size. Thus, while a work-item in SYCL is equivalent to
a thread in CUDA, a work-item in ESIMD differs.

As shown in Fig. 2, when using Intel oneAPI DPC++, an ESIMD kernel is launched
with the "SYCL_ESIMD_KERNEL" property.

⬇

#define SIMD_LEN 16

\parvoid vecAdd(nd_item<1> item, double *a, double *b, double *c) {

const int i = item.get_global_id(0)*SIMD_LEN;

simd<double, SIMD_LEN> va, vb, vc;

va.copy_from(a + i);

vb.copy_from(b + i);

vc = va + vb;

vc.copy_to(c + i);

}

\parint main() {

queue Q(gpu_selector_v);

double *d_a = malloc_device<double>(16*32, Q);

double *d_b = malloc_device<double>(16*32, Q);

double *d_c = malloc_device<double>(16*32, Q);

Q.parallel_for(nd_range<1>(16*32/SIMD_LEN, 32),

[=](nd_item<1> item) SYCL_ESIMD_KERNEL {

vecAdd(item, d_a, d_b, d_c);

});

}

Figure 2: An example of a vector addition written in ESIMD.

##  3 Overview of the Algorithms

In this section, we discuss our approach to parallelizing sparse matrix
operations without giving details on how these algorithms are realized on a
specific architecture. The implementation details specific to a GPU
architecture are covered in the next section, where we account for the
underlying vector architecture, number of registers, cache behavior, etc.

###  3.1 Parallelizing SDDMM

The SDDMM operation multiplies two dense matrices, $\mathbf{C}$ and
$\mathbf{B}^{T}$, and stores the result in a sparse matrix $\mathbf{A}$. The
parallelization is performed along the elements of the resulting matrix
$\mathbf{A}$. A total of $M$ threads are assigned to work concurrently, where
a thread $i$, $1\leq i\leq M$, computes all non-zeros of the row
$\mathbf{A}_{i}$. Figure 3 illustrates this computation for $\mathbf{A}_{i}$
that has four non-zeros at column indices $j_{0},j_{1},j_{2}$, and $j_{3}$.
The non-zero at $j_{0}$ is computed by multiplying the row $\mathbf{C}_{i}$
with a column $j_{0}$ of $\mathbf{B}^{T}$. The rest of the non-zeros are
computed similarly. A high-level description of the algorithm for computing
$\mathbf{A}_{i}$ is outlined in Algorithm 1. In contrast to the above
description, where we implicitly assumed a regular matrix layout for
$\mathbf{A}$, we use the CSR format for $\mathbf{A}$ in the high-level
description of the algorithm. Please note that input to the SDDMM algorithm is
$\mathbf{B}$. However, we are interested in the product of $\mathbf{C}$ with
$\mathbf{B}^{T}$ (cf. Equation (1)).

Algorithm 1 SDDMM-ROW($ia,ja,avalues,B,C,N,\texttt{i}$)

1:Initialize: $\mathbf{C}[i,:]\leftarrow 0$

2:$nnzr\leftarrow ia[i+1]-ia[i]$

3:for $j\leftarrow 0$ to $nnzr-1$ do

4: $k\leftarrow\texttt{ja}[ia[i]+j]$

5: $dp\leftarrow 0$

6: for $l\leftarrow 0$ to $N-1$ do

7: $dp\leftarrow dp+\mathbf{C}[i,l]*\mathbf{B}[k,l]$

8: end for

9: $\texttt{avalues}[ia[i]+j]\leftarrow dp$

10:end for

11:return avalues

The parallelism in SDDMM operation can be increased by assigning more than one
thread to compute non-zeros of $\mathbf{A}_{i}$. For example, we can assign
two threads to compute the non-zeros of a row $\mathbf{A}_{i}$. In this case,
we need a total of $2M$ threads where two threads $2i$ and $2i+1$ collectively
compute all non-zeros of $\mathbf{A}_{i}$, with each thread computing half of
the non-zeros. Figure 4 illustrates this computation for $\mathbf{A}_{i}$ that
has four non-zeros at column indices $j_{0},j_{1},j_{2}$, and $j_{3}$. Thread
$2i$ computes non-zeros at column indices $j_{0}$ and $j_{1}$; and thread
$2i+1$ computes non-zeros at column indices $j_{2}$ and $j_{3}$. The SDDMM
Algorithm 1 for one thread per row of the output matrix can be easily adjusted
to work with two threads per row.

Figure 3: Illustration of SDDMM implementation for computing non-zeros of
$\mathbf{A}_{i}$. A single thread $i$ computes all non-zeros of
$\mathbf{A}_{i}$. Figure 4: SDDMM implementation with two threads assigned to
compute a of $\mathbf{A}$.

###  3.2 Parallelizing SPMM

The SPMM operation multiplies a sparse matrix $\mathbf{A}$ with a dense matrix
$\mathbf{B}$ to compute a dense matrix $\mathbf{C}$. The output matrix
$\mathbf{C}$ is of size $M\times N$. We parallelize the computation over the
rows of the matrix $\mathbf{C}$. A total of $M$ threads are assigned to work
concurrently, where a thread $i$, $1\leq i\leq M$, computes row
$\mathbf{C}_{i}$. This computation requires multiplying a sparse row
$\mathbf{A}_{i}$ with selected elements of the $\mathbf{B}$ matrix to compute
$\mathbf{C}_{i}$. We assume there are three non-zeros in row $\mathbf{A}_{i}$
at column indices $j_{0}$, $j_{1}$ and $j_{2}$, as shown in Figure 5. It
implies that we need rows $\mathbf{B}_{j_{0}}$, $\mathbf{B}_{j_{1}}$, and
$\mathbf{B}_{j_{2}}$ to compute the row $\mathbf{C}_{i}$. We implement this
computation by multiplying $\mathbf{A}_{i,J_{0}}$ with all elements in the row
$\mathbf{B}_{0}$. The result is added to the multiplication of
$\mathbf{A}_{i,j_{1}}$ with $\mathbf{B}_{j_{1}}$, and so on.

A high-level description of the algorithm for computing a row $\mathbf{C}_{i}$
is outlined in Algorithm 5. As before, we use the CSR format for $\mathbf{A}$
to describe the algorithm. The main loop is over the non-zero elements in
$\mathbf{A}_{i}$. In line 4, we obtain the $j$-th non-zero value in row $i$ of
$\mathbf{A}$. The column index, $k$, associated with the $j$-th non-zero value
is fetched in line 5. The for loop starting at line 6 performs a SAXPY
computation [21], which is simply multiplying a scalar ($s$) with a vector
($\mathbf{B}_{k}$) and adding the result to another vector ($\mathbf{C}_{i}$).

Figure 5: A simple example to illustrate multiplication of the row
$\mathbf{A}_{i}$ with $\mathbf{B}$ to compute the row $\mathbf{C}_{i}$ in the
SPMM algorithm. Algorithm 2 SPMM-ROW($ia,ja,avalues,B,C,N,\texttt{i}$)

1:Initialize: $C[i,:]\leftarrow 0$

2:$nnzr\leftarrow ia[i+1]-ia[i]$

3:for $j\leftarrow 0$ to $nnzr-1$ do

4: $s\leftarrow\texttt{avalues}[ia[i]+j]$

5: $k\leftarrow\texttt{ja}[ia[i]+j]$

6: for $l\leftarrow 0$ to $N-1$ do

7: $C[i,l]\leftarrow C[i,l]+s*B[k,l]$

8: end for

9:end for

10:return $C$

###  3.3 Parallelizing FusedMM

The parallel fused version of the SDDMM and SPMM operation is straightforward.
We consider a single thread per rows version of the SDDMM algorithm, see
Figure 3. Each thread is responsible for computing non-zero values for a row
of sparse matrix $\mathbf{A}$. Once a thread has computed all the non-zero
values of a row (note we skip storing these values in the global array
$\mathbf{A}$), we start the SPMM operation as outlined earlier; see Figure 5.
Figure 6 illustrates the FusedMM operation. A high-level description of the
algorithm for computing a row $\mathbf{E}_{i}$ is outlined in Algorithm 3,
which is essentially derived from Algorithm 1 and Algorithm 2.

Figure 6: A simple example to illustrate FusedMM operation. Algorithm 3
FusedMM-ROW($ia,ja,C,B,D,E,N,\texttt{i}$)

1:Initialize: $C[i,:]\leftarrow 0$

2:$nnzr\leftarrow ia[i+1]-ia[i]$

3:$arow[0:nnzr]\leftarrow 0$

4:for $j\leftarrow 0$ to $nnzr-1$ do

5: $k\leftarrow\texttt{ja}[ia[i]+j]$

6: for $l\leftarrow 0$ to $N-1$ do

7: $arow[j]\leftarrow arow[j]+C[i,l]*B[k,l]$

8: end for

9:end for

10:Initialize: $E[i,:]\leftarrow 0$

11:for $j\leftarrow 0$ to $nnzr-1$ do

12: $s\leftarrow\texttt{arow}[j]$

13: $k\leftarrow\texttt{ja}[ia[i]+j]$

14: for $l\leftarrow 0$ to $N-1$ do

15: $E[i,l]\leftarrow E[i,l]+s*B[k,l]$

16: end for

17:end for

18:return $E$

##  4 ESIMD Implementation

###  4.1 Sparse Matrix Collection

The range of matrix sizes in our collection influences the implementation.
Specifically, our implementation is only guaranteed to work for $N=32$ or
$N=128$, which is typically the batch size used in the machine learning
training dataset. For performance comparison, we use the same matrix dataset
that is used in [4]. It consists of 72 tests of varying sizes (cf. Figs. 12\-
17) and randomized input matrices.

###  4.2 ESIMD Implementation of SDDMM

The ESIMD implementation of the SDDMM operation is based on the algorithm
covered in Section 3.1 (Algorithm 1). In the general version of this
algorithm, one or more threads are assigned to compute a row of the output
sparse matrix $\mathbf{A}$. The SDDMM ESIMD kernel is launched with $M\times
NT$ threads 111In our discussion, we use thread and ESIMD work item
interchangeably. with a work-group of size $NT$. A work-group $i$ consisting
of $NT$ threads computes a row $i$ of the sparse matrix of $nnzr$ non-zeros.
Divide the number of non-zeros to be computed, $nnzr$, for a row equally
amongst $NT$ threads, where a thread works on $nnzt=\lceil nnzr/NT\rceil$ non-
zeros. The number of non-zeros assigned to a thread, $nnzt$, are further
partitioned into chunks, where the number of chunks is given by
$nchunks=\lceil nnzt/VLC\rceil$. Here, $VLC$ is a tuning parameter that varies
in the size range of simd objects supported by ESIMD API.

The SDDMM kernel code segment is shown in Figure 7. In line 3, we assign the
work-group id to the variable $i$, and $j$ is assigned the local thread id
that ranges from $0$ to $NT-1$. In line 5, a work-group $i$ load two
consecutive values $ia[i]$ and $ia[i+1]$ into a vector of size $2$. In line
10, we load row $\mathbf{C}_{i}$ into a simd object $regL$ of size $N$. Next,
we setup a nested loop, where the outer loop is over $nchunks$ and the inside
loop is over $VLC$. In the outer loop we load column indices for the current
chunk into a simd object, $ja\\_block$, of size $VLC$. In the inside loop, we
first load the column of $B$ that is needed for the current iteration (line
19-20) into $regR$. Next, we multiply the two simd objects, $regL$ and $regR$
and store the result in a simd object $res$. In line 22, we do the reduction
operation on the simd object $res$, and store the result in position $l0$ of a
simd object $dp$ of size $VLC$. At the end of the inside loop, we have the
$VLC$ non-zero values that are stored in the appropriate location in
$avalues$. The is repeated over $nchunks$, and at the end of the outer loop we
have computed all the non-zero values in row $\mathbf{A}_{i}$ of the sparse
matrix. To keep our code segment simple, we avoid the details for handling the
cases when $nnzr$ is not a multiple of $NT$, and $nnzt$ is not a multiple of
$VLC$.

To hide memory latency and effectively use the vector engine, we add
prefetches and unroll the inside loop as illustrated in the code segment in
Figure 8. The prefetch is added to avoid stalls in loading the $ja\\_block$.
In the outer loop iteration $l$, we prefetch the memory location from where we
will load $ja\\_block$ in $(l+1)$-th iteration (see line 19-20 of Figure 8).
We add pragma to unroll the inside loop (see line 10).

Additionally, we template the kernel on $NT$ (not shown in the code listing).
We implemented a scheme that chooses the value of $NT$ in dependence of the
matrix size such that the occupancy on the device is maximized. To achieve
100% occupancy on a single-stack of the Intel GPU, one needs to launch 4096
esimd work-items or a multiple thereof. We choose $NT$ as the minimum power of
two, not larger than 16, to maximize the occupancy, which is given as follows

$\text{Occupancy}=\frac{\frac{M\times NT}{4096}}{\lceil\frac{M\times
NT}{4096}\rceil}.$

The ceiled value in the denominator ensures that cases where $M\times NT>4096$
are handled appropriately. For example, let $M=1024$, then $NT=4$ is chosen.
For $M=3072$ the choice of $NT=4$ leads to 100% occupancy.

⬇

#define VLC 32

\parconst int i = item.get_group(0);

const int j = item.get_local_id(0);

simd<int, 2> lrowptr2 = lsc_block_load<int, 2>(ia + i);

const int nnzr = lrowptr2[1] - lrowptr2[0];

const int nnzt = (nnzr + NT - 1) / NT;

const int nchunks = (nnzt + VLC - 1) / VLC;

\parsimd<float, N> regL = my_lsc_block_load<float, N>(C + i * N);

simd<int, VLC> ja_block;

simd<float, VLC> a_row;

for (int l = 0; l < nchunks; l++)

{

int idxb = lrowptr2[0] + j * nnzt + l * VLC;

ja_block = lsc_block_load<int, VLC>(ja + idxb);

for (int l0 = 0; l0 < VLC; l0++ )

{

simd<float, N> regR =

my_lsc_block_load<float, N>(B + ja_block[l0] * N);

simd<float, N> res = regL * regR;

a_row[l0] = reduce<float, float, N>(res, std::plus<>());

}

lsc_block_store<float, VLC>(avalues + idxb, a_row);

}

\par

Figure 7: ESIMD implementation of the SDDMM operation.

⬇

\parsimd<float, N> regL = my_lsc_block_load<float, N>(C + i * N);

simd<int, VLC> ja_block;

simd<float, VLC> a_row;

int idxb = lrowptr2[0] + j * nnzt ;

lsc_prefetch<int, VLC, DSZ, L1_C, L3_C>(&ja[idxb]);

for (int l = 0; l < nchunks; l++)

{

ja_block = lsc_block_load<int, VLC>(ja + idxb);

#pragma unroll

for (int l0 = 0; l0 < VLC; l0++ )

{

simd<float, N> regR =

my_lsc_block_load<float, N>(B + ja_block[l0] * N);

simd<float, N> res = regL * regR;

a_row[l0] = reduce<float, float, N>(res, std::plus<>());

}

lsc_block_store<float, VLC>(avalues + idxb, a_row);

idxb = lrowptr2[0] + j * nnzt + (l+1) * VLC;

lsc_prefetch<int, VLC, DSZ, L1_C, L3_C>(&ja[idxb]); }

Figure 8: SDDMM kernel with prefetching and unrolling of the inside loop.

###  4.3 Implementation of SPMM

The ESIMD implementation of SPMM is based on the Algorithm 2. As outlined
there, we parallelize over the rows of matrix $\mathbf{C}$. The SPMM ESIMD
kernel is launched with $M$ threads with work-group of size one. In Algorithm
2, the main loop is over $nnzr$, the number of non-zero elements in
$\mathbf{A}_{i}$. For our implementation, we partition $nnzr$ into $nchunk$
chunks, where $nchunks=\lceil nnzt/VLC\rceil$. Here, $VLC$ is a tuning
parameter that varies in the size range of simd objects supported by ESIMD
API.

The SPMM kernel code segment is shown in Figure 9. We setup a nested loop,
where the outer loop is over $nchunks$ and the inside loop is over $VLC$. In
the outer loop we load column indices for the current chunk into a simd
object, $ja\\_block$, of size $VLC$ (line 18). In the outer loop, in line 20
we load a chunk of size $VLC$ of non-zero values from the row of sparse matrix
into a simd object $arow$. In the inside loop, we load the $j0$-th element
from $arow$, a scalar value into $s$. In line 25, we load the required row of
$\mathbf{B}$ of size $N$ into a simd object $b\\_row$. Next, we multiply the
scalar value $s$ with $b\\_row$ and accumulate the result in a simd object
$c\\_row$ of size $N$. At the end of the outer loop, $c\\_row$ holds the
values of $\mathbf{C}_{i}$, which is then stored in $\mathbf{C}$ in line 30.
To keep our code segment simple, we avoid the details for handling the case
when $nnzr$ is not a multiple of $VLC$. We add prefetching and unrolling of
the inside loop in the SPMM code similar to SDDMM code, see Figure 10. For
improving occupancy, particularly for small-size matrices, having more threads
in the work-group is desirable. The code in Figure 10 can be easily adjusted
to support $NT$ threads in a work-group, which will compute $NT$ rows of the
output matrix $\mathbf{C}$. Note that the kernel will still be launched with a
total of $M$ threads.

⬇

#define VLC 32

\parconst int i = item.get_group(0);

simd<int, 2> lrowptr2 = lsc_block_load<int, 2>(ia + i);

const int nnzr = lrowptr2[1] - lrowptr2[0];

const int nchunks = (nnzr + VLC - 1) / VLC;

\par\parsimd<int, VLC> ja_block;

simd<float, VLC> a_row;

simd<float, N> c_row;

simd<float, N> b_row;

c_row = 0.0;

\parfor (int l = 0; l < nchunks; l++)

{

int idxb = lrowptr2[0] + l * VLC;

ja_block = lsc_block_load<int, VLC>(ja + idxb);

ja_block = ja_block * N;

a_row = lsc_block_load<float, VLC>(avalues + idxb);

\parfor (int j0 = 0; j0 < VLC; j0++)

{

float s = a_row[j0];

int colid = ja_block[j0];

b_row.copy_from(B + colid);

c_row = c_row + s * b_row;

}

}

c_row.copy_to(C + i * N);

\par\par

Figure 9: SPMM Kernel.

⬇

simd<int, VLC> ja_block;

simd<float, VLC> a_row;

simd<float, N> c_row;

simd<float, N> b_row;

c_row = 0.0;

\parint idxb = lrowptr2[0] ;

lsc_prefetch<int, VLC, DSZ, L1_C, L3_C>(&ja[idxb]);

lsc_prefetch<float, VLC, DSZ, L1_C, L3_C>(&avalues[idxb]);

\parfor (int l = 0; l < nchunks; l++)

{

ja_block = lsc_block_load<int, VLC>(ja + idxb);

ja_block = ja_block * N;

a_row = lsc_block_load<float, VLC>(avalues + idxb);

idxb = lrowptr2[0] + (l + 1) * VLC;

lsc_prefetch<int, 32, DSZ, L1_C, L3_C>(&ja[idxb]);

lsc_prefetch<float, 32, DSZ, L1_C, L3_C>(&avalues[idxb]);

#pragma unroll

for (int j0 = 0; j0 < VLC; j0++)

{

float s = a_row[j0];

int colid = ja_block[j0];

b_row.copy_from(B + colid);

c_row = c_row + s * b_row;

}

}

c_row.copy_to(C + i * N);

\par\par\par

Figure 10: SPMM Kernel with prefetching and unrolling of the inside loop.

###  4.4 Implementation of FusedMM

We can easily fuse the SDDMM with SPMM for a work-group with $NT=1$ thread.
The choice of $NT=1$ ensures the two operations SDDMM and SPMM require an
identical number of total threads. Recall that the SDDMM kernel is launched
with $M\times NT$ threads and that the SPMM kernel is launched with $M$
threads (independent of $NT$). The fused kernel for $NT=1$ is shown in Figure
11. The inputs to the kernel are $\mathbf{C}$, $\mathbf{B}$, $\mathbf{D}$
matrices, along with the two arrays, $ia$ and $ja$, that capture the sparsity
of the $\mathbf{A}$ matrix, see Equation 3. The output is the $\mathbf{E}$
matrix. One can observe that the first part of the FusedMM code 11, until line
28, is a copy of the SDDMM kernel code without storing the $a\\_row$ into
$avalues$ (line 18 of Figure 8). Once the data in $a\\_row$ is ready, we start
the SPMM kernel. More specifically, we took the inside loop of the SPMM
kernel, line 20-26 of Figure 10 and used it in the fused kernel at line 31-37
after taking into account that we are working with different matrices. At the
end of the outer loop in 11 we copy $e\\_row$ to $\mathbf{E}_{i}$ similar to
line 28 of the SPMM kernel, Figure 10.

⬇

const int i = item.get_group(0);

simd<int, 2> lrowptr2 = lsc_block_load<int, 2>(ia + i);

const int nnzr = lrowptr2[1] - lrowptr2[0];

const int nchunks = (nnzr + VLC - 1) / VLC;

\parsimd<float, N> regL = my_lsc_block_load<float, N>(C + i * N);

simd<int, VLC> ja_block;

\parint idxb = lrowptr2[0];

lsc_prefetch<int, 32, DSZ, L1_C, L3_C>(&ja[idxb]);

\parsimd<float, VLC> a_row;

simd<float, N> e_row;

simd<float, N> d_row;

e_row = 0.0;

\parfor (int l = 0; l < nchunks - 1; l++)

{

ja_block = lsc_block_load<int, VLC>(ja + idxb);

ja_block = ja_block * N;

\par#pragma unroll

for (int l0 = 0; l0 < VLC; l0++)

{

simd<float, N> regR = my_lsc_block_load<float, N>(B + ja_block[lo]);

simd<float, N> res = regL * regR;

a_row[l0] = reduce<float, float, N>(res, std::plus<>());

}

\par#pragma unroll

for (int j0 = 0; j0 < VLC; j0++)

{

float s = a_row[j0];

int colid = ja_block[j0];

d_row.copy_from(D + colid);

e_row = e_row + s * d_row;

}

idxb = lrowptr2[0] + (l + 1) * VLC;

lsc_prefetch<int, 32, DSZ, L1_C, L3_C>(&ja[idxb]);

}

e_row.copy_to(E + i*N);

\par

Figure 11: FusedMM Kernel with prefetching and unrolling of the inside loop.

##  5 Experiments

In this section we compare our implementations to the implementations
available in MKL. All tests were performed on a single stack of an Intel Data
Center GPU Max 1550 on an Intel-internal test system. The code was compiled
with Intel’s icpx compiler (2023.2.0.20230721), which is included in Intel’s
oneAPI, version 2023.2.1. We used the compile options "-fsycl", "-O3" and, in
the case of the oneMKL tests, "-qmkl". Further, we used an unreleased
engineering GPU driver for the tests 222agama-ci-devel/682.16.

Figures 12 and 13 show the performance of our SDDMM implementation compared to
oneMKL’s (included in oneAPI version 2023.2.1 [22]) dense blas::gemm
implementation [23]. This comparison is highly disadvantageous for oneMKL
since 70%-90% (depending on the sparsity of the workload) of the computed
flops in a dense matrix-matrix multiplication are irrelevant for the output of
SDDMM. Thus, although oneMKL blas::gemm achieves a peak performance of
approximately 25,000 gigaflops per second (Gflops/s) in the dense matrix
multiplication, our specialized SDDMM implementation can outperform it. The
reason we compare to this oneMKL function is that there is not yet a dedicated
SDDMM functionality included in oneMKL. Overall, our implementation shows an
average speedup of 1.84x compared to oneMKL’s dense gemm and a maximum speedup
of 3.4x.

Compared to the implementation shown in [4], we increase the peak performance
from approximately 2,700 Gflops/s (achieved on the 4k/1k/128/70% case) to
close to 10,700 Gflops/s (achieved on the 12k/4k/128/70% case). Further, the
average performance gains achieved compared to [4] in the SDDMM case is a
factor of approximately 4.5x with a maximum relative performance increase of
approximately 10x (achieved on the 32k/8k/32/90% case).

1k,1k,32,70%1k,1k,32,80%1k,1k,32,90%3k,1k,32,70%3k,1k,32,80%3k,1k,32,90%4k,1k,32,70%4k,1k,32,80%4k,1k,32,90%2k,2k,32,70%2k,2k,32,80%2k,2k,32,90%6k,2k,32,70%6k,2k,32,80%6k,2k,32,90%8k,2k,32,70%8k,2k,32,80%8k,2k,32,90%4k,4k,32,70%4k,4k,32,80%4k,4k,32,90%12k,4k,32,70%12k,4k,32,80%12k,4k,32,90%16k,4k,32,70%16k,4k,32,80%16k,4k,32,90%8k,8k,32,70%8k,8k,32,80%8k,8k,32,90%24k,8k,32,70%24k,8k,32,80%24k,8k,32,90%32k,8k,32,70%32k,8k,32,80%32k,8k,32,90%$0$$2{,}000$$4{,}000$Gflops/sSDDMM
Benchmark for $N=32$Our SDDMMMKL GEMM Figure 12: Comparison of our SDDMM
method with MKL blas::gemm. The x-label indicates M/K/N/Sparsity.
1k,1k,128,70%1k,1k,128,80%1k,1k,128,90%3k,1k,128,70%3k,1k,128,80%3k,1k,128,90%4k,1k,128,70%4k,1k,128,80%4k,1k,128,90%2k,2k,128,70%2k,2k,128,80%2k,2k,128,90%6k,2k,128,70%6k,2k,128,80%6k,2k,128,90%8k,2k,128,70%8k,2k,128,80%8k,2k,128,90%4k,4k,128,70%4k,4k,128,80%4k,4k,128,90%12k,4k,128,70%12k,4k,128,80%12k,4k,128,90%16k,4k,128,70%16k,4k,128,80%16k,4k,128,90%8k,8k,128,70%8k,8k,128,80%8k,8k,128,90%24k,8k,128,70%24k,8k,128,80%24k,8k,128,90%32k,8k,128,70%32k,8k,128,80%32k,8k,128,90%$0$$5{,}000$$10{,}000$Gflops/sSDDMM
Benchmark for $N=128$Our SDDMMMKL GEMM Figure 13: Comparison of our SDDMM
method with MKL blas::gemm. The x-label indicates M/K/N/Sparsity.

Figures 14 and 15 show the performance of our SPMM implementation in
comparison to oneMKL’s sparse::gemm [24] implementation. Our implementation
increases the performance on average by 1.37x and at most by 2.16x. These
performance increases may be attributed to the high specialization of our code
to the specific sizes and sparsities relevant in machine learning
applications. The average relative performance gain compared to [4] is
approximately 2.4x with a maxium performance increase of approximately 4x
(achieved on the 32k/8k/128/90% case).

1k,1k,32,70%1k,1k,32,80%1k,1k,32,90%3k,1k,32,70%3k,1k,32,80%3k,1k,32,90%4k,1k,32,70%4k,1k,32,80%4k,1k,32,90%2k,2k,32,70%2k,2k,32,80%2k,2k,32,90%6k,2k,32,70%6k,2k,32,80%6k,2k,32,90%8k,2k,32,70%8k,2k,32,80%8k,2k,32,90%4k,4k,32,70%4k,4k,32,80%4k,4k,32,90%12k,4k,32,70%12k,4k,32,80%12k,4k,32,90%16k,4k,32,70%16k,4k,32,80%16k,4k,32,90%8k,8k,32,70%8k,8k,32,80%8k,8k,32,90%24k,8k,32,70%24k,8k,32,80%24k,8k,32,90%32k,8k,32,70%32k,8k,32,80%32k,8k,32,90%$0$$2{,}000$$4{,}000$$6{,}000$Gflops/sSpMM
Benchmark for $N=32$Our SPMMMKL SPMM Figure 14: Comparison of our SPMM method
with MKL sparse::gemm. The x-label indicates M/K/N/Sparsity.
1k,1k,128,70%1k,1k,128,80%1k,1k,128,90%3k,1k,128,70%3k,1k,128,80%3k,1k,128,90%4k,1k,128,70%4k,1k,128,80%4k,1k,128,90%2k,2k,128,70%2k,2k,128,80%2k,2k,128,90%6k,2k,128,70%6k,2k,128,80%6k,2k,128,90%8k,2k,128,70%8k,2k,128,80%8k,2k,128,90%4k,4k,128,70%4k,4k,128,80%4k,4k,128,90%12k,4k,128,70%12k,4k,128,80%12k,4k,128,90%16k,4k,128,70%16k,4k,128,80%16k,4k,128,90%8k,8k,128,70%8k,8k,128,80%8k,8k,128,90%24k,8k,128,70%24k,8k,128,80%24k,8k,128,90%32k,8k,128,70%32k,8k,128,80%32k,8k,128,90%$0$$5{,}000$$10{,}000$Gflops/sSpMM
Benchmark for $N=128$Our SPMMMKL SPMM Figure 15: Comparison of our SPMM method
with MKL sparse::gemm. The x-label indicates M/K/N/Sparsity.

Figures 16 and 17 compares the performance of our FusedMM operation to
oneMKL’s dense::gemm and sparse::gemm operation. For the oneMKL case, we are
not including the time required for the sampling after the dense matrix
multiplication for the SDDMM operation. Our implementation increases the
performance of MKL by approximately 2x on average and 3.25x at most. The fused
operation achieves on average 1.3x the throughput of our SDDMM implementation
and 1.074x of our SPMM implementation. The average performance gain is thus
higher than in either of the cases before, underlining the importance of
operator fusion for these kind of operations. The peak performance of 12,647
Gflops/s for this operations represents 48% of the theoretical peak
performance of the device.

1k,1k,32,70%1k,1k,32,80%1k,1k,32,90%3k,1k,32,70%3k,1k,32,80%3k,1k,32,90%4k,1k,32,70%4k,1k,32,80%4k,1k,32,90%2k,2k,32,70%2k,2k,32,80%2k,2k,32,90%6k,2k,32,70%6k,2k,32,80%6k,2k,32,90%8k,2k,32,70%8k,2k,32,80%8k,2k,32,90%4k,4k,32,70%4k,4k,32,80%4k,4k,32,90%12k,4k,32,70%12k,4k,32,80%12k,4k,32,90%16k,4k,32,70%16k,4k,32,80%16k,4k,32,90%8k,8k,32,70%8k,8k,32,80%8k,8k,32,90%24k,8k,32,70%24k,8k,32,80%24k,8k,32,90%32k,8k,32,70%32k,8k,32,80%32k,8k,32,90%$0$$2{,}000$$4{,}000$$6{,}000$$8{,}000$Gflops/sFusedMM
Benchmark for $N=32$Our FusedMKL SDDMM+SPMM Figure 16: Comparison of our
FusedMM method with MKL. The x-label indicates M/K/N/Sparsity.
1k,1k,128,70%1k,1k,128,80%1k,1k,128,90%3k,1k,128,70%3k,1k,128,80%3k,1k,128,90%4k,1k,128,70%4k,1k,128,80%4k,1k,128,90%2k,2k,128,70%2k,2k,128,80%2k,2k,128,90%6k,2k,128,70%6k,2k,128,80%6k,2k,128,90%8k,2k,128,70%8k,2k,128,80%8k,2k,128,90%4k,4k,128,70%4k,4k,128,80%4k,4k,128,90%12k,4k,128,70%12k,4k,128,80%12k,4k,128,90%16k,4k,128,70%16k,4k,128,80%16k,4k,128,90%8k,8k,128,70%8k,8k,128,80%8k,8k,128,90%24k,8k,128,70%24k,8k,128,80%24k,8k,128,90%32k,8k,128,70%32k,8k,128,80%32k,8k,128,90%$0$$5{,}000$$10{,}000$Gflops/sFusedMM
Benchmark for $N=128$Our FusedMKL SDDMM+SpMM Figure 17: Comparison of our
FusedMM method with MKL. The x-label indicates M/K/N/Sparsity.

###  5.1 Performance Analysis

In this section, we investigate the efficiency of our implementation using the
roofline model [25]. In particular, we examine the utilization of caches and
how close we are to the theoretical peak of the underlying hardware. We
selected a high-performant case with $M=8192$, $K=8192$, $N=128$, and
$\alpha=0.7$.

For this problem size, we observed a performance of $9,732$ Gflops/s for the
SDDMM operation (cf. Figure 13). The main computation in the SDDMM
implementation is on line 15-16 of Figure 8. Line 15 is a vector
multiplication followed by a reduction. Hence, the performance of SDDMM is
bounded by the "SP Vector Add Peak" of $12,600$ Gflops/s as shown in the
roofline model in Figure 18. The SDDMM performance of $9,732$ Gflops/s is 77%
of the theoretical peak. As indicated earlier, we ran a kernel in a loop over
20 iterations for all our experiments and selected the minimum execution time
for plotting. Since the matrices are sufficiently small, they are cached in
$L3$ after the first iteration. The minimum execution time thus corresponds to
the iteration when most of the data is in L3. This is confirmed with the
roofline model plotted by the Intel Advisor tool shown Figure 18. The red dot
on the right in Figure 18 represents the HBM traffic and it shows an amount of
traffic to HBM which corresponds to only the first iteration. Note that the
numbers shown in Figures 12 and 13 are consistent since we used the same
methodology to generate the MKL data: multiple consecutive runs in a loop and
selected the minimum execution time. In separate experiments, which are not
shown here, we observed a slowdown of close to $10\%$ for the first iteration
compared to the rest of the iterations due to the necessary memory access to
HBM instead of L3 cache. The aggregated performance over 20 iterations shown
by the Advisor tool is $7800$ Gflops/s, which accounts for 10% slow down for
the first iteration and the overhead of the profiler. For the L3 we observed
an arithmetic intensity of 5.6, which is significantly better compared to the
theoretical worst case of the arithmetic intensity of 0.49. For the
calculation of the worst case arithmetic intensity, we assume we load the
columns of the $\mathbf{B}^{T}$ matrix always from L3; that is, there is no
caching of the column in L1. In other words, the total number of bytes loaded
for $\mathbf{B}^{T}$ from L3 is $N\times nnz\times 4$ bytes. On the other
hand, the theoretical peak arithmetic intensity, which would be possible under
the assumption of infinite caches, is close to $57$. This number is calculated
based on the assumption that the required data for the SDDMM operation is only
loaded a single from the L3.

Figure 18: Roofline model illustrating performance of the SDDMM kernel for
$M=8192$, $K=8192$, $N=128$, and $\alpha=0.7$

For our SPMM implementation, we observed a performance of $11,924$ Gflops/s
for the high-performant case with $M=8192$, $K=8192$, $N=128$, and
$\alpha=0.7$. The memory access pattern and operation count for SPMM are
similar to the SDDMM operation; however, we observed a better absolute
performance for SPMM ($10,884$ Gflops/s) compared to SDDMM ($9,732$ Gflops/s).
Thhe main reason for the increased performance compared to SDDMM is that the
main computation in the SPMM implementation on line 25 is a multiply-add (MAD)
operation; see Figure 10. While the absolute performance of SPMM is better, it
achieves a lower relative performance compared to the theoretical peak
performance of the hardware. This is due to the performance of SPMM being
bounded by the "SP Vector MAD Peak" of $26,000$ Gflops/s as shown in the
roofline model in Figure 19 (We note that 2x the single precision vector add
peak would actually amount to $25,200$ Gflops/s. We use $26,000$ Gflops/s to
stay consistent with the Advisor tool). The SPMM performance is thus 37% of
the theoretical peak. Based on this data, we believe there is still room for
improvement of the performance of this operation. As in the case of SDDMM, the
matrices are cached in L3 after the first iteration. The aggregated
performance over 20 iterations shown by the advisor tool is $10,342$ Gflops/s,
which accounts for 10% slow down for the first iteration and the overhead of
the profiler. For the L3 cache we observed an arithmetic intensity of 8.06,
which is better than in the SDDMM case and within the expectation of 0.49
(worst case) to $57$ (best case).

Figure 19: Roofline model illustrating performance of the SPMM kernel for
$M=8192$, $K=8192$, $N=128$, and $\alpha=0.7$

##  6 Conclusion and Future Work

In this contribution, we showed that highly optimized implementations of
sparse matrix operations occurring in machine learning applications outperform
existing approaches. In particular, we showed that the difference between the
general SPMM implementation available in Intel’s oneMKL can be outperformed by
our implementation, which is highly specialized for relatively dense sparse
matrices by up to 2.16x. Further, we showed a gap in Intel’s oneMKL
functionalities by not providing a dedicated SDDMM implementation. While the
dense matrix multiply provided by MKL performs well on Intel’s Max GPU, it
suffers from the sparsity inherent in the SDDMM operation. When fusing SPMM
and SDDMM operations, the performance gains increase further to up to a factor
of 3.25. It is thus highly recommended to implement fused operations wherever
possible in the context of machine learning applications. We further showed
that our sparse matrix operations can be run up to 10x faster on an Intel Data
Center GPU compared to an optimized CUDA code on a NVIDIA V100 GPU. For our
future work, we plan to explore sparse matrix operations utilizing bfloat16
data types and evaluate performance on a larger sparse matrix data set for
different machine learning applications.

## Acknowledgements

This effort has been supported by the Intel oneAPI Center of Excellence at Old
Dominion University. We want to thank Xiao Zhu of Intel, who provided support
throughout this project, making resources available whenever we needed them.

## References

  * [1] PyTorch Project a Series of LF Projects, LLC.  PyTorch on GPUs, 2023.  last accessed 10/27/23.
  * [2] Intel Corporation.  Tensorflow on GPUs, 2023.  last accessed 10/27/23.
  * [3] Trevor Gale, Erich Elsen, and Sara Hooker.  The state of sparsity in deep neural networks.  CoRR, abs/1902.09574, 2019.
  * [4] Trevor Gale, Matei Zaharia, Cliff Young, and Erich Elsen.  Sparse gpu kernels for deep learning.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’20. IEEE Press, 2020.
  * [5] Trevor Gale, Deepak Narayanan, Cliff Young, and Matei Zaharia.  Megablocks: Efficient sparse training with mixture-of-experts, 2022.
  * [6] Changwan Hong, Aravind Sukumaran-Rajam, Israt Nisa, Kunal Singh, and P. Sadayappan.  Adaptive sparse tiling for sparse matrix multiplication.  In Proceedings of the 24th Symposium on Principles and Practice of Parallel Programming, PPoPP ’19, page 300–314, New York, NY, USA, 2019. Association for Computing Machinery.
  * [7] Md. Khaledur Rahman, Majedul Haque Sujon, and Ariful Azad.  Fusedmm: A unified sddmm-spmm kernel for graph embedding and graph neural networks.  In 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 256–266, 2021.
  * [8] V. Bharadwaj, A. Buluc, and J. Demmel.  Distributed-memory sparse kernels for machine learning.  In 2022 IEEE International Parallel and Distributed Processing Symposium (IPDPS), pages 47–58, Los Alamitos, CA, USA, jun 2022. IEEE Computer Society.
  * [9] Y. Saad.  Iterative Methods for Sparse Linear Systems.  Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2nd edition, 2003.
  * [10] I. S. Duff, A. M. Erisman, and J. K. Reid.  Direct Methods for Sparse Matrices.  OxfordPress, 1986.
  * [11] J. A. George and J. W. Liu.  Computer Solution of Large Sparse Positive Definite Systems.  PrenticeHall, Englewood Cliffs, NJ, USA, 1981.
  * [12] Yuwei Hu, Zihao Ye, Minjie Wang, Jiali Yu, Da Zheng, Mu Li, Zheng Zhang, Zhiru Zhang, and Yida Wang.  Featgraph: A flexible and efficient backend for graph neural network systems.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’20. IEEE Press, 2020.
  * [13] NVIDIA Corporation.  CUDA C Programming Guide.  last accessed 7/14/23.
  * [14] Khronos Group.  SYCL, 2020.  last accessed 8/24/20.
  * [15] Intel Corporation.  DPC++ Explicit SIMD API, 2023.  last accessed 5/25/23.
  * [16] Mohammad Zubair, Aaron Walden, Gabriel Nastac, Eric Nielsen, Christoph Bauinger, and Xiao Zhu.  Optimization of Ported CFD Kernels on Intel Data Center GPU Max 1550 using oneAPI ESIMD [manuscript accepted for publication].  In ScalAH23: 14th Workshop on Latest Advances in Scalable Algorithms for Large-Scale Heterogeneous Systems, SC’23, 2023.
  * [17] Intel Corporation.  Intel Data Center GPU Max 1550, 2023.  last accessed 7/14/23.
  * [18] Intel Corporation.  Intel Xe GPU Architecture, 2023.  last accessed 7/17/23.
  * [19] Intel Corporation.  Explicit SIMD SYCL Extension, 2023.  last accessed 5/25/23.
  * [20] Intel Corporation.  Intel Thread Mapping.  last accessed 10/30/23.
  * [21] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen.  LAPACK Users’ Guide.  Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999.
  * [22] Intel Corporation.  Intel oneAPI.  last accessed 10/27/23.
  * [23] Intel Corporation.  oneMKL blas::gemm, 2023.  last accessed 10/27/23.
  * [24] Intel Corporation.  oneMKL blas::gemm, 2023.  last accessed 10/27/23.
  * [25] Intel Corporation.  Intel Advisor GPU roofline.  last accessed 10/30/23.


# Model Predictive Control approach to adaptive messaging intervention for
physical activity
††thanks: This work has been partially supported by National Institutes of
Health (NIH) Grant R01 HL142732 and National Science Foundation (NSF) Grant
1808266.

Ibrahim E. Bardakci  School of Electrical Engineering & Computer Science
The Pennsylvania State University
University Park, USA
<EMAIL_ADDRESS>Sahar Hojjatinia  School of Electrical Engineering &
Computer Science
The Pennsylvania State University
University Park, USA
<EMAIL_ADDRESS>Sarah Hojjatinia  AS & UX
Aptiv
Troy, USA
<EMAIL_ADDRESS>Constantino M. Lagoa  School of Electrical
Engineering & Computer Science
The Pennsylvania State University
University Park, USA
<EMAIL_ADDRESS>David E. Conroy  Department of Kinesiology
The Pennsylvania State University
University Park, USA
<EMAIL_ADDRESS>
###### Abstract

In this work, we have developed a framework for synthesizing data driven
controllers for a class of uncertain switched systems arising in an
application to physical activity interventions. In particular, we present an
application of probabilistic model predictive control to design an efficient,
tractable, and adaptive intervention using behavioral data sets i.e. physical
activity behavior. The models of physical activity for each individual are
provided for the design of controllers that maximize the probability of
achieving a desired physical activity goal subject to intervention
specifications. We have tailored the mixed-integer programming-based approach
for evaluating the Model Predictive Control decision at each time step.

###### Index Terms:

just-in-time adaptive intervention, mixed integer programming, control
systems, behavior modification

##  I Introduction

Physical inactivity is recognized as one of the major risk factors for many
chronic diseases, including diabetes, cardiovascular diseases, stroke, and
different types of cancer. It has also been identified as the fourth leading
risk factor for global mortality ($6\%$ of deaths globally) [2]. While
physical activity is widely recommended and is known to have multi-system
benefits for health, most American adults fail to achieve recommended levels
of physical activity [3]. Given the economic burden of physical inactivity on
one hand and health issues caused by it on the other hand, intervention
programs to increase physical activity and reduce the sedentary behavior are
needed.

With the prevalence of smart phones and activity monitors, an effective method
to increase physical activity has shown to be message interventions [4]. The
fact that physical activity is a dynamic behavior and can be influenced by
external factors such as weather conditions should be considered in
interpreting the outcomes of the intervention [5]. Personalized dynamical
models of physical activity based on message interventions developed by our
group have shown a positive effect on increasing physical activity level [6,
7]. However, the change in physical activity differs for each individual as a
function of the type of intervention message, temperature, precipitation, and
day of week. These heterogeneous responses of people to different types of
interventions a need for designing and implementing more efficient
interventions, adaptive to each individual’s behavior and current context that
balance the potential for achieving behavioral goals against the threat of
disengagement from excessive intervention burden, as well as reducing the
burden treatment.

In this regard, the objective of this paper is to develop a computationally
tractable framework for synthesizing data-driven controllers for a class of
uncertain switched systems arising in an application to increase physical
activity level in individuals. We present a design of personalized behavior
monitoring systems using control systems tools, aimed at maintaining desired
levels of physical activity. In this regard, we use personalized models of
physical activity described by a discrete-time dynamical system of the form

$y_{k}\triangleq f(y_{k-1},y_{k-2},\ldots,y_{k-n},u_{k},\ldots,u_{k-n},w_{k})$
(1)

for $k=1,\ldots,T$ where $y_{k}\in\mathcal{Y}$ is the system output,
$u_{k}\in\mathcal{U}$ is the control input and $w_{k}\in\mathcal{W}$ is the
process noise or disturbance, all at time step $k$. The function
$f:\mathcal{Y}\times\mathcal{U}\times\mathcal{W}\rightarrow\mathcal{Y}$ are
the input-output transition functions. The process noise trajectory
$(w_{1},w_{2},\ldots,w_{T})\in\mathcal{W}$ is random, with a known
distribution. In equation (1), the state $y_{k}$ at time $k$ is a function of
$n$ previous values of the state (order of the model), inputs
$u^{j}_{k}\in\mathcal{U}$, and uncertainties $w_{k}\in\mathcal{W}$ at time k.
These dynamical systems that can model the physical activity behavior of an
individual are provided as a basis for developing controllers to maintain
optimal levels of physical activity volume.

In this paper, we use the developed model representing the dynamical behavior
of a subject’s physical activity as a Picewise Affine (PWA) model. The
dynamical models of physical activity used in this paper were developed based
on the data recorded from the participants recruited for a 6-month messaging
intervention study with the aim to increase their physical activity [6, 7]. A
PWA model is obtained by partitioning the state space into a finite number of
regions and associating each of these regions with an affine sub-model. These
sub-models describe possible continuous dynamics, and are coupled via a
discrete logic that determines which subsystem is active at each time instant.
PWA models are tractable approximations of many non-linear phenomena [8] and
they can be used to approximate the behavior of any complex dynamical system
[9]. The PWA model to be used for representing the physical activity of a
subject has sub-models of the form

$y_{k}=a^{\sigma(k)}_{0}+\sum_{i=1}^{n}a^{\sigma(k)}_{i}y_{k-i}+\sum_{j=1}^{m}\sum_{i=0}^{n}b^{\sigma(k)}_{ij}u^{j}_{k-i}+w_{k},$
(2)

where the index $\sigma(k)$ denotes the sub-model active at time $k$,
$u^{j}_{k}$ indicates what type of message is sent at time k (if any), and
$w_{k}$ is noise. The active sub-model at each time $k$ is determined based on
the current and past states of the subject and external factors such as
weather conditions, and day of week. The larger the number of sub-models and
order of them, the more complex the PWA model becomes. The sub-models provide
information on controller design to deal with non-linear moderators of
intervention effects since they represent different intervention responses
under different conditions.

Drawing on these personalized models of physical activity, the objective is to
design controllers that maximize the probability of achieving a desired
physical activity goal subject to intervention specifications. The effective
operating regions for these controllers are specified by the intervention
window. For the physical activity controller, the intervention window is
defined as the times when participants are likely to respond to intervention,
and it is obtained through the feedback in the modeling phase. The source of
feedback typically comes from one of two sources: self-monitoring or sensors.
Sensors embedded in wearable and portable devices have emerged as an
attractive alternative to self-monitoring because they provide more intensive
data, participant burden, and are not vulnerable to the same memory biases as
self-monitoring [10].

###  I-A Probabilistic model predictive control

Once the identified model is given, a controller repeatedly selects when and
what type of intervention to apply in a given window. Given the estimated
distribution of the noise $w_{k}$ and information on the distributions of the
external factors, at each time $k$ the controller determines a sequence of
interventions that maximizes the probability of achieving the goal while
avoiding undesirable situations such as overburden from applying too many
interventions. Hence, at every time instant, given a look ahead horizon $T$,
the controller solves an optimization problem of the form for
$k=1,2,\ldots,T$,

maximize $\displaystyle\mbox{Prob}\left\\{\mbox{Desired Change in Physical
Activity}\right\\}$ (3) subject to $\displaystyle\mbox{Model }\&\mbox{
Intervention Constraints}$

However, this decision process of repeated interventions is highly complex.
Therefore, in this study, we use a Model Predictive Control (MPC) approach as
a way to design the intervention algorithm since model predictive control can
efficiently handle such complexities [11, 12, 13] and has successfully
implemented in intervention designs [14, 15, 16]. In the MPC problem, the
first step is predicting future behavior of the system to be controlled using
mathematical model of the system along with uncertainty descriptions. The
second step consists of determining the sequence of control actions that is
optimal in terms of specified performance measure (e.g., maximizing the
desired objective or minimizing the cost). Finally, the first value of these
control actions applies, then this whole process repeats with the new
measurements. In the standard MPC, the horizon slides forward and for this
reason MPC is also called receding horizon control. However, because of the
nature of this intervention design, the intervention window is predefined. In
other words, the horizon of MPC controller is reduced by one step at each time
instant and is referred as a shrinking horizon MPC [17].

The approach proposed here is aimed at providing intervention design
procedures that can handle large amounts of stochastic uncertainty. Each step
of MPC algorithm consists of probability maximization problem or chance-
constrained optimization involving integer variables.

In our approach, we consider a sample-average approximation of the probability
maximization problem (3) where the probability function in the objective
approximated by a sample average estimates using Monte Carlo simulation. The
resulting problem involving a discrete distribution can be difficult to solve
[18]. Moreover, the solutions obtained by these methods can be quite
conservative. However, we consider a finite (or shrinking) horizon MPC,
therefore we only have a finite number of randomly chosen instances of
possible outcomes of parameter values. Hence, the integer-programming approach
can be good candidate as a solution of the MPC problem. Following the
procedure proposed in [18, 19] we consider integer-programming based
approaches for evaluating such an MPC control decision at each time step.

##  II Objectives and Constraints of the Intervention Design

In this problem, given the input-output description of the model, we aim at
finding an input policy $(u^{j}_{1},u^{j}_{2},\ldots,u^{j}_{T})$, i.e. when
and which message to send, that maximizes the probability of achieving the
desired number of step count over $T-n$ periods subject to model and input
constraints.

$\displaystyle\\!\max_{u_{k^{*}}^{j},...,u_{T}^{j}}$
$\displaystyle\mathbb{P}\left\\{\sum_{k=1}^{k^{*}-1}y_{k}+\sum_{k^{*}}^{T}y_{k}\geq
Goal\right\\}$ (4a) subject to $\displaystyle
y_{k}=f(y_{k-1},y_{k-2},...,y_{k-n},u^{j}_{k},w_{k})$ (4b) $\displaystyle
u^{j}_{k}\in\mathcal{U}$ (4c)

where $\mathcal{U}$ is the input constraint and the desired step count $Goal$
is determined by the historical data collected. At each time period $k$, the
algorithm determines the control action $u^{j}_{k}$. When the control action
is taken at time $k^{*}$, the previous measurements
$(y_{1},\ldots,y_{k^{*}-1})$ are known; hence, the uncertainties until time
$k^{*}$, $(w_{1},\ldots,w_{k^{*}-1})$, are zero. However, the current period
and future outputs $(y_{k}^{*},\ldots,y_{T})$ are still to be predicted and
the uncertainties $(w_{k}^{*},\ldots,w_{T})$ are present in the model.

###  II-A Constraints on interventions

In this section, we quantify the constraints on the text message interventions
(or control inputs). These constraints are aimed at addressing the
intervention or treatment burden [20, 21]. In the intervention design phase,
we elaborate algorithms that constrain the frequency and number of
interventions. Using the feedback of participants in the study, these
constraints on the text message interventions are specified mathematically as
the input constraint set $\mathcal{U}$

$\displaystyle\mathcal{U}\triangleq\Bigl{\\{}$ $\displaystyle
u^{j}_{k}\in\\{0,1\\}:\sum_{k}^{k+1}\sum_{j=1}^{m}u^{j}_{k}\leq 1,\
\sum_{k=1}^{T}\sum_{j=1}^{3}u^{j}_{k}\leq\alpha,$ (5)
$\displaystyle\sum_{k=1}^{T}\sum_{j=1}^{3}c_{k}u^{j}_{k}\leq\beta,\
j=1,...,m,\ k=1,\ldots,T\Bigr{\\}}$

where input $u^{j}_{k}$ can take values $0$ or $1$ representing either
applying or not applying intervention at each time. Sum of two consecutive
inputs is considered to be less than 1 preventing more than one intervention
(of any types) within 30 minutes (two sampling time of 15 minutes). The number
of interventions applied in the given period $T$ (intervention window or
horizon) is constrained to be less than or equal to $\alpha$, that should be
predefined. Total cost of interventions is less than $\beta$. Here, $c_{k}$
represents the cost of applying intervention $u^{j}_{k}$ at time $k$ for
$j=1,\ldots,3$.

###  II-B Determining the costs on interventions

In this model, we use a progressive penalty function on the text message
interventions (inputs). It is a composite based on the time of day (uniform
and progressively lower cost for messages sent later in the day) and person-
specific temporal patterning of physical activity (greater cost for messages
at times when a person is typically more active). The cost function associated
with the time of day is introduced to prevent pile-up of messages early in the
day (depleting the bank of messages that might be sent later in the day). The
cost function associated with temporal patterning of physical activity
introduced as a constraint based on national data showing hourly patterns in
physical activity.

More specifically we have sequence of cost functions $c^{time}_{k}$ and
$c^{step}_{k}$ on input sequence $u^{t}_{k}$ for $k=n+1,\ldots,T$. The cost
function $c^{step}_{k}$ is obtained from the historical data which is
determined as the normalized hourly averages of step count of the participant
for the specified intervention window (see Figure 1-a). The function
$c^{time}_{k}$ is chosen as monotonically decreasing function of a time to
reflect the cost of applying intervention later in the day is lower (see
Figure 1-b). Then, the composite cost function is obtained as
$c_{k}=c^{time}_{k}\times c^{step}_{k}$ for $k=n+1,\ldots,T$.

(a)

(b)

(c)

Figure 1: (a) Costs as a function of temporal patterning of physical activity
($c^{step}_{k}$) (b) Costs as a function of time of the day ($c^{time}_{k}$).
(c) Composite costs $c_{k}$.

##  III Stochastic Model Predictive Control for Physical Activity Intervention
Design

We now describe a general optimization problem and show how model predictive
control can be implemented on the problem. At every time instant, given the
last $k$ measurements, the controller repeatedly solves the optimization
problem (4). Following the approach in [18], the problem (4) can be formulated
as the following stochastic mixed-integer program over the time interval
$k=n+1,n+2,\ldots,T$

$\displaystyle\\!\max$ $\displaystyle\frac{1}{N}\sum_{s=1}^{N}p_{s}$ (SMIP-
MPC) subject to
$\displaystyle\left(\sum_{k=1}^{k^{*}-1}y_{k}+\sum_{k=k^{*}}^{T}y^{s}_{k}-Goal\right)\geq
M^{s}(p^{s}-1)\qquad\qquad$ $\displaystyle
y^{s}_{k}=a_{0}+\sum_{i=1}^{n}a_{i}y^{s}_{k-i}+\sum_{j=1}^{m}\sum_{i=0}^{n}b_{ij}u^{j}_{k-i}+w^{s}_{k}$
$\displaystyle p_{s}\in\\{0,1\\}$ $\displaystyle u^{j}_{k}\in\mathcal{U}$

for $s=1,\ldots,N$ where $M^{s}$ is a large positive number such that
$M^{s}\geq\max_{u^{j}_{k}\in\mathcal{U}}\left(\sum_{k=1}^{k^{*}-1}y_{k}+\sum_{k=k^{*}}^{T}y^{s}_{k}-Goal\right)$
.

Algorithm 1 Shrinking Horizon MPC algorithm

(0) Given the step count measurements $y_{1},\ldots,y_{k^{*}-1}$, past input
sequence $u^{j}_{1},u^{j}_{2},\ldots,u^{j}_{k^{*}-1}$, dynamical model $f$,
desired step count Goal, the distribution of $w_{k}$
(1) Solve (SMIP-MPC) and obtain a future control sequence
$u_{k^{*}}^{j},u_{k^{*}+1}^{j},\ldots,u_{T}^{j}$
(2) Apply $u_{k^{*}}^{j}$ as the current input.
(3) Set $k^{*}\doteq k^{*}+1$, and return to (1), till $k^{*}=T$

More specifically, at time $k^{*}$, we first collect the measurements
$(y_{1},y_{2},\ldots,y_{k^{*}-1})$, and past input sequence
$(u^{j}_{1},u^{j}_{2},\ldots,u^{j}_{k^{*}-1})$, and generate the current and
future disturbance sequences $(w_{k}^{*},w_{k^{*}+1},\ldots,w_{T})$. Then, we
solve the stochastic mixed-integer MPC problem over the remaining period,
i.e., $k=k^{*},\ldots,T$, using the information given in the first step, and
obtain the optimal current and future input sequence
$(u_{k^{*}}^{j},u_{k^{*}+1}^{j},\ldots,u^{j}_{T})$. Finally, we use
$u_{k}^{j^{*}}$ as the current input.

##  IV Application to behavioral data sets

In this section, we provide an implementation of the data driven approach to
controller design for developing personalized treatments aimed at improving
physical activity volume as measured by step counts. Dynamic model of physical
activity is provided as a 5-th order switched affine system with inputs being
binary variables indicating if a message of any type was received or not, and
output the current step counts. The model consists of two sub-systems based on
weekdays and weekends. The process noise is normally distributed and $\mu_{w}$
and $\sigma_{k}$ are known for each time step $k$. The PWA model of physical
activity of an individual for weekdays is given as

$\displaystyle y_{k}=$ $\displaystyle
80.51-0.0052y_{k-1}+0.0043y_{k-2}+0.0421y_{k-3}$ (7) $\displaystyle-$
$\displaystyle 0.0674y_{k-4}+0.4607y_{k-5}$ $\displaystyle-$ $\displaystyle
20.418u^{1}_{k}+33.621u^{1}_{k-1}-9.370u^{1}_{k-2}+9.534u^{1}_{k-3}$
$\displaystyle+$ $\displaystyle 9.417u^{1}_{k-4}+4.002u^{1}_{k-5}$
$\displaystyle+$ $\displaystyle
2.383u^{2}_{k}-4.976u^{2}_{k-1}+5.695u^{2}_{k-2}+25.936u^{2}_{k-3}$
$\displaystyle-$ $\displaystyle 1.737u^{2}_{k-4}+19.616u^{2}_{k-5}$
$\displaystyle-$ $\displaystyle
14.345u^{3}_{k}+14.103u^{3}_{k-1}+25.980u^{3}_{k-2}+32.900u^{3}_{k-3}$
$\displaystyle-$ $\displaystyle 17.978u^{3}_{k-4}-61.739u^{3}_{k-5}$

where the disturbance $w$ is assumed to be independent and normal with a mean
$\mu_{w}$ and standard deviation $\sigma_{w}$. These parameters are also
provided as a part of model identification phase given as $\mu_{w}=-0.0155$
and $\sigma_{w}=268.679$. Since the model is a 5-th order affine system and
sampling time is 15 minutes, inputs to the model are step count measurements
and intervention history in the last 75 minutes and process noise at the
current time.

The objective is to decide when and what type of messages to be sent (input)
based on the identified switched model in order for the subject’s physical
activity (state) to achieve the desired objectives. Then at each time point
during the implementation of the intervention for a new individual, the
developed algorithm uses the data collected on the individual to recommend an
intervention for this individual at this time. Unlike the model given in (7),
the controller requires step count measurements and intervention history from
the beginning of intervention window and noise trajectory until the end of
intervention window. The controller needs step count measurements as an input
because the objective requires total number of step counts from $k=1$ to $k=T$
where $T$ is the final period of intervention window. Since the interventions
are subject to the constraints such as frequency and maximum number of
interventions in the intervention window, the necessity for intervention
history is clear. Given the model above, we use the Algorithm 1, to solve this
optimization problem.

In this application, the step count goal to be achieved is obtained by
analyzing the physical activity of last 30 days. More specifically, the goal
is set to average number of step counts taken in intervention window plus 500
steps for this period. For example, the intervention window for this subject
is between 9:00am and 07:00pm. The average number of step count for this
window is obtained as $5516$ steps from the available data. Hence, the goal is
set to $6016$ steps. Furthermore, the measured step counts are generated using
normal distribution with the mean $\mu_{window}$ and standard deviation
$\sigma_{window}$ of the step counts in intervention window. These parameters
are also estimated using the last 30 days of activity. In these simulations,
the controller performance is analyzed for three different scenarios as
regular, low and high volume physical activity, shown in figures 2, 3, 4,
respectively.

Scenario 1: Regular physical activity
In this scenario, the measurements of step counts are generated using normal
distribution with the mean $\mu_{window}=137$ and standard deviation
$\sigma_{window}=51$.

(a)

(b)

Figure 2: Text message interventions and the probability of achieving the
desired goal for regular physical activity

Scenario 2 & 3: Low and high physical activity
The measurements for low and high activity are generated using normal
distribution with the mean $0.3\mu_{window}$ and $1.5\mu_{window}$
respectively, and standard deviation $\sigma_{window}$ of the step counts in
intervention window.

(a)

(b)

Figure 3: Text message interventions and the probability of achieving the
desired goal for low physical activity.

(a)

(b)

Figure 4: Text message interventions and the probability of achieving the
desired goal for high physical activity.

Figure 3, that is related to low physical activity, shows that although the
maximum number of messages allowed is set to 6, controller sends 4 messages.
This can be interpreted as controller decides to stop sending messages since
the probability of achieving the desired step count goal is so low later in
the day. Similarly, figure 4, in the case of high physical activity,
controller only sends 3 messages out of allowing 6. Again, this can be
interpreted as controller decides to stop sending messages since the
probability of achieving the desired step count goal is approaching to 1.
Thus, this strategy optimizes treatment to achieve the minimum discrepancy
between future activity and a desired step count (maximizing likelihood of
goal attainment) with the fewest messages possible (minimizing burden).

##  V Conclusion

This study contributes to design a model-based intervention MPC controller and
an adaptation strategy that leads to implementing an efficient personalized
adaptive messaging intervention for increasing physical activity, as well as
reducing the burden of treatment. If successful, this strategy should help to
sustain engagement and improve health and well-being.

## References

  * [1]
  * [2] W. H. Organization, Global recommendations on physical activity for health; 2010\. geneva, World Health Organization, Switzerland (2010).
  * [3] J. A. Bennie, K. De Cocker, M. J. Teychenne, W. J. Brown, S. J. H. Biddle, International Journal of Behavioral Nutrition and Physical Activity 16 (1) (2019) 34.  doi:10.1186/s12966-019-0797-2.
  * [4] D. M. Smith, L. Duque, J. C. Huffman, B. C. Healy, C. M. Celano, Text message interventions for physical activity: a systematic review and meta-analysis, American journal of preventive medicine 58 (1) (2020) 142–151.
  * [5] T. B. Turrisi, K. M. Bittel, A. B. West, S. Hojjatinia, S. Hojjatinia, S. K. Mama, C. M. Lagoa, D. E. Conroy, Seasons, weather, and device-measured movement behaviors: a scoping review from 2006 to 2020, International Journal of Behavioral Nutrition and Physical Activity 18 (1) (2021) 1–26.
  * [6] D. E. Conroy, S. Hojjatinia, C. M. Lagoa, C.-H. Yang, S. T. Lanza, J. M. Smyth, Personalized models of physical activity responses to text message micro-interventions: A proof-of-concept application of control systems engineering methods, Psychology of sport and exercise 41 (2019) 172–180.
  * [7] S. Hojjatinia, S. Hojjatinia, C. M. Lagoa, D. Brunke-Reese, D. E. Conroy, Person-specific dose-finding for a digital messaging intervention to promote physical activity, manuscript submitted for publication (2021).
  * [8] E. Sontag, Nonlinear regulation: The piecewise linear approach, IEEE Transactions on automatic control 26 (2) (1981) 346–358.
  * [9] A. Bemporad, A. Garulli, S. Paoletti, A. Vicino, A bounded-error approach to piecewise affine system identification, IEEE Transactions on Automatic Control 50 (10) (2005) 1567–1580.
  * [10] H. J. de Vries, T. J. Kooiman, M. W. van Ittersum, M. van Brussel, M. de Groot, Do activity monitors increase physical activity in adults with overweight or obesity? a systematic review and meta-analysis, Obesity 24 (10) (2016) 2078–2091.
  * [11] M. Morari, J. H. Lee, Model predictive control: past, present and future, Computers & Chemical Engineering 23 (4-5) (1999) 667–682.
  * [12] D. Q. Mayne, Model predictive control: Recent developments and future promise, Automatica 50 (12) (2014) 2967–2986.
  * [13] A. Mesbah, Stochastic model predictive control: An overview and perspectives for future research, IEEE Control Systems Magazine 36 (6) (2016) 30–44.
  * [14] C. M. Lagoa, K. Bekiroglu, S. T. Lanza, S. A. Murphy, Designing adaptive intensive interventions using methods from engineering., Journal of consulting and clinical psychology 82 (5) (2014) 868.
  * [15] A. Zafra-Cabeza, D. E. Rivera, L. M. Collins, M. A. Ridao, E. F. Camacho, A risk-based model predictive control approach to adaptive interventions in behavioral health, IEEE transactions on control systems technology 19 (4) (2010) 891–901.
  * [16] N. N. Nandola, D. E. Rivera, A novel model predictive control formulation for hybrid systems with application to adaptive behavioral interventions, in: proceedings of the 2010 American Control Conference, IEEE, 2010, pp. 6286–6292.
  * [17] J. Skaf, S. Boyd, A. Zeevi, Shrinking-horizon dynamic programming, International Journal of Robust and Nonlinear Control 20 (17) (2010) 1993–2002.
  * [18] J. Luedtke, S. Ahmed, G. L. Nemhauser, An integer programming approach for linear programs with probabilistic constraints, Mathematical programming 122 (2) (2010) 247–272.
  * [19] B. K. Pagnoncelli, S. Ahmed, A. Shapiro, Sample average approximation method for chance constrained programming: theory and applications, J. Optim. Theory Appl. 142 (2) (2009) 399–416.
  * [20] D. T. Eton, T. A. Elraiyah, K. J. Yost, J. L. Ridgeway, A. Johnson, J. S. Egginton, R. J. Mullan, M. H. Murad, P. J. Erwin, V. M. Montori, A systematic review of patient-reported measures of burden of treatment in three chronic diseases, Patient related outcome measures 4 (2013) 7.
  * [21] B. W. Heckman, A. R. Mathew, M. J. Carpenter, Treatment burden and treatment fatigue as barriers to health, Current opinion in psychology 5 (2015) 31–36.