# A revision of the forces exerted in the Rayleigh regime by a tighlty focused
optical vortex tweezer

Khalid Aloufi, Vasileios E. Lembessis and Omar M. Aldossary  Quantum
Technology Group, Department of Physics and Astronomy, College of Science,
King Saud University, Riyadh 11451, Saudi Arabia

###### Abstract

Optical tweezers use light from a tightly focused laser beam to manipulate the
motion of tiny particles. When the laser light is strongly focused, but still
paraxial, its e/m field is characterized by a longitudinal component which is
of magnitude comparable to the transverse ones and which has been ignored in
the theoretical analysis of the tweezing forces. In our work we revise the
calculations of the various components of the radiation pressure force, within
the limits of Rayleigh regime or dipole approximation, in the case where a
tiny particle interacts, in free space, with a circularly polarized optical
vortex beam, by taking into account this ignored field term. We show that this
term is responsible for considerable modifications in the magnitude of the
various components, moreover and also for the appearance of terms involving
the coupling of the spin angular momentum (SAM) and the orbital angular
momentum (OAM) of the photons of the vortex beam. We compare our findings with
the ones taken ignoring the longitudinal field component.

††preprint: APS/123-QED

##  I  Introduction

The origin of optical tweezers is an experiment conducted in 1970 by Ashkin
where he succeeded in manipulating the motion of microscopic rubber beads by
means of a laser beam [1]. The optical tweezer was invented in 1986 by Ashkin
and Steven Chu when they managed to move, and at the same time trap, small
dielectric beads [2]. Since then we have witnessed a very broad range of
applications of optical tweezers in biology, micromechanics, modelling of
thermodynamic systems, collidal systems, and microchemistry [3].

Optical tweezers are based on the forces that laser light exerts on tiny
particles. Thus, the exact knowledge of these forces is prerequisite for any
study involving the particles’ dynamics in optical tweezers. When the
dimensions of a particle are smaller than the wavelength of laser light then
the motion of it can be described as a result of the action of radiation
pressure force which is comprised of three distinct components namely:
scattering force, dipole force and spin-curl force [3].

Optical tweezers can operate at tight focussing conditions where the involved
laser beams have a small beam waist. In this regime the electric and magnetic
fields of light are characterized by a longitudinal component which has a
magnitude comparable to the transverse electromagnetic fields components [4,
5]. After the seminal work by Bliokh and Nori, [6], there is a growing
interest on important contribution of this term in various physical effects
including the mechanical effects of optical vortices on two-level atoms [5],
the modifications of the dressed states formalism in the case of two-level
atoms [7], the optical properties of structured light beams [8, 9, 10], the
realization of transverse angular momentum of light [11], and its relation to
the optical spin Hall effect [12].

There are two important issues in the optical tweezers forces related with
this field component. First it is usually ignored in the literature from the
theoretical and mostly from numerical analysis. Second, when the light beam is
in the form of a circularly polarized vortex, the longitudinal field component
is responsible for the appearance of terms involving the coupling of the spin
angular momentum (SAM) with the orbital angular momentum (OAM) of the beam’s
photons in the expressions of the radiation pressure force components [5].

This is exactly what we aim to do in the present work. We are going to show
that this, so far ignored, longitudinal field component gives rise to
considerable modifications in the magnitude of the radiation pressure forces
while the coupling of photon SAM with OAM (which also arises from the
logitudinal field component) brings in chiral effects in the forces
expressions.

In our work we will consider a tiny particle (whose dimensions are within the
limits of Rayleigh regime or dipole approximation) irradiated, in free space,
by a tightly focussed (but still paraxial) circularly polarized and we are
going to find out the modified expressions for the various components of the
radiation pressure force and we are going to compare our results with the
”conventional” expressions for the radiation pressure force components which
are in use by the community.

The structure of the paper is as follows. In Section II we give the analytical
expressions for the three components of the radiation pressure force, at the
beam focus $z=0$, as they have been calculated by taking into consideration
the longitudinal field components. We present numerical results where we
compare their magnitude with that of the radiation pressure force components
when we ignore the longitudinal field component. In Section III we present our
conclusions.

##  II Optical tweezers forces

The optical tweezers force is given by the following relation [3]:

$\mathbf{F}=\frac{1}{4}\alpha^{\prime}\nabla|\mathbf{E}\left(\mathbf{r}\right)|^{2}+\frac{1}{c}\sigma_{\text{ext
}}\mathbf{S}\left(\mathbf{r}\right)-\frac{1}{2}\sigma_{\text{ext}}c\nabla\times\mathbf{s},$
(1)

where is
$\mathbf{S}=\frac{1}{2}\operatorname{Re}\left\\{\mathbf{E}\times\mathbf{B}^{*}\right\\}$
the time-averaged Poynting vector and
$\mathbf{s}=i\frac{\varepsilon_{0}}{2\omega}\left\\{\mathbf{E}\times\mathbf{E}^{*}\right\\}$
is the time-averaged spin density of the incoming wave, $c$ is the speed of
light. The quantity
$\sigma_{ext}=\frac{k_{0}\alpha^{\prime\prime}}{\varepsilon_{0}}$ is the
extinction cross section, where $\alpha^{\prime}$ and $\alpha^{\prime\prime}$
are the real and imginary part of the complex polarizability respectively.

As we may see the radiation pressure force has three distinct and physically
different components [3, 13, 14]:
(i) The gradient force
$\mathbf{F}_{\text{grad}}\left(\mathbf{r}\right)=\frac{1}{4}\alpha^{\prime}\nabla|\mathbf{E}\left(\mathbf{r}\right)|^{2}$
which is a conservative force, responsible for confinement in optical tweezers
and it arises from the potential energy of a dipole in the electric field.
(ii) The scattering force
$\quad\mathbf{F}_{\text{scat}}\left(\mathbf{r}\right)=\frac{\sigma_{\mathrm{ext}}}{c}\mathbf{S}\left(\mathbf{r}\right)$
which is non-conservative and is responsible for the transfer of momentum and
OAM from the field to the particle.
(iii) The spin-curl force, which is also non-conservative,
$\mathbf{F}_{\mathrm{spin-
curl}}\left(\mathbf{r}\right)=-\frac{1}{2}\sigma_{\mathrm{ext}}c\nabla\times\mathbf{s}\left(\mathbf{r}\right)$
which is also non-conservative and arises from polarization gradients in the
electromagnetic field.

The two non-conservative components of the radiation pressure force draw their
origin from the time averaged momentum, which in the free space can be
separated into two parts, one which independent of the polarization and one
which depends on circular polarization. The theoretical justification of this
distinction rests Belifante’s symmetrization of the canonincal stress energy
tensor, which makes the results gauge independent. The spin-curl force is zero
in homeogeneous fields and emerges when spatial inhomogeneities of
polarization and intensity are present [15]. This force has been recently
experimentally measured in an experiment involving circularly polarized
evanescent fields[16].

The electric field of a circularly polarized tightly focussed beam can be
written as follows [4] :

$\mathbf{E}=\frac{1}{2}\biggl{[}\alpha E\vec{{x}}+\beta
E\vec{{y}}+\frac{i}{k}\left(\alpha\frac{\partial E}{\partial
x}+\beta\frac{\partial E}{\partial y}\right)\vec{{z}}\biggr{]}e^{(ikz-i\omega
t)}+C.C.,$ (2)

where $\alpha$ and $\beta$ are in general complex constants with modulus
$|\alpha|=|\beta|=1/\sqrt{2}$. Due to the cylindrical symmetry of optical
vortex beams it is convenient to write the function $E$ in cylindrical polar
coordinates as follows: $\ E=u(\rho,\phi,z)\exp(i\Theta(\rho,\phi,z))$

We consider the case where our vortex beam operates in a Laguerre-Gaussian
mode where the functions $u(\rho,\theta,\phi)$ and $\Theta(\rho,\theta,\phi)$
are given by:

$\begin{split}u&=E_{0}\sqrt{\frac{p!}{(p!+|\ell|!)}}\left(\frac{\rho\sqrt{2}}{w(z)}\right)^{|\ell|}\exp\left[-\frac{\rho^{2}}{w^{2}(z)}\right]L_{p}^{|\ell|}\left[\frac{\rho^{2}}{w^{2}(z)}\right],\\\
\Theta&=\ell\phi-(2p+|\ell|+1)\arctan\left(\frac{z}{z_{R}}\right)+\frac{k\rho^{2}z}{2\left(z^{2}+z_{R}^{2}\right)}.\end{split}$
(3)

In the Eq.(3) $\ell$ is the winding number, $p$ the radial index,
$L_{p}^{|l|}$ is the associated Laguerre polynomial,
$w(z)=w_{0}\sqrt{1+z^{2}/z_{R}^{2}}$ ,$w_{0}$ the beam waist, $z_{R}=\pi
w_{0}^{2}/\lambda$ the Rayleigh range of the beam and $E_{0}$ the amplitude
associated with the mode and is related to the beam power via the equation
$E_{0}=\sqrt{\frac{4P}{\pi\epsilon_{0}c^{3}k^{2}w_{0}^{2}}}$ with
$k=2\pi/\lambda$ being the wavenumber [8]. With this notation the vortex
electric field can be written as

$\displaystyle\begin{split}\vec{E}=\frac{1}{2}\biggl{\\{}\alpha u\vec{x}+\beta
u\vec{y}+\frac{1}{k}\biggl{[}i\left(\alpha\frac{\partial u}{\partial
x}+\beta\frac{\partial u}{\partial y}\right)\\\
-u\left(\alpha\frac{\partial\Theta}{\partial
x}+\beta\frac{\partial\Theta}{\partial
y}\right)\biggr{]}\vec{z}\biggr{\\}}e^{ikz+i\Theta-i\omega t}+C.C.\end{split}$
(4)

Figure 1: The radial component of the dipole force, at the beam focus $(z=0)$,
(a) up when the beam is tightly focused $(w_{0}=0.6\lambda)$ (b) down when the
beam is less focused to a beam waist $(w_{0}=2\lambda)$.In both (a) and (b)
the radial distances are in beam waist $w_{0}$ units. The solid line
represents the conventional result where the longitudinal field components
have not been taken into account, while the dashed line represents the case
where $\ell=+5$ and the dashed-dotted line represents the case $\ell=-5$

###  II.1  Dipole force

Figure 2: The axial component of the dipole force, at the radial distance
$(\rho={w_{0}}\sqrt{\frac{|\ell|}{2}})$ (a) up when the beam is tightly
focused $(w_{0}=0.6\lambda)$ (b) down when the beam is not tightly focused
$(w_{0}=2\lambda)$. In both (a) and (b) axial distances are in Rayleigh range
units $(z_{R})$ . the sold line represents the conventional result where the
longitudinal field components have not been taken into account, dash line
represents the case $\ell=+5$ and dash dot line represents the case $\ell=-5$

For the calculation of the dipole force we need the quantity
$\nabla|\mathbf{E}\left(\mathbf{r}\right)|^{2}$. After calculations the
expression for the dipole force, for a circularly polarized beam, is given by:

$\begin{split}&\mathbf{F}_{\text{grad}}\left(\mathbf{r}\right)=\frac{1}{4}\alpha^{\prime}\nabla|\mathbf{E}\left(\mathbf{r}\right)|^{2}\\\
&=\frac{1}{4}\alpha^{\prime}\biggl{\\{}\frac{1}{2}u\nabla u-\frac{i\sigma
l}{4k^{2}}\nabla\left(\frac{u}{\rho}\frac{\partial
u}{\partial\rho}\right)+\frac{1}{8k^{2}}\biggl{[}2\left(\frac{\partial
u}{\partial\rho}\right)\nabla\frac{\partial u}{\partial\rho}\\\
&+\nabla\left(u\left(\frac{\partial\Theta}{\partial\rho}\right)^{2}\right)+\frac{ul^{2}}{\rho^{2}}\left(\nabla
u-\frac{u}{\rho}\widehat{\rho}\right)\biggr{]}\biggr{\\}},\end{split}$ (5)

