11institutetext: German Research Center for Artificial Intelligence (DFKI) and
Saarland University, Saarbrücken, Germany
11email: {peter.pfeiffer, johannes.lahann<EMAIL_ADDRESS>
# Multivariate Business Process Representation Learning utilizing Gramian
Angular Fields and Convolutional Neural Networks

Peter Pfeiffer  Johannes Lahann  Peter Fettke

###### Abstract

Learning meaningful representations of data is an important aspect of machine
learning and has recently been successfully applied to many domains like
language understanding or computer vision. Instead of training a model for one
specific task, representation learning is about training a model to capture
all useful information in the underlying data and make it accessible for a
predictor. For predictive process analytics, it is essential to have all
explanatory characteristics of a process instance available when making
predictions about the future, as well as for clustering and anomaly detection.
Due to the large variety of perspectives and types within business process
data, generating a good representation is a challenging task. In this paper,
we propose a novel approach for representation learning of business process
instances which can process and combine most perspectives in an event log. In
conjunction with a self-supervised pre-training method, we show the
capabilities of the approach through a visualization of the representation
space and case retrieval. Furthermore, the pre-trained model is fine-tuned to
multiple process prediction tasks and demonstrates its effectiveness in
comparison with existing approaches.

###### Keywords:

Predictive Process Analytics Representation Learning Multi-view Learning

##  1 Introduction

Current machine-learning-based methods for predictive problems on business
process data, e.g., neural-network-based methods like LSTMs or CNNs, achieve
high accuracies in many tasks such as next activity prediction or remaining
time prediction on many publicly available datasets [15]. In recent time, a
large variety of new architectures for next step and outcome prediction have
been proposed and evaluated [19, 17, 24]. These machine-learning-based methods
are mostly task-specific and not generic, i.e., they are designed and tested
on predictive process analytics tasks like next step prediction, outcome
prediction, anomaly detection, or clustering. Moreover, most of the proposed
approaches process only one or a limited set of attributes including activity,
resource [6] or timestamp [23]. On one hand, predicting the next activity in
an ongoing case using the control-flow information is very similar to
predicting the next word in a sentence. On the other hand, data in an event
log recorded from business process executions is very rich in information.
Usually, there are various attributes with different types, scales, and
granularity, which makes generating good representations containing all
characteristics a difficult tasks. Characteristics of a process instance that
are embedded in the data attributes, can improve the performance of predictive
models [13]. For example, the next activity in a business process can depend
on multiple attributes in previous events and complex dependencies between
these attributes [2].

Inspired by recent research in the language modeling domain [5], we propose
and evaluate a novel and generic network architecture and training method –
the Multi-Perspective Process Network (MPPN) – which learns a meaningful,
multi-variate representation of a case in a general-purpose feature vector
that can be used for a variety of tasks. The research contributions of this
paper is threefold:

  1. 1.

We introduce a novel neural-network-based architecture to process a flexible
number of process perspectives of a process instance that examines all of its
characteristics using gramian angular fields.

  2. 2.

For this architecture, we propose an self-supervised pre-training method to
generate a feature-vector representation of a process instance that can be
fine-tuned to various tasks, thus making a contribution to representation
learning for predictive process analytics.

  3. 3.

We show the effectiveness of this approach by analyzing the representation
space and performing an unsupervised case-retrieval task. Furthermore, we
fine-tune and compare the model on a variety of predictive process analytic
tasks such as next step and outcome prediction against existing approaches.

The structure of the remaining chapters unfolds as follows: Section 2
introduces the reader to preliminary concepts. Section 3 discusses related
work on the use of machine learning in predictive process analytics and
representation learning for business process data. Section 4 and 5 present the
proposed approach and the evaluation on a variety of predictive process
analytic tasks. Section 6 closes the paper with a summary of the main
contributions and findings as well as an outline of future work.

##  2 Foundations

###  2.1 Business Process Event Log and Perspectives

Event logs contain records from process-aware information systems in a
structured format. These recordings contain information about what activities
have been conducted by whom at what time as well as additional contextual
data. The following definitions will be used in later sections of the paper.

###### Definition 1

Event Log
An event log is a tupel $L=(E,<,V_{1},\dots,V_{n},a_{1},\dots,a_{n})$, where

  * •

E is the set of events

  * •

$<$ is a total order on E

  * •

$V_{1},\dots,V_{n}$ are the sets of the attribute values

  * •

attributes $a_{i}:E\rightarrow V_{i}$, maps an event to an attribute value.

In the following, we expect an event log to have at least the attributes case-
id, activity, resource, and timestamp.

###### Definition 2

Case
Let $a_{j}$ be the attribute case-id and $v\in V_{j}$. C is the set of events
of one case, iff

  1. 1.

$C\subseteq E$

  2. 2.

For each $e\in C:a_{j}(e)=v$

  3. 3.

For each $e\in E\setminus C:a_{j}(e)\neq v$

Let $C=\\{e_{1},\dots,e_{n}\\}$ be the events of a case $c$. These events
follow the order $<$. We use the notation: $c=<e_{1},\dots,e_{n}>$, where
$e_{i}<e_{j}$ if $i<j$.

###### Definition 3

Business Process Perspective
Given $L=(E,<,V_{1},\dots,V_{n},a_{1},\dots,a_{n})$ and a sequence of events
$\langle e_{1},\dots,e_{n}\rangle$, we define a perspective on each attribute
$a_{i}$ as $\Pi_{a_{i}}:=\langle a_{i}(e_{1}),\dots,a_{i}(e_{n})\rangle$

In the next sections, we frequently use the control-flow perspective
$\Pi_{control-flow}$, the resource perspective $\Pi_{resource}$ and the
temporal perspective $\Pi_{time}$.

Last, we distinguish between event attributes and case attributes. A case
attribute returns the same value for all events in a case, i.e. it fulfils
equation 2 in definition 2. Otherwise, it is an event attribute.

###  2.2 Business Process Data and Representation Learning

Representation learning is the task of learning ”representations of data that
make it easier to extract useful information when building classifiers or
other predictors” [1]. For example, embeddings are utilized in natural
language processing to learn a vectorized representation of words. In recent
times, attention-based networks [25] have shown superior performance in many
language tasks such as machine translation, mainly due to their ability to
generate meaningful representations. Usually, these types of networks are pre-
trained in an unsupervised fashion on extensive datasets. For business process
data, embeddings are commonly used in predictive process analytics to
represent the control-flow [3], or certain perspectives [6] in a vector space.
Thus, it should allow the model to exploit the vector representation more
effectively than a one-hot or integer encoding. When learning good
representations of cases, one tries to represent all relevant characteristics
within the representation. For predictive tasks, this includes the underlying
distribution of the exploratory factors, i.e., attributes that influence the
prediction [2]. Usually, there are several attributes with different types of
data – categorical, numerical, and temporal ones. Within each type, the
dimensions, scales, and variabilities can be different. For instance, there
can be a numerical attribute cost ranging from [0, 1,000,000] and another one,
e.g., discount that is in range [0, 30]. Temporal attributes also have
different scales (daily, weekly, monthly, etc.) and high variability.
Categorical attributes often vary in their dimensionality. While the activity
has few different values, the resource, e.g., persons involved in a process,
often has much more distinct values. Some perspectives are very spare while
others are rich. Furthermore, different from event-related attributes, there
are also case-related attributes. Some attributes change their values only in
certain events of a case, while others have a different value for each event.
This, in turn, means that perspectives have different levels of granularity.
When learning representations of cases, the multivariate, multi-scalar, and
multi-granular nature of the data must be considered and depicted.

##  3 Related Work

In [4], the authors introduced methods to learn representations of activities,
cases, and process models. They trained a model on the next step prediction
task to learn representations similar to obtaining embeddings for words,
sentences, and documents. Although they explained that other attributes
besides activity are important, they only consider the control-flow
perspective. Furthermore, they did not include an extensive evaluation to
elaborate on the effectiveness of the proposed representations. Apart from
that, representation learning is not explicitly tackled in existing predictive
analytic approaches. However, these approaches learn a representation
alongside a specific prediction task. [6] was the first to introduce neural
networks to the field of process prediction. They applied recurrent neural
networks to next activity and remaining time prediction. They trained separate
neural networks considering the control-flow, resource, and time perspectives.
[23] examined the next step and remaining time prediction task. They used an
LSTM-based approach that optimized both tasks simultaneously and elaborated on
the effect of separated or shared LSTM network layers. [3] elaborated three
different LSTM architectures for predicting the next activity, resource, and
timestamp. In their first architecture, they used specialized layers for each
attribute to predict. The second version combines the categorical attributes
in a shares layer, while the third version shares categorical attributes and
the timestamp. In [17], the authors proposed an LSTM-based method that can
combine multiple attributes. In their model, the authors use embeddings for
each categorical attribute while non-categorical attributes are concatenated
with the categorical ones’ embedded representation. Other approaches for next
step prediction used different techniques such as decision trees (DT) [9],
autoencoder (AE) with n-grams [10], attention networks [12], CNNs [18] or
generative adversarial networks (GAN) [24]. Similarly, CNNs [19], LSTMs [14]
or autoencoder [11] are used for outcome prediction. In order to detect
anomalies in business process data, LSTMs [16] or Bayesian neural networks
[20] were applied. Table 1 gives an overview of existing predictive approaches
and categorizes them by prediction task, examined perspective, encoding
technique per perspective, as well as used machine learning method. Also, it
clarifies what information is available to which model in what form. While all
predictive approaches use at least information from two perspectives, only a
few approaches are able to encode and process all types of attributes. We
differentiated between the most commonly used perspectives and other attribute
types to delimit generic and non-generic approaches. Only [17] used a generic
encoding approach that can process and represent all types of attributes.
However, the approach is tailored towards next step prediction and does not
focus on the learned representation within the model. Thus, we propose a
generic multi-attribute representation learning approach that is not tailored
to a specific prediction task.

Table 1: Overview of encoding techniques used in literature. duration:
timestamps to duration after a certain timestamp, IMG: one or multiple
perspectives $\Pi$ to a single or multiple matrices $M^{h\times w\times c}$ of
size $height\times width\times\\#channels$.

Task | Approach |  ML Method | Encoding Techniques for event attributes
---|---|---|---
Most common perspectives | Other event log attributes
$\Pi_{control-flow}$ |  $\Pi_{resource}$ |  $\Pi_{timestamp}$ |  Categorical |  Numerical |  Temporal
Next Step Prediction | [6] |  LSTM | embedding |  - |  - |  - |  -
[23] |  LSTM |  one-hot |  - |  custom |  - |  - |  -
[10] |  AE + FF |  N-gram |  one-hot |  - |  one-hot |  as-is |  -
[9] |  DT |  integer |  |  |  integer |  |
[18] |  CNN | IMG |  - |  - |  -
[12] |  Attention |  one-hot |  one-hot |  - |  one-hot |  as-is |  -
[17] |  LSTM |  embedding |  embedding |  duration |  integer |  as-is |  duration
[3] |  LSTM |  embedding |  embedding |  duration |  - |  - |  -
[24] |  GAN |  one-hot |  - |  duration |  - |  - |  -
Outcome Prediction | [19] |  CNN | IMG |  - |  - |  -
[11] |  FF |  custom |  custom |  custom |  custom |  custom |  custom
| [14] |  LSTM |  one-hot |  one-hot |  custom |  one-hot |  - |  -
Anomaly Detection | [16] |  LSTM |  integer |  integer |  - |  integer |  - |  -
[20] |  Bayesian NN |  probability |  probability |  probability |  probability |  probability |  probability
All | MPPN |  CNN |  IMG |  IMG |  IMG |  IMG |  IMG |  IMG

##  4 Multi-Perspective Process Network (MPPN)

The MPPN approach for representation learning is mainly built on two concepts
– graphical event log encoding and neural-network-based processing. The first
part is to encode the perspectives of interest $\hat{\Pi}$ in the event log
$L$ uniformly as 2D images by transforming them to distinct gramian angular
fields (GAF) – no matter if the perspective contains categorical, numerical or
temporal information. The second part is a convolutional neural network
architecture and training method that learns representations of cases using
the GAF-encoded perspectives. Figure 1 shows the architecture and processing
pipeline of the MPPN approach. In this example, the six perspectives of
interest $\hat{\Pi}=\\{\Pi_{control-flow},\\\
\Pi_{timestamp},\Pi_{type},\Pi_{resource},\Pi_{travel\\_start}\,\Pi_{cost}$ of
the case with case-id 1565 are first encoded as six individual GAFs, which can
then be processed by the MPPN. CNN1 extracts the features of each perspective
$\hat{\Pi}$ which are combined in the perspective pooling layer. The combined
features capture all characteristics of interest of a particular case.

Figure 1: Architecture and processing pipeline of MPPN using a case from the
MobIS event log [7].

The forwards pass in MPPN of a single case $c$ is as follows: For each GAF-
encoded perspective, CNN1 extracts a feature vector that is pooled before
being passed to NN2. NN2 then takes the pooled features from all perspectives,
processes them, and produces a single feature vector $FV$. This two-stage
architecture allows CNN1 to focus on the features within each perspective
while NN2 captures and models the dependencies between perspectives. By
transforming all perspectives uniformly to GAFs, all attributes lie within the
same range, no matter what scale or variability they had before. At the same
time, NN2 can learn what features from what perspectives are important. Unlike
RNN-based models, MPPN consumes the whole case at once instead of being fed
with cases event by event.
Inspired by multi-view learning, e.g., Multi-View Convolutional Neural
Networks (MVCNN) for 3D object detection using 2D renderings [22], the Multi-
Perspective Process Network creates a feature vector $FV$ for a case using the
GAF-encoded perspectives. Analogous to using multiple renders from different
views to represent 3D structures, we use different perspectives of a process
to represent a case. Another important aspect of MPPN is its ability to be
used for several tasks instead of being task-specific. This is achieved by an
self-supervised pre-training phase, as also done in[10], that learns a
representation in the form of a feature vector $FV$. Afterwards, a task-
specific layer can be added to the pre-trained model allowing to fine-tune the
model on different tasks. The learned representation $FV$ thus serves as the
basis for any downstream task.

###  4.1 Graphical Representation of Event Log Data

In order to encode all types of attributes into a single representation in a
generic way, we decided to choose a graphical encoding instead of the methods
used in related work shown in table 1. We see a strong similarity in the
characteristics of time-series and a single perspective of a process case.
Furthermore, all types of attributes in a case can easily be transformed to
time-series. For this reason, we treat the perspectives $\hat{\Pi}$ as
multivariate time series. A naive way to get a both machine-readable and
visualizable representations of perspectives $\Pi$ is to represent and plot
them as time series. Thus, each value of a perspective is encoded as a real
number and visualized in a 2D representation. The y-coordinate corresponds to
the value $v$, and the x-coordinate to $t$. In figure 1, one can see the 6
perspectives $\hat{\Pi}$ of case 1565 encoded and plotted as 6 distinct time
series. Although this representation is a nice visualization for humans,
presenting the perspectives $\Pi$ as a time series plot is a very naive way.
Such a plot is very sparse, i.e., most of the plot is empty with just a fine
line drawn, containing only little information for convolutional neural
networks.

Gramian angular fields (GAF), originally proposed for time-series
classification, transform sequential data to 2D images, which contain more
information for machine learning methods as time-series plots [26]. For a
sequence $\langle v_{1},...,v_{n}\rangle$, a gramian angular field is a matrix
of size $n\times n$ where each entry is the cosine of the sum of two polar
coordinates in time – the polar coordinate of $v_{i}$ plus $v_{j}$. This
projection is bijective and preserves temporal relations. To transform event
log data, i.e., all perspectives $\hat{\Pi}$ of categorical, numerical, and
temporal event log attributes and case attributes into gramian angular fields,
they must be treated and transformed to distinct sequences of numerical
values. In order to get numerical sequences from each type, the following
transformations and encodings are performed. Other types of attributes can
also be used (e.g., textual data) if they are encoded as numerical sequences.

  1. 1.

For categorical attributes $a_{i}$, we applied an integer encoding
$integer:\mathcal{V}_{i}\rightarrow int$ where
$int\in[0,1,2,..,|\mathcal{V}_{i}|-1]$.

  2. 2.

Timestamps are transformed into the duration in seconds from the earliest
timestamp.

  3. 3.

Numerical attributes are used unchanged.

Case attributes are first duplicated to the case length before being encoded
in the same way as event attributes. Once encoded as numerical sequences, each
perspective can easily be encoded as a gramian angular field after scaling
them to a [-1, 1] range. To ensure equal size images where characteristics are
equally represented, the sequences are adjusted to equal length, either by
padding or truncating. This results in distinct GAF-representations for each
perspective $\hat{\Pi}$ of a case as shown in figure 1.
By using graphical encodings we can transform attributes of different types to
images and use state-of-the-art image processing neural networks. This way we
avoid building networks that process cases with customized architectures for
specific attribute types as well as the complexity of training embeddings.

###  4.2 Architecture

The MPPN architecture consists of three parts as shown in figure 1: CNN1 for
feature extraction, NN2 for modeling dependencies and relations between
perspectives, and one or multiple task-specific layers called HEAD. Between
CNN1 and NN2, a pooling layer combines the features produced by CNN1 for each
GAF-encoded perspective to a single vector. The weights in CNN1 are shared
between all perspectives, i.e. the same CNN1 is applied on all perspectives.
For CNN1, we use Alexnet [8]. However, as gramian angular fields are different
from natural images, pre-training CNN1 on GAFs significantly reduces the later
training time. NN2 is a fully-connected neural network. Together, CNN1 and NN2
form the model used for representation learning that produces $FV$. $FV$ can
either be used directly or by any other task-specific layer; e.g., a fully-
connected HEAD with softmax for next step prediction or a HEAD for remaining
time prediction.

###  4.3 Training Method

One integral part of MPPN is its ability to learn representations of all
perspectives $\hat{\Pi}$ of cases in an event log. In order to obtain good
feature vectors $FV$, one must ensure that all relevant characteristics are
fully captured in the model. We distinguish three stages of training that
should be performed successively.

####  4.3.1 Pre-Training CNN1 on GAFs

As GAFs are very different from natural images, pre-trained CNNs like Alexnet
need to be fine-tuned. While lower-level features like edges and corners are
present in GAFs too, higher-level features differ. In order to make the MPPN
sensitive to GAF-specific feature, we fine-tuned the CNN1 once by classifying
cases according to their variant. This task has been chosen as the process
variant is always directly derivable from the sequence of events and the model
can focus on learning features from the single GAF-encoded perspective. All
relevant information for this task is entailed in the GAF image which is what
we want the model to focus on. However, many other tasks are also possible. In
detail, we build a MPPN with a pre-trained Alexnet consuming only
$\Pi_{control-flow}$ and predicting the variant on the MobIS dataset. For each
case $c$ the whole sequence of $activity$ was used as input and the variant
used as the target. Afterwards, the weights of CNN1 are saved on disk and can
be used on any dataset, any perspective, and any task for MPPN in the future.

####  4.3.2 Representation Learning

To obtain meaningful feature vectors $FV$ of business process cases, MPPN must
be trained to hold all characteristics of a case in it. One can train MPPN on
next step prediction tasks, e.g., to predict the next activity in an ongoing
case given $\hat{\Pi}$. This works fine, but the model will learn the relation
in the data, which are important for the next activity. This leads to a
feature vector $FV$ that by design holds features that are important to
predict the next activity. Attributes that do not have relevance for the next
activity will be less present in $FV$.
To obtain more generic feature vectors of cases, a self-supervised multi-task
next event prediction training method is applied that trains the network to
predict $a_{i}(e_{t+1})$ for each attribute in $\hat{\Pi}$. For this task, the
MPPN architecture is extended by small networks $HEAD_{a_{i}}$ – one for each
attribute $a_{i}$ to predict. Each $HEAD$ is a task-specific layer that
consumes $FV$ and predicts $a_{i}(e_{t+1})$. During representation learning,
the task’s criterion is to minimize the sum of all losses of all predictions,
measured as mean absolute error (for numerical and temporal attributes) and
cross-entropy (for categorical attributes). During training, all $HEADs$ are
trained in parallel and in conjunction with the rest of MPPN. Thereby, the
MPPN and especially NN2 learns to focus on important features in all
perspectives $\hat{\Pi}$ and produces a $FV$ that holds information relevant
for the attributes in the next event. Using this method, a representation can
be learned without the need for manual labeled data. However, depending on the
final task to be solved, other training methods are also possible. As long as
all relevant characteristics of the case are enclosed in $FV$, any training
method is appropriate. We chose the multi-task next event prediction task as
it allows the model to incorporate all attributes for each prediction. While
making predictions for each attribute the model is forced to not drop relevant
characteristics of a case. Afterwards, the weights of MPPN (without the heads)
are saved on disk. The $FV$ produced in this state can directly be used for
tasks where additional labels are hard to obtain or unavailable, such as
clustering, retrieval or anomaly detection using the same dataset.

####  4.3.3 Fine-tuning on Specific Tasks

After being trained to learn good representations, MPPN can also be fine-tuned
on other tasks using the same event log and given appropriate labels.
Therefore, one or multiple $HEADs$ are added that consume $FV$. With each
HEAD, the model and especially the HEAD can be trained on a large variety of
tasks, e.g., outcome prediction, next step prediction or (supervised) anomaly
detection. Thereby, the model makes use of the representation in $FV$ to solve
a certain problem.

###  4.4 Implementation Details

We implemented MPPN with the following hyperparameter choices: We padded or
truncated all cases $c$ to length 64 which results in GAF images of size
$64\times 64$ pixel. CNN1 consists of four CNN layers with max-pooling and
dropout. NN2 is a two-layer fully-connected network with dropout. We pooled
the perspectives behind CNN1 by concatenation. The HEADs consist of shallow
fully-connected networks with a softmax or regression layer. More details can
be found in the implementation.

##  5 Evaluation

This section elaborates on two experiments. The first experiment visualizes
the learned representations during the self-supervised pre-training phase and
demonstrates a contextual retrieval task. In the second experiment, we compare
the MPPN model to existing approaches on next event and outcome prediction
tasks by fine-tuning the pre-trained model.

###  5.1 Representation Visualization and Retrieval

In the following, we demonstrate how MPPN’s internal representations $FV$ can
be used for case-based case retrieval. Figure 3 visualizes $FVs$ of each cases
after they were reduced to a two-dimensional representation space using PCA.
The training of the MPPN was performed analog to section 4.3.2 using the same
input attributes as described in table 2111We added $travel\\_start$ as
another attribute but complete cases $c$ instead of prefixes. Note that the
feature vectors hold information of all perspectives. Therefore, the clusters
do not solely depend on the control-flow.

Figure 2: Visualization of the representation space learned by the MPPN on all
MobIS cases. Different colors indicate different control-flow variants.

$c_{query}$ | $\hat{\mathcal{C}}$ ID | $FV$ distance | DLD | MAE
---|---|---|---|---
timestamp | cost | travel_start
5523 | 5511 | 0.00411 | 0 | 101.66 | 240 | 0.53
5613 | 0.00479 | 0 | 154.82 | 154 | 6.47
5036 | 0.01665 | 5 | 1307.83 | 244 | 28.53
5911 | 0.01755 | 8 | 1004.02 | 203 | 24.47
6034 | 0.01937 | 8 | 917.40 | 253 | 31.93
5980 | 0.02088 | 8 | 933.96 | 237 | 28.55
| 5868 | 0.02115 | 8 | 1045.91 | 69 | 21.55
2056 | 4819 | 0.01388 | 0 | 2066.35 | 49 | 174.00
4960 | 0.02068 | 5 | 3587.52 | 218 | 181.33
4765 | 0.02295 | 5 | 729.69 | 253 | 169.15
4497 | 0.02340 | 5 | 632.51 | 217 | 153.00
4715 | 0.02428 | 5 | 717.51 | 263 | 167.00
5044 | 0.02453 | 5 | 847.96 | 375 | 188.00
| 4657 | 0.02465 | 5 | 689.45 | 233 | 162.39
7222 | 7109 | 0.00006 | 0 | 14.71 | 97 | 7.35
7092 | 0.00006 | 0 | 16.86 | 94 | 8.42
7073 | 0.00012 | 0 | 18.01 | 77 | 9.00
7090 | 0.00015 | 0 | 17.10 | 24 | 8.54
7133 | 0.00016 | 0 | 11.72 | 32 | 5.86
7231 | 0.00017 | 0 | 0.54 | 100 | 0.27
| 7052 | 0.00021 | 0 | 18.01 | 41 | 9.00
3006 | 3227 | 0.00048 | 0 | 55.97 | 392 | 153.00
2403 | 0.00105 | 0 | 164.74 | 1123 | 0.00
2624 | 0.00118 | 0 | 54.72 | 501 | 30.00
3748 | 0.0012 | 0 | 206.65 | 662 | 153.00
2859 | 0.00123 | 0 | 103.17 | 629 | 38.00
2861 | 0.0014 | 0 | 89.77 | 474 | 52.00
| 2116 | 0.00153 | 0 | 287.39 | 250 | 40.00

Figure 3: Similarities in the perspectives of the retrieved cases

Figure 3 shows that some clusters consist of cases with the same process
variant. Other clusters are formed based on specific attribute combinations.
For example, the biggest bulk shows all finished cases, i.e., complete cases
from start to end containing the most common variant, represented by case
3006. One can make use of this representation for case-based case retrieval.
Given $L$ and a query case $c_{query}$, the task is to generate an ordered set
of cases $\hat{\mathcal{C}}$ such that all cases in $\hat{\mathcal{C}}$ have
similar characteristics as $c_{query}$. Instead of applying different filters
on an event log to retrieve cases with particular characteristics, one can
also retrieve cases starting with a specific case of interest. For this, the
same feature vectors $FVs$ can now be used for retrieving such cases that
share similar characteristics as a query case. First, the feature vector of
the query case $FV_{query}$ is computed and compared to all other $FV$ of
cases in $L$ using the cosine similarity. Next, the cases are sorted by their
similarity, and those with the highest similarity are returned. We picked four
cases as shown in table 3 for retrieval and marked the retrieved cases with
bold symbols in figure 3. We see that the control-flow still is the deciding
feature for the model as most of the retrieved cases have the same sequence of
activities. Additionally, the retrieved cases have similar other
characteristics as the query case:

  * •

5523: Different process variants starting and ending with the same activities
performed around the same date with cost below 1000.

  * •

2056: Cases that looped through the same activities with various number of
this loop.

  * •

7222: Cases consisting of the first two events in the process performed around
the same date with costs around 200.

  * •

3006: Complete cases from start to end of the most common variant.

From table 3 one can see that the retrieved cases $\hat{\mathcal{C}}$ are
similar in all perspectives to the query cases. We calculate the cosine
distance of the $FV$, the Damerau-Levenshtein distance (DLD) and the mean
absolute error for the three perspectives $\Pi_{cost}$, $\Pi_{travel\\_start}$
and $\Pi_{timestamp}$ (the MAE is computed after transforming the timestamps
to durations).

###  5.2 Next Step and Outcome Prediction

This experiment evaluates the performance of the fine-tuned MPPN model in
comparison to four baselines on the tasks next activity, last activity, next
resource, last resource, event duration, and remaining time prediction.

####  5.2.1 Datasets

In this experiment, we consider seven event logs from different application
domains. The Helpdesk222http://doi.org/10.17632/39bp3vv62t.1 event log
contains events from a ticketing management process of the help desk of an
Italian software company. Five event logs from the BPI Challenge
2012333https://doi.org/10.4121/uuid:3926db30-f712-4394-aebc-75976070e91f. The
original event log is taken from a Dutch Financial Institute and represents
the application process for a personal loan or overdraft within a global
financing organization. We included the original log as well as each sub-
process individually. The event log within BPI Challenge
2013444https://doi.org/10.4121/uuid:a7ce5c55-03a7-4583-b855-98b86e1a2b07 is an
export from Volvo IT Belgium and contains events from an incident and problem
management system called VINST. The event log within BPI Challenge
2017555https://doi.org/10.4121/uuid:5f3067df-f10b-45da-b98b-86ae4c7a310b is an
updated, richer version of BPI Challenge 2012. The event log from BPI
Challenge
2020666https://doi.org/10.4121/uuid:52fb97d4-4588-43c9-9d04-3604d4613b51 was
collected data from the reimbursement process at TU/e. We only included the
request-for-payment log. The MobIS event
log777http://dx.doi.org/10.13140/RG.2.2.11870.28487 was elaborated in the
MobIS Challenge [7]. It describes the execution of a business travel
management process in a medium-sized consulting company. We chose Helpdesk,
BPIC 2012, and BPIC 2013 to achieve high comparability with existing
approaches. BPIC 2017 and BPIC 2020 are selected as significantly more complex
event logs that pose new challenges to prediction approaches while also
revealing weaknesses of current approaches. MobIS contains several attributes
and relationships, making it well-suited to demonstrate MPPN’s multi-
perspective approach’s benefits. Table 2 lists characteristics of each log and
presents the attributes used as inputs for the process prediction tasks.

Table 2: Event logs, statistics and attributes used

| #Traces | #Events |  Avg. trace length |  Avg. trace duration | Input Attributes |
---|---|---|---|---|---|---
|  |  |  categorical |  numerical |  temporal
Helpdesk | 4580 | 21348 |  4.66 |  62.9 days |  activity, resource |  |  timestamp
BPIC12 | 13087 | 262200 |  20.04 |  150.2 days |  activity, resource |  AMOUNT_REQ |  timestamp
BPIC12_Wc | 9658 | 72413 |  7.50 |  95.6 days |  activity, resource |  AMOUNT_REQ |  timestamp
BPIC13_CP | 1487 | 6660 |  4.48 |  426.5 days |  activity, resource, resource country, organization country, organization involved, impact, product, org:role |  |  timestamp
BPIC17_O | 42995 | 193849 |  4.51 |  23.9 days |  activity, Action, NumberOfTerms, resource |  FirstWithdrawalAmount, MonthlyCost, OfferedAmount, CreditScore |  timestamp
BPIC20_RFP | 6886 | 36796 |  5.34 |  31.6 days |  org:role, activity, resource, Project, Task, OrganizationalEntity |  RequestedAmount |  timestamp
MobIS | 6555 | 166512 |  25.40 |  1194.4 days |  activity, resource, type |  cost |  timestamp

####  5.2.2 Experimental Setup

We compare MPPN with four different approaches [6, 23, 3, 17]. For each task,
the models receive as input case prefixes of increasing length, starting with
the prefix that contains only the first event of a case up to the prefix that
omits just the last event; i.e., for each case $\langle
e_{1},...,e_{n}\rangle$ we create $n$ prefixes $\langle
e_{1},\dots,e_{t}\rangle$ with $0<=t<n$. In addition, we front-padded all
prefixes to equal length. To make the results reproducible, we apply a random
split between training, validation, and test cases for each data set,
utilizing $20\%$ of the cases as test and $10\%$ of the cases as validation
data. While the test set is fixed through all experiments and runs, the split
between train and validation is performed randomly from run to run. All models
were trained and validated with the same sets in each run. Each model was
trained in the same fashion with a batch size of 512 while utilizing cyclical
learning rates and early stopping [21]. The learning rate was picked with the
learning rate finder algorithm as defined in [21]. Other than that, we picked
the hyper-parameters of the baselines as mentioned in the corresponding
papers. While [6, 23, 3] only considered control flow, resource and timestamp
perspectives, the MiDA and the MPPN model is fed with all attributes listed in
table 2. We only removed attributes that contained duplicated information.
Last, we decided to remove all cases that are longer than 64 events since
these are mostly outliers that falsify the prediction results and
significantly increase training time. Each model was trained and tested ten
times on all datasets and tasks.

####  5.2.3 Prediction Tasks and Evaluation Metrics

For this experiment, we formalize the prediction tasks and evaluation metrics
as follows:
Given a prefix $p_{t}=\langle e_{1},...,e_{t}\rangle$ of a case $c=\langle
e_{1},...,e_{n}\rangle$ with $0<=t<n;t,n\in\mathbb{N}$, we define next step
prediction of an attribute $a$ as the task $NSP_{a}(p_{t})$ that predicts
$a(e_{t+1})$ based of the prefix $p_{t}$. We define outcome prediction
analogously to next step prediction as the task $OUT_{a}(p_{t})$ that predicts
$a(e_{n})$ based on a prefix $p_{t}$. We measure the prediction performance of
a model through a metric function, which is pairwise applied to all
predictions and ground truth values for all prefixes over all cases and
afterward combined to a final score. According to the type of the predicted
attribute, it is necessary, to use different metric functions. In this
experiment, we predict activity, resource, and timestamp. For activity and
resource, we select the metric function accuracy. For timestamp, we convert it
in the duration in days and then compute the mean absolute error.

####  5.2.4 Baselines

We re-implemented eight models as baselines from [6, 23, 3, 17] based on the
original papers and the corresponding source code. Our main objective is to
reproduce the different network architectures, to be able to compare them in a
fair and unified test setting with our MPPN.

Table 3: Process prediction results

Dataset |  Model |  $NSP_{activity}$ |  $NSP_{resource}$ |  $OUT_{activity}$ |  $OUT_{resource}$ |  $NSP_{timestamp}$ |  $OUT_{timestamp}$
---|---|---|---|---|---|---|---
Helpdesk |  Evermann[6] |  0.651+-0.128 |  0.222+-0.005 |  0.994+-0.000 |  0.811+-0.000 |  — |  —
|  Ca_Spez.[3] |  0.693+-0.168 |  0.289+-0.071 |  0.994+-0.000 |  0.811+-0.000 |  7.95+-0.576 |  6.654+-0.101
|  Ca_concat[3] |  0.696+-0.116 |  0.421+-0.035 |  0.994+-0.000 |  0.811+-0.000 |  7.63+-0.052 |  6.739+-0.253
|  Ca_full[3] |  0.774+-0.077 |  0.432+-0.000 |  0.994+-0.000 |  0.811+-0.000 |  5.308+-0.288 |  7.018+-0.225
|  Tax_Spez.[23] |  0.763+-0.082 |  — |  0.994+-0.000 |  — |  7.777+-0.526 |  6.895+-0.253
|  Tax_Mixed[23] |  0.3+-0.003 |  — |  0.994+-0.000 |  — |  14.849+-0.034 |  7.197+-0.101
|  Tax_Shared[23] |  0.793+-0.004 |  — |  0.994+-0.000 |  — |  5.088+-0.129 |  6.67+-0.100
|  MiDA[17] |  0.693+-0.120 |  0.263+-0.089 |  0.994+-0.000 |  0.811+-0.000 |  4.898+-0.043 |  6.629+-0.166
|  MPPN |  0.805+-0.003 |  0.691+-0.006 |  0.994+-0.000 |  0.847+-0.008 |  5.197+-0.126 |  6.691+-0.089
BPIC12 |  Evermann[6] |  0.595+-0.107 |  0.149+-0.000 |  0.417+-0.000 |  0.172+-0.000 |  — |  —
|  Ca_Spez.[3] |  0.795+-0.030 |  0.333+-0.282 |  0.417+-0.000 |  0.177+-0.015 |  0.693+-0.208 |  7.82+-0.033
|  Ca_concat[3] |  0.74+-0.071 |  0.426+-0.164 |  0.417+-0.000 |  0.172+-0.000 |  0.722+-0.206 |  7.849+-0.076
|  Ca_full[3] |  0.756+-0.064 |  0.283+-0.197 |  0.417+-0.000 |  0.184+-0.025 |  0.687+-0.226 |  6.649+-0.084
|  Tax_Spez.[23] |  0.585+-0.194 |  — |  0.417+-0.000 |  — |  0.734+-0.155 |  7.477+-0.127
|  Tax_Mixed[23] |  0.615+-0.182 |  — |  0.417+-0.000 |  — |  0.544+-0.205 |  6.678+-0.101
|  Tax_Shared[23] |  0.824+-0.008 |  — |  0.487+-0.019 |  — |  0.542+-0.167 |  6.693+-0.080
|  MiDA[17] |  0.565+-0.123 |  0.149+-0.000 |  0.417+-0.000 |  0.172+-0.000 |  0.625+-0.041 |  6.587+-0.047
|  MPPN |  0.846+-0.006 |  0.775+-0.002 |  0.53+-0.005 |  0.316+-0.004 |  0.82+-0.079 |  6.694+-0.066
BPIC12_Wc |  Evermann[6] |  0.774+-0.000 |  0.104+-0.000 |  0.435+-0.000 |  0.11+-0.000 |  — |  —
|  Ca_Spez.[3] |  0.775+-0.002 |  0.104+-0.000 |  0.435+-0.000 |  0.113+-0.005 |  1.799+-0.088 |  8.31+-0.058
|  Ca_concat[3] |  0.794+-0.027 |  0.104+-0.000 |  0.435+-0.000 |  0.115+-0.013 |  1.843+-0.169 |  8.333+-0.041
|  Ca_full[3] |  0.792+-0.026 |  0.104+-0.000 |  0.443+-0.026 |  0.112+-0.005 |  1.81+-0.125 |  7.455+-0.063
|  Tax_Spez.[23] |  0.713+-0.081 |  — |  0.435+-0.000 |  — |  1.765+-0.098 |  7.932+-0.086
|  Tax_Mixed[23] |  0.774+-0.000 |  — |  0.435+-0.000 |  — |  1.595+-0.064 |  7.51+-0.106
|  Tax_Shared[23] |  0.773+-0.001 |  — |  0.537+-0.057 |  — |  1.645+-0.070 |  7.409+-0.110
|  MiDA[17] |  0.805+-0.022 |  0.104+-0.000 |  0.435+-0.000 |  0.155+-0.028 |  1.767+-0.115 |  7.424+-0.103
|  MPPN |  0.815+-0.006 |  0.237+-0.011 |  0.558+-0.01 |  0.147+-0.009 |  1.761+-0.061 |  7.528+-0.072
BPIC13_CP |  Evermann[6] |  0.417+-0.113 |  0.082+-0.005 |  1.0+-0.000 |  0.211+-0.000 |  — |  —
|  Ca_Spez.[3] |  0.481+-0.090 |  0.086+-0.000 |  1.0+-0.000 |  0.211+-0.000 |  50.927+-2.339 |  137.718+-3.125
|  Ca_concat[3] |  0.524+-0.004 |  0.086+-0.000 |  1.0+-0.000 |  0.211+-0.000 |  51.672+-3.62 |  139.412+-5.032
|  Ca_full[3] |  0.493+-0.065 |  0.106+-0.035 |  1.0+-0.000 |  0.211+-0.000 |  67.168+-7.233 |  137.193+-6.280
|  Tax_Spez.[23] |  0.502+-0.067 |  — |  1.0+-0.000 |  — |  50.785+-4.395 |  140.481+-4.913
|  Tax_Mixed[23] |  0.309+-0.003 |  — |  1.0+-0.000 |  — |  112.867+-0.279 |  176.167+-0.930
|  Tax_Shared[23] |  0.51+-0.011 |  — |  1.0+-0.000 |  — |  47.741+-1.217 |  144.528+-21.964
|  MiDA[17] |  0.434+-0.110 |  0.083+-0.005 |  1.0+-0.000 |  0.211+-0.000 |  54.949+-4.044 |  128.185+-10.555
|  MPPN |  0.562+-0.009 |  0.178+-0.024 |  1.0+-0.000 |  0.216+-0.008 |  54.922+-3.948 |  127.824+-3.806
BPIC17_O |  Evermann[6] |  0.818+-0.000 |  0.067+-0.005 |  0.509+-0.032 |  0.186+-0.041 |  — |  —
|  Ca_Spez.[3] |  0.818+-0.000 |  0.064+-0.000 |  0.513+-0.018 |  0.192+-0.000 |  3.628+-0.057 |  9.604+-0.017
|  Ca_concat[3] |  0.818+-0.000 |  0.226+-0.261 |  0.501+-0.027 |  0.192+-0.000 |  3.611+-0.082 |  9.606+-0.014
|  Ca_full[3] |  0.818+-0.000 |  0.081+-0.048 |  0.52+-0.001 |  0.192+-0.000 |  3.627+-0.105 |  9.519+-0.025
|  Tax_Spez.[23] |  0.67+-0.065 |  — |  0.454+-0.019 |  — |  3.529+-0.019 |  9.688+-0.145
|  Tax_Mixed[23] |  0.726+-0.178 |  — |  0.458+-0.014 |  — |  3.999+-0.503 |  9.768+-0.184
|  Tax_Shared[23] |  0.818+-0.000 |  — |  0.519+-0.000 |  — |  3.531+-0.037 |  9.47+-0.021
|  MiDA[17] |  0.836+-0.030 |  0.064+-0.000 |  0.828+-0.002 |  0.192+-0.000 |  3.297+-0.037 |  8.946+-0.059
|  MPPN |  0.818+-0.000 |  0.553+-0.061 |  0.518+-0.001 |  0.208+-0.001 |  3.567+-0.068 |  9.534+-0.016
BPIC20_RFP |  Evermann[6] |  0.699+-0.099 |  0.817+-0.084 |  0.957+-0.000 |  0.958+-0.000 |  — |  —
|  Ca_Spez.[3] |  0.756+-0.087 |  0.841+-0.020 |  0.957+-0.000 |  0.958+-0.000 |  2.556+-0.142 |  6.068+-0.185
|  Ca_concat[3] |  0.704+-0.09 |  0.997+-0.000 |  0.957+-0.000 |  0.958+-0.000 |  2.631+-0.199 |  6.062+-0.079
|  Ca_full[3] |  0.804+-0.025 |  0.997+-0.001 |  0.957+-0.000 |  0.958+-0.000 |  2.634+-0.252 |  5.931+-0.117
|  Tax_Spez.[23] |  0.791+-0.085 |  — |  0.957+-0.000 |  — |  2.269+-0.085 |  5.933+-0.087
|  Tax_Mixed[23] |  0.431+-0.252 |  — |  0.957+-0.000 |  — |  3.827+-2.194 |  8.55+-0.058
|  Tax_Shared[23] |  0.849+-0.001 |  — |  0.957+-0.000 |  — |  2.12+-0.095 |  5.468+-0.181
|  MiDA[17] |  0.55+-0.109 |  0.997+-0.001 |  0.957+-0.000 |  0.958+-0.000 |  2.673+-0.173 |  5.842+-0.086
|  MPPN |  0.849+-0.001 |  0.997+-0.000 |  0.957+-0.000 |  0.958+-0.000 |  3.018+-0.849 |  6.495+-0.909
MobIS |  Evermann[6] |  0.767+-0.140 |  0.163+-0.000 |  0.798+-0.000 |  0.075+-0.000 |  — |  —
|  Ca_Spez.[3] |  0.87+-0.040 |  0.163+-0.000 |  0.798+-0.000 |  0.075+-0.000 |  4.648+-0.560 |  30.106+-0.814
|  Ca_concat[3] |  0.836+-0.034 |  0.163+-0.000 |  0.798+-0.000 |  0.075+-0.000 |  4.801+-0.525 |  30.133+-0.526
|  Ca_full[3] |  0.838+-0.038 |  0.163+-0.000 |  0.798+-0.000 |  0.075+-0.000 |  3.966+-0.922 |  24.449+-0.354
|  Tax_Spez.[23] |  0.85+-0.079 |  — |  0.798+-0.000 |  — |  3.919+-0.968 |  28.236+-1.569
|  Tax_Mixed[23] |  0.545+-0.188 |  — |  0.798+-0.000 |  — |  2.333+-0.602 |  21.384+-0.977
|  Tax_Shared[23] |  0.926+-0.008 |  — |  0.805+-0.009 |  — |  2.323+-0.638 |  20.963+-0.420
|  MiDA[17] |  0.7+-0.154 |  0.163+-0.000 |  0.798+-0.000 |  0.075+-0.003 |  2.992+-0.372 |  24.498+-0.405
|  MPPN |  0.934+-0.003 |  0.536+-0.026 |  0.812+-0.002 |  0.121+-0.023 |  4.827+-0.420 |  22.454+-1.011

Thus, we deviate from the original work in some aspects regarding train-test
splitting, sequence generation, and pre-processing, which also leads to
different prediction results. Unfortunately, we cannot guarantee that we have
correctly reproduced all the details of the specifications of the models, due
to missing source code, documentation or test data.

####  5.2.5 Interpretation of Results

For the final comparison, we averaged the prediction scores over ten runs.
Table 3 presents the final results. There is no superior model that performs
best in all tasks on all datasets. However, the results suggest the
effectiveness of the MPPN. This yields in particular for the $NSP_{activity}$
and the $NSP_{resource}$ tasks, where it achieves the highest scores on nearly
all datasets. The MPPN also performs well on the $OUT_{activity}$ and the
$OUT_{resource}$ tasks. However, there is not such a wide performance variety
between the models. Most of the examined processes only have a few outcome
classes. Therefore, the tasks are supposed to be simpler and lead to similar
results. At the same time, the available information in the prefixes may not
always allow for a adequate prediction. For the two regression tasks, the MPPN
achieves solid but no outstanding results. Overall, the results suggest that
the MPPN model is more robust than the other models and does not require
extensive hyperparameter tuning. One explanation might be that the MPPN
utilizes gramian angular fields in combination with CNNs instead of embeddings
and recurrent layers. Also, the CNN in MPPN is based on the Alexnet
architecture, which has been carefully optimized for image recognition tasks.
[23, 3, 17] utilize multi-task learning without fine-tuning, which seem to
fail occasionally to optimize one particular task fully. In contrast, through
the fine-tuning step of the MPPN, it can focus on one task at a time.
Additionally, the MPPN performs reasonable overall tasks and datasets which is
a strong indicator that it can learn effective, general representations of the
underlying process. Another interesting aspect is the influence of the
different perspectives on the process predictions. The MPPN and the MiDA model
utilized almost all available perspectives, while the other models only
examined activity, resource, and timestamp. In the datasets containing
contextual attributes, the MPPN can often outperform other methods indicating
that the model can make use of the additional information and embed them into
the representation. In the future, we plan to further investigate the
influence of different datasets and subsets of perspectives. For example, in
the case of BPI17, we expect that contextual information such as application
type and event origin can positively affect the prediction quality.

####  5.2.6 Reproducibility

All code used for this paper, including the implementation of MPPN as well as
the case retrieval and the prediction experiments, can be found in our git
repository888http://bit.do/fQRbF.

##  6 Conclusion

In this work, we have proposed a novel approach for multivariate business
process representation learning utilizing gramian angular fields and
convolutional neural networks. MPPN is a generic method that generates multi-
purpose vector representations by exposing important characteristics of
process instances without the need for manually labeled data. We showed how
these representations can be exploited for analytics tasks such as clustering
and case retrieval. Furthermore, our work demonstrated the advantages of
meaningful, general representations for later downstream tasks such as next
step and outcome prediction. In the performed experiments, we were able to
outperform existing approaches and generate robust results over several
datasets and tasks. This demonstrates that representation learning can
successfully be applied on business process data. Furthermore, the self-
supervised pre-training makes the model robust and helps in cases where
contextual information is given. Additionally, in spite of recent advances in
NLP, our result indicate that a non-recurrent neural network outperforms other
architectures that use recurrent layers.

One limitation of this paper is a missing systematic hyper-parameter tuning.
In this paper, we investigated the robustness of the models on multiple
datasets and tasks making it a generic approach. In the future, we want to
elaborate on how hyper-parameter tuning can improve the performance of a
specific model on a given dataset. Furthermore, we plan to investigate how the
approach can explain the impact of certain attributes on other events in a
process. The ”black box” nature of deep learning models is still a major issue
in the context of predictive process analysis. Last, we want to elaborate more
approaches and ideas from other domains such as natural language processing
and computer vision to learn richer representations capturing more and finer
characteristics.

## References

  * [1] Bengio, Y., Courville, A., Vincent, P.: Representation learning: A review and new perspectives. IEEE transactions on pattern analysis and machine intelligence pp. 1798–1828 (2013)
  * [2] Brunk, J., Stottmeister, J., Weinzierl, S., Matzner, M., Becker, J.: Exploring the effect of context information on deep learning business process predictions. Journal of Decision Systems (2020)
  * [3] Camargo, M., Dumas, M., Gonzá, lez-Rojas, O.: Learning accurate lstm models of business processes. In: Business Process Management. pp. 286–302. Springer International Publishing (2019)
  * [4] De Koninck, P., vanden Broucke, S., De Weerdt, J.: act2vec, trace2vec, log2vec, and model2vec: Representation learning for business processes. In: Business Process Management. pp. 305–321. Springer International Publishing (2018)
  * [5] Devlin, J., Chang, M.W., Lee, K., Toutanova, K.: Bert: Pre-training of deep bidirectional transformers for language understanding. arXiv preprint arXiv:1810.04805 (2018)
  * [6] Evermann, J., Rehse, J.R., Fettke, P.: Predicting process behaviour using deep learning. Decision Support Systems 100, 129–140 (2017)
  * [7] Houy, C., Rehse, J.R., Scheid, M., Fettke, P.: Model-based compliance in information systems - foundations, case description and data set of the mobis-challenge for students and doctoral candidates. In: Tagungsband der Internationalen Tagung Wirtschaftsinformatik 2019. February 24-27, Siegen, Germany. pp. 2026–2039. Universität Siegen
  * [8] Krizhevsky, A., Sutskever, I., Hinton, G.E.: Imagenet classification with deep convolutional neural networks. NIPS 25, 1097–1105 (2012)
  * [9] Maggi, F.M., Di Francescomarino, C., Dumas, M., Ghidini, C.: Predictive monitoring of business processes. In: CAiSE. pp. 457–472. Springer (2014)
  * [10] Mehdiyev, N., Evermann, J., Fettke, P.: A novel business process prediction model using a deep learning method. Business and Information Systems Engineering 62(2), 143–157 (2020)
  * [11] Mehdiyev, N., Fettke, P.: Local post-hoc explanations for predictive process monitoring in manufacturing. arXiv preprint arXiv:2009.10513 (2020)
  * [12] Moon, J., Park, G., Jeong, J.: Pop-on: Prediction of process using one-way language model based on nlp approach. Applied Sciences 11(2) (2021)
  * [13] Márquez-Chamorro, A.E., Resinas, M., Ruiz-Cortés, A.: Predictive monitoring of business processes: A survey. IEEE Transactions on Services Computing pp. 962–977 (2018)
  * [14] Navarin, N., Vincenzi, B., Polato, M., Sperduti, A.: Lstm networks for data-aware remaining time prediction of business process instances. In: 2017 IEEE Symposium Series on Computational Intelligence (SSCI). pp. 1–7. IEEE (2017)
  * [15] Neu, D.A., Lahann, J., Fettke, P.: A systematic literature review on state-of-the-art deep learning methods for process prediction. Artificial Intelligence Review (2021)
  * [16] Nolle, T., Seeliger, A., Mühlhäuser, M.: Binet: Multivariate business process anomaly detection using deep learning. In: Business Process Management. pp. 271–287 (2018)
  * [17] Pasquadibisceglie, V., Appice, A., Castellano, G., Malerba, D.: A multi-view deep learning approach for predictive business process monitoring. IEEE Transactions on Services Computing (2021)
  * [18] Pasquadibisceglie, V., Appice, A., Castellano, G., Malerba, D.: Predictive process mining meets computer vision. In: International Conference on Business Process Management. pp. 176–192. Springer (2020)
  * [19] Pasquadibisceglie, V., Appice, A., Castellano, G., Malerba, D., Modugno, G.: Orange: Outcome-oriented predictive process monitoring based on image encoding and cnns. IEEE Access 8, 184073–184086 (2020)
  * [20] Pauwels, S., Calders, T.: Detecting anomalies in hybrid business process logs. ACM SIGAPP Applied Computing Review 19, 18–30 (2019)
  * [21] Smith, L.N.: A disciplined approach to neural network hyper-parameters: Part 1–learning rate, batch size, momentum, and weight decay. arXiv preprint arXiv:1803.09820 (2018)
  * [22] Su, H., Maji, S., Kalogerakis, E., Learned-Miller, E.: Multi-view convolutional neural networks for 3d shape recognition. In: 2015 IEEE International Conference on Computer Vision (ICCV). pp. 945–953 (2015)
  * [23] Tax, N., Verenich, I., La Rosa, M., Dumas, M.: Predictive business process monitoring with lstm neural networks. In: International Conference on Advanced Information Systems Engineering. pp. 477–492. Springer (2017)
  * [24] Taymouri, F., La Rosa, M., Erfani, S., Bozorgi, Z.D., Verenich, I.: Predictive business process monitoring via generative adversarial nets: The case of next event prediction. In: International Conference on Business Process Management. pp. 237–256. Springer (2020)
  * [25] Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A.N., Kaiser, L., Polosukhin, I.: Attention is all you need. arXiv preprint arXiv:1706.03762 (2017)
  * [26] Wang, Z., Oates, T.: Encoding time series as images for visual inspection and classification using tiled convolutional neural networks. In: Workshop at the Twenty-Ninth AAAI conference on Artificial Intelligence (2014)


# On the Penalty term for the Mixed Discontinuous Galerkin Finite Element
Method for the Biharmonic Equation

Balaje Kalyanaraman
School of Mathematical and Physical Sciences,
University of Newcastle, Callaghan, NSW 2300, Australia.
and
Danumjaya P.
Birla Institute of Technology and Science, Pilani, K.K.Birla Goa Campus,
NH17B, Zuarinagar, South Goa, Goa 403726, India Corresponding Author. Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we present a study on the effect of penalty term in the mixed
Discontinuous Galerkin Finite Element Method for the biharmonic equation
proposed by Gudi et al. (2008). The proposed mixed Discontinuous Galerkin
Method showed sub-optimal rates of convergence for piecewise quadratic
elements and no significant convergence rates for piecewise linear elements.
We show that by choosing the penalty term proportional to $|e_{k}|^{-1}$
instead of $|e_{k}|^{-3}$, ensures an optimal rate of convergence for the
approximation, including for piecewise linear elements. Finally, we present
numerical experiments to validate our theoretical results.
Keywords:
Finite Elements; Discontinuous Galerkin Finite Element Method; Biharmonic
problem; Optimal Error Estimates.

##  1 Introduction

Discontinuous Galerkin methods are popular finite element techniques which use
discontinuous polynomials to construct approximate solutions. The local nature
of approximation offers flexibility in using higher order polynomials with
non–uniformity in the degree of approximation and in adaptive methods. For a
review of the various discontinuous Galerkin methods we refer the reader to
Arnold et al. (2002); Cockburn et al. (2012); Riviére (2008). For fourth order
problems, conforming methods require imposition of $C^{1}$–continuity across
the inter element boundaries and are computationally expensive. Several finite
element methods have been proposed including the mixed finite element method
(Ciarlet & Raviart, 1974; Monk, 1987; Danumjaya & Pani, 2012),
$C^{0}$–interior penalty methods (Brenner & Sung, 2005; Engel et al., 2002)
and so on to relax the continuity requirements. The idea behind the mixed
finite element method is to split the fourth order problem into two second
order problems by introducing an auxiliary variable $v=-\Delta u$. The
$C^{1}$–continuity requirement is relaxed and the system is then approximated
using a $C^{0}$-finite element method.

An hp mixed discontinuous Galerkin Finite Element Method was proposed by Gudi
et al. (2008). Using a primal formulation leads to integrals involving higher
order derivatives which can be avoided by using a mixed formulation. However,
the method yielded sub-optimal convergence rates for piecewise quadratic
elements and no significant convergence rates for piecewise linear elements.
This is due to the choice of the penalty parameter in the weak formulation
which was taken proportional to the inverse cube $(|e_{k}|^{-3})$ of edge/face
diameter $e_{k}$, a common choice in super-penalization in non-symmetric
interior penalty Galerkin methods (Riviére, 2008). However, modifying the
penalty parameter by taking it proportional to the inverse of the edge length
$(|e_{k}|^{-1})$, essentially reducing the size of the penalty term, yields
optimal convergence.

In this paper, we consider the mixed formulation of the biharmonic equation

$\Delta^{2}u=f\quad\text{in}\;\;\Omega,$ (1)

subject to the clamped boundary conditions

$u=\frac{\partial u}{\partial n}=0\quad\text{on}\;\;\partial\Omega.$

where $\Omega\subset\mathbb{R}^{n}$, $n=2,3$ is a bounded and convex domain
with smooth boundary $\partial\Omega$ and $n$ is the outward unit normal to
$\partial\Omega$. We assume that the data $f$ is sufficiently smooth so that
there exists a unique solution to the problem in $H^{4}(\Omega)$. We introduce
a new variable $v=-\Delta u$ and split the biharmonic equation (1) into two
equations as

$\displaystyle-\Delta v$ $\displaystyle=f\quad\text{in}\;\;\Omega,$ (2)
$\displaystyle-\Delta u$ $\displaystyle=v\quad\text{in}\;\;\Omega,$
$\displaystyle u$ $\displaystyle=\frac{\partial u}{\partial
n}=0\quad\text{on}\;\;\Gamma.$

We then discretize the problem using the mixed discontiunous Galerkin method
and prove some error estimates for the finite element solution to study the
convergence of the method with respect to the mesh size $h$. We then perform
some numerical experiments to validate the theoretical results.

The paper is organized as follows. In Section 2 we derive the weak formulation
of the biharmonic problem. The proof of well–posedness of the problem is
identical to Gudi et al. (2008) and hence we omit the same. In Section 3, we
discuss the error estimates for $h-$refinement and we derive $L^{2}$ and
energy estimates for the primal variable and an error estimate in the
$L^{2}$-norm for the auxiliary variable. Finally, in Section 4 we perform
numerical experiments to validate our theoretical results established in
Section 3.

##  2 Weak Formulation

Let $\bar{\Omega}=\underset{K\in\mathcal{T}_{h}}{\cup}\bar{K}$ be a uniform
partition of $\Omega$. Let us denote the edges of $\mathcal{T}_{h}$ by $e_{k}$
where $k=1,\cdots,M_{h}$. Let $\Gamma_{\text{I}}=\left\\{e_{1},e_{2},\cdots
e_{P_{h}}\right\\}$ denote the set of all interior edges and
$\Gamma_{D}=\left\\{e_{P_{h}+1},\cdots,e_{M_{h}}\right\\}$, the boundary
edges. Let $\Gamma=\Gamma_{\text{I}}\cup\Gamma_{\text{D}}$. For each edge
$e_{k}=K_{i}\cap K_{j}$, associate a unit normal $n_{k}$ outward from $K_{i}$,
so that $n|_{K_{i}}=-n|_{K_{j}}$. Let $h_{K}$ denote the diameter of each
$K\in\mathcal{T}_{h}$ and $h=\underset{K\in\mathcal{T}_{h}}{\max}\,h_{K}$. We
define the broken Sobolev space

$H^{s}(\mathcal{T}_{h})=\left\\{v\in L^{2}(\Omega)\;:\;v_{K_{i}}\in
H^{s}(K)\;\,\forall K\in\mathcal{T}_{h}\right\\}.$ (3)

where $H^{s}(K)$ denotes the standard Sobolev space of order $s$. The
associated broken norm and semi–norm are defined respectively, by

$\|v\|_{s,\mathcal{T}_{h}}=\left(\sum_{K\in\mathcal{T}_{h}}\|v\|_{s,K}^{2}\right)^{1/2}\quad\text{and}\quad|v|_{s,\mathcal{T}_{h}}=\left(\sum_{K\in\mathcal{T}_{h}}|v|_{s,K}^{2}\right)^{1/2}$

where $\|v\|_{s,K}$ and $|v|_{s,K}$ denotes the standard Sobolev norm and
seminorm on $K$, respectively. We denote the $L^{2}$ norm by $\|\cdot\|$ and
its inner–product by $(\cdot,\cdot)$. We set $V=H^{2}(\mathcal{T}_{h})$.

Let $u,v$ be sufficiently smooth functions. We consider the first equation in
(2), multiply by some $\varphi\in V$, and integrate over the domain to obtain

$-\int_{\Omega}(\Delta v)\varphi\,dx=\int_{\Omega}f\varphi\,dx$ (4)

Applying integration by parts and using the fact that $[v]=0$ in
$\Gamma_{\text{I}}$, we obtain

$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla
v\cdot\nabla\varphi\,dx-\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\nabla v\cdot
n_{k}\right\\}[\varphi]\,ds$
$\displaystyle-\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\nabla\varphi\cdot
n_{k}\right\\}[v]\,ds$ (5)
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\int_{K}f\varphi\,dx.$

