# A Spatiotemporal Illumination Model for
3D Image Fusion in Optical Coherence Tomography

###### Abstract

Optical coherence tomography (OCT) is a non-invasive, micrometer-scale imaging
modality that has become a clinical standard in ophthalmology. By raster-
scanning the retina, sequential cross-sectional image slices are acquired to
generate volumetric data. In-vivo imaging suffers from discontinuities between
slices that show up as motion and illumination artifacts. We present a new
illumination model that exploits continuity in orthogonally raster-scanned
volume data to accurately compensate illumination artifacts. Our novel
spatiotemporal parametrization adheres to illumination continuity both
temporally, along the imaged slices, as well as spatially, in the transverse
directions. Yet, our formulation does not make inter-slice assumptions, which
could have discontinuities. This is the first optimization of a 3D inverse
model in an image reconstruction context in OCT. Evaluation in 68 volumes from
eyes with pathology showed reduction of illumination artifacts in 88% of the
data, and only 6% showed moderate residual illumination artifacts. The method
enables the use of forward-warped motion corrected data, which is more
accurate, and enables supersampling and advanced 3D image reconstruction in
OCT.

Index Terms—  Optical coherence tomography, illumination correction, inverse
model, ophthalmology, image reconstruction

##  1 Introduction & problem statement

Optical Coherence Tomography (OCT) is a non-invasive 3D optical imaging
modality and standard of care in ophthalmology [1, 2], where it is used for
retinal imaging. In OCT angiography (OCTA), scans are repeated and, besides
the standard OCT signal, a metric of its variance is computed [3]. This
displays vasculature, because moving blood cells cause varying backscattering,
while it is constant for tissue. In eyes with pathology, a recent study found
severe image artifacts that are associated with the reliability of
quantitative biomarkers in 53.5% of the data [4]. In this study, motion-
related artifacts (16.0%) were among the 3 most prevalent types of artifacts.
Motion artifacts arise because OCT images are acquired by raster-scanning a
laser beam across the retinal surface for a multi-second acquisition duration
and apply to OCT and OCTA data in the same way. In their review on image
artifacts in OCTA, Anvari et al. list 6 types of scanning-related artifacts
[5], which we categorize in dominantly distortion-related (displacement,
ghosting/doubling, stretch and quilting/crisscross) and illumination-related
(blink and banding). Again, illumination-related artifacts exist in OCT and
OCTA data, and originate from signal variation in the OCT signal. A
particularity of OCT illumination artifacts is that they can be highly
anisotropic due to raster scanning, with slow variation along a scan line, but
sharp changes between them (Figure 1a).

Besides eye tracking-based approaches, image distortions are compensated via
retrospective software motion correction methods by aligning multiple scans
[6, 7], which can subsequently be merged to improve image quality. The
combined approach was shown to be superior to tracking-based approaches alone
[8], which are fundamentally limited by latency. Recently, we demonstrated
forward warping motion correction in OCT [9], which achieves a general
improvement in accuracy, reducing displacement and ghosting/doubling
artifacts. More specifically, forward warping correctly localizes gaps, which
solves stretch artifacts, and is able to correctly model eye motion in the
slowly-scanned direction, which enables correction of distortions that were
responsible for quilting artifacts in prior methods. Despite the clear need
for forward warping motion correction, its utilization is currently hindered
in clinical data because it introduces a new challenge for data fusion. It
results from illumination differences between the coregistered scans, which is
closely related to the second category of scanning artifacts.

Illumination / signal level can vary during acquisition of OCT volumes for
various reasons, including: Eye blinks; attenuation when the beam passes
through “floaters”, which are non-stationary opacities in the vitreous;
attenuation from corneal drying; vignetting on the iris caused by insufficient
instrument alignment and eye motion; and depth-dependent attenuation from
intensity roll-off. Many of these effects worsen with age and in pathology.

Eye blinks block the imaging beam entirely, which reduces signal levels to
background noise and can not be corrected. All other forms of change in
illumination are spatially continuous and, normally, vary only slowly in time,
which is tolerated in OCT imaging. Banding artifacts emerge when the temporal
continuity is interrupted by scan interruption (which can occur when scanners
reset the position in eye tracking), or when the spatial continuity is
interrupted between successive B-scans due to eye motion.

a

b

Fig. 1: Global illumination variation (ellipse) and banding artifacts (arrows)
in enface images (depth-averaged volumes). Left: Input with vertical B-scan
orientation. Right: forward-warped & merged.

Retrospective motion correction generally reduces varying illumination via
averaging, but can also introduce additional artifacts: When image
registration maps non-subsequently acquired B-scans next to each other to
correct distortions, only spatial continuity can be restored, while temporal
continuity is not. Similar to scan discontinuities, this can appear as a
banding-like artifact. This especially applies to larger, saccadic motion that
is responsible for quilting artifacts. In these cases, forward warping motion
correction alone, which only corrects distortion, is insufficient to remove
the artifact. Lastly, the coregistered volumes usually differ in illumination,
and the merged data displays the average illumination. However, if one of the
registered volumes has a gap, merged image intensity abruptly changes to the
illumination level of the remaining volume(s), which can again result in a
banding-like artifact. The same effect occurs when eye-blinks are removed from
a scan before merging. Similarly, in forward warping, the merged volume is not
computed by uniform averaging between the coregistered volumes. In contrast,
the distance from the merged voxel to each individual A-scan can be taken into
account by assigning closer A-scans a higher weight. While this can be
advantageous for spatial accuracy, it makes the contribution from each volume
change on a per-voxel basis, and with it, interpolates between the
illumination levels of the merged volumes.

While banding artifacts are distracting for human readers (Figure 1a), they
can confound algorithmic quantitative disease markers: sharp changes can cause
false edges in segmentation of vessels or pathologic features, and change in
intensity can adversely affect thresholding as used in the segmentation of
retinal layers, or hyperreflective foci. The illumination variation in
forward-warped and merged images appears as intensity deviation that follows
the B-scan acquisition directions, which can make them look similar to
quilting artifacts (Figure 1b). However, given the continuity of the
underlying illumination, many artifact lines appear over larger areas, which
can be prohibitive, especially in clinical use of forward warping motion
correction, for both human and algorithmic image analysis.

There are only few studies on the correction of artifacts from inter-scan
illumination differences in OCT, and general illumination normalization
approaches can not be transferred because they do not model the anisotropy in
OCT. In the work of Wang et al., stitched OCTA volumes were linearly
interpolated over a wide range to achieve smooth transitions between
illumination levels [10]. However, wide range interpolation cannot compensate
sharp illumination changes between scan lines or per voxel weight changes in
forward warping. Kraus et al. reported an unsharp mask-like illumination
normalization method as preprocessing step for motion correction [11].
However, by normalization, the approach alters morphology of image features
like the darker appearing optic nerve head (ONH), or geographic atrophy (GA),
a pathology appearing in late-stage age-related macular degeneration. Removal
of such features would be detrimental to clinical image analysis. Furthermore,
the approach does not take the anisotropy of illumination into account,
leaving residual artifacts at sharp changes.

In this paper, we propose an inverse model formulation of illumination bias
which not only models the spatial, but, for the first time, also temporal
(dis)continuity of illumination on a per-A-scan level. By fitting the model
via minimization of inter-scan illumination differences in registered OCT
volumes, strong changes of illumination, as they can appear between B-scans of
a single scan, and merging artifacts from illumination differences between
scans, are resolved. At the same time, contrast of morphology like the foveal
avascular zone or the ONH, and pathologies like GA, are preserved. By using
scans with orthogonal B-scan directions, sharp anisotropic changes in
illumination are accurately estimated and corrected.

##  2 Methods

We will first define the illumination model and proceed with its optimization,
which is performed by iterative image reconstruction. For prior displacement
field computation, and efficient warping of OCT images, we refer to our recent
work on motion correction [9].

###  2.1 Illumination Model and Parameterization

We model the measured signal $t$ as the product of illumination strength $a$
with sample backscattering $t^{*}$, and added noise $n$.

$t=a\cdot t^{*}+n$ (1)

For strong backscattering, the term is dominated by the multiplicative part
which can, given knowledge of $a$, be inverted to obtain an estimate
$t^{*}\approx\frac{1}{a}\cdot t$. For weak backscattering (in the background),
the noise term dominates and $t^{*}$ cannot be estimted. In these cases, a
correction with $\frac{1}{a}$ should not be applied to avoid noise
amplification. More importantly, these signals must be excluded from the
illumination estimation process to avoid bias. To maintain the appearance of
OCT data and because the goal of this paper is limited to illumination
artifact correction, we simply copy $t$, i.e., the noise-measurement, which is
independent from $a$, and reconstruct the image

$\tilde{t}=\begin{cases}\frac{1}{a}\cdot t&t^{*}\text{ in foreground}\\\
t&t^{*}\text{ in background}\end{cases}=\left(\frac{1}{a}\right)^{m}\cdot t$
(2)

To decide between the cases, we perform a first estimate of $t^{*}$ by median-
filtering the input B-scans, to reduce noise, and apply a small threshold
$t_{\text{min}}$. The result is a binary foreground mask-value $m$. Since
deviations in illumination are ignored, this estimation is inaccurate for
large $t^{*}$, but the threshold value is chosen at the transition to the
background noise level, where the multiplicative part has small influence
(because $t^{*}$ is small and $a$ is generally close to $1$).

To simplify computations and increase numeric stability, we evaluate Eq. (2)
in logarithmic scale and define $s=\log t$. Because A-scans are acquired in a
single acquisition, there is no illumination variation within them. Since each
B-scan is acquired spatially and temporally continuous, we assume continuous
illumination within them. Due to temporal discontinuity and potential motion,
we refrain from inter-B-scan assumptions. Therefore, we model the correction
coefficient $\log\frac{1}{a}$ by a Hermite spline $c_{i}^{V}(j)$ for each
B-scan $i$ and volume $V$ and index with the A-scan index $j$ over the spline
parameters $\boldsymbol{c}$. We use a parameter density of 1 per millimeter.
Given a voxel index $k$, the corrected log-signal $\tilde{s}$ can then be
computed via

$\tilde{s}_{i,j,k}^{V}(\boldsymbol{c})=s_{i,j,k}^{V}+m_{i,j,k}^{V}\cdot
c_{i}^{V}(j).$ (3)

###  2.2 Model Optimization

Figure 2 sketches the optimization goal, which is a key to understanding the
model optimization. We minimize the objective