where $\sigma=(a*b^{*}-a^{*}b)$ is the spin of the light beam. The first term
in the RHS is the ”conventional” result which has been so far considered by
the community. All the rest terms owe their existence to the longitudinal
field component. The second term in the curly brackets, clearly shows a
coupling between the SAM and OAM carried by the beam photons and vanishes when
the beam is linearly polarized or when the beam photons carry no OAM (like an
ordinary Gaussian beam). This terms exhibits chirality, in the sense that it
reverses sign when either $\sigma$ or $\ell$ do. Similarly the last term
vanishes when there is no OAM in the beam. It is clear that the analytical
work shows a considerable deviation from the conventional results. Let us see
what the numerical work shows.

In all the numerical examples in this work we consider the following values
for the involve parameters: an LG beam of wavelength
$\lambda=633\mathrm{~{}nm}$, a beam waist $w_{0}=0.65\mu\mathrm{m}$, a power
$P=10\mathrm{~{}mW}$, a radial number $p=0$, a winding number $\ell=5$, and a
counterclockwise circular polarization ($\sigma=i$). The beam irradiates a
glass particle of relative dielectric permittivity $\varepsilon=2.25$ [3].

For the components of the dipole force, at beam focus $(z=0)$, we get the
plots shown in Fig.1 and in Fig.2 . In each plot we show three curves: a curve
representing the ”conventional” result where the longitudinal field components
have not been taken into account (solid line), a curve which represents the
case of positive winding number $+\ell$ (dashed line) and a curve for a
negative winding number $-\ell$ (dashed-dotted line) .

In Fig.1(a) we see that the magnitude of the maximum force increases
considerably when we take into account the longitudinal component of the field
with a slight shift of the maximum along the radial direction to the left in
the case of +$\ell$ and to the right in the case -$\ell$. The origin of this
shifts in the maxima positions has to do with the fact that the the radial
position at which the intensity maximizes is different when we take into
account the longitudinal field components. Indeed if we ignore them we know
that the radial distance at which the intensity maximizes is
$w_{0}\sqrt{|\ell|/2}$ while in the present case this distance is modified by
the existence of new terms in tightly focusing conditions. In Fig.1(b) we see
how the new results approximate the conventional results when we choose a less
tighter focusing with a beam waist $w_{0}=2\lambda$. The same behavior is
repeated when we calculate the axial component of the dipole force as it is
shown in Figs.2(a)-(b). As we may see the plot shows the axial component as a
function of the distance along the $z-$axis, keeping the radial distance equal
to $w_{0}\sqrt{|\ell|/2}$. So the axial force is not plotted at the focus
($z=0$) as the radial component. The reason is that at the focus the axial
force component is everywhere equal to zero along the radial direction. Our
calculations have also shown that there is no component for dipole force along
$\phi$ direction.

###  II.2  Scattering force

The scattering force is proportional to the time averaged Poynting vector and
is given by the relation

$\begin{split}&\mathbf{F}_{\text{scat}}\left(\mathbf{r}\right)=\frac{\sigma_{\mathrm{ext}}}{c}\mathbf{S}\left(\mathbf{r}\right)\\\
&=\frac{\varepsilon_{0}\sigma_{\mathrm{ext}}}{2c}\biggl{[}\mathrm{i}\sigma\left(E\nabla
E^{*}-E^{*}\nabla E\right)+2k|E|^{2}\hat{z}-\frac{\partial|E|^{2}}{\partial
r}\hat{\bm{\phi}}\biggr{]}.\end{split}$ (6)

where, as we have said, $E=ue^{i\Theta}$. The first term in the square bracket
of the second line has a zero contribution for linearly polarized light and is
entirely due to the existence of the longitudinal field component. However new
terms do exist, even in the case of a linearly polarized beam, in the axial
and azimuthal scattering force components due to the new terms that exist in
the expression of $|E|^{2}$.

Figure 3: The azimuthal component of the scattering force, at the beam focus
$(z=0)$, (a) up when the beam is tightly focused $(w_{0}=0.6\lambda)$ (b) down
when the beam is not tightly focused $(w_{0}=2\lambda)$. In both (a) and (b)
the radial distances are in beam waist $w_{0}$ units. The sold line represents
the conventional result where the longitudinal field components have not been
taken into account and dash line represents the case $\ell=+5$

We proceed to plot the scattering force components at $z=0$. At this position
the radial force component is zero and the scattering force has only axial and
azimuthal components. We use the same parameter values as in the case of the
dipole force. In Figs.3(a)-(b) we present the azimuthal component of the
scattering force. The solid one represents the force without taking into
account the longitudinal component of the field, the dashed one represents the
same force but taking into account the longitudinal component of the electric
field with positive $\ell$ and the dashed-dotted the magnitude of the same
force when we consider a negative winding number $-\ell$. In the later case
the sign of the force component becomes negative and this is the reason that
we plot its magnitude. In Fig.3(b) we see the same forces in the case of a
relatively larger beam waist equal to $w_{0}=2\lambda$. As the beam waist
increases (weaker focussing) the azimuthal scattering force component becomes
weaker because the relevant beam intensity is smaller for given power and is
similar in magnitude to the result obtained if we ignore the longitudinal
field component.

In Figs.4(a)-(b) we present the axial component of the scattering force. The
solid one represents the force without taking into account the longitudinal
component of the field, the dashed one represents the same force but taking
into account the longitudinal component of the electric field with positive
$\ell$ and the dashed-dotted the magnitude of the same force when we consider
a negative winding number $-\ell$. In Fig.4(b) we see the same forces in the
case of a relatively larger beam waist equal to $w_{0}=2\lambda$. As the beam
waist increases (weaker focusing) the azimuthal scattering force component
becomes weaker because the relevant beam intensity is smaller for given power.
The relevant relation between the conventional result (shown in solid line)
and the new results not only remains but it becomes stronger (we need
explanation here). In both Fig.4(a)-(b) we do not see the dashed-dotted curve
which corresponds to the case of $-\ell$. The reason is that the axial
scattering force does depend on $|\ell|$ and thus the two plots for $\pm\ell$
coincide at $z=0$ as it can be seen from the Eq.LABEL:Eq:6.

Figure 4: The axial component of the scattering force at the beam focus
$(z=0)$ (a) up when the beam is tightly focused $(w_{0}=0.6\lambda)$ (b) down
when the beam is not tightly focused $(w_{0}=2\lambda)$. In both (a) and (b)
the radial distances are in beam waist $w_{0}$ units. the sold line represents
the conventional result where the longitudinal field components have not been
taken into account and dash line represents the case $\ell=+5$

###  II.3  Spin curl force

Figure 5: The azimuthal component for the spin curl force, at the beam focus
$(z=0)$, (a) up when the beam is tightly focused $(w_{0}=0.6\lambda)$ (b) down
when the beam is not tightly focused $(w_{0}=2\lambda)$. In both (a) and (b)
the radial distances are in beam waist $w_{0}$ units the sold line represents
the conventional result where the longitudinal field components have not been
taken into account also it represents the case $\ell=+5$ and the case
$\ell=-5$

.

The spin curl force, as we have said, is given by
$\mathbf{F}_{\mathrm{spin}}\left(\mathbf{r}\right)=-\frac{1}{2}\sigma_{\mathrm{ext}}c\nabla\times\mathbf{s}\left(\mathbf{r}\right)$.
The expression of the spin curl force involves the calculation of the
following quantity:

$\begin{split}\mathbf{F}_{\mathrm{spin}}\left(\mathbf{r}\right)=\frac{1}{16k}&\alpha^{\prime\prime}\bigg{\\{}\left[-i\frac{\partial}{\partial
z}\left(u\frac{\partial
u}{\partial\rho}\right)-2u\sigma\frac{\ell}{\rho}\frac{\partial u}{\partial
z}\right]\widehat{\rho}\\\ &+\sigma\left[\frac{\partial}{\partial
z}\left(u^{2}\frac{\partial\Theta}{\partial\rho}\right)-2uk\frac{\partial
u}{\partial\rho}\right]\widehat{\phi}\\\
&+{\frac{1}{\rho}}\left[i\frac{\partial}{\partial\rho}\left(u\rho\frac{\partial
u}{\partial\rho}\right)+2u\sigma l\frac{\partial
u}{\partial\rho}\right]\widehat{{z}}\biggr{\\}}\end{split}$ (7)

The magnitude of this force is very small compared to magnitude of dipole and
scattering force and is normally ignored from the analysis. From Eq.7 we see
that the radial component contains partial derivatives with respect to $z$
which means that the radial component vanishes at the focal plane $z=0$. The
azimuthal component is proportional to $\sigma$ and thus vanishes for a
linealry polarized beam. The axial component contains a term which couples the
photon SAM with OAM, as well as the radial term but in the later case it
vanishes on the focus plane. In Figs.5-6 we present the azimuthal and axial
components of the spin-curl force. It can be seen from the diagram that the
inclusion of the longitudinal component of the field in the calculations does
not affect the magnitude of the azimuthal component of the spin curl force.
The force expression in beam foucs plane $(z=0)$ contains the absolute value
of $\ell$ and therefore positive or negative has the same curve.

If the beam is not tightly focused we can see a significant decrease in this
force with same behaviour. Along the axial z-direction our calculations show
that at $z=0$ there is no component in the conventional spin-curl force. Such
a component arises only when we take into account the longitudinal field
component since this field component actually creates a polarization gradient
in the light field. Our calculations have also shown that the spin-curl force
has no radial at the beam focus plane.

Figure 6: The axial component of the spin curl force,at the beam focus
$(z=0)$, (a) up when the beam is tightly focused $(w_{0}=0.6\lambda)$ (b) down
when the beam is not tightly focused $(w_{0}=2\lambda)$. In both (a) and (b)
axial distances are in Rayleigh range units $(z_{R})$. The dash line
represents the case $\ell=+5$ and the dash dot line represents the case
$\ell=-5$ where the longitudinal field components have been taken into
account.

##  III conclusions

In our work we considered the interaction of a dielectric particle with a
tightly focused circularly polarized LG beam. The size of the particle is
smaller than the light wavelength and its interaction with the beam occurs in
the free space. Due to strong focusing the electric and magnetic field of the
beam have a longitudinal component with a magnitude comparable with the
transverse ones. These components have been ignored from the analysis of the
tweezing forces on the particle. We calculated the three components of the
radiation pressure force and we showed that there are considerable
modifications in their magnitude when compared with the results obtained
ignoring the longitudinal field component. Moreover, in the case of a
circularly polarized vortex beam, we get new terms which arise from a coupling
of the photon SAM and OAM. The longitudinal term is responsible for the
appearance of a polarization gradient in the light field. As a result it is
responsible for the appearance of a spin-curl force of magnitude comparable to
the gradient and scattering forces which cannot be ignored, as it has so far
been done in the literature where it was considered negligible. The spin curl
force has a significant effect on the dynamics of a particle trapped in
optical tweezers as has been demonstrated [16]. The longitudinal term is
present in any coherent tightly focused laser beam even in a Gaussian one. We
have performed numerical and analytical work for such a beam and we have seen
that the modifications are negligible. This shows clearly that the
modifications in the expressions of the forces are a result of the coupling of
the photon SAM and OAM which is present in the case of a tightly focused
optical vortex. Our analysis concerns only the the calculation of forces in
the Rayleigh regime. It is of great interest to investigate how the inclusion
of the longitudinal field component will modify the expressions of the forces
when the particle size is far larger than the optical wavelength (ray optics
regime), or between these two regimes where we need to use the Lorentz-Mie
theory or advanced numerical techniques.