Since $[u]=0$ in $\Gamma$, we see that

$J(u,\varphi)=\sum_{e_{k}\in\Gamma}\int_{e_{k}}\alpha_{k}[u][\varphi]\,ds=0$
(6)

where $\alpha_{k}$ is a positive real constant. Adding $J(u,\varphi)$ to (5)
we obtain

$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla
v\cdot\nabla\varphi\,dx$
$\displaystyle-\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\nabla v\cdot
n_{k}\right\\}[\varphi]\,ds-\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\nabla\varphi\cdot
n_{k}\right\\}[v]\,ds$ (7)
$\displaystyle+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\alpha_{k}[u][\varphi]\,ds=\sum_{K\in\mathcal{T}_{h}}\int_{K}f\varphi\,dx,$

which yields the weak formulation of the first equation. Similarly for the
second equation, multiplying some $\chi\in V$, integrating over $\Omega$ and
using the fact that $\nabla u\cdot n=0$ on $\Gamma_{\text{D}}$, $[u]=0$ in
$\Gamma$, we obtain

$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla u\cdot\nabla\chi\,dx$
$\displaystyle-\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\nabla
u\cdot
n_{k}\right\\}[\chi]\,ds-\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\nabla\chi\cdot
n_{k}\right\\}[u]\,ds=\sum_{K\in\mathcal{T}_{h}}\int_{K}v\chi\,dx,$ (8)