$\mathcal{J}(\boldsymbol{c})=\sum\nolimits_{M\in\mathcal{V}}\mathcal{D}^{M}(\boldsymbol{c})+\lambda\cdot\mathcal{R}(\boldsymbol{c})\quad\text{s.t.}\quad\mathcal{C}(\boldsymbol{c})=0,$
(4)

X-fast B-scan signal $s_{i,j}^{X}$Correction coefficient
$c_{i}^{X}(j)$Average($s_{i,j}^{X},s^{Y}$)Correction coefficient
$c^{Y}$$s_{i,j}^{X}+c_{i}^{X}(j)\stackrel{{\scriptstyle!}}{{=}}s^{Y}+c^{Y}$Y-fast
signal $s^{Y}$ Fig. 2: Compensation of illumination artifacts. Top left:
Schematic of X-fast B-scans (blue) and Y-fast B-scans (red) after motion
correction. Bottom: Visualization of the optimization goal in logarithmic
scale, with legend / zoom in the top right. Blue curve: X-fast B-scan signal
corresponding to the dark blue line in the top left. For clarity, we assume a
sample with uniform (foreground) backscattering such that changes in measured
signal only originate from illumination variance. Therefore, given the
continuous illumination assumption in B-scans, the blue curve is continuous.
The red circles represent the intensities of the orthogonal scan, which are
interpolated and therefore written without indices. Because they originate
from different B-scans, their illumination is not necessarily continuous. The
dotted black curve, corresponding to the average of the uncorrected
intensities, is discontinuous where $s^{Y}$ is discontinuous or has gaps. In
forward warping, the averaging weights of the red intensities vary between
points, which would distribute the black dots discontinuously somewhere
between the blue and red inputs. The blue arrows represent the illumination
correction of the X-fast intensities, which are continuous due to their spline
parametrization. Added to the blue curve, they form the dashed green curve,
the corrected X-fast intensities. Correction coefficients of the orthogonal
scan are optimized such that the corrected intensities match the green line up
to noise (in a different addend of the outer sum in Eq. 5).

where the constraint $\mathcal{C}(\boldsymbol{c})=\sum c_{i}$ enforces global
intensity to match the inputs. L2-norm penalization
$\mathcal{R}(\boldsymbol{c})=||\boldsymbol{c}||_{2}^{2}$ of illumination
change preserves the original contrast, because it centers the illumination
corrected result (green line) between the inputs. The dataterm

$\mathcal{D}^{M}(\boldsymbol{c})=\sum_{\begin{subarray}{c}T\in\mathcal{V}\\\
\operatorname{dir}(T)\\\
\neq\operatorname{dir}(M)\end{subarray}}\sum_{i,j,k}\left(\tilde{s}_{i,j,k}^{M}(\boldsymbol{c})-\mathcal{W}\big{(}\tilde{S}^{T},\boldsymbol{\tilde{x}}_{i,j,k}^{M}\big{)}\right)^{2}$
(5)

penalizes the squared difference between each voxel (inner sum) of the current
volume $M$ and each orthogonally scanned volume $T$ (outer sum). For this
purpose, illumination corrected target volumes $\tilde{S}^{T}$ are Hermite-
interpolated ($\mathcal{W}$) to the current voxel location
$\tilde{x}_{i,j,k}^{M}$. If there is a gap in the target volume, the
difference is excluded.111In contrast to Figure 2, backscattering of non-
constant samples varies arbitrarily depending on location. However, given
correct co-registration, this location-dependent effect cancels out in the
difference computation, making it irrelevant for optimization. Therefore,
Figure 2 is a valid visualization of the optimization; and the continuity of
correction coefficients along B-scans generalizes beyond the constant value
assumption. We start the optimization by assuming uniform illumination
($\boldsymbol{c}:=\boldsymbol{0}$) and minimize via iterative momentum
gradient descent.

The crux of the approach is that the inter-B-scan illumination continuity is
enforced by the orthogonal target B-scans (whose spline-corrected illumination
is guaranteed to remain continuous, Figure 2). This also utilizes the
registration to enforce illumination continuity in motion corrected space, as
opposed to scanner space as in [11], achieving continuous illumination in the
merged image. Furthermore, gaps and varying merging weights are resolved:
Because the corrected local intensities differ only by noise, weight changes /
exclusion of a scan no longer impact the merged average.

##  3 Evaluation & Results

###  3.1 Data