## References

  * Ashkin [1970] A. Ashkin, Acceleration and trapping of particles by radiation pressure, Physical review letters 24, 156 (1970).
  * Ashkin  _et al._ [1986] A. Ashkin, J. M. Dziedzic, J. E. Bjorkholm, and S. Chu, Observation of a single-beam gradient force optical trap for dielectric particles, Optics letters 11, 288 (1986).
  * Jones  _et al._ [2015] P. H. Jones, O. M. Maragò, and G. Volpe,  _Optical tweezers: Principles and applications_ (Cambridge University Press, 2015).
  * Allen  _et al._ [1999] L. Allen, M. Padgett, and M. Babiker, Iv the orbital angular momentum of light, in  _Progress in optics_ , Vol. 39 (Elsevier, 1999) pp. 291–372.
  * Lembessis  _et al._ [2021] V. E. Lembessis, K. Koksal, J. Yuan, and M. Babiker, Chirality-enabled optical dipole potential energy for two-level atoms, Physical Review A 103, 013106 (2021).
  * Bliokh and Nori [2015] K. Y. Bliokh and F. Nori, Transverse and longitudinal angular momenta of light, Physics Reports 592, 1 (2015).
  * Lembessis  _et al._ [2023] V. Lembessis, A. Lyras, and O. Aldossary, Two-level atom dynamics induced by a spin-orbit coupled optical vortex: dressed states formulation, Optics Continuum 2, 1256 (2023).
  * Koksal  _et al._ [2023] K. Koksal, M. Babiker, and V. Lembessis, Optical characteristics of radially-polarised twisted light, Journal of Optics 25, 065501 (2023).
  * Koksal  _et al._ [2022a] K. Koksal, M. Babiker, V. Lembessis, and J. Yuan, Truncated optical bessel modes, Physical Review A 105, 063512 (2022a).
  * Koksal  _et al._ [2022b] K. Koksal, M. Babiker, V. Lembessis, and J. Yuan, Hopf index and the helicity of elliptically polarized twisted light, JOSA B 39, 459 (2022b).
  * Aiello  _et al._ [2015] A. Aiello, P. Banzer, M. Neugebauer, and G. Leuchs, From transverse angular momentum to photonic wheels, Nature Photonics 9, 789 (2015).
  * Banzer  _et al._ [2013] P. Banzer, M. Neugebauer, A. Aiello, C. Marquardt, N. Lindlein, T. Bauer, and G. Leuchs, Demonstration of a state of light with purely transverse angular momentum, in  _Conference on Coherence and Quantum Optics_ (Optica Publishing Group, 2013) pp. T2A–1.
  * Ruffner and Grier [2013] D. B. Ruffner and D. G. Grier, Comment on “scattering forces from the curl of the spin angular momentum of a light field”, Physical Review Letters 111, 059301 (2013).
  * Chen  _et al._ [2011] J. Chen, J. Ng, Z. Lin, and C. T. Chan, Optical pulling force, Nature photonics 5, 531 (2011).
  * Svak  _et al._ [2018] V. Svak, O. Brzobohatỳ, M. Šiler, P. Jákl, J. Kaňka, P. Zemánek, and S. Simpson, Transverse spin forces and non-equilibrium particle dynamics in a circularly polarized vacuum optical trap, Nature communications 9, 5453 (2018).
  * Antognozzi  _et al._ [2016] M. Antognozzi, C. Bermingham, R. Harniman, S. Simpson, J. Senior, R. Hayward, H. Hoerber, M. Dennis, A. Bekshaev, K. Bliokh, _et al._ , Direct measurements of the extraordinary optical momentum and transverse spin-dependent force using a nano-cantilever, Nature Physics 12, 731 (2016).


# Hybrid Fuzzy-Crisp Clustering Algorithm: Theory and Experiments

Akira R. Kinjo1 and Daphne Teck Ching Lai2
1Department of Mathematics and 2School of Digital Science,
University Brunei Darussalam

###### Abstract

With the membership function being strictly positive, the conventional fuzzy
c-means clustering method sometimes causes imbalanced influence when clusters
of vastly different sizes exist. That is, an outstandingly large cluster drags
to its center all the other clusters, however far they are separated. To solve
this problem, we propose a hybrid fuzzy-crisp clustering algorithm based on a
target function combining linear and quadratic terms of the membership
function. In this algorithm, the membership of a data point to a cluster is
automatically set to exactly zero if the data point is “sufficiently” far from
the cluster center. In this paper, we present a new algorithm for hybrid
fuzzy-crisp clustering along with its geometric interpretation. The algorithm
is tested on twenty simulated data generated and five real-world datasets from
the UCI repository and compared with conventional fuzzy and crisp clustering
methods. The proposed algorithm is demonstrated to outperform the conventional
methods on imbalanced datasets and can be competitive on more balanced
datasets.

##  1 Introduction

Fuzzy clustering is a clustering method in which each data point may belong to
multiple clusters [11]. The membership of a data point to clusters is
expressed in terms of a membership function. In the conventional fuzzy c-means
clustering (FCM) method [3, 1], the membership function is strictly positive.
This sometimes causes a problem when there are clusters of widely different
sizes. For example, consider a set of samples drawn from the following
Gaussian mixture (Fig. 1):

Figure 1: An example of two clusters of very different sizes (c.f. Eqs. 1 and
2).

$\rho(x)\propto 10^{3}\phi(x;-5,1)+\phi(x;+5,1)$ (1)

where the Gaussian density is given as

$\phi(x;\mu,\sigma)=\frac{1}{\sqrt{2\pi\sigma^{2}}}e^{-\frac{(x-\mu)^{2}}{2\sigma^{2}}}.$
(2)

There are two clearly separated clusters (i.e., $c=2$), and one (Cluster 1) of
them is 1,000 times more populated than the other (Cluster 2). Now, let’s
apply the conventional fuzzy c-means clustering to this data set. That is, we
try to minimize the following objective function:

$J(\\{u_{i}(x)\\},\\{v_{i}\\})=\sum_{i=1}^{c}\int_{V}dx\rho(x)[u_{i}(x)]^{m}\|x-v_{i}\|^{2}$
(3)

where $V$ is the range of data point $x$ (in the present case,
$V=(-\infty,+\infty)$), $m>1$ is a constant, $v_{i}$ is the coordinate of the
$i$-th cluster center, and $u_{i}(x)$ is the membership of the point $x$ to
the $i$-th cluster. This objective function can be optimized by iteratively
applying the following mutually recursive formulae:

$\displaystyle u_{i}(x)$ $\displaystyle\leftarrow$
$\displaystyle\left(\sum_{j=1}^{c}\left(\frac{\|x-v_{i}\|^{2}}{\|x-v_{j}\|^{2}}\right)^{\frac{1}{m-1}}\right)^{-1},$
(4) $\displaystyle v_{i}$ $\displaystyle\leftarrow$
$\displaystyle\frac{\int_{V}dx\rho(x)[u_{i}(x)]^{m}x}{\int_{V}dx\rho(x)[u_{i}(x)]^{m}}.$
(5)

It is noted that the fuzzy c-means clustering reduces to the crisp c-means
(k-means) clustering in the limit $m\to 1+0$. With the data set following
$\rho(x)$ in Eq. (1) and the initial cluster centers given as $v_{1}=-5$ and
$v_{2}=+5$, the conventional fuzzy c-means clustering yields the result shown
in Fig. 2A and B. We can see that the center of Cluster 2 is attracted towards
the center of Cluster 1 on each update, and at convergence, we have
$v_{1}\approx-5.80$ and $v_{2}\approx-4.17$. It is clear that the points with
large membership values to Cluster 2 are buried within Cluster 1, and the
membership at around $x=+5$ (the center of the smaller Gaussian) is very fuzzy
with $u_{1}(5)\approx 0.42$ and $u_{2}(5)\approx 0.58$.

Figure 2: Applications of fuzzy clustering to the data generated by the
Gaussian mixture (Eq. 1, Fig. 1). A and B are based on the conventional FCM
method. C and D are based on the Hybrid Fuzzy-Crisp (HFC) method. A:
Trajectory of cluster centers according to FCM. B: Membership functions
obtained by FCM with initial cluster centers being the means of the mixed
Gaussians. C: Trajectory of cluster centers according to HFC. D: The
membership functions obtained by HFC with initial cluster centers being the
means of the mixed Gaussians.

As we have already mentioned, the source of the above problem is that the
membership functions assume strictly positive values. Thus, when one cluster
is outstandingly large, it will drag the centers of the other clusters however
far they are separated.

One way to alleviate this problem is to set the membership of a data point to
a cluster to exactly zero if the data point is “sufficiently” far from the
cluster center. Kinjo [7] has shown that this is indeed possible by using an
objective function that combines linear and quadratic terms of the membership
functions. However, his proposed algorithm was rather primitive and extremely
inefficient (exponential in the number of clusters). Later, Klawonn and
Höppner [8] independently proposed essentially the same objective function,
but with a more efficient, practical algorithm. In the simplest version of
this hybrid fuzzy-crisp clustering method (HFC), the objective function to be
minimized is given as

$J(\\{u_{i}(x)\\},\\{v_{i}\\})=\sum_{i=1}^{c}\int_{V}dx\rho(x)\left(\frac{u_{i}^{2}(x)}{2}+u_{i}(x)\right)\|x-v_{i}\|^{2}$
(6)

where $u_{i}^{2}(x)=[u_{i}(x)]^{2}$. We will provide a more general objective
function below. If we apply HFC to the distribution given by Eq. (1), we
obtain the result shown in Fig. 2C, D. The center of the smaller cluster is no
longer attracted to the larger cluster. In fact, the centers converge at
$v_{1}\approx-5.00$ and $v_{2}\approx+4.98$ (Fig. 2C). Furthermore, the
membership functions can take values exactly equal to 1 or 0 near cluster
centers (Fig. 2D). As we can observe in this example, HFC combines
characteristics of fuzzy and crisp (hard) clustering methods thereby solving
some difficulties in the conventional clustering methods.

In this paper, we present a new algorithm for the fuzzy-crisp clustering based
on a geometric interpretation of the objective function minimization problem.

##  2 Theory

###  2.1 Hybrid Fuzzy-Crisp c-means clustering (HFC)

Suppose we have the data distribution $\rho(x)$ with $x\in\mathbb{R}^{d},x\in
V$ where $V$ is the range of the data points. When the data consist of $n$
discrete points, $x_{1},x_{2},\cdots,x_{n}$, the distribution can be written
as $\rho(x)=\frac{1}{n}\sum_{k=1}^{n}\delta(x-x_{k})$ where $\delta$ is the
Dirac’s delta function. We assume that there are $c$ clusters whose centers
are $v_{1},v_{2},\cdots,v_{c}\in\mathbb{R}^{d}$. The membership of the data
point $x$ to cluster $i$ is represented by the membership functions
$u_{i}(x)$. We define the distance function between data point $x$ and cluster
center $v_{i}$ by the Euclidean norm $d_{i}(x)=\|x-v_{i}\|$. In the following,
we denote the square of the function value $f(x)$ by $f^{2}(x)$ (e.g.,
$u_{i}^{2}(x)$ indicates $[u_{i}(x)]^{2}$). Our task is to find the optimal
cluster centers $\\{v_{i}\\}$ and membership functions $\\{u_{i}(x)\\}$ that
minimize the following objective function

$J(\\{u_{i}(x)\\},\\{v_{i}\\})=\sum_{i=1}^{c}\int_{V}dx\rho(x)\left(\frac{u_{i}^{2}(x)}{2\alpha_{i}}+\frac{u_{i}(x)}{\beta}\right)d_{i}^{2}(x)$
(7)