We define the bilinear forms

$\displaystyle B(w,z)$ $\displaystyle=\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla
w\cdot\nabla z\,dx-\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\nabla w\cdot
n_{k}\right\\}[z]\,ds-\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\nabla
z\cdot n_{k}\right\\}[w]\,ds$ $\displaystyle B_{\text{I}}(w,z)$
$\displaystyle=\sum_{K\in\mathcal{T}_{h}}\int_{K}\nabla w\cdot\nabla
z\,dx-\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\nabla w\cdot
n_{k}\right\\}[z]\,ds-\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\nabla
z\cdot n_{k}\right\\}[w]\,ds$

Now the weak formulation of the problem is as follows: Find $(u,v)\in V\times
V$, such that

$\displaystyle B(v,\varphi)+J(u,\varphi)$ $\displaystyle=(f,\varphi),$ (9)
$\displaystyle B(\chi,u)$ $\displaystyle=(v,\chi),$

for all $(\varphi,\chi)\in V\times V$. Define a finite dimensional subspace
$V_{h}$ of $V$ as

$V_{h}=\left\\{v\in
H^{s}(\mathcal{T}_{h})\;:\;v|_{K}\in\mathbb{P}_{p}(K)\;\;\text{for
all}\;\;K\in\mathcal{T}_{h}\right\\},$