This research study was conducted retrospectively using human subject data.
OCT imaging was performed at New England Eye Center (NEEC) in line with the
principles of the Declaration of Helsinki. Approval was granted by Committee
on the Use of Humans as Experimental Subjects at MIT (#1901669061) and Tufts
Health Sciences Institutional Review Board at NEEC (#MODCR-03-13240). The
dataset contains 18 subjects (70 $\pm$ 11 years old) with ocular pathologies,
described in detail in [9]. Four $6\times 6$ mm fovea-centered scans were
acquired with alternating B-scan orientation. Illumination correction was
performed on all 4 combinations of orthogonal scan pairs, totalling to 72
merged volumes. One subject was retrospectively excluded because illumination
in the same peripheral region was so low in all 4 scans, that all retinal
layers except the RPE vanished in background noise. A correction at this
location is therefore impossible. Yet, other image areas in these scans were
well corrected. The prototype spectral domain OCT instrument images at 840 nm
wavelength with 128 kHz A-scan rate, and acquired $500\times 500$ A-scans with
3 µm resolution (FWHM) [12]. A-scans were reconstructed with 775 pixels and
1.78 µm spacing.

###  3.2 Qualitative Evaluation

1a

1b

1c

1d

1e

2a

2b

2c

2d

2e

Fig. 3: Representative enface images (depth-averaged volumes) from volumes
merged without (row 1) and with illumination correction (row 2). (a) Good
results after large and (b) small improvement. The zooms in (a) show a
checkerboard overlay of the registered volumes before merging. Illumination
difference makes the individual tiles in (1a) clearly distinguishable, while
they are almost uniform after correction in (2a). (c) Minor residual banding
artifacts (the black squares are gaps which originate from large saccadic
motion). (d) Worst result (besides the residual banding artifact, slight lines
with residual illumination deviation are visible throughout on a monitor). (e)
The volume with most worsening (the other two worsened much less). In B-scans
(which are not depth-averaged), the same change of illumination would be less
visible due to higher noise levels. Contrast is adjusted jointly per column,
corruptions in dark areas arise from the publication process.

We categorized the result images by (a) severity of residual illumination
artifacts and (b) change of illumination artifacts, and present the results in
Table 1. Minor artifacts are roughly evenly split among typically small-extent
border artifacts, and more globally appearing but weak artifacts. The latter
ones are barely visible in the enface images, and vanish in the higher noise
present in non-averaged volumetric data. Representative examples are presented
in Figure 3.

$\downarrow$Result \ change$\rightarrow$ | + | o | - |  |
---|---|---|---|---|---
no artifacts | 45 | 2 | 0 | 47 | 69.1%
minor artifacts | 14 | 2 | 1 | 17 | 25.0%
moderate artifacts | 1 | 1 | 2 | 4 | 5.9%
| 60 | 5 | 3 | 68 |
| 88.2% | 7.4% | 4.4% |  |
Table 1: Contingency table of categories no residual artifacts, minor residual
artifacts and moderate residual artifacts; and improvement (+), inconclusive
change (o), and worsening (-).

###  3.3 Quantitative Evaluation

While the strong banding artifacts often result from gaps in one of the
volumes, the line artifacts that show up in forward-warped data (Figure 1b)
originate from illumination difference between the registered volumes. To
confirm our previous observations, we evaluated, if the illumination
correction reduces the signal difference between the registered (but not
merged) volumes. To reduce the impact of noise, we averaged the volumes along
depth, where we only included registered voxels that were present in both
volumes to avoid bias. Complying with the predominantly additive noise
assumption, resampling and averaging was performed in linear scale. (Images in
Figure 3 were computed in the same way). When computing the differences, we
found a reduction in mean absolute difference of 22.5%. Note that additional
effects, like change in focus, misregistration, noise, and others, also
increase the difference, so illumination correction alone can only reduce the
difference, but not reach zero. 66 differences decreased, while 2 increased
(which coincided with 2 of the 3 cases observed to have worsened in Table 1,
while the third one had the second-smallest decrease and was close to zero).

##  4 Discussion

We performed an in-depth analysis on the formation of image artifacts in OCT
and designed a physically-motivated model that can compensate for illumination
variance. We parameterized the model using only the relevant degrees of
freedom to achieve most-accurate compensation of temporal changes via
overdetermined estimation, while the spline-based formulation ensures that
clinically relevant features are maintained. By enforcing the corrected
intensities to fit the average intensities, spatial illumination effects are
averaged and therefore converge to zero with increasing number of scans.
Optimization was performed using the first 3D inverse model in OCT.

In our challenging clinical dataset, our approach reduces illumination-related
imaging artifacts in 88% of cases, while only 4.4% worsened. After correction,
only 6% of cases had moderate residual illumination artifacts, which were
easily identifiable in B-scan slices. We reconfirmed these numbers with a
quantitative evaluation of residual difference between registered data after
illumination correction. Weak artifacts only become visible after axial
averaging, but vanish in noise when viewing B-scans, which is the typical
clinical use-case. Small residual effects can sometimes be seen at the image
boundary. While gaps inside volumes are not an issue for curve estimation,
extrapolation at the image borders is more error-prone. More sophisticated
extrapolation could potentially further reduce artifacts in these areas. On a
second level, this work enables use of forward warping-registration for
accurate motion correction and merging of clinical data. Therefore, largely
improved accuracy, resolution of stretch artifacts, and more reliable
correction of doubling and quilting artifacts are enabled [9].

Despite the promising results, our method has limitations: Orthogonal scanning
and accurate registration is a prerequisite. Furthermore, image continuity and
smallest correctable scale of illumination is a trade-off, therefore shadows
from small-scale floaters are not corrected. However, since floaters are
moving, they would be reduced by volume averaging.

By largely resolving the multiplicative illumination effect, this work enables
high-quality 3D denoising via iterative reconstruction methods. Furthermore,
super-resolution reconstruction of volumes and enface images is enabled
including the transverse directions. Here, due to the denser sampling in the
merged volume, individual voxels are only merged from a subset of the input
scans, while neighboring pixels are filled from other scans that differ in
illumination. Scan merging weights and, consequently, averaged brightness
would vary from pixel to pixel. As illumination inconsistency increases with
scan time, its correction will remain relevant in widefield imaging and in OCT
angiography scan protocols [3], especially if multiple repeats are required
for flow speed analysis [13].

##  5 Author Contributions

SP is the main author. Algorithm design, evaluation, implementations, initial
writing: SP. Data: JW, JG, KL, NW. Funding: SP, AM, JF. All authors discussed
the results and their presentation, and contributed to the final manuscript.

##  6 Acknowledgment

We acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation), project 508075009, and National Institutes of Health,
project 5-R01-EY011289-36.

## References

  * [1] David Huang, Eric A. Swanson, and Charles P. Lin,  “Optical coherence tomography,”  Science, vol. 254, no. 5035, pp. 1178–1181, 11 1991.
  * [2] James Fujimoto and Eric Swanson,  “The development, commercialization, and impact of optical coherence tomography,”  Invest. Ophthalmol. Vis. Sci., vol. 57, no. 9, pp. OCT1, 2016.
  * [3] RF Spaide, JG Fujimoto, NK Waheed, SR Sadda, and G Staurenghi,  “Optical coherence tomography angiography,”  Progress in Retinal and Eye Research, vol. 64, pp. 1–55, 2018.
  * [4] Ian C. Holmen, Sri Meghana Konda, and Jeong W. Pak,  “Prevalence and Severity of Artifacts in Optical Coherence Tomographic Angiograms,”  JAMA Ophthalmology, vol. 138, no. 2, pp. 119–126, 02 2020.
  * [5] Pasha Anvari, Maryam Ashrafkhorasani, and Abbas Habibi,  “Artifacts in optical coherence tomography angiography,”  J. Ophthal. & Vis. Research, vol. 16, no. 2, pp. 271–286, 2021.
  * [6] Luisa Sánchez Brea, Danilo Andrade De Jesus, and Muhammad Faizan Shirazi,  “Review on retrospective procedures to correct retinal motion artefacts in OCT imaging,”  Applied Sciences, vol. 9, no. 13, 2019.
  * [7] Stefan B. Ploner, Martin F. Kraus, and Eric M. Moult,  “Efficient and high accuracy 3-d oct angiography motion correction in pathology,”  Biomed. Opt. Express, vol. 12, no. 1, pp. 125–146, Jan 2021.
  * [8] Acner Camino, Miao Zhang, and Simon S. Gao,  “Evaluation of artifact reduction in oct angiography with real-time tracking and motion correction technology,”  Biomed. Opt. Express, vol. 7, no. 10, pp. 3905–3915, Oct 2016.
  * [9] Stefan Ploner, Siyu Chen, and Jungeun Won,  “A spatiotemporal model for precise and efficient fully-automatic 3d motion correction in oct,”  in Medical Image Computing and Computer Assisted Intervention – MICCAI 2022, 2022, pp. 517–527.
  * [10] Jie Wang, Acner Camino, Xiaohui Hua, Liang Liu, David Huang, Thomas S. Hwang, and Yali Jia,  “Invariant features-based automated registration and montage for wide-field oct angiography,”  Biomed. Opt. Express, vol. 10, no. 1, pp. 120–136, Jan 2019.
  * [11] Martin Kraus, Jonathan Liu, and Julia Schottenhamml,  “Quantitative 3D-OCT motion correction with tilt and illumination correction, robust similarity measure and regularization,”  Biomed. Opt. Express, vol. 5, no. 8, pp. 2591–2613, 7 2014.
  * [12] ByungKun Lee, Siyu Chen, and Eric M. Moult,  “High-Speed, Ultrahigh-Resolution Spectral-Domain OCT with Extended Imaging Range Using Reference Arm Length Matching,”  Trans. Vis. Sci. & Tech., vol. 9, no. 7, pp. 12–12, 06 2020.
  * [13] Stefan B. Ploner, Eric M. Moult, and WooJhon Choi,  “Toward quantitative optical coherence tomography angiography,”  Retina, vol. 36, pp. S118–S126, 2016.


CLUSTERING VIETNAMESE CONVERSATIONS FROM FACEBOOK PAGE TO BUILD TRAINING
DATASET FOR CHATBOT

Trieu Hai Nguyen1, Thi-Kim-Ngoan Pham1, Thi-Hong-Minh Bui1,

Thanh-Quynh-Chau Nguyen1

1 Faculty of Information Technology, Nha Trang University,

02 Nguyen Dinh Chieu Street, Nha Trang City, Vietnam

e-mails<EMAIL_ADDRESS>
###### Abstract

The biggest challenge of building chatbots is training data. The required data
must be realistic and large enough to train chatbots. We create a tool to get
actual training data from Facebook messenger of a Facebook page. After text
preprocessing steps, the newly obtained dataset generates FVnC and Sample
dataset. We use the Retraining of BERT for Vietnamese (PhoBERT) to extract
features of our text data. K-Means and DBSCAN clustering algorithms are used
for clustering tasks based on output embeddings from PhoBERTbase. We apply
V-measure score and Silhouette score to evaluate the performance of clustering
algorithms. We also demonstrate the efficiency of PhoBERT compared to other
models in feature extraction on the Sample dataset and wiki dataset. A
GridSearch algorithm that combines both clustering evaluations is also
proposed to find optimal parameters. Thanks to clustering such a number of
conversations, we save a lot of time and effort to build data and storylines
for training chatbot. ††footnotetext: Corresponding author: Trieu Hai Nguyen
<EMAIL_ADDRESS>
Preprint submitted to JJCIT December 10, 2021

Keywords: BERT, Clustering, Language models, Feature Extraction, Word
Embeddings

##  1 Introduction

Chatbot is certainly not an unfamiliar name in the field of Natural Language
Processing (NLP). Previously, traditional chatbot simply interacts with the
user via predefined rules, which means the user is only allowed to enter by
these rules to get answers. However, NLP chatbot is not only a word
recognition algorithm but it can also understand what the user is saying. It
is one of the pioneering applications using Artificial Intelligence (AI)
namely, NLP to help humans interact with machine like humans with humans via
Virtual Assistant autoresponder. Currently, there are many parties developing
NLP chatbot, which can be mentioned as Google’s DialogFlow, Watson of IBM and
Rasa.

In the process of building NLP chatbots, all chatbots require real datasets
for training bot. The training datasets can be large or small depending on the
size and intelligence level of the chatbots. Raw training data can be
collected from past conversations through social media, archived user chats,
previous questions, email chains, or live telephone transcripts. But these
data are messy, not in any structure or order, and come from various sources
collected with huge amounts of raw data. Thus the first priority when
constructing a chatbot is to transform those raw data into useful data for the
purpose of training bot.

In order to fit chatbot building orientation, that raw dataset needs to be
divided into specific intents, which serves to build conversations to train a
chatbot. There are many ways to process that raw dataset into specific intents
(topics, conversations). The first method to be mentioned is using Supervised
Learning task [1] to classify intents. In particular, this method requires
labeling for the input examples and then it predicts labels for remaining data
in the raw dataset. In this case, label prediction corresponds to classified
raw data into the intents that have been labeled previously. However, building
and labeling manual intents on large datasets lead to big challenges for
chatbot developers. With a simple raw dataset of about 8000 conversations,
analyzing how many intents are created is a conundrum.

Instead, we can approach the above problem by using the second method, which
is clustering similar raw data together into corresponding intents. The
advantage of this approach is that it uses Unsupervised Learning technical
[1], which is only based on the features of the data to perform specific tasks
such as Clustering. As expected, this method has a significant effect on the
analysis of raw data. It saves us a lot of time and effort to make a training
dataset for chatbot. There are many clustering algorithms such as K-Means,
DBSCAN, BIRCH, Spectral clustering [2, 3]. In this article, we use K-Means [4]
and DBSCAN [5] techniques to cluster our dataset and consider that clustering
is a downstream NLP task. Each technique has its own advantages and
disadvantages. K-Means algorithm is a simple and fast implementation but it
requires knowing the number of clusters to perform clustering whereas DBSCAN
does not. Nevertheless, DBSCAN technique is more difficult to implement and
requires finding the optimal parameters [5, 6, 7], which leads to drastically
increased costs, especially for large datasets. Thus, we can combine the
advantages of both techniques to serve the purpose of efficient clustering.

The input of clustering algorithms in particular and downstream NLP tasks in
general is document embeddings extracted from the dataset. There are many ways
to extract information from text datasets, for example, we can use traditional
machine learning algorithm like TF-IDF [8], proposed word embedding models in
recent years such as Word2Vec, GloVe [9, 10], FastText[11] or popular language
models like GPT-2 [12], BERT model and its variations [13, 14, 15]. In this
work, we use BERT (Bidirectional Encoder Representations from Transformers),
which is state-of-the-art embeddings [13] to extract features of documents.
Recently a clustering approach with the BERT model has been proposed by O.
Gencoglu [16]. As suggested in [17] clustering techniques using pre-trained
transformer language models are applied to short text clustering. The
combination of word embeddings using BERT models and clustering algorithms to
obtain topics was presented in [18]. In distinctive, PhoBERT is Pre-trained
language models for Vietnamese, which is used to embed our Vietnamese dataset
[15]. V-measure score [19] and Silhouette score [20] are used to evaluate the
performance of a clustering algorithms as well as the feature extraction
efficiency of the language models.

The aim of the present paper is to study and apply PhoBERT model to our
Facebook Vietnamese conversations dataset, thereby deriving document
embeddings in order to serve the clustering task. The combination of both
K-Means and DBSCAN clustering algorithms is proposed by us to achieve the best
clustering results on the actual dataset. The finding of these data clusters
allows us to simplify and accelerate the building of a training dataset for
chatbot. In section 2, we recall some theories of Transformer and BERT
architecture proposed by Vaswani et al. in [21] and Devlin et al. in [13]
respectively. In section 3, we offer an approach to apply PhoBERT to the
clustering task from the idea of classification task [15, 22]. Next, we also
recall clustering algorithms in machine learning and evaluation metrics for
unsupervised learning algorithms in Section 4. Some experiments on our
Facebook Vietnamese conversations dataset (includes FVnC and Sample dataset)
and wiki dataset such as searching optimal parameters, clustering performance
evaluations as well as clustering results are considered in Section 5. In
particular, we show that among the models that support Vietnamese, PhoBERT’s
feature extraction efficiency is the best based on V-measure score.
Ultimately, we give some conclusions in Section 6. The code, datasets and pre-
trained models are available at
https://github.com/trieuntu/conversation_clustering.

##  2 Related Work

We provide some background knowledge about Transformer architecture, Pre-
Trained Language Models especially BERT. From these theoretical constructs, we
apply them to solve our NLP tasks.

###  2.1 Transformer

Transformer architecture was first introduced in paper “Attention Is All You
Need” by [21]. At the time of launch, this architecture is a new breakthrough
in the field of natural language processing and related tasks. Currently, when
dealing with sequence-to-sequence models in NLP, the transformer is still one
of the state-of-the-art (SOTA) types of model and completely replaces RNN/LSTM
[23]. Transformer architecture overcomes the disadvantages of RNN and its
variations. For instance, it doesn’t take advantage of GPU parallelism because
it has to process input word by word sequentially into encoder/decoder, and
the information is easily lost during propagation through hidden layers for
long input sentences.

Transformer architecture contains two parts, Encoder attention and Decoder
attention. According to the original article of [21], the encoder part has 6
layers, each of which has two sublayers, which are multi-head self-attention
and fully connected feed-forward. Decoder part is similar to the encoder part
but it adds a masked multi-head attention sublayer and the last layer of the
encoder part will be passed to the multi-head attention sublayer in the
decoder part. Note that, the input of both parts is the sum of positional
encoding vector and word vector embedding.

The attention mechanism is the most important component of transformer
architecture. Self-attention sublayer is an attention mechanism, which
contains the weight sets of the model $W_{q},W_{k},W_{v}$ to be trained. The
attention mechanism presents the relation of a word to all its related words
in the sequence based on the adjustment of the above sets of weights. The
product of the input embedding layer and $W_{q},W_{k},W_{v}$ is matrixes Query
$Q$, Value $V$ and Key $K$. In order to calculate Attention vector of word $i$
to the rest of the words, Vaswani et al. [21] have given the formula

$Attention_{i}(Q_{i},K_{i},V_{i})=softmax\left(\frac{Q_{i}K_{i}^{T}}{\sqrt{d_{K_{i}}}}\right).V_{i}$

where $d_{K}$ is the dimension of $K$. Each computed $Attention$ obtains a
head-attention. We can compute the $Attention$ in parallel, which leads to the
multi-head attention mechanism by concatenating head-attentions

$MultiHead(Q,K,V)=Concatenate(head_{1},head_{2},\ldots,head_{n}).W_{O}$

where $head_{1}$ corresponds to $Attention_{1}$. Matrix $W_{O}$ has the same
number of columns as the input matrix.

###  2.2 Pre-Trained Language Models: BERT

Training models from scratch on large datasets is impossible for most people.
Thus using pre-trained models is an inevitable trend in the development of
Artificial Intelligence. Taking into account the advantage of the weights that
can be learned from trained models, we just need to fine-tune them to suit
specific purposes. Formerly, pre-trained models in NLP have been mentioned in
many studies [9, 10, 24, 25, 26]. One of the great advantages of the
transformer’s architecture is that it allows the creation of NLP models
trained, which can be reused in downstream NLP tasks. Some of the pre-trained
language models based on transformer’s architecture have achieved state-of-
the-art results like BERT of [13] from Google, GPT of Radford and Narasimhan
[27] from Open AI and their variations. These new models can do things that
the old models can’t such as that it allows transfer learning in NLP with both
low and high-level features. Transfer learning is a combination of reusing the
architecture of pre-trained model and fine-tune parameters of the original
layers to accommodate downstream tasks.

Specifically, BERT is an easily fine-tuned pre-train word embedding on a large
unlabelled text corpus (unsupervised) which is trained based on Masked
Language Model Task and Next Sentence Prediction Task. BERT’s architecture is
built only on the Encoder part of the Transformer. The input text before
applying fine-tuning for Vietnamese in particular and other languages in
general is a combination of Token Embeddings, Segment Embedding and Position
Embeddings. If the input text consists of two or more sentences (pair-
sequence), we must add token [CLS] at the beginning of the sentence and token
[SEP] to separate the sentences.

Masked Language Model task allows us to fine-tune word representations on any
unsupervised text corpus. This task creates embeddings for the above
Vietnamese dataset. The principle of operation of model training can be
understood by predicting a missing word in the sequence instead of trying to
predict the next word in the sequence itself. A missing word is equivalent to
[MASK] token. We randomly mask 15% of the total tokens in the sequence and
predict these [MASK] tokens. Note that a missing word can be replaced by
[MASK] token 80% of the time, 10% of the time for a random token and 10% of
the time for the unchanged token.

Next Sentence Prediction (NSP) is a binary classification task applied
practically to the Question Answering (QA) task. NSP helps us to understand
the relationship between sentences. The input of the model is a pair-sequence,
which has been added tokens [CLS], [SEP]. During model training, we select 50%
of the time of the second sentence, which is the next sentence of the first
one and labeled as IsNext, while the remaining 50% of the second sentence is
randomly chosen from unrelated sentences in the dataset and labeled as
NotNext.

There are many versions of BERT with different parameters on transformer
architectures. The two most basic models are BERTBASE and BERTLARGE. In
essence, both models are the same but they are different in size.
Specifically, according to Devlin et al. [13] these models have the following
sizes

$\text{BERT}_{BASE}(L=12,H=768,A=12,Total\text{ }Parameters=110M),$
$\text{BERT}_{LARGE}(L=24,H=1024,A=16,Total\text{ }Parameters=340M)$

where $L$ is the number of layers in the Encoder part of transformer
architecture, $H$ is the hidden size and $A$ is the number of heads in multi-
head self-attention.

##  3 PhoBERT for Text Clustering

PhoBERT is Pre-trained language models for Vietnamese proposed by Nguyen and
Nguyen [15]. At the time of launch, pre-trained PhoBERT models established
state-of-the-art results in most tasks related to Vietnamese NLP. Although
BERT can be applied to many tasks like Classification, Clustering, Dependency
parsing, Sentiment analysis, Summarization text, Part-of-speech tagging,
Question Answering, Named-entity recognition and Machine translation, in this
work we only focus on clustering task to analyze our Vietnamese conversations
dataset.

PhoBERTbase and PhoBERTlarge are two versions of PhoBERT, whose architectures
are similar to the BERTBASE and BERTLARGE above. PhoBERT uses RoBERTa , which
is based on pytorch framework [28] to retrain the BERT models on new 20GB pre-
training Vietnamese dataset. Since PhoBERT architecture is based on RoBERTa,
it only trains BERT model with Masked Language Model task. Another difference
between PhoBERT and RoBERTa is fastBPE used to tokenize input sentences.
Currently there are many methods to tokenize such as Word Level Tokenizer,
Multi-Word-Level Tokenizer, Character Level Tokenizer, Subword Units Level
(BPE algorithm) Tokenizer but only BPE (Byte-Pair Encoding proposed by
Sennrich et al. [29]) achieves SOTA and is applied to most modern NLP models.

BPE is a compression technique and is adapted for word segmentation tasks.
Most words can be represented by subwords using the BPE method. It overcomes
the disadvantages of Word and Character Tokenizers, for instance, words that
do not appear in the dictionary can be represented in these subwords and the
index length of sequence output is significantly shorter than Character
Tokenizers. Code111Scripts are available at
https://github.com/rsennrich/subword-nmt of BPE algorithm to segment word into
subword units was published by [29]. For example, assume that the given
Vietnamese vocabulary is

vocab = {’x i n h ¡/w¿’: 10, ’đ ẹ p ¡/w¿’: 20, ’x i n h _ đ ẹ p ¡/w¿’: 10,

’x i n h _ x ắ n ¡/w¿’: 15, ’x ắ n ¡/w¿’: 8}

Notice that unlike English, the Vietnamese language does not use white space
to separate words because Vietnamese words can have more than one syllable.
For illustration take a simple Vietnamese sentence “Cô ấy rất xinh đẹp”
(English version is “She is very beautiful”), which can be rewritten in the
monosyllable form “Cô_ấy${}_{\texttt{ She }}$ rất${}_{\texttt{ very }}$
xinh_đẹp${}_{\texttt{ beautiful }}$”. Therefore, we can apply a Multi-Word-
Level Tokenizer on the pre-training Vietnamese dataset before going into BPE.
There are many toolkits to support word segmentation based on Multi-Word-Level
Tokenization like RDRSegmenter from VnCoreNLP [30], pyvi [31], underthesea
[32]. In the example above, tokens </w> are appended to the end of the words
to mark the end of a word in Vietnamese vocabulary. After merging the most
frequent pair at the 9th iteration, we obtain a new vocabulary as follows

vocabnew={’xinh’: 10, ’¡/w¿’: 10, ’đẹp¡/w¿’: 30, ’xinh_’: 25, ’xắn¡/w¿’: 23}

It is clear that the word ’xinh_đẹp¡/w¿’ can be represented by subwords
’xinh_’ and ’đẹp¡/w¿’ from the above vocabnew. Especially, word ’xinh_xinh’
(English meaning is pretty) is out of vocabulary words, which can also be
represented by the word pair ’xinh_’ and ’xinh’.

Figure 1: The first token [CLS] as feature of the input sentence

In order to fine-tune PhoBERT for downstream tasks, we can use library
packages such as Transformers of Hugging Face [33] and FAIRSeq of Facebook
[34] to implement. We use PhoBERTbase with 12 block sub-layers of the Encode
part to obtain the Embedding vectors as features of input sequences. More
specifically, this Embedding vector is an output vector of the first token
[CLS] from the final hidden state h (Fig.1). According to the idea of [22,
35], the vector of token [CLS] is the feature of the whole sentence for
Classification task. To verify the idea just mentioned earlier, let’s consider
the following three sentences in Table 1.

Table 1: Three Vietnamese sentences are used as examples to extract Embedding vectors with PhoBERT model. Embedding vectors of sentences A,B and C are E${}_{\texttt{[CLS]}}^{\text{A}}$, E${}_{\texttt{[CLS]}}^{\text{B}}$ and E${}_{\texttt{[CLS]}}^{\text{C}}$, respectively Sentence |  Embedding Vector | Vietnamese | English
---|---|---|---
A |  E${}_{\texttt{[CLS]}}^{\text{A}}$ |  hà_nội là thủ_đô của việt_nam |  hanoi is the capital of vietnam
B |  E${}_{\texttt{[CLS]}}^{\text{B}}$ |  thủ_đô của nước chxhcn việt_nam có tên gọi là hà_nội |  the capital of the socialist republic of vietnam is called hanoi
C |  E${}_{\texttt{[CLS]}}^{\text{C}}$ |  hôm_nay trời sẽ có mưa dông, gió mạnh |  today there will be thunderstorms and strong winds

######  Statement 1.

Assume that Embedding vectors $\text{{E}}_{\texttt{[CLS]}}^{i}$ and
$\text{{E}}_{\texttt{[CLS]}}^{j}$ represent the whole sentences $i$ and $j$,
respectively. If sentences $i$ and $j$ are similar, then the Cosine Similarity
between $\text{{E}}_{\texttt{[CLS]}}^{i}$ and
$\text{{E}}_{\texttt{[CLS]}}^{j}$ must be sufficiently larger than a certain
threshold and converges to 1 with identical sentences.

###### Proof.

The Cosine Similarity Formula is

$Cosine\\_Similarity\left(\text{{E}}_{\texttt{[CLS]}}^{i},\text{{E}}_{\texttt{[CLS]}}^{j}\right)=\frac{\text{{E}}_{\texttt{[CLS]}}^{i}\cdot\text{{E}}_{\texttt{[CLS]}}^{j}}{\parallel\text{{E}}_{\texttt{[CLS]}}^{i}\parallel\parallel\text{{E}}_{\texttt{[CLS]}}^{j}\parallel}.$
(1)

The statement above can be easily demonstrated through the example in Table 1.
As observed, sentences A, B are almost similar and have higher similarity over
sentence C. The computation of cosine similarity between Embedding vectors
E${}_{\texttt{[CLS]}}^{\text{A}}$ , E${}_{\texttt{[CLS]}}^{\text{B}}$ and
E${}_{\texttt{[CLS]}}^{\text{C}}$ is shown in Table 2. Since A, B are almost
alike, their similarity metric will be high and converge to 1 and vice versa
for C. ∎

Table 2: Compute Cosine Similarity between Embedding vectors for Table 1 with
PhoBERT Cosine Similarity(E${}_{\texttt{[CLS]}}^{i}$,
E${}_{\texttt{[CLS]}}^{j}$)
---
E${}_{\texttt{[CLS]}}^{\text{A}}$, E${}_{\texttt{[CLS]}}^{\text{A}}$ | E${}_{\texttt{[CLS]}}^{\text{A}}$, E${}_{\texttt{[CLS]}}^{\text{B}}$ | E${}_{\texttt{[CLS]}}^{\text{A}}$, E${}_{\texttt{[CLS]}}^{\text{C}}$ | E${}_{\texttt{[CLS]}}^{\text{B}}$, E${}_{\texttt{[CLS]}}^{\text{C}}$
1.0 | 0.8519334 | 0.4461445 | 0.43029878

Figure 2: The overall flow for Clustering task with PhoBERTbase model.
Starting from Vietnamese pre-training data, passing the layers Encoder0→11 to
obtain embedding vectors from the final hidden state $h$ through the Linear
projection layer. Finally using unsupervised Learning algorithm K-mean and
DBSCAN to cluster the text data.

Using the Embedding vector of token [CLS] as a feature of the whole sentence,
we adapt this idea to our Clustering task. The Clustering implementation
process with PhoBERTbase model is shown in Figure 2. After obtaining the
output embeddings of sentences with PhoBERTbase model, we use the algorithm
K-mean and DBSCAN to cluster our text data. The output embeddings of sentences
have the form as follows

$\textbf{E}_{\texttt{[CLS]}}^{i}=\textbf{h}\textbf{W}$ (2)

where $\textbf{W}\in\Re^{d,H}$ and h are projection matrices at the linear
projection layer and the final hidden state, respectively.

##  4 Clustering algorithm

K-Means is one of the most basic algorithms in unsupervised learning [3].
According to the K-Means algorithm, a set of $N$ samples
$\textbf{E}_{\texttt{[CLS]}}$ is divided into $K$ disjoint clusters ($K<N$).
Let $Y$ is the set of all label vectors for $N$ samples, i.e., each sample
$\textbf{E}_{\texttt{[CLS]}}^{i}$ has a label vector
$y_{i}=[y_{i1},y_{i2},\ldots,y_{iK}]\in Y$. If vector
$\textbf{E}_{\texttt{[CLS]}}^{i}$ belongs to cluster $k$, then $y_{ik}=1$ and
$y_{ij}=0,\forall i\neq k$. Each cluster is characterized by a cluster
“centroids”. A set of centroids is denoted
$M=[\textbf{m}_{1},\textbf{m}_{2},\ldots,\textbf{m}_{K}]$. In K-means
algorithm, the clustering problem will be reduced to the optimization for loss
function $\mathcal{L}\left(\textbf{Y},\textbf{M}\right)$.

$\mathcal{L}\left(\textbf{Y},\textbf{M}\right)=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{K}y_{ij}\left|\left|\textbf{E}_{\texttt{[CLS]}}^{i}-\textbf{m}_{j}\right|\right|^{2}$
(3)

Data: $\textbf{E}_{\texttt{[CLS]}}$

Input : $\varepsilon$, MinPts

1 Input : $metric$ // calculating distance between data points

Input : cluster

2 foreach  _$p$ in $\textbf{E}_{\texttt{[CLS]}}$_ do // linear scan for all
data points in $\textbf{E}_{\texttt{[CLS]}}$

3 if  _cluster( $p$) $!=$ unassigned_ then continue;

4 // p is unassigned to any cluster or noise Neighbors $N\leftarrow$
find_neighbors($\textbf{E}_{\texttt{[CLS]}}$, $metric$, $p$, $\varepsilon$)

5 if  _$|N| <minPts$_ then

6 cluster($p$) $\leftarrow$ noise

7 continue

8 end if

9 $c\leftarrow$ create a new cluster

10 cluster($p$) $\leftarrow$ $c$

11 S $\leftarrow$ $N\setminus\\{p\\}$

12 foreach  _$q$ in S_ do

13 if  _cluster($q$)$==$noise_ then cluster($q$) $\leftarrow$ $c$;

14 if  _cluster( $q$) $!=$ unassigned_ then continue;

15 Neighbors $N\leftarrow$ find_neighbors($\textbf{E}_{\texttt{[CLS]}}$,
$metric$, $q$, $\varepsilon$)

16 cluster($q$) $\leftarrow$ $c$

17 if  _$|N| <minPts$_ then continue;

18 S $\leftarrow$ $N\cup\textbf{{S}}$

19 end foreach

20

21 end foreach

Algorithm 1 Pseudocode of original DBSCAN algorithm for our data

Along with that, we also apply DBSCAN technique to cluster data points [5]. In
real-life data, DBSCAN can work well for nonconvex clusters with arbitrary
shapes and noises. DBSCAN algorithm focuses on radius eps–$\varepsilon$ and
the minimum number of neighbors required to create a cluster $minPts$. Radius
eps defines a circle for each point to determine its neighbors. A point
becomes core point if the circle surrounding this point with radius eps
contains more than MinPts neighbors. In case the number of neighbor points is
less than MinPts, the core point is the border points. On the other hand, a
point without any neighbors within radius eps is called noise. The
relationship state of two points in DBSCAN can be direct density reachable,
density reachable or density connected. A point is called direct density
reachable for $C_{i}$ point if and only if it lies within the circle centered
core point $C_{i}$. If a core point is connected unidirectionally to any other
core point through a chain of core points, there is a density reachable state
between them. In case there are two points, which are density reachable from
the same point, they are density connected states. Pseudocode describing
DBSCAN clustering algorithm [36, 5] is shown in Algorithm 1.

Unlike the evaluation metrics for supervised learning algorithms, the
evaluation of clustering performance can be applied to datasets with known or
unknown ground truth labels. If the ground truth class assignment of dataset
is known, we use entropy-based measure, V-measure proposed in [19] to evaluate
clustering performance for our sample dataset. The sample dataset will be
described in detail in Section 5. Based on the conditional entropy analysis of
two terms of homogeneity and completeness, V-measure is a harmonic mean
function of those terms and can be calculated as follows

$v(\beta,h,c)=\frac{(1+\beta)\times h\times c}{(\beta\times h)+c},$ (4)

where $h=\frac{1-H(C|K)}{H(C)}$ and $c=\frac{1-H(K|C)}{H(K)}$ are homogeneity
and completeness respectively. The conditional entropy $H(K|C)$ and entropy
$H(K)$ are symmetric. In formula (4), $\beta$ weight represents the
contributions of homogeneity or completeness and the default value of $\beta$
is equal to 1.

Unfortunately, in fact, we don’t know anything about the ground truth classes
for document clustering task. Thus, we can evaluate clustering performance
based on the partition obtained from clustering techniques and two types of
proximities, which are similarity or dissimilarity between objects. Suggested
in [20] Silhouette is a typical evaluation for this case. Besides providing a
graphical overview of the partitional clustering (silhouette plot), Silhouette
also allows evaluating clustering validity based on the average silhouette
width. From those analyses, we can obtain a suitable number of clusters for
the K-means algorithm. The Silhouette Coefficient $s(i)$ for object $i$ has
the form

$s(i)=\frac{b(i)-a(i)}{max\\{a(i),b(i)\\}}.$ (5)

Here $a(i)$ is the average dissimilarity of object $i$ to the remaining
objects in the same cluster and $b(i)$ is the average dissimilarity of $i$ to
all objects of the next nearest cluster. The value of Silhouette Coefficient
is in the range $[-1,+1]$, where near -1 indicates that the object for
incorrect clustering and vice versa for +1. The value around 0 represents
overlapping clusters.

##  5 Experiment

We apply K-means and DBSCAN algorithm with PhoBERTbase to cluster our Facebook
Vietnamese conversations dataset (FVnC) and Sample dataset. In order to
implement clustering task, we use PhoBERTbase with the Transformers package of
Hugging Face and Scikit-learn library [33, 37]. Furthermore, we also search
the optimal parameters for the clustering algorithms in this article. The
clustering results will be used to build Intents for chatbot later.

### Clustering task dataset

We evaluate our approach on Facebook Vietnamese conversations dataset. There
are plenty of free tools or extensions to download conversations from a
personal page because it is quite simple. However collecting conversations
from public page is more difficult so most tools or extensions to carry out
this task are paid. In order to acquire this dataset, we created a tool named
NTUCrawler222Tool is available at https://archive.org/download/NTUCrawler for
scraping conversations from a Facebook messenger page of our University. This
tool is written in Python language and based on Facebook’s Graph API platform
to get messages. It has two versions, one is linux executable (run on Ubuntu
distribution) and the other is Windows executable. The UI of the Windows
version is shown in Figure 3. NTUCrawler requires users to provide four
parameters, start time and finish time to get data, PageID and Token of page.
Downloaded dataset contains 8000 conversations with more than 150 thousand raw
text sentences of clients and admins of a Facebook page in the six-month
period of the year 2020. The contents of the conversations are FAQ (frequently
asked questions), which are related to information already published on the
university website, for example, tuition fees, insurance, dormitory, English
language test, course registration, etc.

Figure 3: User interface of NTUCrawler on Windows

### Data preprocessing

The raw downloaded dataset above must be preprocessed. Some punctuations “!”,
“"”, “?”, “’”, “:” are removed from the dataset. Further, we eliminate stop
words that have no value or no meaning to the NLP model such as “dạ”, “vâng”,
“chào ad”, “vâng ạ”, “alo”, “ừ”, “vậy”, “ok”, “nhé”. We also filter out
duplicate sentences or those that contain less than 4 words. After data
preprocessing, 44846 text sentences are obtained and it is our desired FVnC
dataset. From FVnC dataset, we randomly selected 95 sample text sentences
(0.2% of FVnC dataset size) to form the Sample dataset. Besides that, we also
use another sub dataset which is called the wiki dataset. This dataset
contains 396 text sentences of articles on 5 topics collected from Wikipedia.
The Wikipedia library was applied to access and parse data from Wikipedia. The
reason we use those Sample dataset and wiki sub dataset is that we can
evaluate the effectiveness of applying PhoBERTbase for downstream task
(clustering). Label assignment to sample text sentences and analyzing the
number of clusters will cost less in our task. In addition, clustering on the
small Sample dataset not only helps to easily evaluate clustering performance
but also saves time compared to the original FVnC dataset. These sample text
sentences were completed manual labeling by us and divided into 3 classes,
which describe questions between users (students) and admin (university) about
information related to insurance, dormitory and English language test.
Specifically, in order to specify the lables of classes, we rely on the
experience of the specialists of the training department, who are responsible
for answering students’ questions directly or via social platforms. In their
opinion, first-year students of our university who are often concerned with
insurance, dormitory and English language test. During data collection for
each label, we carefully selected the sentences in the dataset that matched
the recommendations of the specialists. Those 3 labels are one of the intents
used to train the chatbot. Table 3 shows the details of the two datasets.

Table 3: Brief description of the datasets used for PhoBERTbase

Dataset | Label | Description | Tasks |  $\textbf{E}_{\texttt{[CLS]}}^{i}$ size
---|---|---|---|---
FVnC | unknown |  clean downloaded dataset |  Silhouette evaluation, clustering | 44846 $\times$ 768
| class 1 |  features relating to insurance |  | 31 $\times$ 768
Sample | class 2 |  features relating to english language test |  clustering, clustering performance evaluation | 38 $\times$ 768
| class 3 |  features relating to dormitory |  | 26 $\times$ 768

According to [15], input text must be already word-segmented before going
through the BPE algorithm. We use “pyvi” toolkit of Tran [31] to perform word
segmentation in our datasets. After passing the fastBPE step, we have the
index of tokens and attention masks for the text data. Taking them through
PhoBERT’s architecture leads to output embeddings of
$\textbf{E}_{\texttt{[CLS]}}$. From this step, we use these embeddings as
feature vectors to cluster text data.

### Dealing with varying length

Because the length of the sentences in FVnC dataset is different, we have to
use padding to make sure the input texts have the same length.

Figure 4: The distribution of the sentence lengths in FVnC dataset

In particular, the maximum sequence length of PhoBERT is 256. We truncate
sentences with padding length less than 256 tokens. To choose the optimal
padding length for all sentences, we can analyze the distribution of sentence
lengths in Figure 4. Based on the above distribution, most sentences have
lengths of less than 33 words. Thus, we decide that the padding length is
equal to 33 in all datasets.

### Parameter optimization

When using the DBSCAN technique, we need to take care of two parameters MinPts
and $\varepsilon$. Choosing these two parameters is not easy. Their influence
is very large on the clustering results. There is no way to accurately
determine the parameter MinPts. However, there are several ways of choosing
MinPts which have been proposed in [5, 6]. Besides, the value of MinPts must
also depend on domain knowledge and the data distribution observation. So we
derive that MinPts should be greater than the number of dimensionality of
feature vector $\textbf{E}_{\texttt{[CLS]}}\left(44846\times 768\right)$. As
we can see, the number of dimensions of $\textbf{E}_{\texttt{[CLS]}}$ is too
large (768), which affects computation time and cost for large data sets.
Therefore, we use the dimensionality reduction method to reduce
$\textbf{E}_{\texttt{[CLS]}}$ to lower-dimensional while it retains most of
the original information. Principal Component Analysis (PCA) and t-distributed
Stochastic Neighbour Embedding (t-SNE) are common techniques for data
dimensionality reduction. PCA relies on eigenvalues and eigenvectors of
$\textbf{E}_{\texttt{[CLS]}}$ to reduce the original data to a specific number
of dimensions (commonly known as principal components) but it still ensures a
threshold of allowable variance. If we use PCA then MinPts can be selected as
follows

$MinPts\geq\textit{principal components}+1$

Figure 5: Data Points are sorted ascending by the average distance to MinPts
nearest neighbors. (a) Calculating on Sample dataset; (b) Calculating on FVnC
dataset

Parameter $\varepsilon$ can be found from a K-Distance graph, which is based
on the average distance between objects and its MinPts nearest neighbors [7].
The K-Distance graph with MinPts=3 for FVnC and Sample dataset is shown in
Figure 5. The blue solid curve and red dashed curve correspond to the average
distance of objects to MinPts nearest neighbors which are sorted in ascending
order for FVnC and Sample dataset respectively. Usually, a point at the
position with the largest slope change in K-Distance graph or what we
popularly call the “knee/elbow” of the graph is the optimal value of parameter
$\varepsilon$ [38]. Especially, the greatest slope change zones are
highlighted in Fig.5(a) and Fig.5(b) for specific datasets. In order to take
exactly the point mentioned above or the “knee point” of the graph, the
kneedle algorithm is considered in our work [39]. The knee point obtained from
the kneedle algorithm is determined by the intersection of the specific data
curve with the vertical straight line in Figure 6. The optimum values for
parameter $\varepsilon$ are 0.57 and 0.12 in the case of Fig.6(a) and
Fig.6(b), respectively.

Figure 6: Determining knee points at the greatest slope change zones using the
kneedle algorithm for Figure 5. (a) Solid blue curves – the greatest slope
change zone of Sample dataset; (b) Dashed red curve – the greatest slope
change zone of FVnC dataset.

Data : $\textbf{E}_{\texttt{[CLS]}}$, label ;

// label: points labels (unassigned or assigned)

1

Input : K-Distance

Input : $n\\_components$, step ;

// step: minimum distance between points

2

3 Output : Index // the appropriate $\varepsilon$ and $MinPts$ values for
DBSCAN

Initial : $\text{Max}_{Silhouette}\leftarrow-1$ ;
$\text{Max}_{Vmeasure}\leftarrow 0$ ; $\text{Max}\leftarrow-0.5$

4 slope $\leftarrow$ CalculatingSlope(K-Distance)// the greatest slope change
zone

5 nearest_neighbors $\leftarrow$ arange($n\\_components+1$, $2\times
n\\_components+1$, step)

6 foreach  _$\varepsilon$ in slope_ do

7 foreach  _MinPts in nearest_neighbors_ do

8
$p\leftarrow\texttt{PCA}\left(\textbf{E}_{\texttt{[CLS]}},n\\_components\right)$

9
$ClusterAssignment\leftarrow\texttt{DBSCAN}\left(p,\textit{MinPts},\varepsilon\right)$

10 $\text{SilCoeff}\leftarrow\texttt{SilhouetteScore}(p,ClusterAssignment)$

11 if  _label $=$ unassigned_ then

12 if  _$\text{SilCoeff} >\text{Max}_{Silhouette}$_ then

13 $\text{Max}_{Silhouette}\leftarrow\text{SilCoeff}$

14 $\text{Index}\leftarrow(\varepsilon,MinPts)$

15

16 else

17
$\text{VScore}\leftarrow\texttt{VMeasureScore}\left(label,ClusterAssignment\right)$

18 if  _$\left(\text{SilCoeff}+\text{VScore}\right)/2 >\text{Max}$_ then

19 $\text{Max}\leftarrow\left(\text{SilCoeff}+\text{VScore}\right)/2$

20 $\text{Index}\leftarrow(\varepsilon,MinPts)$

21

22 end if

23

24 end foreach

25

26 end foreach

Algorithm 2 Pseudocode of the proposed technique to find appropriate
$\varepsilon$ and MinPts

However, optimizing parameter $\varepsilon$ by choosing a fixed value of knee
point in some cases does not lead to good clustering efficiency. Specifically
in the test of the good separability between clusters in the subsection below,
Silhouette score is quite low. This could lead to objects being assigned to
the wrong clusters. With the expectation of improved clustering performance,
in this work, we propose another technique based on the combination of
K-Distance graph and clustering performance evaluations to find the right
optimal value $\varepsilon$. Besides observing the greatest slope change zone
of the line from a pair of points on the K-Distance graph, the maximum value
of the Silhouette coefficient and V-measure score are also considered. The
pseudocode for our technique is given in Algorithm 2. The sequence of MinPts
is taken from $\textit{principal components}+1$ to $2\times\textit{principal
components}+1$ and incremented by a step of the minimum distance between data
points. Gridsearch technique is applied in algorithm 2 for the sequence of
MinPts and the greatest slope change zone in K-Distance graph. At the position
where the Silhouette Coefficient is maximum, we obtain the optimal pair of
values $\left(\varepsilon,MinPts\right)$ for unlabeled data. For the labeled
Sample dataset, the search of the optimal values
$\left(\varepsilon,MinPts\right)$ is based on the greatest mean of V-measure
and Silhouette evaluation.

### Result

Figure 7: The graphical overview of the partitional clustering for 2, 3, 4
clusters using K-Means algorithm on the Sample dataset. (a) The silhouette
plot; (b) The visualization of the clustered data for corresponding silhouette
plots.

In this part, our first task is to cluster text documents and evaluate
clustering performance on the Sample dataset. As a consequence, analyzing on
the Sample dataset will be generalized to the general dataset FVnC such as
choosing the number of clusters using silhouette analysis for K-Means
algorithm. As mentioned above, the Sample dataset has three labeled clusters
(see Table 3). We use Silhouette evaluation to confirm that the Sample dataset
has exactly three clusters and the way to choose the right number of clusters
when using it. Besides considering average silhouette scores, the silhouette
plot is also an important factor in determining the number of clusters. Figure
7 represents the graphical overview of the partitional clustering using
K-Means algorithm on Sample dataset. X-axis of Fig.7(a) corresponds to the
silhouette coefficient values and the vertical dashed line is the average of
the silhouette coefficients of data points. Clusters corresponding to
silhouette plots in Fig.7(b) are visualized in two dimension space using PCA.
The maximum value of average silhouette scores is close to 0.481 for 3
clusters. On the other hand, the thickness of the silhouette plot for clusters
is similar. The analysis outlined above fits the facts in Table 3 for the
Sample dataset.

Table 4: Evaluating feature extraction efficiency of models through V-measure scores for clustering task on the Sample dataset and wiki dataset. Approach | Dataset
---|---
Sample | wiki
FastText (cc.vi.300) | 0.64 | 0.27
GloVe (glove.6B.100d) | 0.49 | 0.26
GloVe (vncorpus.3B.100d) | 0.61 | 0.27
BERTbase Uncased | 0.11 | 0.19
BERTbase Multilingual Uncased | 0.36 | 0.43
DistilBERTbase Multilingual Cased | 0.37 | 0.12
GPT-2 | 0.04 | 0.04
PhoBERTbase | 0.76 | 0.62

Moreover, in order to evaluate the efficiency of the PhoBERTbase model in
feature extraction for clustering tasks, we compared it with the other models
such as BERTbase Uncased, BERTbase Multilingual Uncased [13], DistilBERTbase
Multilingual Cased [14], GPT-2 [12] based on V-measure score, among which
GPT-2 and BERTbase Uncased models do not support Vietnamese. We also use
traditional approaches like FastText, GloVe to compare with transformer
models. Both models are applicable to Vietnamese language. However, both of
these models are commonly used for learning word representations. Thus, to
obtain sentence embedding, our method is averaging the word embeddings of all
the words in the sentence. Especially, in order to make the comparison better,
we also retrained the GloVe model on a new Vietnamese 20GB dataset based on
the improvement proposed in [40]. Pre-trained word vectors vncorpus.3B.100d of
new GloVe model corresponds to the corpus 3B tokens, 1.3M vocab, and 100d
vectors and 1.17 GB download. Comparison results on the Sample dataset and
wiki dataset using K-Means algorithm with 3 clusters and 5 clusters
respectively are presented in Table 4. Through the obtained results, the
models that support Vietnamese achieve significantly better V-measure scores
than the rest of the models, especially the PhoBERTbase model, which obtained
the highest V-measure score on both the Sample dataset and wiki dataset (0.76
and 0.62 respectively). The comparisons obtained are in good agreement with
the pointed theoretical and experimental works. Therefore, we use PhoBERTbase
model to extract features for the FVnC dataset, which leads to the output
embeddings serving the clustering task.

Table 5: Experiment to find the optimal values of parameters $\varepsilon$,
MinPts for DBSCAN algorithm on the Sample dataset.

Approach | eps ($\varepsilon$) | MinPts | V-measure score | Silhouette score | Average | N-Clusters
---|---|---|---|---|---|---
Algorithm 2 | 0.84 | 4 | 0.56 | 0.28 | 0.42 | 5
0.78 | 4 | 0.52 | 0.27 | 0.4 | 5
0.72 | 4 | 0.5 | 0.22 | 0.36 | 6
0.76 | 3 | 0.47 | 0.22 | 0.34 | 4
0.54 | 5 | 0.37 | 0.0 | 0.18 | 8
Knee point | 0.57 | 3 | 0.42 | 0.15 | 0.29 | 9
0.57 | 4 | 0.41 | 0.1 | 0.26 | 8
0.57 | 5 | 0.41 | 0.05 | 0.23 | 7

Figure 8: DBSCAN clustering on the Sample dataset

For DBSCAN clustering method on the Sample dataset, we need to find the
optimal parameters $\varepsilon$, MinPts using Algorithm 2 and kneedle
algorithm. Based on K-Distance graph of Sample dataset in Fig.5(a) and knee
visualization in Fig.6(a), the greatest slope zone is selected in the range of
$[0.5,0.85]$ and the knee point value is 0.57. Due to “principal components”
equal to 2 in PCA dimensionality reduction, the value of MinPts will be
selected from 3 to 5. Some experimental results in finding the optimal
parameters are shown in Table 5. As shown in the table, $\varepsilon$, MinPts
and the average of V-measure and Silhouette score equal to 0.84, 4 and 0.42
respectively are the best choice to cluster data points. Results in Table 5
also show that our approach – Algorithm 2 gives better clustering performance
than kneedle algorithm on the Sample dataset. Using DBSCAN algorithm with the
obtained parameters, 5 clusters are formed in Figure 8. In which cluster “-1”
contains noisy objects in Fig.8(a), noisy objects are removed in Fig.8(b) and
Silhouette plot of denoised Sample dataset in Fig.8(c). After removing noise,
our clustering result has 4 clusters while the Sample dataset has only 3
clusters. Based on actual observation in Fig.8(b) and Silhouette plot in
Fig.8(c), cluster 0 and 1 must be merged into one cluster.

######  Remark 1.

DBSCAN is a density-based spatial clustering algorithm. The biggest
disadvantage of DBSCAN is working in cases of varying density clusters. From
the Silhouette plot in Fig.8(c), there are two separate clusters of very
different density comparing to the other two clusters. Two low-density
clusters (cluster 0 and 1) are a matter of concern to us for predicting the
number of clusters. Obviously, we hope that our clusters will be more equally
distributed to choose the number of clusters more precisely. In order to solve
this problem, we propose some solutions as follows
(i) Based on the thickness of the silhouette plot, the silhouette plots for
clusters 0 and 1 are much smaller than clusters 2, 3. In particular, clusters
0 and 1 are close neighbors. Therefore, we have reasons to merge these small
two clusters into a larger cluster. As a result, 3 clusters are the correct
number of clusters when using DBSCAN clustering on the Sample dataset. After
all, we should combine the selection of the optimal parameters and the width
of Silhouette plots for clusters when using DBSCAN algorithm to obtain the
most accurate number of clustering.
(ii) On the other hand, we can use a criterion for determining the minimum
number of core points in a cluster. If at least two clusters are close to each
other and have a smaller the minimum number of core points than the specified
criteria, we can merge them into a larger cluster.

Table 6: Silhouette scores using K-Means algorithm on original FVnC and
denoised FVnC datasets for several different cluster numbers were found from
DBSCAN algorithm.

DBSCAN |  K-Means Silhouette Scores
---|---
eps $\varepsilon$ | MinPts | N-Clusters | Original FVnC | Denoised FVnC
0.14 | 22 | 10 | 0.3359 | 0.3436
0.14 | 25 | 12 | 0.3359 | 0.3460
0.13 | 12 | 15 | 0.3310 | 0.3368
0.14 | 19 | 16 | 0.3294 | 0.3326
0.13 | 16 | 20 | 0.3267 | 0.3314
0.12 | 30 | 23 | 0.3252 | 0.3418
0.11 | 26 | 29 | 0.3225 | 0.3343
0.1 | 14 | 53 | 0.3196 | 0.3257
0.09 | 12 | 75 | 0.3196 | 0.3275
0.1 | 3 | 81 | 0.3217 | 0.3205
0.09 | 4 | 102 | 0.3220 | 0.3229
0.09 | 3 | 128 | 0.3207 | 0.3228
0.08 | 4 | 136 | 0.3212 | 0.3235

Finally, we perform the main task, which is text clustering for the FVnC
dataset. Note that, the ground truth classes of the FVnC dataset are unknown.
In order to estimate the number of specific clusters for K-Means clustering
method, we take full advantage of DBSCAN algorithm. In the same way, we get
the number of clusters from Algorithm 2. Specifically, we choose the number of
clusters through Silhouette scores with large values corresponding to each
pair of parameters $\varepsilon$ and MinPts. For FVnC dataset, the optimal
values of $\varepsilon$ can be found in Fig.5(b). After applying DBSCAN
algorithm, the number of clusters actually obtained must be subtracted by 1
for noisy objects (points labeled -1). Before using K-Means algorithm to
perform clustering, these noisy objects are removed from FVnC dataset, which
leads to the form of the corresponding denoised FVnC datasets. Lastly, we
obtain clustering results of FVnC dataset without noise from the K-Means
method based on the Silhouette evaluation. Silhouette scores test on original
FVnC and denoised FVnC datasets for several different cluster numbers can be
examined in Table 6. In case of using $\varepsilon=0.14$ and
$\textit{MinPts}=25$, the best obtained Silhouette scores for the original
FVnC and denoised FVnC dataset are 0.3359 and 0.3460 respectively.

Figure 9: Silhouette plots and its clustering visualization

The Silhouette plots and its clustering visualization of the best Silhouette
scores can be seen in Figure 9. By comparing the Silhouette plots, we believe
that the clustering result on the denoised FVnC dataset (see Fig.9(b)) is
better than the original FVnC dataset (see Fig.9(a)), because there are some
noise points outside the clusters and the size of the silhouette plots of
clusters 3, 4 are wide fluctuation in the visualization of Fig.9(a). With the
number of clusters received from the DBSCAN algorithm, the K-Means algorithm
gives very good clustering results for denoised FVnC datasets. After all, the
clusters obtained from the clustering process are considered as big intents to
help us build data for chatbot. For this reason, we can save time and effort
and build chatbot faster.

##  6 Conclusion

In this work, we research transformer architecture as well as pre-trained
language models such as BERT and PhoBERT. We also cover how to apply PhoBERT
to our Facebook Vietnamese conversations dataset (FVnC). Furthermore, we have
also built a tool to crawl conversations from a Facebook Messenger Page. After
extracting embedding vectors at the final hidden layer, we use the
unsupervised learning algorithm K-means and DBSCAN to cluster text data.
V-measure score and Silhouette score are used to evaluate the performance of
clustering algorithms. A GridSearch algorithm that combines these two
clustering evaluations is also proposed to find optimal parameters for the
DBSCAN algorithm. The algorithm proposed by us obtained better clustering
performance than kneedle algorithm through experimentations based on V-measure
scores and Silhouette score on the Sample and FVnC datasets. In addition, we
compare the efficiency of the PhoBERTbase model in feature extraction for
clustering tasks with the other models. PhoBERTbase model achieves the best
V-measure score on the Sample dataset and wiki dataset. We apply the K-Means
clustering method with the number of clusters received from the DBSCAN
algorithm to cluster the FVnC dataset. Topics obtained from clustering are
similar to intents in building chatbot. From a pre-analysis data screening
perspective, clustering results are valuable for building stories in our
chatbot. Thanks to the implementation of this clustering, we save a lot of
time and effort to build data and storylines for training chatbot.

## Acknowledgments

This work was partially supported by Nha Trang University (project
TR2020-13-42). The authors thank Hien Thao Le for proofreading our manuscript
and fruitful discussions. We also address special thanks to the reviewers for
their helpful comments and suggestions.

## References

  * Jung [2018] A. Jung, “A gentle introduction to supervised machine learning,” _CoRR_ , vol. abs/1805.05052, pp. 6–7, 2018.
  * Kanungo et al. [2002] T. Kanungo, D. Mount, N. Netanyahu, C. Piatko, R. Silverman, and A. Wu, “An efficient k-means clustering algorithm: analysis and implementation,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , vol. 24, no. 7, pp. 881–892, 2002.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011. [Online]. Available: https://scikit-learn.org/stable/modules/clustering.html#k-means
  * MacQueen [1967] J. B. MacQueen, “Some methods for classification and analysis of multivariate observations,” in _Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability_ , L. M. L. Cam and J. Neyman, Eds., vol. 1. University of California Press, 1967, pp. 281–297.
  * Ester et al. [1996] M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “A density-based algorithm for discovering clusters in large spatial databases with noise,” in _Proc. of 2nd International Conference on Knowledge Discovery and_ , 1996, pp. 226–231.
  * Sander et al. [1998] J. Sander, M. Ester, H. Kriegel, and X. Xu, “Density-based clustering in spatial databases: The algorithm gdbscan and its applications,” _Data Mining and Knowledge Discovery_ , vol. 2, pp. 169–194, 1998.
  * Gaonkar and Sawant [2013] M. Gaonkar and K. Sawant, “Autoepsdbscan : Dbscan with eps automatic for large dataset,” vol. 2, 2013, pp. 11–16.
  * Salton and McGill [1986] G. Salton and M. J. McGill, “Introduction to modern information retrieval,” 1986\.
  * Mikolov et al. [2013] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” in _1st International Conference on Learning Representations, ICLR 2013, Scottsdale, Arizona, USA, May 2-4, 2013, Workshop Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2013.
  * Pennington et al. [2014] J. Pennington, R. Socher, and C. Manning, “GloVe: Global vectors for word representation,” in _Proceedings of the 2014 Conference on Empirical Methods in Natural Language Processing (EMNLP)_. Doha, Qatar: Association for Computational Linguistics, Oct. 2014, pp. 1532–1543.
  * Bojanowski et al. [2016] P. Bojanowski, E. Grave, A. Joulin, and T. Mikolov, “Enriching word vectors with subword information,” _arXiv preprint arXiv:1607.04606_ , 2016.
  * Radford et al. [2019] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, and I. Sutskever, “Language models are unsupervised multitask learners,” 2019.
  * Devlin et al. [2019] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “BERT: Pre-training of deep bidirectional transformers for language understanding,” in _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_. Minneapolis, Minnesota: Association for Computational Linguistics, Jun. 2019, pp. 4171–4186.
  * Sanh et al. [2019] V. Sanh, L. Debut, J. Chaumond, and T. Wolf, “Distilbert, a distilled version of bert: smaller, faster, cheaper and lighter,” _ArXiv_ , vol. abs/1910.01108, 2019.
  * Nguyen and Nguyen [2020] D. Q. Nguyen and A. T. Nguyen, “PhoBERT: Pre-trained language models for Vietnamese,” in _Findings of the Association for Computational Linguistics: EMNLP 2020_ , 2020, pp. 1037–1042.
  * Gencoglu [2019] O. Gencoglu, “Deep representation learning for clustering of health tweets,” _CoRR_ , vol. abs/1901.00439, 2019.
  * Pugachev and Burtsev [2021] L. Pugachev and M. Burtsev, “Short text clustering with transformers,” _arXiv preprint arXiv:2102.00541_ , 2021.
  * Sia et al. [2020] S. Sia, A. Dalmia, and S. J. Mielke, “Tired of topic models? clusters of pretrained word embeddings make for fast and good topics too!” _arXiv preprint arXiv:2004.14914_ , 2020.
  * Rosenberg and Hirschberg [2007] A. Rosenberg and J. Hirschberg, “V-measure: A conditional entropy-based external cluster evaluation measure,” in _Proceedings of the 2007 Joint Conference on Empirical Methods in Natural Language Processing and Computational Natural Language Learning (EMNLP-CoNLL)_. Prague, Czech Republic: Association for Computational Linguistics, Jun. 2007, pp. 410–420.
  * Rousseeuw [1987] P. J. Rousseeuw, “Silhouettes: A graphical aid to the interpretation and validation of cluster analysis,” _Journal of Computational and Applied Mathematics_ , vol. 20, pp. 53–65, 1987.
  * Vaswani et al. [2017] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. u. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in Neural Information Processing Systems_ , I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, Eds., vol. 30. Curran Associates, Inc., 2017\.
  * Sun et al. [2019] C. Sun, X. Qiu, Y. Xu, and X. Huang, “How to fine-tune bert for text classification?” in _Chinese Computational Linguistics_ , M. Sun, X. Huang, H. Ji, Z. Liu, and Y. Liu, Eds. Cham: Springer International Publishing, 2019, pp. 194–206.
  * Chung et al. [2014] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” in _NIPS 2014 Workshop on Deep Learning, December 2014_ , 2014.
  * Le and Mikolov [2014] Q. Le and T. Mikolov, “Distributed representations of sentences and documents,” in _Proceedings of the 31st International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, E. P. Xing and T. Jebara, Eds., vol. 32, no. 2. Bejing, China: PMLR, 22–24 Jun 2014, pp. 1188–1196.
  * Joulin et al. [2016] A. Joulin, E. Grave, P. Bojanowski, and T. Mikolov, “Bag of tricks for efficient text classification,” _CoRR_ , vol. abs/1607.01759, 2016.
  * Peters et al. [2018] M. E. Peters, M. Neumann, M. Iyyer, M. Gardner, C. Clark, K. Lee, and L. Zettlemoyer, “Deep contextualized word representations,” 2018.
  * Radford and Narasimhan [2018] A. Radford and K. Narasimhan, “Improving language understanding by generative pre-training,” 2018.
  * Liu et al. [2019] Y. Liu, M. Ott, N. Goyal, J. Du, M. Joshi, D. Chen, O. Levy, M. Lewis, L. Zettlemoyer, and V. Stoyanov, “Roberta: A robustly optimized bert pretraining approach,” 2019.
  * Sennrich et al. [2016] R. Sennrich, B. Haddow, and A. Birch, “Neural machine translation of rare words with subword units,” in _Proceedings of the 54th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_. Berlin, Germany: Association for Computational Linguistics, Aug. 2016, pp. 1715–1725.
  * Vu et al. [2018] T. Vu, D. Q. Nguyen, D. Q. Nguyen, M. Dras, and M. Johnson, “VnCoreNLP: A Vietnamese natural language processing toolkit,” in _Proceedings of the 2018 Conference of the North American Chapter of the Association for Computational Linguistics: Demonstrations_. New Orleans, Louisiana: Association for Computational Linguistics, Jun. 2018, pp. 56–60.
  * Tran [2020] V.-T. Tran, “Python vietnamese toolkit, https://pypi.org/project/pyvi/,” 2020. [Online]. Available: https://pypi.org/project/pyvi/
  * Anh et al. [2018] V. Anh, B. N. Anh, and D. V. Dung, “Open-source vietnamese natural language process toolkit,” 2018. [Online]. Available: https://github.com/undertheseanlp/word_tokenize
  * Wolf et al. [2020] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz, J. Davison, S. Shleifer, P. von Platen, C. Ma, Y. Jernite, J. Plu, C. Xu, T. L. Scao, S. Gugger, M. Drame, Q. Lhoest, and A. M. Rush, “Transformers: State-of-the-art natural language processing,” in _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_. Online: Association for Computational Linguistics, Oct. 2020, pp. 38–45.
  * Ott et al. [2019] M. Ott, S. Edunov, A. Baevski, A. Fan, S. Gross, N. Ng, D. Grangier, and M. Auli, “fairseq: A fast, extensible toolkit for sequence modeling,” in _Proceedings of NAACL-HLT 2019: Demonstrations_ , 2019.
  * Neeraj [2020] T. Neeraj, “Feature-based approach with bert,” _trishalaneeraj.github.io_ , 2020.
  * Schubert et al. [2017] E. Schubert, J. Sander, M. Ester, H. Kriegel, and X. Xu, “DBSCAN revisited, revisited: Why and how you should (still) use DBSCAN,” _ACM Trans. Database Syst._ , vol. 42, no. 3, pp. 19:1–19:21, 2017.
  * Scikit-learn [2011] Scikit-learn, “Clustering,” 2011. [Online]. Available: https://scikit-learn.org/stable/modules/clustering.html#k-means
  * Rahmah and Sitanggang [2016] N. Rahmah and I. S. Sitanggang, “Determination of optimal epsilon (eps) value on DBSCAN algorithm to clustering data on peatland hotspots in sumatra,” _IOP Conference Series: Earth and Environmental Science_ , vol. 31, p. 012012, jan 2016.
  * Satopaa et al. [2011] V. Satopaa, J. Albrecht, D. Irwin, and B. Raghavan, “Finding a" kneedle" in a haystack: Detecting knee points in system behavior,” in _2011 31st international conference on distributed computing systems workshops_. IEEE, 2011, pp. 166–171.
  * Nguyen [2020] T. H. Nguyen, “Analyze the effects of weighting functions on cost function in the glove model,” _arXiv preprint arXiv:2009.04732_ , 2020.