with the constraints

$u_{i}(x)\geq 0,i=1,\cdots,c$ (8)

and

$\sum_{i=1}^{c}u_{i}(x)=1$ (9)

for all $x\in V$. $\alpha_{i}$’s and $\beta$ are positive parameters, and we
impose the constraint

$\displaystyle\sum_{i=1}^{c}\alpha_{i}$ $\displaystyle=$ $\displaystyle 1.$
(10)

We assume that $\beta$ is given and fixed. $\alpha_{i}$’s may be either fixed
or optimized. We have found that the weight of the linear term ($\beta$)
cannot be cluster-dependent (such as $\\{\beta_{i}\\}$) as it may cause the
membership functions to explode.

Given the membership functions $\\{u_{i}(x)\\}$, the optimal cluster centers
$\\{v_{i}\\}$ are obtained by solving

$\frac{\partial J}{\partial v_{i}}=0,$ (11)

which gives

$v_{i}=\frac{\int_{V}dx\rho(x)({u_{i}^{2}(x)}/{2\alpha_{i}}+u_{i}(x)/\beta)x}{\int_{V}dx\rho(x)({u_{i}^{2}(x)}/{2\alpha_{i}}+u_{i}(x)/\beta)}.$
(12)

Optimizing the objective function with respect to $\alpha_{i}$’s with the
constraint (Eq. 10), with all the other variables being fixed, gives

$\alpha_{i}=\frac{\sqrt{\int_{V}dx\rho(x)u_{i}^{2}(x)d_{i}^{2}(x)}}{\sum_{j=1}^{c}\sqrt{\int_{V}dx\rho(x)u_{j}^{2}(x)d_{j}^{2}(x)}}.$
(13)

Given the cluster centers $\\{v_{i}\\}$, the objective function $J$ is
minimized with constraints Eqs. (8) and (9). This minimization problem can be
solved with a geometric method (see the next section). First, we need the
following result:

######  Proposition 1.

For each $x\in V$, there exists a unique subset of clusters $H\subset
C=\\{1,2,\cdots,c\\}$ such that

$\displaystyle{d_{i}^{2}(x)}$ $\displaystyle<$ $\displaystyle
d_{H}^{2}(x)~{}~{}\text{if $i\in H$},$ (14) $\displaystyle{d_{i}^{2}(x)}$
$\displaystyle\geq$ $\displaystyle d_{H}^{2}(x)~{}~{}\text{if $i\notin H$}$
(15)

where

$d_{H}^{2}(x)=\frac{\beta+\sum_{j\in H}\alpha_{j}}{\sum_{j\in
H}\alpha_{j}/d_{j}^{2}(x)}.$ (16)

Using this subset $H$ of clusters, the optimal membership functions are given
for each $x\in V$ as