where $\mathbb{P}_{p}(K)$ denotes the space of polynomials of degree $\leq p$.
Now the discontinuous Galerkin finite element weak formulation reads: Find
$(u_{h},v_{h})\in V_{h}\times V_{h}$ such that

$\displaystyle B(v_{h},\varphi)+J(u_{h},\varphi)$ $\displaystyle=(f,\varphi),$
(10) $\displaystyle B(\chi,u_{h})$ $\displaystyle=(v_{h},\chi),$

for all $(\varphi,\chi)\in V_{h}\times V_{h}$. For the purpose of error
analysis, we define the following mesh dependent energy norms,

$\displaystyle|||w|||^{2}$
$\displaystyle=\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla
w|^{2}\,dx+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{|e_{k}|}{\sigma_{1}}\left\\{\frac{\partial
w}{\partial
n_{k}}\right\\}^{2}\,ds+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[w]^{2}\,ds\right),$
$\displaystyle|||w|||_{\text{I}}^{2}$
$\displaystyle=\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla
w|^{2}\,dx+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{|e_{k}|}{\sigma_{1}}\left\\{\frac{\partial
w}{\partial
n_{k}}\right\\}^{2}\,ds+\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[w]^{2}\,ds\right).$

Below, we state some properties of the bilinear form $B(\cdot,\cdot)$ without
proof.