$u_{i}(x)=\left\\{\begin{array}[]{cc}\frac{\alpha_{i}}{\beta}\left(\frac{d_{H}^{2}(x)}{d_{i}^{2}(x)}-1\right)&\text{if
$i\in H$},\\\ 0&\text{if $i\notin H$}.\end{array}\right.$ (17)

The proof is given in the subsection 2.5 (Geometric interpretation) below.

###  2.2 Algorithm for finding the membership function of a given data point

Here is an algorithm to compute the optimal membership function
$\\{u_{i}(x)\\}_{i=1,\cdots,c}$ for a given data point $x\in V$.

  1. 1.

Let $\nu=0$ and $H^{(0)}=\\{1,2,\cdots,c\\}$.

  2. 2.

Calculate $d^{2}_{H^{(\nu)}}(x)$ by Eq. (16).

  3. 3.

If $d_{i}^{2}(x)<d^{2}_{H^{(\nu)}}(x)$ for all $i\in H^{(\nu)}$, then set

$u_{i}(x)=\left\\{\begin{array}[]{cc}\frac{\alpha_{i}}{\beta}\left(\frac{d_{H^{(\nu)}}^{2}(x)}{d_{i}^{2}(x)}-1\right)&\text{if
$i\in H^{(\nu)}$},\\\ 0&\text{if $i\notin H^{(\nu)}$}.\end{array}\right.$ (18)

and we are done. Otherwise, go to the next step.

  4. 4.

Set

$H^{(\nu+1)}=\\{i\in H^{(\nu)}\mid d_{i}^{2}(x)<d^{2}_{H^{(\nu)}}(x)\\}.$ (19)

  5. 5.

Update $\nu:=\nu+1$ and go to Step 2.

In Step 4, multiple clusters may be eliminated at once. This is in contrast to
the algorithm by Klawonn and Höppner [8] where only one cluster is eliminated
at a time.

###  2.3 Relationship to fuzzy c-means clustering

We can show that HFC can recover FCM (with $m=2$) in the limit
$\beta\to\infty$. In this limit, we have

$\lim_{\beta\to\infty}d_{H}^{2}=\infty$ (20)

so that Eq. (14) is trivially satisfied for any subset $H\subset C$, in
particular, for $H=C$. Furthermore, the membership function becomes

$u_{i}(x)=\lim_{\beta\to\infty}\frac{\alpha_{i}}{\beta}\left(\frac{d_{H}^{2}(x)}{d_{i}^{2}(x)}-1\right)=\frac{\alpha_{i}/d_{i}^{2}}{\sum_{j\in
H}\alpha_{j}/d_{j}^{2}(x)}$ (21)

which is the membership function of the FCM with $m=2$ (with the weighting
factors $\\{\alpha_{i}\\}$).

###  2.4 Analytical results with $c=2$

Figure 3: The case for $c=2$ with $v_{1}=(-1,0),v_{2}=(1,0)$ and
$\alpha_{1}=0.2,\alpha_{2}=0.8$. Two cluster centers are marked with black
dots labeled “Cluster 1” and “Cluster 2”. The green circle indicates the area
where $d_{2}\geq d_{C}$ so that $u_{1}=1,u_{2}=0$; the magenta circle
indicates the area where $d_{1}\geq d_{C}$ so that $u_{1}=0,u_{2}=1$; the area
outside the two circles indicates that with $u_{1}>0,u_{2}>0$.

We have already given an example with $c=2$ (two clusters) in Fig. 2 (C and
D). Here, we illustrate in detail the case with $c=2$ with fixed cluster
centers.

Suppose that there are two cluster centers in $\mathbb{R}^{2}$, one at
$v_{1}=(-x_{0},0)$ and the other at $v_{2}=(x_{0},0)$ with $x_{0}>0$. Given an
arbitrary data point $(x,y)$, what are the membership functions?

To simplify the notation, let us define the distances

$\displaystyle d_{1}^{2}$ $\displaystyle=$ $\displaystyle(x+x_{0})^{2}+y^{2},$
(22) $\displaystyle d_{2}^{2}$ $\displaystyle=$
$\displaystyle(x-x_{0})^{2}+y^{2}.$ (23)

First, consider the case where $H=C=\\{1,2\\}$. The threshold distance is

$d_{C}^{2}=\frac{\beta+\alpha_{1}+\alpha_{2}}{\alpha_{1}/d_{1}^{2}+\alpha_{2}/d_{2}^{2}}.$
(24)

#### Case 1: $d_{1}^{2}<d_{C}^{2}$ and $d_{2}^{2}<d_{C}^{2}$.

Then, the two inequalities

$\left\\{\begin{array}[]{c}d_{1}^{2}<d_{C}^{2}\\\
d_{2}^{2}<d_{C}^{2}\end{array}\right.$ (25)

are equivalent to

$\displaystyle\left[x-(1+2\alpha_{2}/\beta)x_{0}\right]^{2}+y^{2}$
$\displaystyle>$ $\displaystyle
4(\alpha_{2}/\beta)(1+\alpha_{2}/\beta)x_{0}^{2}$ (26)
$\displaystyle\left[x+(1+2\alpha_{1}/\beta)x_{0}\right]^{2}+y^{2}$
$\displaystyle>$ $\displaystyle
4(\alpha_{1}/\beta)(1+\alpha_{1}/\beta)x_{0}^{2}.$ (27)

Eq. (26) corresponds to the exterior of the circle centered at
$((1+2\alpha_{2}/\beta)x_{0},0)$ with radius
$2x_{0}\sqrt{(\alpha_{2}/\beta)(1+\alpha_{2}/\beta)}$, and Eq. (27) to
exterior of the circle centered at $(-(1+2\alpha_{1}/\beta)x_{0},0)$ with
radius $2x_{0}\sqrt{(\alpha_{1}/\beta)(1+\alpha_{1}/\beta)}$. Thus, when the
data point $(x,y)$ is in the intersection of these regions, the membership
functions $u_{1}$ and $u_{2}$ are both non-zero and are given by

$\displaystyle u_{1}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{1}}{\beta}\left(d_{C}^{2}/d_{1}^{2}-1\right)=\alpha_{1}\left(\frac{\beta+\alpha_{1}+\alpha_{2}}{\alpha_{1}+\alpha_{2}d_{1}^{2}/d_{2}^{2}}-1\right),$
(28) $\displaystyle u_{2}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{2}}{\beta}\left(d_{C}^{2}/d_{2}^{2}-1\right)=\alpha_{2}\left(\frac{\beta+\alpha_{1}+\alpha_{2}}{\alpha_{1}d_{2}^{2}/d_{1}^{2}+\alpha_{2}}-1\right).$
(29)

It is noted that we indeed have $u_{1}>0$, $u_{2}>0$ and, using
$\alpha_{1}+\alpha_{2}=1$,

$u_{1}+u_{2}=1.$ (30)

#### Case 2: $d_{1}^{2}<d_{C}^{2}$ and $d_{2}^{2}\geq d_{C}^{2}$.

In this case, we set $H_{1}=\\{i\in C|d_{i}^{2}<d_{C}^{2}\\}=\\{1\\}$ and

$d_{H_{1}}^{2}=\frac{\beta+\alpha_{1}}{\alpha_{1}/d_{1}^{2}}=(1+\beta/\alpha_{1})d_{1}^{2}.$
(31)

Let us see if we indeed have

$\displaystyle d_{1}^{2}$ $\displaystyle<$ $\displaystyle d_{H_{1}}^{2},$ (32)
$\displaystyle d_{2}^{2}$ $\displaystyle\geq$ $\displaystyle d_{H_{1}}^{2}.$
(33)

These inequalities are equivalent, respectively, to

$\displaystyle d_{1}^{2}$ $\displaystyle>$ $\displaystyle 0,$ (34)
$\displaystyle\left[x+(1+2\alpha_{1}/\beta)x_{0}\right]^{2}+y^{2}$
$\displaystyle\leq$ $\displaystyle
4(\alpha_{1}/\beta)(1+\alpha_{1}/\beta)x_{0}^{2}.$ (35)

Eq. (34) is trivially satisfied. Eq. (35) is the complement of the region
defined by Eq. (27) which, in turn, follows from $d_{2}^{2}<d_{C}^{2}$. Eq.
(35) follows from $d_{2}^{2}\geq d_{C}^{2}$, and hence is satisfied. The
membership function is given by

$\displaystyle u_{1}$ $\displaystyle=$
$\displaystyle\alpha_{1}(d_{H_{1}}^{2}/d_{1}^{2}-1)=1,$ (36) $\displaystyle
u_{2}$ $\displaystyle=$ $\displaystyle 0.$ (37)

#### Case 3: $d_{1}^{2}\geq d_{C}^{2}$ and $d_{2}^{2}<d_{C}^{2}$.

By symmetry, we set $H_{2}=\\{2\\}$ and define

$d_{H_{2}}^{2}=\frac{\beta+\alpha_{2}}{\alpha_{2}/d_{2}^{2}}=(1+\beta/\alpha_{2})d_{2}^{2}.$
(38)

Thus, we have

$\displaystyle d_{1}^{2}$ $\displaystyle\geq$ $\displaystyle d_{H_{2}}^{2},$
(39) $\displaystyle d_{2}^{2}$ $\displaystyle<$ $\displaystyle d_{H_{2}}^{2}$
(40)

and

$\displaystyle u_{1}$ $\displaystyle=$ $\displaystyle 0,$ (41) $\displaystyle
u_{2}$ $\displaystyle=$ $\displaystyle 1.$ (42)

#### Case 4: $d_{1}^{2}\geq d_{C}^{2}$ and $d_{2}^{2}\geq d_{C}^{2}$.

This case cannot happen. The distance between the centers of the circles
defined by $d_{1}^{2}=d_{C}^{2}$ and $d_{2}^{2}=d_{C}^{2}$ is

$(1+2\alpha_{2}/\beta)x_{0}+(1+2\alpha_{1}/\beta)x_{0}=2x_{0}(1+\alpha_{1}/\beta+\alpha_{2}/\beta)=2x_{0}(1+1/\beta).$
(43)

For any positive real numbers $a$ and $b$, we have the inequality

$\sqrt{ab}\leq\frac{a+b}{2}.$ (44)

Using this, the sum of the radii of the two circles, divided by $2x_{0}$, is

$\sqrt{(\alpha_{2}/\beta)(1+\alpha_{2}/\beta)}+\sqrt{(\alpha_{1}/\beta)(1+\alpha_{1}/\beta)}\leq
1+1/\beta.$ (45)

Thus the sum of the radii is smaller than the distance between the centers,
and hence the two circles cannot intersect.

An example is shown in Fig. 3 where we set $\alpha_{1}=0.2$ and
$\alpha_{2}=0.8$. The colored areas indicate the regions where one of the
membership functions is exactly 0.

###  2.5 Geometric interpretation

Here, we show that the problem of finding the optimal membership function
$u_{i}(x)$ for a single data point $x$ can be interpreted as finding the point
in a $(c-1)$-simplex in the $c$-dimensional space that is closest to the
origin. By solving the latter problem, we prove Proposition 1. We will omit
the argument “$x$” henceforth. Thus the objective function is

$J(\\{u_{i}\\})=\sum_{i=1}^{c}\left(\frac{u_{i}^{2}}{2\alpha_{i}}+\frac{u_{i}}{\beta}\right)d_{i}^{2}$
(46)

with constraints

$u_{i}\geq 0$ (47)

and

$\sum_{i=1}^{c}u_{i}=1.$ (48)

Now define

$\displaystyle a_{i}$ $\displaystyle=$
$\displaystyle\frac{d_{i}}{\sqrt{2\alpha_{i}}},$ (49) $\displaystyle b_{i}$
$\displaystyle=$
$\displaystyle\sqrt{\frac{\alpha_{i}}{2}}\frac{d_{i}}{\beta},$ (50)
$\displaystyle\xi_{i}$ $\displaystyle=$ $\displaystyle a_{i}u_{i}+b_{i}.$ (51)

Then the objective function is represented as

$J(\\{\xi_{i}\\})=\sum_{i=1}^{c}(\xi_{i}^{2}-b_{i}^{2})$ (52)

with the constraints

$\xi_{i}\geq b_{i},i=1,2,\cdots,c,$ (53)

and

$\sum_{i=1}^{c}\frac{\xi_{i}-b_{i}}{a_{i}}=1.$ (54)

Eq. (54) defines a hyperplane in $\mathbb{R}^{c}$ that contains points
$P_{1},P_{2},\cdots,P_{c}$ where the $j$-th coordinate of point $P_{i}$ is

$(P_{i})_{j}=a_{j}\delta_{i,j}+b_{j}.$ (55)

Here, $\delta_{i,j}$ is Kronecker’s delta: $\delta_{i,j}=1$ if $i=j$, or
$\delta_{i,j}=0$, otherwise. Together with the restraints Eq. (53), it defines
the $(c-1)$-simplex consisting of vertices $P_{1},P_{2},\cdots,P_{c}$. Since
$b_{i}$’s are constants, minimizing the objective function $J(\\{\xi_{i}\\})$
is equivalent to finding the point in the $(c-1)$-simplex that is closest to
the origin.

Let $C=\\{1,2,\cdots,c\\}$ be the set of indices of vertices. For a subset
$S\subset C$, let $\pi(S)$ denote the hyperplane defined by the points
$\\{P_{i}\\}_{i\in S}$. Similarly, let $\mathcal{S}(S)$ denote the
$(|S|-1)$-simplex with vertices $(P_{i})_{i\in S}$. $\mathcal{S}(S)$ is a
facet of the simplex $\mathcal{S}(C)$. It is clear that the closest point from
the origin to the simplex $\mathcal{S}(C)$ resides in one of the facets of the
simplex (including $\mathcal{S}(C)$).

######  Definition 1 (Optimal subset).

We say that a subset $H\subset C$ is _optimal_ if the closest point from the
origin to the simplex $\mathcal{S}(C)$ exists in the interior of the facet
$\mathcal{S}(H)$.

Note, in particular, that $H$ is _not_ optimal if the closest point from the
origin exists on the boundary of $\mathcal{S}(H)$.

Now let $G$ be the foot of the perpendicular to $\pi(C)$ from the origin $O$.
Then, the coordinates of $G=(\xi_{i}^{G})$ are given as

$\xi_{i}^{G}=\frac{1+\sum_{j\in C}\frac{b_{j}}{a_{j}}}{a_{i}\sum_{j\in
C}\frac{1}{a_{j}^{2}}}.$ (56)

If $\xi_{i}^{G}>b_{i}$ for all $i\in C$, then this $G$ is in the _interior_ of
$\mathcal{S}(C)$ and it is the closest point from the origin. If
$\xi_{i}^{G}\leq b_{i}$ for some $i\in C$ and $\xi_{i}^{G}>0$ for all the
other $i\in C$, then the point $G$ is outside or on one of the facets of the
simplex $\mathcal{S}(C)$, and the closest point $G^{*}$ must exist on one (and
only one) of the facets of $\mathcal{S}(C)$. Let the the set of vertices of
that facet be $H\subset C$. A facet of a simplex is again a simplex (in a
lower dimension). Thus, the point $G^{*}$ is an interior point of simplex
$\mathcal{S}(H)$. We define the following quantities for convenience:

$\displaystyle B$ $\displaystyle=$ $\displaystyle(b_{1},b_{2},\cdots,b_{c}),$
(57) $\displaystyle\mathbf{p}_{i}$ $\displaystyle=$
$\displaystyle\overrightarrow{BP}_{i}=(a_{i}\delta_{i,j})_{j=1,\cdots,c},$
(58) $\displaystyle A_{S}^{2}$ $\displaystyle=$ $\displaystyle\left(\sum_{i\in
S}\frac{1}{a_{i}^{2}}\right)^{-1}\text{ for subset $S\subset C$},$ (59)
$\displaystyle D_{S}$ $\displaystyle=$ $\displaystyle\sum_{i\in
S}\frac{b_{i}}{a_{i}}\text{ for subset $S\subset C$},$ (60)
$\displaystyle\mathbf{n}$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{c}\frac{A_{C}^{2}}{a_{i}^{2}}\mathbf{p}_{i},$ (61)
$\displaystyle\mathbf{q}_{i}$ $\displaystyle=$
$\displaystyle\mathbf{n}-\mathbf{p}_{i}.$ (62)

For any $i,j,k\in C$, we have

$\overrightarrow{P_{j}P_{k}}\cdot\mathbf{q}_{i}=a_{i}^{2}(\delta_{ij}-\delta_{ik}).$
(63)

Thus, for any subset $H\subset C$, the subspaces spanned by
$\\{\mathbf{p}_{i}\\}_{i\in H}$ and by $\\{\mathbf{q}_{i}\\}_{i\in H^{c}}$ are
orthogonal to each other ($H^{c}=C\setminus H$ is the complement of $H$).

######  Lemma 1.

Let $H$ be the optimal subset of $C$, and $G$ be the foot of the perpendicular
to the hyperplane $\pi(C)$ from the origin (Eq. 56). Then, we have

$\overrightarrow{BG}=\sum_{i\in H}s_{i}\mathbf{p}_{i}+\sum_{i\in
H^{c}}s_{i}^{\prime}\mathbf{q}_{i}$ (64)

where $\sum_{i\in H}s_{i}=1$, $s_{i}>0$ for all $i\in H$, $s_{i}^{\prime}\geq
0$ for all $i\in H^{c}$.

###### Proof.

There are three cases.

#### Case 1: $G$ is in the interior of $\mathcal{S}(C)$

In this case, $H=C$ and $H^{c}=\emptyset$. Eq. (64) clearly holds with
$\sum_{i\in C}s_{i}=1$, and $s_{i}>0$ for all $i\in C$ (the second sum over
$i\in H^{c}(=\emptyset)$ is absent).

#### Case 2: $G$ is in the interior of the facet $\mathcal{S}(H)$ and
$H\subsetneq C$

In this case, Eq. (64) holds with $\sum_{i\in H}s_{i}=1$, $s_{i}>0$ for all
$i\in H$, and $s_{i}^{\prime}=0$ for all $i\in H^{c}$.

#### Case 3: $G$ is in the exterior of $\mathcal{S}(C)$

Let $G^{*}$ be the closest point in $\mathcal{S}(C)$ from the origin. Then,
$G^{*}$ is in the interior of some facet $\mathcal{S}(H)$ where $H\subsetneq
C$. We have

$\overrightarrow{BG}=\overrightarrow{BG^{*}}+\overrightarrow{G^{*}G}.$

Since $G^{*}$ is in the interior of $\mathcal{S}(H)$, we have

$\overrightarrow{BG^{*}}=\sum_{i\in H}s_{i}\mathbf{p}_{i}$ (65)

where $\sum_{i\in H}s_{i}=1$ and $s_{i}>0$ for all $i\in H$. Furthermore, the
vector $\overrightarrow{G^{*}G}$ is perpendicular to the facet
$\mathcal{S}(H)$ (i.e., $G^{*}$ is the point on the hyperplane $\pi(H)$ that
is closest to $G$). Noting that the subspace spanned by
$\\{\mathbf{q}_{i}\\}_{i\in H^{c}}$ is orthogonal to the subspace $\pi(H)$
(Eq. 63), the vector $\overrightarrow{G^{*}G}$ can be represented as

$\overrightarrow{G^{*}G}=\sum_{i\in H^{c}}s_{i}^{\prime}\mathbf{q}_{i}.$ (66)

We show $s_{i}^{\prime}>0$ for all $i\in H^{c}$, then we are done. Using Eq.
(63), we have, for each $i\in H^{c}$,

$\overrightarrow{G^{*}P_{i}}\cdot\mathbf{q}_{i}=-a_{i}^{2}<0.$ (67)

This indicates that the vector $\mathbf{q}_{i}$ is pointing “outward” from the
simplex $\mathcal{S}(C)$. Since $P_{i}$ for any $i\in H^{c}$ and $G$ are on
opposite sides of $\mathcal{S}(H)$, it follows that $s_{i}^{\prime}>0$. ∎ By
substituting $\mathbf{n}$ and $\mathbf{q}_{i}$’s, we have

$\displaystyle\overrightarrow{BG}$ $\displaystyle=$ $\displaystyle\sum_{i\in
H}\left(s_{i}+S^{\prime}\frac{A_{C}^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}+\sum_{i\in
H^{c}}\left(S^{\prime}\frac{A_{C}^{2}}{a_{i}^{2}}-s_{i}^{\prime}\right)\mathbf{p}_{i}$
(68)

where $S^{\prime}=\sum_{i\in H^{c}}s_{i}^{\prime}$. From the coordinates of
$G$ (Eq. 56), we have

$\overrightarrow{BG}=\sum_{i=1}^{c}\left(\frac{(1+D_{C})A_{C}^{2}}{a_{i}^{2}}-\frac{b_{i}}{a_{i}}\right)\mathbf{p}_{i}.$
(69)

Comparing this equation with Eq. (68) and noting that the vectors
$(\mathbf{p}_{i})_{i\in C}$ comprise an orthogonal basis set, we have

$\displaystyle s_{i}+S^{\prime}\frac{A_{C}^{2}}{a_{i}^{2}}$ $\displaystyle=$
$\displaystyle\frac{(1+D_{C})A_{C}^{2}}{a_{i}^{2}}-\frac{b_{i}}{a_{i}},~{}i\in
H,$ (70) $\displaystyle S^{\prime}\frac{A_{C}^{2}}{a_{i}^{2}}-s_{i}^{\prime}$
$\displaystyle=$
$\displaystyle\frac{(1+D_{C})A_{C}^{2}}{a_{i}^{2}}-\frac{b_{i}}{a_{i}},~{}i\in
H^{c}.$ (71)

By summing Eq. (71) over $i\in H^{c}$ and recalling $S^{\prime}=\sum_{i\in
H^{c}}s_{i}^{\prime}$, we have

$S^{\prime}=1+D_{C}-(1+D_{H})\frac{A_{H}^{2}}{A_{C}^{2}}.$ (72)

Substituting this into Eqs. (70) and (71), we have the coefficients
$\\{s_{i}\\}$ and $\\{s_{i}^{\prime}\\}$:

$\displaystyle s_{i}$ $\displaystyle=$
$\displaystyle\frac{(1+D_{H})A_{H}^{2}}{a_{i}^{2}}-\frac{b_{i}}{a_{i}},i\in
H;$ (73) $\displaystyle s_{i}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{(1+D_{H})A_{H}^{2}}{a_{i}^{2}}+\frac{b_{i}}{a_{i}},i\in
H^{c}.$ (74)

Thus, $s_{i}>0$ and $s_{i}^{\prime}\geq 0$ translate into

$\displaystyle a_{i}b_{i}$ $\displaystyle<$
$\displaystyle(1+D_{H})A_{H}^{2},i\in H,$ (75) $\displaystyle a_{i}b_{i}$
$\displaystyle\geq$ $\displaystyle(1+D_{H})A_{H}^{2},i\in H^{c}.$ (76)

Eqs. (75) and (76) correspond to Eqs. (14) and (15) in Proposition 1,
respectively. Note that, given that $H$ is the optimal subset, Eq. (65) holds
for all the above cases. The coordinates of
$G^{*}=(\xi_{1}^{*},\xi_{2}^{*},\cdots,\xi_{c}^{*})$ are given by

$\xi_{i}^{*}=\left\\{\begin{array}[]{cl}a_{i}s_{i}+b_{i}&(i\in H),\\\
b_{i}&(i\in H^{c}).\end{array}\right.$ (77)

Thus, the optimal membership function is given by

$u_{i}^{*}=\left\\{\begin{array}[]{cl}s_{i}&(i\in H),\\\ 0&(i\in
H^{c}).\end{array}\right.$ (78)

(c.f. Eq. 51). Substituting Eq. (73), we obtain Eq. (17) in Proposition 1. It
is noted from Eq. (68) that

$\xi_{i}^{G}-b_{i}=(s_{i}+S^{\prime}A_{C}^{2}/a_{i}^{2})a_{i}>0$ (79)

for $i\in H$ because $s_{i},S^{\prime},A_{C}^{2}$ and $a_{i}$ are all
positive. Thus, we have

######  Lemma 2.

If a vertex $i\in C$ is in the optimal facet $H$, i.e., $i\in H$, then
$\xi_{i}^{G}>b_{i}$.

Its contrapositive is a useful corollary:

######  Corollary 1.

If $\xi_{i}^{G}\leq b_{i}$, then $i\notin H$.

Thus, we can eliminate the vertices that do not belong to the optimal subset
based on the sign of $\xi_{i}^{G}-b_{i}$ (or $u_{i}$).

In the above lemmas, the optimal facet $H$ is assumed to be known. To find the
optimal facet, we note the representation of the vector $\overrightarrow{BG}$
in terms of $\\{\mathbf{p}_{i}\\}_{i\in H}$ and $\\{\mathbf{q}_{i}\\}_{i\in
H^{c}}$ (Eq. 64) . The following lemma shows that the subset of $C$ admitting
such a representation is unique.

######  Lemma 3.

The subset $H\subset C$ with the following properties is unique:

$\overrightarrow{BG}=\sum_{i\in H}s_{i}\mathbf{p}_{i}+\sum_{i\in
H^{c}}s_{i}^{\prime}\mathbf{q}_{i}$ (80)

with

$\displaystyle s_{i}$ $\displaystyle>$ $\displaystyle 0,i\in H;$ (81)
$\displaystyle\sum_{i\in H}s_{i}$ $\displaystyle=$ $\displaystyle 1;$ (82)
$\displaystyle s^{\prime}_{i}$ $\displaystyle\geq$ $\displaystyle 0,i\in
H^{c}$ (83)

###### Proof.

Suppose there are two such subsets $H_{1}$ and $H_{2}$. Then, we can represent
the vector $\overrightarrow{BG}$ in two different ways:

$\displaystyle\overrightarrow{BG}$ $\displaystyle=$ $\displaystyle\sum_{i\in
H_{1}}\left(s_{i}+S^{\prime}\frac{A^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}+\sum_{i\in
H_{1}^{c}}\left(S^{\prime}\frac{A^{2}}{a_{i}^{2}}-s_{i}^{\prime}\right)\mathbf{p}_{i}$
(84) $\displaystyle=$ $\displaystyle\sum_{i\in
H_{2}}\left(t_{i}+T^{\prime}\frac{A^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}+\sum_{i\in
H_{2}^{c}}\left(T^{\prime}\frac{A^{2}}{a_{i}^{2}}-t_{i}^{\prime}\right)\mathbf{p}_{i}$
(85)

where $s_{i}>0,s_{i}^{\prime}\geq 0,S^{\prime}=\sum_{i\in
H_{1}^{c}}s_{i}^{\prime}\geq 0$ and $t_{i}>0,t_{i}^{\prime}\geq
0,T^{\prime}=\sum_{i\in H_{2}^{c}}t_{i}^{\prime}\geq 0$. Rearrange the second
equality into:

$\displaystyle\sum_{i\in H_{1}\cap
H_{2}}\left(s_{i}-t_{i}+[S^{\prime}-T^{\prime}]\frac{A^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}+\sum_{i\in
H_{1}\cap
H_{2}^{c}}\left(s_{i}+t_{i}^{\prime}+[S^{\prime}-T^{\prime}]\frac{A^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}$
$\displaystyle+\sum_{i\in H_{1}^{c}\cap
H_{2}}\left(-s_{i}^{\prime}-t_{i}+[S^{\prime}-T^{\prime}]\frac{A^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}+\sum_{i\in
H_{1}^{c}\cap
H_{2}^{c}}\left(-s_{i}^{\prime}+t_{i}^{\prime}+[S^{\prime}-T^{\prime}]\frac{A^{2}}{a_{i}^{2}}\right)\mathbf{p}_{i}$
$\displaystyle=$ $\displaystyle 0$

Because $(\mathbf{p}_{i})$ are an orthogonal basis set, all the coefficients
of $\mathbf{p}_{i}$ must be zero. In particular,

$\displaystyle
s_{i}+t_{i}^{\prime}+[S^{\prime}-T^{\prime}]\frac{A^{2}}{a_{i}^{2}}$
$\displaystyle=$ $\displaystyle 0,~{}~{}i\in H_{1}\cap H_{2}^{c},$ (87)
$\displaystyle-
s_{i}^{\prime}-t_{i}+[S^{\prime}-T^{\prime}]\frac{A^{2}}{a_{i}^{2}}$
$\displaystyle=$ $\displaystyle 0,~{}~{}i\in H_{1}^{c}\cap H_{2}.$ (88)

Since $s_{i}>0$, we have for $i\in H_{1}\cap H_{2}^{c}$

$s_{i}=-t_{i}^{\prime}-[S^{\prime}-T^{\prime}]A^{2}/a_{i}^{2}>0,$ (89)

hence

$[T^{\prime}-S^{\prime}]A^{2}/a_{i}^{2}>t_{i}^{\prime}.$ (90)

However, we also have $t_{i}^{\prime}\geq 0$, which indicates that the left-
hand side is positive so that

$T^{\prime}>S^{\prime}.$ (91)

Similarly, for $i\in H_{1}^{c}\cap H_{2}$, we have $t_{i}>0$, that is,

$t_{i}=-s_{i}^{\prime}+[S^{\prime}-T^{\prime}]A^{2}/a_{i}^{2}>0,$ (92)

and hence

$[S^{\prime}-T^{\prime}]A^{2}/a_{i}^{2}>s_{i}^{\prime}.$ (93)

But we also have $s_{i}^{\prime}\geq 0$, which leads to

$S^{\prime}>T^{\prime}.$ (94)

This is a contradiction. ∎

######  Corollary 2.

The subset $H\subset C$ with the properties in Lemma 3 is the optimal subset.

##  3 Experimental Results

In this section we provide the results of the proposed HFC method applied to
simulated and real-world data, and compare them to conventional methods.

###  3.1 Datasets

Twenty simulated data generated using the MixSim R package [10] and five real-
world datasets from the UCI repository [2] were used. The dataset
specifications detailing the number of clusters $k$, cluster distribution
$Pi$, the dimension of data points $p$, the number of data points $N$, and the
average overlap of MixSim datasets are listed in Table 1 and the UCI datasets
Appendicitis, Iris, Seeds, WDBC and Wine in Table2.

Table 1: MixSim generated dataset specifications.

Dataset | ka | Pib | pc | Nd | overlape
---|---|---|---|---|---
Sim2 | 2 | 0.5, 0.5 | 2 | 200 | 0.001
Sim3 | 2 | 0.5, 0.5 | 2 | 200 | 0.01
Sim6 | 2 | 0.5, 0.5 | 2 | 200 | 0.05
Sim9 | 2 | 0.9, 0.1 | 2 | 200 | 0.0005
Sim14 | 2 | 0.9, 0.1 | 2 | 30 | 0.001
Sim8 | 2 | 0.9, 0.1 | 2 | 200 | 0.01
Sim7 | 2 | 0.9, 0.1 | 2 | 200 | 0.05
Sim10 | 2 | 0.8, 0.2 | 2 | 200 | 0.0001
Sim17 | 2 | 0.8, 0.2 | 2 | 50 | 0.005
Sim1 | 2 | 0.2, 0.8 | 2 | 200 | 0.001
Sim4 | 2 | 0.2, 0.8 | 2 | 200 | 0.01
Sim5 | 2 | 0.2, 0.8 | 2 | 200 | 0.05
Sim13 | 2 | 0.7, 0.3 | 2 | 30 | 0.001
Sim18 | 2 | 0.7, 0.3 | 2 | 50 | 0.005
Sim11 | 2 | 0.7, 0.3 | 2 | 200 | 0.01
Sim12 | 2 | 0.7, 0.3 | 2 | 200 | 0.05
Sim19 | 3 | 0.8, 0.1, 0.1 | 2 | 50 | 0.005
Sim15 | 3 | 0.7, 0.2, 0.1 | 2 | 50 | 0
Sim20 | 3 | 0.7, 0.2, 0.1 | 2 | 50 | 0.005
Sim16 | 4 | 0.6, 0.2, 0.1, 0.1 | 2 | 50 | 0

aNumber of clusters. bCluster distribution. cData dimension. dNumber of data
points. eAverage overlap (BarOmega).

Table 2: UCI dataset specifications.

Dataset | ka | Pib | pc | Nd
---|---|---|---|---
Appendicitis | 2 | 0.198, 0.802 | 9 | 106
Iris | 3 | 0.333, 0.333, 0.333 | 4 | 150
Seeds | 3 | 0.333, 0.333, 0.333 | 7 | 210
WDBC | 2 | 0.373, 0.627 | 32 | 569
Wine | 3 | 0.331, 0.399, 0.270 | 13 | 178

aNumber of clusters. bCluster distribution. cData dimension. dNumber of data
points.

###  3.2 Evaluation measures

Four evaluation measures were recorded, external measure Corrected Rand Index
or Adjusted Rand Index (ARI) [4] and internal measures entropy of the
distribution of cluster memberships (Entropy) [9], separation index (SIndex)
and Dunn2 index [5] based on the fpc R package [6].

Corrected Rand Index assessed the similarity between the clustering result and
the ground truth of the same set of objects, adjusted for chance. Values of
near 1 indicate high similarity. Entropy distribution of cluster memberships
is calculated based on the negative of the sum of product of probability and
log probability of the data points belonging to a cluster, across all
clusters. A near zero value indicate little uncertainty in the structure of
the cluster. The separation index is based on the distances of every point to
the closest point not in the same cluster. Larger values indicate
distinctiveness of cluster. The results of the separation index are found in
the supplementary material. Dunn2 calculates the ratio of minimum average
dissimilarity between two clusters and the maximum average within cluster
dissimilarity, where a larger value indicates more distinct clusters.

We also applied the plotcluster function from fpc R package to plot the data
points and the cluster they belong to (indicated by the numbered labels while
the colors ground truth).

###  3.3 Clustering algorithms

We compared the performance of our proposed hybrid fuzzy-crisp HFC algorithm
with two conventional clustering algorithms, K-means (KM) and Fuzzy c-means
(FCM) to investigate in the behaviors of the three algorithms in handling
datasets with imbalanced clusters.

###  3.4 Experiment settings

The experiments were repeated 30 times with each run having a maximum number
of iterations set at 50 and the algorithm stops when the performance by a
threshold of less than 0.0001. $\alpha$ is initialized as $1/c$ and updated at
each iteration but $\beta$ is kept fixed at $1/c$ throughout all iterations.

The purpose of this study is not to investigate in cluster initialization, but
the behavior of the clustering algorithms. To initialize the cluster centers
and avoid bad start, two data points from the same clusters are chosen for
each cluster and averaged. The same cluster centers are used to initialize the
three algorithms HFC, KM and FCM.

For the synthetic datasets generated using MixSim, a different dataset is
generated for each of the 30 runs, each with the same settings as detailed in
Table 1. The reason we did this is to see how well the algorithms perform on
different datasets with same settings since on the same dataset, the
performance usually do not vary much across different runs. This was also
observed on UCI datasets by their low standard deviation in Table 10. The
average ARI, Entropy and Dunn2 results and their standard deviation are
presented. For full results with Sindex measure and their best result, please
refer to the supplementary materials.

###  3.5 MixSim datasets

MixSim is an R package [10] that generates random points from Gaussian
mixtures. We generated various datasets with two to four clusters with
different cluster sizes (Table 1).

#### Balanced datasets

Table 3 shows the performance of HFC, KM and FCM on balanced MixSim datasets
with cluster distribution of (0.5, 0.5). HFC outperforms the other algorithms
in terms of ARI measure on the Sim2 data which has the least overlap while KM
performs best on Sim3 and Sim6. HFC is observed to produce clusters with the
best Entropy in all three datasets. We observed that for balanced dataset
where there is little or no overlap, HFC can be as competitive as KM and FCM,
if not outperforming them as shown in Figure 4 (1a) where more data points at
the edge of the green cluster are correctly assigned by HFC than KM (in 1b)
and FCM (in 1c). Where there is more overlap in Figure 4 (3a), HFC is unable
to distinguish clusters for data points sitting at the boundaries where the
overlap is occurring.

#### Imbalanced datasets

Table 4 shows the performance of HFC, KM and FCM on imbalanced MixSim datasets
with cluster distribution of (0.9, 0.1). HFC performs the best out of all
algorithms in terms of ARI, Entropy and Sindex measures on the Sim9, Sim14,
Sim8 and Sim7 data which have increasing amount of overlap from 0.0005 for
Sim9 to 0.05 for Sim7. As higher degree of overlap is introduced, the
performance decrease for all algorithms, with HFC performing the best. This
demonstrates that HFC outperforms K and FCM when the dataset is highly
imbalanced. Figure 5 shows the clustering results on datasets of the two
degrees of overlapping, Sim9 with very slight overlap of 0.0005 and Sim7 with
slight overlap of 0.05. In both datasets, we observed visually that HFC were
able to cluster slightly better than KM while FCM performed much worse than
the other two algorithms. Similar trends were found using simulated datasets
with cluster distribution of (0.8, 0.2) in Table 5 and Figure 6.

In Table 6 and Figure 7, we observed that HFC can still outperform KM and FCM
in datasets with cluster distribution (0.7, 0.3). However, on Sim12 data with
slightly more overlap, KM and FCM outperforms HFC. This demonstrate that HFC
can handle imbalanced dataset with very slight overlapping ($<0.01$).

In Table 7 and Figure 8, we observed the results on Sim19 with cluster
distribution of (0.8, 0.1, 0.1) and slight overlap where HFC outperformed KM
and FCM on average across 30 such datasets on all measures, achieving the
maximum average ARI, Sindex and Dunn2 measures and minimum average Entropy. As
there is great disparity of cluster ratios of 0.8 and 0.1, we observed in 8
that KM and FCM wrongly assign labels 2 and 3 to the blue cluster while HFC
wrongly labels a few data points of the blue cluster to 2.

Table 8 and Figure 9 show results on imbalanced datasets with 3 clusters and
cluster distribution of (0.7, 0.2, 0.1). KM and FCM outperformed HFC in Sim15
which has no overlap while KM outperformed HFC in Sim20 which has 0.005
overlap. This may indicate that as disparity of cluster distribution reduces,
KM and FCM will outperform HFC. As the number of cluster increases, there may
be a tendency for cluster distribution differences to be smaller such that KM
and FCM will outperform HFC, as shown in Table 9 and Figure 10.

Table 3: Evaluation on unique MixSim datasets with Pi = (0.5, 0.5) (the average overlap value indicated in parentheses). See also Figure 4.  |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim2 | Mean | 0.831 | 0.772 | 0.73 | 0.538 | 0.561 | 0.573 | 2.14 | 2.099 | 2.079
(0.001) | StDev | 0.264 | 0.285 | 0.296 | 0.068 | 0.076 | 0.075 | 0.421 | 0.428 | 0.443
Sim3 | Mean | 0.798 | 0.804 | 0.724 | 0.516 | 0.527 | 0.554 | 2.104 | 2.077 | 1.983
(0.01) | StDev | 0.257 | 0.244 | 0.28 | 0.087 | 0.076 | 0.084 | 0.549 | 0.548 | 0.56
Sim6 | Mean | 0.673 | 0.684 | 0.682 | 0.675 | 0.678 | 0.681 | 1.747 | 1.756 | 1.749
(0.05) | StDev | 0.18 | 0.173 | 0.149 | 0.023 | 0.017 | 0.014 | 0.276 | 0.268 | 0.27
Figure 4: Cluster plots of one run on 1) Sim2 (overlap of 0.001) 2) Sim3 (overlap of 0.01) and 3) Sim6 (overlap of 0.05) with cluster distribution (0.5, 0.5). The ARI values are reported in parentheses. See also Table 3.  Table 4: Evaluation on unique MixSim datasets with Pi = (0.9, 0.1) (the average overlap value indicated in parentheses). See also Figure 5 |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim9 | Mean | 0.963 | 0.959 | 0.849 | 0.342 | 0.347 | 0.391 | 2.764 | 2.721 | 2.632
(0.0005) | StDev | 0.131 | 0.132 | 0.318 | 0.063 | 0.061 | 0.118 | 0.817 | 0.847 | 0.894
Sim14 | Mean | 0.851 | 0.831 | 0.681 | 0.433 | 0.438 | 0.502 | 2.302 | 2.306 | 2.106
(0.001) | StDev | 0.319 | 0.323 | 0.42 | 0.122 | 0.126 | 0.137 | 0.745 | 0.742 | 0.798
Sim8 | Mean | 0.797 | 0.715 | 0.57 | 0.376 | 0.419 | 0.48 | 2.358 | 2.215 | 2.033
(0.01) | StDev | 0.329 | 0.371 | 0.435 | 0.139 | 0.156 | 0.172 | 0.568 | 0.611 | 0.593
Sim7 | Mean | 0.465 | 0.403 | 0.3 | 0.509 | 0.541 | 0.594 | 1.828 | 1.765 | 1.656
(0.05) | StDev | 0.391 | 0.394 | 0.362 | 0.173 | 0.163 | 0.14 | 0.384 | 0.304 | 0.278
Figure 5: Cluster plots of one run on 1) Sim9 (overlap of 0.0005) and 2) Sim7 (overlap of 0.05) with cluster distribution (0.9, 0.1). The ARI values are reported in parentheses. See also Table 4. Table 5: Evaluation on unique MixSim datasets with Pi = (0.8, 0.2) (the average overlap value indicated in parentheses). See also Figure 6. |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim10 | Mean | 0.967 | 0.952 | 0.943 | 0.504 | 0.51 | 0.512 | 2.581 | 2.546 | 2.542
(0.0001) | StDev | 0.091 | 0.106 | 0.105 | 0.048 | 0.048 | 0.049 | 0.642 | 0.652 | 0.656
Sim17 | Mean | 0.94 | 0.903 | 0.869 | 0.506 | 0.516 | 0.53 | 2.206 | 2.204 | 2.149
(0.005) | StDev | 0.105 | 0.135 | 0.179 | 0.063 | 0.07 | 0.074 | 0.495 | 0.489 | 0.49
Sim1 | Mean | 0.87 | 0.816 | 0.797 | 0.517 | 0.536 | 0.547 | 2.236 | 2.189 | 2.157
(0.001) | StDev | 0.17 | 0.227 | 0.21 | 0.053 | 0.064 | 0.063 | 0.591 | 0.605 | 0.62
Sim4 | Mean | 0.871 | 0.841 | 0.808 | 0.512 | 0.523 | 0.545 | 2.131 | 2.092 | 2.005
(0.01) | StDev | 0.188 | 0.198 | 0.231 | 0.064 | 0.068 | 0.062 | 0.512 | 0.497 | 0.526
Sim5 | Mean | 0.645 | 0.616 | 0.611 | 0.521 | 0.562 | 0.582 | 1.837 | 1.727 | 1.672
(0.05) | StDev | 0.285 | 0.299 | 0.301 | 0.095 | 0.086 | 0.078 | 0.291 | 0.314 | 0.3
Figure 6: Cluster plots of one run on 1) Sim10 (overlap of 0.0001) and 2) Sim5 (overlap of 0.05) with cluster distribution (0.8, 0.2). The ARI values are reported in parentheses. See also Table 5. Table 6: Evaluation on unique MixSim datasets with Pi = (0.7, 0.3) (the average overlap value indicated in parentheses). See also Figure 7. |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim13 | Mean | 0.91 | 0.89 | 0.886 | 0.602 | 0.614 | 0.613 | 2.284 | 2.289 | 2.291
(0.001) | StDev | 0.194 | 0.231 | 0.214 | 0.088 | 0.081 | 0.083 | 0.661 | 0.662 | 0.656
Sim18 | Mean | 0.862 | 0.859 | 0.817 | 0.593 | 0.605 | 0.607 | 2.084 | 2.033 | 2.055
(0.005) | StDev | 0.154 | 0.164 | 0.208 | 0.059 | 0.053 | 0.054 | 0.557 | 0.546 | 0.552
Sim11 | Mean | 0.803 | 0.785 | 0.768 | 0.615 | 0.619 | 0.633 | 1.927 | 1.957 | 1.917
(0.01) | StDev | 0.237 | 0.24 | 0.223 | 0.047 | 0.048 | 0.045 | 0.409 | 0.38 | 0.41
Sim12 | Mean | 0.631 | 0.673 | 0.658 | 0.598 | 0.61 | 0.624 | 1.934 | 1.913 | 1.877
(0.05) | StDev | 0.264 | 0.235 | 0.231 | 0.078 | 0.073 | 0.057 | 0.339 | 0.33 | 0.347
Figure 7: Cluster plots of one run on 1) Sim13 (overlap of 0.001) and 2) Sim12 (overlap of 0.05) with cluster distribution (0.7, 0.3). The ARI values are reported in parentheses. See also Table 6.  Table 7: Evaluation on unique MixSim datasets with Pi = (0.8, 0.1, 0.1) (the average overlap value indicated in parentheses). See also Figure 8. |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim19 | Mean | 0.879 | 0.872 | 0.817 | 0.7 | 0.708 | 0.741 | 2.348 | 2.342 | 2.246
(0.005) | StDev | 0.21 | 0.215 | 0.273 | 0.141 | 0.142 | 0.162 | 0.694 | 0.679 | 0.734
Table 8: Evaluation on unique MixSim datasets with Pi = (0.7, 0.2, 0.1) (the average overlap value indicated in parentheses). See also Figure 9. |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim15 | Mean | 0.929 | 0.932 | 0.933 | 0.849 | 0.854 | 0.855 | 4.412 | 4.409 | 4.41
(0) | StDev | 0.211 | 0.202 | 0.199 | 0.09 | 0.1 | 0.101 | 1.781 | 1.785 | 1.784
Sim20 | Mean | 0.824 | 0.839 | 0.788 | 0.846 | 0.869 | 0.882 | 1.997 | 1.972 | 1.958
(0.005) | StDev | 0.218 | 0.194 | 0.242 | 0.12 | 0.111 | 0.117 | 0.681 | 0.675 | 0.661
Figure 8: Cluster plots of one run on 1) Sim19 (overlap of 0.005) with cluster distribution (0.8, 0.1, 0.1). The ARI values are reported in parentheses. See also Table 7. Figure 9: Cluster plots of one run on 1) Sim15 (overlap of 0) and 2) Sim20 (overlap of 0.005) with cluster distribution (0.7, 0.2, 0.1). The ARI values are reported in parentheses. See also Table 8. Table 9: Evaluation on unique MixSim datasets with Pi = (0.6, 0.2, 0.1, 0.1) (the average overlap value indicated in parentheses). See also Figure 10. |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Sim16 | Mean | 0.996 | 1 | 1 | 1.077 | 1.076 | 1.076 | 3.956 | 4.054 | 4.054
(0) | StDev | 0.018 | 0 | 0 | 0.104 | 0.104 | 0.104 | 1.71 | 1.606 | 1.606
Figure 10: Cluster plots of one run on 1) Sim16 (overlap of 0) with cluster
distribution (0.6, 0.2, 0.1, 0.1). The ARI values are reported in parentheses.
See also Table 9.