######  Lemma 2.1.

For sufficiently large constant $\sigma_{1}$, it can be shown that for any
$w_{h}\in V_{h}$,

$C|||w_{h}|||_{\text{I}}^{2}\leq
B_{\text{I}}(w_{h},w_{h})+\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[w_{h}]^{2}\,ds,$
(11)

where $C$ is a constant independent of $h$ and $p$. It can also be shown that
for all $w,q\in V$ there exists a constant independent of $h$ such that

$|B(w,q)|\leq C|||w|||_{\text{I}}\;|||q|||.$ (12)

The inequalities (11) and (12) refer to the coercivity and the boundedness
property of the bilinear forms. For proofs, we refer the reader to Riviére
(2008); Gudi et al. (2008) and the references therein. Below we define an
interpolation operator for functions on the Broken Sobolev space.

######  Definition 2.1.

For given $\phi\in H^{s}(\mathcal{T}_{h})$, define $I_{h}\phi\in V_{h}$ be an
interpolation operator satisfying optimal approximation properties, i.e.,

$|||\phi-I_{h}\phi|||\leq Ch^{\mu-1}\|\phi\|_{s,\mathcal{T}_{h}}$

where $\mu=\min\\{p+1,s\\}$.

We recall the following trace inequality on the finite element space.

######  Lemma 2.2.

Let $v_{h}\in V_{h}$. There exists a constant $C>0$ such that

$\|\nabla^{l}v_{h}\|_{L^{2}(e_{k})}\leq
Ch_{K}^{-1/2}\|\nabla^{l}v_{h}\|_{L^{2}(K)},\quad l=0,1.$ (13)

Below, we state the inverse inequality without proof.

######  Lemma 2.3.

Let $v_{h}\in V_{h}$. There exists a constant $C>0$ such that

$|v_{h}|_{H^{1}(K)}\leq Ch_{K}^{-1}\|v_{h}\|_{L^{2}(K)},$ (14)

For proofs of the trace inequality and the inverse inequality, we refer the
reader to Brenner & Scott (2007) and Riviére (2008). We state and prove the
following Lemma.

######  Lemma 2.4.

For $w_{h}\in V_{h}$, there exists a positive constant $C>0$, such that

$|||w_{h}|||^{2}\leq C_{T}\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla
w_{h}|^{2}+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[w_{h}]^{2}\,ds\right).$

###### Proof.

From the definition of the energy norm and the trace inequality, we have for
$w_{h}\in V_{h}$

$\displaystyle|||w_{h}|||^{2}$
$\displaystyle=\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla
w_{h}|^{2}\,dx+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{|e_{k}|}{\sigma_{1}}\left\\{\frac{\partial
w_{h}}{\partial
n_{k}}\right\\}^{2}\,ds+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[w_{h}]^{2}\,ds\right),$
$\displaystyle\leq C_{T}\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla
w_{h}|^{2}\,dx+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[w_{h}]^{2}\,ds\right),$

which yields the desired result. ∎

We are interested to study the effect of the penalty parameter on the
convergence of the discrete solution. Let us set

$\alpha_{k}=\sigma_{0}|e_{k}|^{-i}p^{2}.$ (15)

where $p$ is the degree of approximation of the polynomial, $\sigma_{0}$ is a
constant which will be defined later and $i$ is an integer which describes the
degree of penalization. It is well known in literature that for the case of
the Non–symmetric interior penalty Galerkin (NIPG) method for second order
elliptic problems, the choice of the penalty parameter plays a crucial role on
the convergence of the solution. If $i=1$, under normal penalization, the
convergence in the $L^{2}-$norm for the piecewise quadratic case is
suboptimal. However, for the case when $i=3$, under super penalization, the
convergence is optimal. For the mixed Discontinuous Galerkin Method considered
here, we observe that the penalization term with $i=1$ produces optimal
convergence in all cases, as opposed to $i=3$ where optimal convergence rates
are observed only for piecewise cubic elements. The latter case was well
studied by Gudi et al. (2008). In the next section, we derive optimal error
estimates with $i=1$.

##  3 Error Analysis

Define the auxiliary projection $\Pi_{h}:H^{s}(\mathcal{T}_{h})\to V_{h}$ by

$\displaystyle B_{\text{I}}(\phi-\Pi_{h}\phi,\chi)$
$\displaystyle+\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\phi-\Pi_{h}\phi][\chi]\,ds=0,$
(16) $\displaystyle\left(\phi-\Pi_{h}\phi,1\right)$ $\displaystyle=0.$

The auxiliary projection can be shown to satisfy optimal error estimates in
the $L^{2}$ and energy norms. We define $\mu=\min\\{p+1,s\\}$,
$\bar{\mu}=\min\\{p+1,\bar{s}\\}$, $\theta=\min\\{p+1,4\\}$. We state and
prove the following Lemma.

######  Lemma 3.1.

For a given $\phi\in V$, there exists a unique $\Pi_{h}\phi\in V_{h}$
satisfying (16). Moreover, there exists positive constants independent of $h$
such that

$\displaystyle|||\phi-\Pi_{h}\phi|||$ $\displaystyle\leq
Ch^{\mu-1}\,\|\phi\|_{s,\mathcal{T}_{h}}$ (17)
$\displaystyle\|\phi-\Pi_{h}\phi\|$ $\displaystyle\leq
Ch^{\mu}\,\|\phi\|_{s,\mathcal{T}_{h}}.$ (18)

provided $\sigma_{1}$ is a sufficiently large positive constant, i.e.,
$\sigma_{1}\geq\sigma_{1}^{*}$ for some $\sigma_{1}^{*}\in\mathbb{R}^{+}$.

###### Proof.

For the proof, we refer the reader to Theorem 5.1 in Gudi et al. (2008). ∎

Subtracting (10) from (9), we obtain the error equations

$\displaystyle B(v-v_{h},\varphi)+J(u-u_{h},\varphi)$ $\displaystyle=0,$ (19)
$\displaystyle B(\chi,u-u_{h})-(v-v_{h},\chi)$ $\displaystyle=0.$ (20)

for $(\varphi,\chi)\in V_{h}\times V_{h}$. Set $e_{u}:=u-u_{h}$ and
$e_{v}:=v-v_{h}$. The following Lemma is useful in proving the error
estimates.

######  Lemma 3.2.

Let $\xi_{v}=v_{h}-\Pi_{h}v$. Provided $\sigma_{0}>\sigma_{*}$, there exists a
positive constant $C_{*}$ such that

$\displaystyle C_{*}|||\xi_{v}|||^{2}\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}^{2}+\frac{\sigma_{0}\epsilon^{2}}{2C_{1}}J(\xi_{u},\xi_{u})$

where $C_{1}$ is a positive constant independent of $\sigma_{0}$ and
$0<\epsilon\ll 1$.

###### Proof.

Let $\eta_{v}=v-\Pi_{h}v$. From equation (19) we obtain,

$B(\xi_{v},\varphi)=B(\eta_{v},\varphi)+J(\eta_{u},\varphi)-J(\xi_{u},\varphi),$

where $\eta_{u}=u-I_{h}u$ and $\xi_{u}=u_{h}-I_{h}u$. Setting
$\varphi=\xi_{v}$ and expanding the bilinear form $B(\cdot,\cdot)$, we obtain

$\displaystyle\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla\xi_{v}|^{2}\,dx$
$\displaystyle=\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\frac{\partial\xi_{v}}{\partial
n_{k}}\right\\}[\xi_{v}]\,ds+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\frac{\partial\xi_{v}}{\partial
n_{k}}\right\\}[\xi_{v}]\,ds$
$\displaystyle+B(\eta_{v},\xi_{v})+J(\eta_{u},\xi_{v})-J(\xi_{u},\xi_{v}).$
(21)

Using Lemma 2.4 with $w_{h}=\xi_{v}$, and substituting the result of equation
(21) onto the Lemma, we obtain the inequality

$\displaystyle c|||\xi_{v}|||^{2}$
$\displaystyle\leq\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\frac{\partial\xi_{v}}{\partial
n_{k}}\right\\}[\xi_{v}]\,ds+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\frac{\partial\xi_{v}}{\partial
n_{k}}\right\\}[\xi_{v}]\,ds$
$\displaystyle+B(\eta_{v},\xi_{v})+J(\eta_{u},\xi_{v})-J(\xi_{u},\xi_{v})+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\xi_{v}]^{2}\,ds,$

with $c=\frac{1}{C_{T}}$ being a constant which arises due to the trace
inequality and is dependent on the triangulation. We observe that using the
boundedness properties of the bilinear forms together with the trace and
Young’s inequalities, we have

$\displaystyle B(\eta_{v},\xi_{v})$ $\displaystyle\leq
C|||\eta_{v}|||\;|||\xi_{v}|||\leq
Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}^{2}+\frac{\epsilon^{2}}{4}|||\xi_{v}|||^{2},$
$\displaystyle J(\eta_{u},\xi_{v})$
$\displaystyle\leq(J(\eta_{u},\eta_{u}))^{1/2}(J(\xi_{v},\xi_{v}))^{1/2}\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+\frac{\epsilon^{2}}{4}|||\xi_{v}|||^{2},$

Using the Cauchy-Schwarz’s inequality and subsequently the Young’s inequality,
the trace inequality and the inverse inequality, the term,

$\displaystyle\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\frac{\partial\xi_{v}}{\partial
n_{k}}\right\\}[\xi_{v}]\,ds$
$\displaystyle\leq\left(\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{|e_{k}|}{\sigma_{0}}\left\\{\frac{\partial\xi_{v}}{\partial
n_{k}}\right\\}\,ds\right)^{1/2}\left(\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{0}}{|e_{k}|}[\xi_{v}]^{2}\,ds\right)^{1/2},$
$\displaystyle\leq\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}\frac{1}{\sigma_{0}}|\nabla\xi_{v}|^{2}\,dx\right)^{1/2}\left(\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{0}}{|e_{k}|}[\xi_{v}]^{2}\,ds\right)^{1/2},$
$\displaystyle\leq\frac{C_{1}}{\sigma_{0}}|||\xi_{v}|||\;|||\xi_{v}|||,$
$\displaystyle\leq\frac{C_{1}}{2\sigma_{0}\epsilon^{2}}|||\xi_{v}|||^{2}+\frac{\epsilon^{2}}{2}|||\xi_{v}|||^{2},$

where $C_{1}$ is a constant independent of $\sigma_{0}$ and depends on the
triangulation. Similarly, we have

$J(\xi_{u},\xi_{v})\leq\frac{\sigma_{0}\epsilon^{2}}{2C_{1}}J(\xi_{u},\xi_{u})+\frac{C_{1}}{2\sigma_{0}\epsilon^{2}}|||\xi_{v}|||^{2}.$

Finally we have the bound,

$\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\xi_{v}]^{2}\,ds\leq\frac{\sigma_{1}}{\sigma_{0}}|||\xi_{v}|||^{2}.$

Combining all the bounds for a small $\epsilon$, we obtain the intermediate
inequality,

$\displaystyle\left(c-\frac{C_{1}}{\sigma_{0}\epsilon^{2}}-\frac{\sigma_{1}}{\sigma_{0}}\right)|||\xi_{v}|||^{2}\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}+\frac{\sigma_{0}\epsilon^{2}}{2C_{1}}J(\xi_{u},\xi_{u})$

A condition that

$\displaystyle
C_{*}=\left(c-\frac{C_{1}}{\sigma_{0}\epsilon^{2}}-\frac{\sigma_{1}}{\sigma_{0}}\right)>0\quad\text{or}\quad\sigma_{0}>\frac{C_{1}+\sigma_{1}\epsilon^{2}}{c\epsilon^{2}},$

yields a lower bound on $\sigma_{0}$. Thus,

$\displaystyle C_{*}|||\xi_{v}|||^{2}\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}+\frac{\sigma_{0}\epsilon^{2}}{2C_{1}}J(\xi_{u},\xi_{u}),$

which is the desired estimate with
$\sigma_{*}=\frac{C_{1}+\sigma_{1}\epsilon^{2}}{c\epsilon^{2}}$ for a
sufficiently small $\epsilon$. ∎

Lemma 3.2 relates the constant $\sigma_{0}$ with $\sigma_{1}$ by providing a
bound on $\sigma_{0}$ in terms of $\sigma_{1}$ which in turn depends on the
triangulation and is bounded below by $\sigma_{1}^{*}$, which is independently
known from $\sigma_{0}$. Now we state and prove the following theorems which
yields the error estimates for the DG method.

######  Theorem 3.1.

Let $u_{h},v_{h}$ satisfy (10). Then there exists some positive constant $C$
such that

$\displaystyle\|e_{v}\|^{2}+J(e_{u},e_{u})\leq
C\left(h^{2\mu-2}\|u\|_{s,\mathcal{T}}+h^{2\bar{\mu}-2}\|u\|_{\bar{s},\mathcal{T}}\right)$

provided $\sigma_{*}<\sigma_{0}<\sigma^{*}$ for some
$\sigma_{*},\sigma^{*}\in\mathbb{R}^{+}$.

###### Proof.

Let $\delta\ll 1$. First, we split

$\displaystyle u-u_{h}$
$\displaystyle=(u-I_{h}u)-(u_{h}-I_{h}u):=\eta_{u}-\xi_{u},$ $\displaystyle
v-v_{h}$ $\displaystyle=(v-\Pi_{h}v)-(v_{h}-\Pi_{h}v):=\eta_{v}-\xi_{v}.$

Substituting the above definitions into the error equations (19) and (20)
setting $\varphi=\xi_{u}$ in (19) and $\chi=\xi_{v}$ in (20) and subtracting
the resulting equations, we obtain

$\displaystyle\|\xi_{v}\|^{2}+J(\xi_{u},\xi_{u})=B(\eta_{v},\xi_{u})+J(\eta_{u},\xi_{u})-B(\xi_{v},\eta_{u})+(\eta_{v},\xi_{v}).$

Using the definition of the auxiliary projection we obtain,

$\displaystyle B(\eta_{v},\xi_{u})$
$\displaystyle=B_{\text{I}}(\eta_{v},\xi_{u})-\sum_{e_{k}\in\Gamma_{\text{D}}}\int_{e_{k}}\frac{\partial\eta_{v}}{\partial
n_{k}}\xi_{u}\,ds,$
$\displaystyle=-\sum_{e_{k}\in\Gamma_{\text{I}}}\int\frac{\sigma_{1}}{|e_{k}|}[\eta_{v}][\xi_{u}]\,ds-\sum_{e_{k}\in\Gamma_{\text{D}}}\int_{e_{k}}\frac{\partial\eta_{v}}{\partial
n_{k}}\xi_{u}\,ds.$

We observe that using the trace inequality, inverse inequality and the Young’s
inequality,

$\displaystyle B(\eta_{v},\xi_{u})\leq
Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}^{2}+\frac{\delta^{2}}{2}J(\xi_{u},\xi_{u}).$

for some $\delta\ll 1$. Using the Cauchy-Schwarz inequality and the Young’s
inequality, we obtain

$\displaystyle J(\eta_{u},\xi_{u})\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+\frac{\delta^{2}}{2}J(\xi_{u},\xi_{u}).$

We observe that using the Young’s inequality and the result from Lemma 3.2 for
the estimate of $|||\xi_{v}|||$ we obtain,

$\displaystyle B(\xi_{v},\eta_{u})$ $\displaystyle\leq
C|||\xi_{v}|||\;|||\eta_{u}|||,$ $\displaystyle\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+\frac{\delta^{2}C_{*}}{2}|||\xi_{v}|||^{2},$
$\displaystyle\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}^{2}+\frac{\sigma_{0}\epsilon^{2}\delta^{2}}{4C_{1}}J(\xi_{u},\xi_{u}),$

provided $\sigma_{0}>\sigma_{*}$. We note that the constant $C_{*}$ which
appears in front of $|||\eta_{v}|||$ is absorbed into the estimate involving
$h$. Since $C_{*}$ does not depend on $h$, the estimate is unaffected.
Finally, we observe that

$(\eta_{v},\xi_{v})\leq
Ch^{2\bar{\mu}}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}+\frac{\delta^{2}}{2}\|\xi_{v}\|^{2}.$

Combining the bounds and using the triangle inequality, we obtain

$\|\xi_{v}\|^{2}+J(\xi_{u},\xi_{u})\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}+Ch^{2\bar{\mu}-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}+\frac{\sigma_{0}\epsilon^{2}\delta^{2}}{4C_{1}}J(\xi_{u},\xi_{u}).$
(22)

Setting

$\displaystyle
1-\frac{\sigma_{0}\epsilon^{2}\delta^{2}}{4C_{1}}>0\quad\text{or}\quad\sigma_{0}<\frac{4C_{1}}{\epsilon^{2}\delta^{2}},$