###  3.6 Real-world data

In this section, we present and analyze the results from testing selected
real-world datasets. In Table 10, HFC outperformed KM and FCM on the
Appendicitis, WDBC and Seeds datasets in terms of ARI measures.

The Appendicitis dataset is highly imbalanced and HFC outperformed the others
in all four measures. Studying Figure 11 (1a) to (1c), we observed that many
that belong to the green class are incorrectly assigned to cluster 1 for KM
(1b) and FCM (1c).

For WDBC, there is considerable overlapping between the two clusters as shown
in Figure 11 (2a) to (2c). Despite this, HFC was able to overcome the
overlapping while handling the imbalanced dataset, demonstrating HFC’s ability
to handle imbalanced datasets.

Despite Seeds being a balanced dataset, HFC outperformed KM and FCM. Studying
the cluster plots on one run in Figure 12 (2a) to (2c), those that belong to
the purple group are more often incorrectly assigned to cluster 1 for KM and
FCM than HFC.

From Table 10, HFC did not outperform KM and FCM in Iris and Wine datasets.
While ARI values showed competitive performance, with the near balanced and
overlapping nature of the clusters, KM and FCM performed better on such
dataset, particularly for Iris (see Figure 12 (1a) to (1c)) where there is
more overlapping between the green and purple groups.