yields an upper bound for $\sigma_{0}$ in terms of the constant $C_{1}$. The
desired estimate follows from the triangle inequality and with
$\sigma_{*}=\frac{C_{1}+\sigma_{1}\epsilon^{2}}{c\epsilon^{2}}$ and
$\sigma^{*}=\frac{4C_{1}}{\epsilon^{2}\delta^{2}}$ for sufficiently small
$\delta$. ∎

The consequence of having the condition that the constant $\sigma_{0}$ to be
bounded between $[\sigma_{*},\sigma^{*}]$ may explain why the convergence rate
deteriorates for large values of penalty parameters. Choosing the value of
$\sigma_{0}$ too large tends the value of the penalty terms to be proportional
to $|e_{k}|^{-3}$ and thusumay erode the estimates derived above. With the
above conditions, the estimate in Theorem 3.1 predicts a higher rate of
convergence, $\mu-1$ for $\|e_{v}\|$ than $\mu-2$ which was predicted by Gudi
et al. (2008). This will lead to better estimates for $|||e_{u}|||$ and
$\|e_{v}\|$. We state and prove the following theorems.

######  Theorem 3.2.

Let $u_{h},v_{h}$ satisfy (10). Then there exists some positive constant $C$
such that

$\displaystyle\|e_{u}\|\leq
C\left(h|||e_{u}|||+h^{2}\|e_{v}\|+h^{\theta-1}|||v-\Pi_{h}v|||+h^{\theta-2}\|\Pi_{h}v-v_{h}\|+h^{\theta-1}(J(e_{u},e_{u}))^{1/2}\right).$

###### Proof.

For the construction of the proof, we use the Aubin-Nitsche duality argument.
Consider the dual problem,

$\displaystyle-\Delta\phi$ $\displaystyle=z\quad\text{in}\quad\Omega,$
$\displaystyle-\Delta\psi$ $\displaystyle=\phi\quad\text{in}\quad\Omega,$
$\displaystyle\psi$ $\displaystyle=\frac{\partial\psi}{\partial
n}=0\quad\text{on}\quad\partial\Omega,$

where the functions $\phi$ and $\psi$ satisfy the regularity result

$\displaystyle\|\psi\|_{H^{4}}+\|\phi\|_{H^{2}}\leq C\|z\|.$

Take $z=e_{u}$, multiply the first equation of the dual problem by $e_{u}$ and
integrate over the domain to obtain

$B(\phi,e_{u})+J(e_{u},\psi)=\|e_{u}\|^{2}.$ (23)

where we have used the fact that $[\psi]=0$ on $\Gamma$ and $[\phi]=0$ on
$\Gamma_{\text{I}}$. Next multiply the second equation of the dual problem by
$e_{v}$ and integrate to obtain

$B(e_{v},\psi)-(e_{v},\phi)=0.$ (24)

where we have used $\frac{\partial\psi}{\partial n}=0$ on $\Gamma_{\text{D}}$
and $[\psi]=0$ on $\Gamma$. Adding the equations (23) and (24) we obtain

$\|e_{u}\|^{2}=B(\phi,e_{u})+J(e_{u},\psi)+B(e_{v},\psi)-(e_{v},\phi).$ (25)

Denoting $\eta_{\phi}=\phi-I_{h}\phi$ and $\eta_{\psi}=\psi-I_{h}\psi$, from
the orthogonality result (19) and (20) we have

$\displaystyle B(e_{v},I_{h}\psi)+J(e_{u},I_{h}\psi)$ $\displaystyle=0,$
$\displaystyle B(I_{h}\phi,e_{u})-(e_{v},I_{h}\phi)$ $\displaystyle=0.$

Subtracting the above equations from (25), we obtain the following equation

$\|e_{u}\|^{2}=B(\eta_{\phi},e_{u})+J(e_{u},\eta_{\psi})+B(e_{v},\eta_{\psi})-(e_{v},\eta_{\phi}).$
(26)

Constructing the error bounds and using the regularity of the functions $\phi$
and $\psi$ we obtain,

$\displaystyle B(\eta_{\phi},e_{u})$ $\displaystyle\leq
C|||\eta_{\phi}|||\,|||e_{u}|||\leq Ch|||e_{u}|||\;\|\phi\|_{H^{2}},$
$\displaystyle(e_{v},\eta_{\phi})$ $\displaystyle\leq
Ch^{2}\|e_{v}\|\;\|\phi\|_{H^{2}},$ $\displaystyle B(e_{v},\eta_{\psi})$
$\displaystyle\leq C|||e_{v}|||\;|||\eta_{\psi}|||\leq
Ch^{\theta-1}|||e_{v}|||\;\|\psi\|_{H^{4}},$ $\displaystyle
J(e_{u},\eta_{\psi})$
$\displaystyle\leq(J(\eta_{\psi},\eta_{\psi}))^{1/2}\,(J(e_{u},e_{u}))^{1/2}\leq
Ch^{\theta-1}\;(J(e_{u},e_{u}))^{1/2}\|\psi\|_{H^{4}}.$

Combining all the bounds,

$\displaystyle\|e_{u}\|$ $\displaystyle\leq
C\left(h|||e_{u}|||+h^{2}\|e_{v}\|+h^{\theta-1}|||e_{v}|||+h^{\theta-1}(J(e_{u},e_{u}))^{1/2}\right),$
$\displaystyle\leq
C\left(h|||e_{u}|||+h^{2}\|e_{v}\|+h^{\theta-1}|||v-\Pi_{h}v|||+h^{\theta-2}\|\Pi_{h}v-v_{h}\|+h^{\theta-1}(J(e_{u},e_{u}))^{1/2}\right),$
(27)

which yields the desired estimate. ∎

######  Theorem 3.3.

Let $u_{h},v_{h}$ satisfy (10). There exists some positive constant $C$ such
that

$\displaystyle|||e_{u}|||$ $\displaystyle\leq
Ch^{\min\\{\mu+\theta-3,\mu-1\\}}\|u\|_{s,\mathcal{T}_{h}}+Ch^{\min\\{\bar{\mu}+\theta-3,\bar{\mu}-1\\}}\|u\|_{\bar{s}+2,\mathcal{T}_{h}},$

provided the conditions on $\sigma_{0}$ in Theorem 3.1 holds.

###### Proof.

Let the constant $\epsilon\ll 1$ and $\xi_{u}=u_{h}-I_{h}u$. We split the
error $u-u_{h}=(u-I_{h}u)-(u_{h}-I_{h}u):=\eta_{u}-\xi_{u}$. From Lemma 2.4,
we note that

$\displaystyle c|||\xi_{u}|||^{2}$
$\displaystyle\leq\left(\sum_{K\in\mathcal{T}_{h}}\int_{K}|\nabla\xi_{u}|^{2}\,dx+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\xi_{u}]^{2}\,ds\right)$
$\displaystyle\leq
B(\xi_{u},\xi_{u})+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\frac{\partial\xi_{u}}{\partial
n_{k}}\right\\}[\xi_{u}]\,ds+\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\frac{\partial\xi_{u}}{\partial
n_{k}}\right\\}[\xi_{u}]\,ds+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\xi_{u}]^{2}\,ds$

From (20), we obtain

$B(\chi,\xi_{u})=B(\chi,\eta_{u})-(e_{v},\chi).$

Thus, with $\chi=\xi_{u}$ we obtain the inequality

$\displaystyle c|||\xi_{u}|||^{2}\leq B(\xi_{u},\eta_{u})-(e_{v},\xi_{u})$
$\displaystyle+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\frac{\partial\xi_{u}}{\partial
n_{k}}\right\\}[\xi_{u}]\,ds+\sum_{e_{k}\in\Gamma_{\text{I}}}\int_{e_{k}}\left\\{\frac{\partial\xi_{u}}{\partial
n_{k}}\right\\}[\xi_{u}]\,ds$
$\displaystyle+\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\xi_{u}]^{2}\,ds.$

We construct the following upper bounds for the terms on the right hand side,

$\displaystyle B(\xi_{u},\eta_{u})$ $\displaystyle\leq
C|||\xi_{u}|||\;|||\eta_{u}||||\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+\frac{\epsilon^{2}}{2}|||\xi_{u}|||^{2},$
$\displaystyle(e_{v},\xi_{u})$
$\displaystyle\leq\frac{C}{2\epsilon^{2}}\|e_{v}\|^{2}+\frac{\epsilon^{2}}{2}\|\xi_{u}\|^{2},$
$\displaystyle\sum_{e_{k}\in\Gamma}\int_{e_{k}}\left\\{\frac{\partial\xi_{u}}{\partial
n_{k}}\right\\}[\xi_{u}]\,ds$
$\displaystyle\leq\left(\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{|e_{k}|}{\sigma_{0}}\left\\{\frac{\partial\xi_{u}}{\partial
n_{k}}\right\\}^{2}\,ds\right)^{1/2}\,\left(\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{0}}{|e_{k}|}[\xi_{u}]^{2}\,ds\right)^{1/2},$
$\displaystyle\leq\frac{C}{2\epsilon^{2}}J(\xi_{u},\xi_{u})+\frac{\epsilon^{2}}{2}|||\xi_{u}|||^{2}.$

Fixing $\sigma_{0}$ and noting that $k_{1}\sigma_{0}\leq\sigma_{1}\leq
k_{2}\sigma_{0}$, we have the final upper bound,

$\sum_{e_{k}\in\Gamma}\int_{e_{k}}\frac{\sigma_{1}}{|e_{k}|}[\xi_{u}]^{2}\,ds\leq
k_{2}J(\xi_{u},\xi_{u}).$

Combining all the bounds and applying the triangle inequality, we obtain

$\displaystyle c|||\xi_{u}|||^{2}\leq
Ch^{2\mu-2}\|u\|_{s,\mathcal{T}_{h}}^{2}+C\left(\|e_{v}\|^{2}+J(\xi_{u},\xi_{u})\right)+\frac{\epsilon^{2}}{2}\|\xi_{u}\|^{2}.$
(28)

Writing

$\displaystyle e_{v}=v-v_{h}=v-\Pi_{h}v-(v_{h}-\Pi_{h}v),$

we obtain $\Pi_{h}v-v_{h}=(v-\Pi_{h}v)-e_{v}$. From estimate (27), we have

$\displaystyle\epsilon\|e_{u}\|$ $\displaystyle\leq
C\epsilon\left(h|||e_{u}|||+h^{2}\|e_{v}\|+h^{\theta-1}|||\Pi_{h}v-v|||+h^{\theta-2}\|\Pi_{h}v-v_{h}\|+h^{\theta-1}(J(e_{u},e_{u}))^{1/2}\right),$
$\displaystyle\leq
C\epsilon\left(h|||e_{u}|||+h^{2}\|e_{v}\|+h^{\theta-1}|||v-\Pi_{h}v|||+h^{\theta-2}\|v-\Pi_{h}v\|+h^{\theta-2}\left(\|e_{v}\|+J(e_{u},e_{u})^{1/2}\right)\right).$

Now we have from the results of Theorem 3.1 and Lemma 3.1 the estimates,

$\displaystyle
h^{2}\|e_{v}\|+h^{\theta-2}\left(\|e_{v}\|+J(e_{u},e_{u})^{1/2}\right)$
$\displaystyle\leq
Ch^{\min\\{\theta-2,2\\}}\left(h^{\mu-1}\|u\|_{s,\mathcal{T}_{h}}+h^{\bar{\mu}-1}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}\right),$
$\displaystyle h^{\theta-2}\|v-\Pi_{h}v\|$ $\displaystyle\leq
Ch^{\bar{\mu}+\theta-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}},$ $\displaystyle
h^{\theta-1}|||v-\Pi_{h}v|||$ $\displaystyle\leq
Ch^{\bar{\mu}+\theta-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}},$

respectively. Thus the estimate for $e_{u}$ reads,

$\displaystyle\epsilon\|e_{u}\|\leq
C\epsilon\left(h|||e_{u}|||+h^{\bar{\mu}+\theta-2}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}+h^{\min\\{\theta-2,2\\}}\left[h^{\mu-1}\|u\|_{s,\mathcal{T}_{h}}+h^{\bar{\mu}-1}\|u\|_{\bar{s}+2,\mathcal{T}_{h}}\right]\right).$
(29)

Combining estimates (28) and (29) and using the triangle inequality, we obtain
the desired result for $|||e_{u}|||$. ∎

Finally using Theorem 3.2 and Theorem 3.3, we can prove the following estimate
for $\|e_{u}\|$.

######  Theorem 3.4.

Let $u_{h},v_{h}$ satisfy (10). There exists some positive constant $C$ such
that

$\displaystyle\|e_{u}\|\leq
Ch^{\min\\{\mu+\theta-3,\mu\\}}\|u\|_{s,\mathcal{T}_{h}}+Ch^{\min\\{\bar{\mu}+\theta-3,\bar{\mu}\\}}\|u\|_{\bar{s}+2,\mathcal{T}_{h}},$

provided the conditions on $\sigma_{0}$ in Theorem 3.1 holds.

##  4 Numerical Experiments

In this section, we present some numerical experiments to validate the
theoretical results. To perform the grid refinement analysis, we choose the
exact solution $u(x,y)=1000\,x^{4}y^{4}(1-x)^{4}(1-y)^{4}$ and calculate the
right hand side. We then calculate the DGFEM solution $u_{h}(x,y)$ and compute
the error in the $L^{2}$ and the energy norms. The numerical experiments were
conducted using FreeFem++. We consider two cases with different choices of
penalty parameters and compare the results of the numerical scheme. In all the
numerical results, we set the parameter $\sigma_{0}=1$.

$p$ |  $h$ |  $\|u-u_{h}\|$ |  $|||u-u_{h}|||$ |  O($\|u-u_{h}\|$) |  O($|||u-u_{h}|||$)
---|---|---|---|---|---
1 |  0.141421 |  0.000419 |  0.011011 |  - |  -
|  0.094281 |  0.000375 |  0.007711 |  0.270249 |  0.878539
|  0.070711 |  0.000371 |  0.006441 |  0.036169 |  0.625625
|  0.056569 |  0.000372 |  0.005790 |  -0.009047 |  0.477075
|  0.047140 |  0.000373 |  0.005412 |  -0.017329 |  0.371163
|  0.040406 |  0.000374 |  0.005172 |  -0.017346 |  0.293645
2 |  0.282843 |  0.000402 |  0.017695 |  - |  -
|  0.141421 |  0.000027 |  0.002197 |  3.887976 |  3.009653
|  0.094281 |  0.000007 |  0.000960 |  3.183386 |  2.043090
|  0.070711 |  0.000003 |  0.000544 |  2.876605 |  1.973371
|  0.056569 |  0.000002 |  0.000350 |  2.695863 |  1.977692
|  0.047140 |  0.000001 |  0.000244 |  2.565206 |  1.983976
3 |  0.282843 |  0.000034 |  0.002809 |  - |  -
|  0.141421 |  0.000002 |  0.000240 |  4.168647 |  3.550827
|  0.094281 |  0.000000 |  0.000068 |  3.915067 |  3.120461
|  0.070711 |  0.000000 |  0.000028 |  3.929946 |  3.049651
|  0.056569 |  0.000000 |  0.000014 |  3.947376 |  3.031438
|  0.047140 |  0.000000 |  0.000008 |  3.960072 |  3.024891
Table 1: Error values and rates of convergence of the Mixed DGFEM for the
biharmonic equation. The penalty parameter is chosen to be
$\alpha_{k}=\sigma_{0}|e_{k}|^{-3}p^{2}$ with $\sigma_{0}=1$. The rate of
convergence for the linear DGFEM is not significant and a sub-optimal
convergence for the piecewise quadratic elements is observed which is an
observation made by Gudi et al. (2008).

First, we compute the rates of convergence for the DGFEM considered by Gudi et
al. (2008). From Table 1 we observe that the rates of convergence for the
linear DGFEM is not significant and a sub-optimal convergence (rate of
convergence decreases rapidly) in the $L^{2}$–norm was observed for piecewise
quadratic elements. The penalty parameter was chosen as
$\alpha_{k}=\sigma_{0}|e_{k}|^{-3}p^{2}$ for the mixed DGFEM. However, the
choice of the penalty parameter is crucial for the mixed method as choosing a
higher value may result in sub-optimal convergence.

$p$ |  $h$ |  $\|u-u_{h}\|$ |  $|||u-u_{h}|||$ |  O($\|u-u_{h}\|$) |  O($|||u-u_{h}|||$)
---|---|---|---|---|---
1 |  0.141421 |  0.004153 |  0.135755 |  - |  -
|  0.094281 |  0.001258 |  0.060027 |  2.944717 |  2.012619
|  0.070711 |  0.000549 |  0.034655 |  2.882794 |  1.909588
|  0.056569 |  0.000291 |  0.022947 |  2.839182 |  1.847442
|  0.047140 |  0.000174 |  0.016491 |  2.814693 |  1.811970
|  0.040406 |  0.000113 |  0.012517 |  2.799510 |  1.788672
2 |  0.282843 |  0.006098 |  0.361025 |  - |  -
|  0.141421 |  0.000578 |  0.071355 |  3.400331 |  2.339021
|  0.094281 |  0.000142 |  0.026470 |  3.452994 |  2.445662
|  0.070711 |  0.000053 |  0.013065 |  3.462819 |  2.454307
|  0.056569 |  0.000024 |  0.007532 |  3.473787 |  2.468250
|  0.047140 |  0.000013 |  0.004798 |  3.478131 |  2.473979
3 |  0.282843 |  0.000159 |  0.018969 |  - |  -
|  0.141421 |  0.000013 |  0.003115 |  3.562213 |  2.606456
|  0.094281 |  0.000003 |  0.001173 |  3.409106 |  2.409145
|  0.070711 |  0.000001 |  0.000522 |  3.817537 |  2.811388
|  0.056569 |  0.000000 |  0.000267 |  4.012594 |  3.007459
|  0.047140 |  0.000000 |  0.000151 |  4.124334 |  3.121269
Table 2: Error values and rates of convergence for the mixed DGFEM for the
biharmonic equation. The rate of convergence is close to $p+1$ in the
$L^{2}$-norm and close to $p$ in the energy norm indicating optimal rates of
convergence. The penalty parameter is chosen to be
$\alpha_{k}=\sigma_{0}|e_{k}|^{-1}p^{2}$ with $\sigma_{0}=1$.

In Table 2, we present the convergence results by choosing a lower value of
the penalty parameter with $\alpha_{k}=\sigma_{0}|e_{k}|^{-1}p^{2}$. The
constant $\sigma_{0}$ in all the cases was set to be equal to $1$.

Figure 1: Rates of convergence of various errors ($y$-axes) with mesh size
($x$-axes) for piecewise linear Discontinuous Galerkin Method. The black solid
line shows the expected rate of convergence predicted by the theoretical
results.

We observe a significant improvement in the convergence rates for piecewise
linear elements and for piecewise quadratic elements. We observe that on
choosing a lesser value of penalty, the solution becomes more accurate with
refinement. This is especially strong in the linear case, where the solution
converges rapidly (rate close to $\approx 2.5$) to the exact solution and the
magnitude of the $L^{2}-$error $\|u-u_{h}\|$ is significantly lower in the
last iteration. This is illustrated in Figure 4 where convergence is observed
only when $\alpha_{k}=\sigma_{0}|e_{k}|^{-1}$. However, the theoretical
results in Theorem 3.3 and 3.4 predict optimal estimates for the energy norm
$|||e_{u}|||$ $(\approx o(1))$ and sub-optimal estimates for $\|e_{u}\|$
$(\approx o(1))$, repectively, despite the observed higher convergence rates
in the case of $\|e_{u}\|$.

Figure 2: Same as Figure 1 but for piecewise quadratic Discontinuous Galerkin
Method.

We observe a higher convergence rate for piecewise quadratic elements in the
$L^{2}$ and energy norms while the error magnitude is comparable to the
previous case in Table 1. The higher convergence rates are also predicted by
the theoretical estimates for piecewise quadratic case. We observe that
optimal convergence rates are preserved for the piecewise cubic case, although
the error magnitudes are higher than that observed from Table 1. The
theoretical error estimates predict optimal convergence rates for the errors
$|||e_{u}|||$ and $\|e_{u}\|$ for both piecewise quadratic and cubic cases.

Figure 3: Same as Figure 1 but for piecewise cubic Discontinuous Galerkin
Method.

(a) $\sigma_{0}=1$

(b) $\sigma_{0}=10$

Figure 4: Convergence of the solution in terms of the maximum value of the
approximate solution ($y$-axis) with the mesh parameter ($x$-axis) for
piecewise linear polynomials. The red curve (+) denotes maximum value of the
DG solution with $\alpha_{k}=\sigma_{0}|e_{k}|^{-1}$ and the blue curve (O)
denotes the maximum value of the DG solution with
$\alpha_{k}=\sigma_{0}|e_{k}|^{-3}$ considered by Gudi et al. (2008). The
horizontal black line denotes the maximum value of the exact solution. We
observe that the convergence of the solution is observed only in the red curve
(lower penalty term), whereas convergence is not observed in the blue curve
(higher penalty term).

Similar observations were made for the convergence in the auxiliary variable
$v$ which is summarized in Table 3. Gudi et al. (2008) observed a sub–optimal
convergence rate ($\approx p-1$) for the auxiliary variable in the $L^{2}$
norm. This was not observed in the case when the penalty parameter was chosen
as $\alpha_{k}=\sigma_{0}|e_{k}|^{-1}p^{2}$ and the convergence rates are
optimal $(\approx p+1)$ in the $L^{2}$ norm. The observed rate of convergence
is better than the theoretically established result in Theorem 3.1, which
predicts a sub–optimal convergence in $v$ $(\approx p)$. We observe that the
current choice of penalty term works well to approximate the auxiliary
variable even for the piecewise linear case. Figures 1, 2 and 3 shows the
decay in the approximation error with respect to the mesh parameter $h$ for
piecewise linear, quadratic and cubic polynomial approximations, respectively.

| $\alpha_{k}=\sigma_{0}|e_{k}|^{-1}p^{2}$ | $\alpha_{k}=\sigma_{0}|e_{k}|^{-3}p^{2}$
---|---|---
$p$ |  $N$ |  $\|v-v_{h}\|$ |  O($\|v-v_{h}\|$) |  $\|v-v_{h}\|$ |  O($\|v-v_{h}\|$)
1 |  0.141421 |  0.062298 |  - |  0.030978 |  -
|  0.094281 |  0.033646 |  1.519314 |  0.022990 |  0.735562
|  0.070711 |  0.020921 |  1.651590 |  0.024028 |  -0.153535
|  0.056569 |  0.014233 |  1.726174 |  0.025236 |  -0.219860
|  0.047140 |  0.010340 |  1.752746 |  0.026058 |  -0.175885
|  0.040406 |  0.007888 |  1.755539 |  0.026604 |  -0.134564
2 |  0.282843 |  0.069201 |  - |  0.038757 |  -
|  0.141421 |  0.006641 |  3.381222 |  0.005352 |  2.856424
|  0.094281 |  0.001834 |  3.173283 |  0.002105 |  2.301325
|  0.070711 |  0.000736 |  3.176304 |  0.001346 |  1.555525
|  0.056569 |  0.000362 |  3.171694 |  0.001035 |  1.176948
|  0.047140 |  0.000205 |  3.136815 |  0.000855 |  1.044048
3 |  0.282843 |  0.013709 |  - |  0.006130 |  -
|  0.141421 |  0.001251 |  3.454083 |  0.000502 |  3.609629
|  0.094281 |  0.000250 |  3.972908 |  0.000288 |  1.370246
|  0.070711 |  0.000076 |  4.143599 |  0.000188 |  1.479595
|  0.056569 |  0.000030 |  4.227323 |  0.000128 |  1.722811
|  0.047140 |  0.000014 |  4.278536 |  0.000091 |  1.869907
Table 3: The rates of convergence for the auxiliary variable for the DGFEM
with respect to the mesh size and penalty parameter. We observe optimal
convergence rate for the auxiliary variable in the $L^{2}$ and the energy norm
for the current choice of penalty parameter and suboptimal rate of $p-1$ for
the choice in Gudi et al. (2008). We set $\sigma_{0}=1.$

##  5 Conclusion

In this work, we have considered a mixed Discontinuous Galerkin Finite Element
Method to solve the biharmonic equation subject to clamped boundary
conditions. We have derived the weak formulation of the problem and
established error estimates for $h$-refinement. The theoretical results
predict optimal convergence rates in the energy norm for the primal variable
$u$, whereas optimal convergence is predicted for the piecewise quadratic and
cubic case in the $L^{2}$-norm. We performed a series of numerical experiments
using FreeFem++, and verified the theoretical results. We observed that the
choice of the penalty term is crucial and must be chosen to be of the form
$\alpha_{k}=\sigma_{0}|e_{k}|^{-1}p^{2}$ to obtain optimal error estimates in
$h$-refinement. Significant improvements in convergence rates for the
piecewise linear and quadratic elements were observed as a result.

## References

  * Arnold et al. (2002) Arnold, D. N., Brezzi, F., Cockburn, B. & Marini, L. D. (2002)  Unified Analysis of Discontinuous Galerkin Methods for Elliptic Problems.  SIAM Journal on Numerical Analysis, 39, 1749–1779.
  * Brenner & Scott (2007) Brenner, S. C. & Scott, R. (2007)  The Mathematical Theory of Finite Element Methods, vol. 15.  Springer Science & Business Media.
  * Brenner & Sung (2005) Brenner, S. C. & Sung, L.-Y. (2005)  $C^{0}$ Interior Penalty Methods for Fourth Order Elliptic Boundary Value Problems on Polygonal Domains.  Journal of Scientific Computing, 22, 83–118.
  * Ciarlet & Raviart (1974) Ciarlet, P. G. & Raviart, P.-A. (1974)  A Mixed Finite Element Method for the Biharmonic Equation.  Mathematical Aspects of Finite Elements in Partial Differential Equations.  Elsevier, pp. 125–145.
  * Cockburn et al. (2012) Cockburn, B., Karniadakis, G. E. & Shu, C.-W. (2012)  Discontinuous Galerkin Methods: Theory, Computation and Applications, vol. 11.  Springer Science & Business Media.
  * Danumjaya & Pani (2012) Danumjaya, P. & Pani, A. K. (2012)  Mixed Finite Element Methods for a Fourth Order Reaction Diffusion Equation.  Numerical Methods for Partial Differential Equations, 28, 1227–1251.
  * Engel et al. (2002) Engel, G., Garikipati, K., Hughes, T., Larson, M., Mazzei, L. & Taylor, R. L. (2002)  Continuous/Discontinuous Finite Element Approximations of Fourth-Order Elliptic Problems in Structural and Continuum Mechanics with Applications to Thin Beams and Plates, and Strain Gradient Elasticity.  Computer Methods in Applied Mechanics and Engineering, 191, 3669–3750.
  * Gudi et al. (2008) Gudi, T., Nataraj, N. & Pani, A. K. (2008)  Mixed Discontinuous Galerkin Finite Element Method for the Biharmonic Equation.  Journal of Scientific Computing, 37, 139–161.
  * Monk (1987) Monk, P. (1987)  A mixed finite element method for the biharmonic equation.  SIAM Journal on Numerical Analysis, 24, 737–749.
  * Riviére (2008) Riviére, B. (2008)  Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations: Theory and Implementation.  SIAM.