Table 10: Evaluation of algorithms on UCI datasets. Their cluster distribution are indicated in parentheses.  |  | ARI | Entropy | Dunn2
---|---|---|---|---
Dataset |  | HFC | KM | FCM | HFC | KM | FCM | HFC | KM | FCM
Appendicitis | Mean | 0.506 | 0.3 | 0.294 | 0.511 | 0.607 | 0.604 | 1.731 | 1.47 | 1.475
(0.2, 0.8) | StDev | 0 | 0.05 | 0 | 0 | 0.029 | 0 | 0 | 0.054 | 0
Iris | Mean | 0.608 | 0.62 | 0.63 | 1.046 | 1.097 | 1.098 | 1.478 | 1.6 | 1.578
(0.33, 0.33, 0.33) | StDev | 0.016 | 0 | 0 | 0.018 | 0 | 0 | 0.06 | 0 | 0
Seeds | Mean | 0.793 | 0.784 | 0.772 | 1.098 | 1.098 | 1.098 | 1.611 | 1.635 | 1.641
(0.33, 0.33, 0.33) | StDev | 0.024 | 0.012 | 0 | 0.001 | 0 | 0 | 0.048 | 0.019 | 0
WDBC | Mean | 0.694 | 0.663 | 0.683 | 0.631 | 0.638 | 0.647 | 1.212 | 1.196 | 1.184
(0.37, 0.63) | StDev | 0 | 0.009 | 0 | 0 | 0.003 | 0 | 0 | 0.007 | 0
Wine | Mean | 0.894 | 0.897 | 0.897 | 1.095 | 1.093 | 1.093 | 1.255 | 1.287 | 1.287
(0.33, 0.4, 0.27) | StDev | 0.023 | 0 | 0 | 0.001 | 0 | 0 | 0.004 | 0 | 0
Figure 11: Cluster plots of one run on imbalanced UCI datasets 1) Appendicitis
with cluster distribution (0.2, 0.8) 2) WDBC with cluster distribution (0.37,
0.63) and 3) Wine with cluster distribution (0.33, 0.4, 0.27). The ARI values
are reported in parentheses.  Figure 12: Cluster plots of one run on balanced
UCI datasets 1) Iris and 2) Seeds both with cluster distribution (0.33, 0.33,
0.33). The ARI values are reported in parentheses.

###  3.7 Convergence of the HFC algorithm

In Figure 13, we presented the fitness value calculated by the objective
function against the iterations for selected datasets which HFC perform both
well and poorly in. The datasets are Sim7, Sim9, Appendicitis, Iris, Seeds and
WDBC. We observed that HFC convergences at certain iterations except for HFC
on Iris which has reached the maximum iteration of 50.

Figure 13: Convergence plot of one run on selected datasets.

##  4 Conclusion

We introduced a hybrid fuzzy-crisp (HFC) clustering approach to solve the
problem of centers of other small clusters being pulled towards outstandingly
larger clusters. In the HFC algorithm, the membership of a data point to a
cluster becomes exactly zero when the data point is deemed “far” from the
center. We showed that the HFC clustering can be interpreted as the problem of
findind the closest point in a $(c-1)$-simplex from the origin in the
$c$-dimensional Euclidean space, and proved the correctness of the algorithm.
The HFC algorithm was then tested on simulated and real data of varying
cluster sizes and degrees of overlapping between clusters. Based on the
results and analysis, HFC has been demonstrated to perform better on
imbalanced datasets and can be competitive with KM and FCM on more balanced
datasets. This was particularly evident on the two-cluster dataset with the
mix ratios of (0.9, 0.1) and (0.8,0.2) where HFC outperformed KM and FCM. On
dataset with mix ratios of (0.7,0.3), HFC can outperform KM and FCM where
degrees of overlapping is small (average overlap $\leq 0.01$). Through
evaluation measures such as ARI and Entropy as well as visual inspection via
t-SNE biplots, we observed such performance from HFC.

The experiments revealed the relative weakness of HFC in handling datasets
that are well-balanced and/or largely overlapping, compared to other methods.
One possible way to overcome this weakness may be to optimize the relative
weights between the linear and quadratic terms ($\beta$ and
$\\{\alpha_{i}\\}$, respectively) in the objective function (Eq. 7). This is
left for future studies.

## Acknowledgments

A. R. K. thanks Nobuhiro Go and Yusuke Umatani for suggesting this problem.
Software used in this study is available at https://github.com/arkinjo/fccm.

## CRediT author statement

AR Kinjo: Conceptualization, Methodology, Formal Analysis, Validation, Writing
- Original draft preparation, Writing - Reviewing and Editing, Visualization,
Supervision

DTC Lai: Software, Validation, Investigation, Data Curation, Writing- Original
draft preparation, Writing - Reviewing and Editing

## References

  * Bezdek et al. [1984] J. C. Bezdek, R. Ehrlich, and W. Full.  FCM: The fuzzy c-means clustering algorithm.  _Comput. geosci._ , 10(2-3):191–203, 1984.
  * Dua and Graff [2017] D. Dua and C. Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Dunn [1973] J. C. Dunn.  A fuzzy relative of the ISODATA process and its use in detecting compact well-separated clusters.  _J. Cybern._ , 3:32–57, 1973.
  * Gordon [1999] A. D. Gordon.  _Classification_.  CRC Press, 1999.
  * Halkidi et al. [2001] M. Halkidi, Y. Batistakis, and M. Vazirgiannis.  On clustering validation techniques.  _Journal of Intelligent Information Systems_ , 17(2-3):107–145, 2001.
  * Hennig [2018] C. Hennig.  fpc: Flexible procedures for clustering. 2015.  _R package version_ , pages 2–1, 2018.
  * Kinjo [1998] A. Kinjo.  A fuzzy description of local backbone structures of proteins.  Master’s thesis, Kyoto University, 1998.
  * Klawonn and Höppner [2003] F. Klawonn and F. Höppner.  What is fuzzy about fuzzy clustering? understanding and improving the concept of the fuzzifier.  In _International symposium on intelligent data analysis_ , pages 254–264. Springer, 2003.
  * Meilă [2007] M. Meilă.  Comparing clusterings—an information based distance.  _Journal of multivariate analysis_ , 98(5):873–895, 2007.
  * Melnykov et al. [2012] V. Melnykov, W.-C. Chen, and R. Maitra.  MixSim: An R package for simulating data to study performance of clustering algorithms.  _Journal of Statistical Software_ , 51(12):1–25, 2012.  URL https://doi.org/10.18637/jss.v051.i12.
  * Ruspini et al. [2019] E. H. Ruspini, J. C. Bezdek, and J. M. Keller.  Fuzzy clustering: A historical perspective.  _IEEE Comput. Intell. Mag._ , 14(1):45–55, 2019\.