# The extraordinary frequency pattern variation in $\delta$ Scuti stars

L. A. Balona  South African Astronomical Observatory, P.O. Box 9, Observatory
7935, Cape Town, South Africa

###### Abstract

Inspection of the periodograms of TESS $\delta$ Scuti stars indicates that
there is little, if any, similarity between the frequencies of stars in the
same region of the H–R diagram. This is difficult to understand because
pulsation models predict that stars with similar physical parameters should
have similar frequencies. To investigate the problem, a quantitative measure
of similarity between frequency patterns is described. When applied to non-
adiabatic pulsation models with similar temperatures and luminosities, a
strong correlation is found between the frequency patterns, as expected. The
correlation increases when rotational frequency splitting is included. When
applied to observations, no correlation can be found, confirming the
impression from visual inspection. It seems that each star has its own unique
frequency pattern, unrelated to its position in the H–R diagram. This presents
a problem in our understanding of stellar pulsation. The existence of a
period-luminosity law and the effect of combination frequencies in $\delta$
Scuti stars is briefly discussed.

###### keywords:

stars:oscillations — stars:variables:$\delta$ Scuti

##  1 Introduction

Pulsational driving in models of A and F stars depends on the effective
temperature, $T_{\rm eff}$, luminosity, $L/L_{\odot}$, and chemical
composition. Rotation has the effect of introducing additional frequencies.
For stars located close together within a small region in the H–R diagram, the
frequency spectra should be similar for stars with the same rotation rate.
Differences in chemical composition will also affect the frequencies, but to a
lesser extent. Furthermore, all stars within the instability strip should
pulsate.

Inspection of the periodograms of several thousand $\delta$ Sct stars observed
by TESS gives a strong impression that the above conclusions are not correct.
Bedding et al. (2023) has recently remarked that the characteristics of the
pulsation spectra in TESS $\delta$ Sct stars are varied and do not correlate
with stellar temperature. More than half of the stars within the instability
strip do not seem to pulsate at all (Balona, 2018). The frequency patterns in
$\gamma$ Dor stars are very different from $\delta$ Sct stars, even though
they might have the same effective temperature and luminosity. This cannot be
attributed to rotation because the distribution of projected rotational
velocities for $\gamma$ Dor and $\delta$ Sct stars are the same.

Pulsation models cannot account for these observations. A better understanding
of why certain pulsation modes are selected in preference to others is
required. The study of mode selection in pulsating stars involves non-
adiabatic, non-linear calculations for nonradial modes in rotating stars,
which is beyond current capabilities. It is no surprise that scant attention
has been payed to this problem. The review by Dziembowski (1993) is still
mostly relevant. For more recent work, see the excellent review by Smolec
(2014).

The aim of this paper is to present a sample of periodograms to illustrate the
problem. A method which measures the difference in frequency patterns in a
quantitative way, while allowing for effects of rotation and chemical
abundance, is described. This is applied to TESS $\delta$ Sct stars in several
regions of the H–R diagram. The size of each region is sufficiently small that
stars within it may be considered to have the same effective temperature and
luminosity.

The calculations indicate that frequency patterns of stars within each region
do not match at all. It is shown that rotation cannot account for the
disparity. This means that there cannot be a general period - luminosity law
for $\delta$ Sct stars (Poretti et al., 2008; McNamara, 2011; Garg et al.,
2010; Poleski et al., 2010; Cohen and Sarajedini, 2012; Ziaali et al., 2019).
The idea that combination frequencies may contribute to the disparity in
frequency spectra can also be discarded.

##  2 Data and variability classification

Each TESS sector consists of about 27 d of continuous photometry with 2-min
cadence. While stars near the ecliptic equator are only observed for one
sector, stars near the ecliptic poles are observed in every sector. The light
curves are corrected for instrumental signatures and long-term drift. These
are called pre-search data conditioning (PDC) light curves. The data are from
sectors 1–66 of the TESS mission.

As each sector becomes available, the light curves and periodograms of stars
with $T_{\rm eff}>6000$ K are used to classify them according to variability
type. The classification scheme follows that of the General Catalogue of
Variable Stars (GCVS, Samus et al. 2017). In this way, many thousands of stars
have been classified. A full description of the identification and
classification of TESS variables is given in Balona (2022). As expected, most
of the variables are of the $\delta$ Sct type. Among the 125000 stars already
classified, there are 14000 $\delta$ Sct stars.

The TESS data present several problems in the classification of light
variability. It seems that distinct regions of instability, as predicted by
the models, do not actually exist. There is no region in the H–R diagram
entirely free of pulsating stars (Balona and Ozuyar, 2020). In order to assign
a meaningful variability class, arbitrary boundaries in effective temperature,
$T_{\rm eff}$, and pulsation frequency, $\nu$, must be defined. For example,
the boundary between the $\delta$ Sct and Maia stars is set at $T_{\rm
eff}=10000$ K. The periodograms of these two classes of variable look the
same: the Maia stars seem to be just an extension of the $\delta$ Sct
instability strip to the early B stars (Balona, 2023). The arbitrary
boundaries in $T_{\rm eff}$ and $\nu$ are forced by the lack of a suitable
theory, but are required to define the traditional classification scheme.

##  3 Variety of frequency patterns in $\delta$ Sct stars

To illustrate the variety of frequency patterns in $\delta$ Sct stars with
similar stellar parameters, the instability strip was divided into 12 regions
parallel to the zero-age main sequence (ZAMS), each region having a size
$(\Delta\log T_{\rm eff},\Delta\log L/L_{\odot})=(0.02,0.25)$. This is
sufficiently small that one may regard all stars within a region to have
approximately the same effective temperature and luminosity, while at the same
time sufficiently large to include a suitable number of stars. The region size
also matches the uncertainties: $\sigma(\log T_{\rm eff})\approx 0.2$ and
$\sigma(\log L/L_{\odot})\approx 0.1$. Luminosities are from GAIA DR3
parallaxes (Gaia Collaboration et al., 2016, 2021). Parameters of each region
are listed in Table 1. Fig. 1 shows the names and locations of the regions in
the H–R diagram.

Figure 1: The top panel shows the location of the regions in the H–R diagram and the zero-age main sequence. The bottom panel shows, in addition, the $\delta$ Sct stars detected by TESS. Table 1: Data for rotating and non-rotating models located within a small region in the H–R diagram. The first column is the region number followed by the central $\log T_{\rm eff}$ and $\log L/L_{\odot}$ for the region. The last column is the number of models, $N_{\rm mod}$, or $N_{\rm obs}$, the number of stars, within each region. Region | $\log T_{\rm eff}$ | $\log\tfrac{L}{L_{\odot}}$ | $N_{\rm mod}$ | $N_{\rm obs}$
---|---|---|---|---
(1,1) | 3.85 | 0.695 | 13 | 129
(1,2) | 3.85 | 0.945 | 5 | 747
(1,3) | 3.85 | 1.195 | 10 | 468
(2,1) | 3.87 | 0.821 | 7 | 726
(2,2) | 3.87 | 1.071 | 4 | 1048
(2,3) | 3.87 | 1.321 | 14 | 459
(3,1) | 3.89 | 0.943 | 8 | 1099
(3,2) | 3.89 | 1.193 | 5 | 748
(3,3) | 3.89 | 1.443 | 14 | 356
(4,1) | 3.91 | 1.063 | 6 | 899
(4,2) | 3.91 | 1.313 | 5 | 468
(4,3) | 3.91 | 1.563 | 15 | 206

Evolutionary stellar models were computed using the Warsaw - New Jersey
evolution code (Paczyński, 1970), assuming an initial hydrogen fraction,
$X_{0}=0.70$ and metal abundance, $Z=0.020$ and using the chemical element
mixture of Asplund et al. (2009) and OPAL opacities (Rogers and Iglesias,
1992) Overshooting from the convective core was not included. A mixing length
parameter $\alpha_{\rm MLT}=1.0$ was adopted for the convective scale height.
All models are non-rotating. The non-adiabatic code developed by Dziembowski
(1977) was used to obtain the pulsation frequencies and growth rates.
Pulsations in these models are purely driven by the opacity $\kappa$ mechanism
operating in the He II partial ionization region.

Figure 2: Illustration of the variety of frequency patterns in $\delta$ Sct
stars compared with theoretical models. Each box consists of three panels. The
top panel shows stars with rich frequency spectra, while stars with just a
single dominant frequency peak are shown in the bottom panel. The middle panel
(green) shows corresponding unstable frequencies from models. The label on the
top of each box refers to a small region in the H–R diagram (Table 1). Figure
3: Examples of TESS $\delta$ Scuti periodograms for stars within a small
region of the H–R diagram centered on $\log T_{\rm eff}=3.87$, $\log
L/L_{\odot}=1.07$ and with box size $0.02\times 0.025$ dex.

Examples of the disparity of frequency patterns within a given region are
illustrated in Fig. 2. In each region, examples of three stars with rich
frequency spectra and three stars with just a single dominant frequency peak
are shown. All stars have good spectroscopic estimates of $T_{\rm eff}$. The
middle row in Fig. 2 shows unstable modes from non-rotating models. The
amplitudes were made to be proportional to the growth rate and corrected
according to the visibility factor for the given spherical harmonic degree,
$l$. Modes with $l\leq 8$ were selected.

The figure shows that stars with the same physical parameters can have
completely different frequency patterns. Frequencies of single-mode stars may
be very different from star to star within the same region. This indicates
that the spherical harmonic degree in these single-mode stars must differ from
star to star. It is evident that the pattern in stars with rich spectra cannot
be obtained from rotational frequency splitting of the single-mode pulsators.

While Fig.2 shows stars across the H–R diagram, Fig. 3 shows a larger sample
of randomly chosen stars with similar, well determined, effective temperatures
and luminosities in region (2,2). Each star seems to have its own unique
frequency pattern.

##  4 Correlation

While a subjective impression of the disparity of frequency patterns can be
obtained by inspection of Figs. 2 and 3, a method which quantitatively
measures the degree of similarity between a pair of frequency patterns is
required. One possibility is to calculate the correlation coefficient.

Suppose there are two frequency spectra, $i$ and $j$, each sampled with equal
frequency spacing, $\delta\nu$, between frequencies $\nu_{1}$ and $\nu_{2}$.
Let $y_{i}(\nu)$ and $y_{j}(\nu)$ be the amplitudes at frequency $\nu$ of
spectra $i$ and $j$. The correlation coefficient, $r$, between the two spectra
is given by

$\displaystyle
r=\frac{1}{A_{i}A_{j}}\sum_{\nu=\nu_{1}}^{\nu_{2}}y_{i}(\nu)y_{j}(\nu)\delta\nu,$

where $A_{i}$, $A_{j}$ are constants. The values of these constants are chosen
so that correlating a periodogram with itself gives unity. Then

$\displaystyle
A_{i}^{2}=\sum_{\nu=\nu_{1}}^{\nu_{2}}y_{i}^{2}(\nu)\delta\nu,~{}~{}A_{j}^{2}=\sum_{\nu=\nu_{1}}^{\nu_{2}}y_{j}^{2}(\nu)\delta\nu.$

The value of $r$ lies between $r=0$ (no similarity) and $r=1$ (periodograms
are identical).

Given any two stars, the extracted pulsation frequencies will never be exactly
the same. This means that the value of $r$ will always be near zero if the
sampling frequency is very small. This problem is avoided by constructing a
histogram of the frequency distribution with a suitable frequency bin size.
Each extracted frequency peak can either be given the same weight or weighted
according to its amplitude.

Examination of pulsation models within a range of 0.02 dex in $\log T_{\rm
eff}$ and 0.25 dex in $\log L/L_{\odot}$ shows that for a fixed spherical
harmonic degree and radial order $(l,n)$, the rms scatter in the frequencies
is typically 2–6 d-1. This can be higher for g modes or mixed modes. A change
in abundance leads to a somewhat different convective core size. Since g modes
mostly sample the interior of a star near the convective core, there can be
changes in frequencies of some g modes, but these are only of the order of
several microhertz (Guzik et al., 1998). This is smaller than the effect of
rotation. Other effects, such as mass loss and the presence of a magnetic
field, do not alter the mean density and have no measurable effect on the
pulsation frequencies. These considerations suggest that a bin size of about
$\delta\nu=5$ d-1 will be an appropriate choice.

##  5 Correlation from pulsation models

A small change in effective temperature and luminosity in a model of a non-
rotating star will lead to small changes in the predicted frequencies. It is
therefore expected that the frequency distributions derived from models with
similar values of ($\log T_{\rm eff}$, $\log L/L_{\odot}$) will be strongly
correlated. The aim of this section is to show that this is indeed the case.
It is also important to study how the correlation coefficient might change
between models with different rotation rates.

Figure 4: Distribution of correlation coefficients for models in each region.
The top half are for non-rotating models and the bottom half for models with a
distribution of rotation rates typical of $\delta$ Sct stars.

For a model in a particular region, the correlation coefficient between the
frequency pattern of unstable modes was correlated with the frequency pattern
for all other models in the region. The number of model for each region,
$N_{\rm mod}$, is shown in Table 1. Each pair of models gives rise to one
correlation coefficient, $r$. The resulting distribution of $r$ for all pairs
of models is shown in the top panel of Fig. 4. As expected, the distribution
of $r$ tends towards $r=1$.

To test the effect of rotation, frequency splitting was applied to each
unstable mode with spherical harmonic degree $l$, so that a frequency in the
non-rotating model is replaced by $2l+1$ frequencies of equal amplitude. Reese
et al. (2006) provides a formula for third-order rotational perturbation as
well as the necessary coefficients for $l\leq 3$. For higher values of $l$,
coefficients for $l=3$ were used. The distribution of rotational equatorial
velocities for mid-A stars was used to derive the distribution of rotational
frequencies required in the formula.

The bottom panel of Fig. 4 shows the distribution of $r$ within each region
when rotational splitting is applied. Rotation has the effect of introducing
additional frequency peaks. These extra peaks tend to broaden the overall
frequency pattern by an amount roughly proportional to the rotational
frequency. As a result, the overlap in the frequency distributions increases,
leading to higher correlation.

The models indicate that one should expect a high level of correlation of
frequency patterns for $\delta$ Sct stars within a given region. Furthermore,
it suggests that differing rotation rates among the stars actually contribute
to enhance correlation.

##  6 Correlation from observations

Frequencies of $\delta$ Sct stars were extracted from the periodograms. The
extracted frequency is deemed significant if the signal-to-noise ratio in the
amplitude exceeding 4.7. Each significant extracted frequency peak was
assigned a weight proportional to its amplitude when calculating the histogram
of the frequency distribution. As before, the correlation coefficients, $r$,
were obtained for all possible pairs of stars within a given region. Finally,
the distribution of $r$ for each region was determined.

The calculation was performed for 1673 stars with precise measurements of
effective temperature as well as 7353 stars for which any estimate of
effective temperature was accepted. The results are very similar. In Table 1,
$N_{\rm obs}$ is the number of stars in each region for the latter case. The
distributions of correlation coefficients are shown in Fig. 5.

Figure 5: Distribution of correlation coefficients for TESS $\delta$ Sct stars
within a small region in the H–R diagram.

Fig. 5 shows that, contrary to the model predictions, there is no correlation
between the frequency patterns in $\delta$ Sct stars with similar effective
temperatures and luminosities. Moreover, the variety of frequency patterns
cannot be explained by rotation.

##  7 Mode coupling

Mode coupling involving nonlinear mode interactions between two or more modes
may result in a complex frequency spectrum which could account for the wide
disparity in the frequency patterns in $\delta$ Sct stars. Mourabit and
Weinberg (2023) used the theory of three-mode coupling to study the strength
and prevalence of nonlinear mode interactions in $\delta$ Sct models across
the instability strip. They conclude that resonant mode interactions can be
signiﬁcant. Mode coupling may explain the rapid changes in amplitude and
frequency seen in some $\delta$ Sct stars (Bowman et al., 2016, 2021).

A combination frequency involving two parent modes, $\nu_{1},\nu_{2}$ is given
by $\nu=n_{1}\nu_{1}+n_{2}\nu_{2}$ where $n_{1}$ and $n_{2}$ are positive or
negative integers. If the frequency $\nu$ matches an observed frequency peak
to within the observational error, then $\nu$ can be regarded as a combination
frequency. To test the degree of mode coupling that may be present in $\delta$
Sct stars, all possible combination frequencies with $|n_{1}|\leq 7$ and
$|n_{2}|\leq 7$ were identified among 6840 $\delta$ Sct stars.

It was found that at least one or more combination frequencies are present in
4151 stars (60 percent of the total), but in most of these stars only a small
fraction of the frequency peaks are identified as combination frequencies. The
distribution of stars with combination frequencies is shown in Fig. 6.

Figure 6: Fraction of $\delta$ Sct stars as a function of percentage of
frequencies that are combinations.

Higher order combinations involving three or even more modes could occur, but
these would be difficult to identify. Considering that very few combination
frequency peaks were identified in most stars, it is very unlikely that mode
interaction may be the explanation for the diversity of frequency patterns.

##  8 Period–luminosity relation

Considering the large disparity in frequency patterns at constant temperature
and luminosity, it would be surprising if a period–luminosity law, similar to
the Cepheid P–L law, can be found for $\delta$ Sct stars. Nevertheless, the
existence of a P–L relation in $\delta$ Sct stars has been reported from time
to time. Since there are usually many periods to be chosen in any given
$\delta$ Sct star, the period corresponding to the highest amplitude is used.
Barac et al. (2022) found that stars lie on a ridge that corresponds to
pulsation in the fundamental radial mode, while others pulsate in shorter
period overtones. The correlation is best seen in the period - mean density
relationship.

Figure 7: Effective temperature and luminosity as a function of pulsation
period, $P$ (days), using the period which corresponds to the peak with
highest amplitude in TESS $\delta$ Sct stars. The left panels include all
stars, while the panels on the right are restricted to stars with one dominant
frequency.

It should be borne in mind is that there are several correlations between the
pulsation period, $P$, and the physical parameters. Pulsations in $\delta$ Sct
stars are mainly driven by the $\kappa$ opacity mechanism operating in the He
II partial ionization zone. In cool stars, the zone is deep, the resonant
cavity is large and the radial order of the pulsations is low, so the periods
are typically long. In hot stars, the zone is shallow, the resonant cavity is
small and the radial order is high, so the periods are shorter. This effect
can be seen in the model frequencies in Fig. 2.

From the observations of 11833 $\delta$ Sct stars, the frequency of maximum
amplitude was extracted. The left panels of Fig. 7 shows the relationship
between $\log T_{\rm eff}$ and $\log P$ (top panel) and between $\log
L/L_{\odot}$ and $\log P$ (bottom panel). As predicted by the models, the
pulsation period decreases as the temperature increases, though the
correlation is not particularly strong. However, there is no significant
correlation between $\log P$ and $\log L/L_{\odot}$.

It is possible that stars with a single dominant mode may be best suited to
test the idea of a P–L relationship because these could be radial modes (or at
least nonradial modes of the same degree). The corresponding relationships for
these stars are shown in the right panels of Fig. 7. Once again, the only
correlation is between $\log T_{\rm eff}$ and $\log P$. This indicates that
modes in stars with only one dominant frequency do not have the same value of
$l$ or radial order.

In conclusion, there seems to be a significant correlation between the
dominant pulsation period and effective temperature in $\delta$ Sct stars, in
agreement with model predictions. However, a period - luminosity law similar
to that in Cepheids does not exist.

##  9 Discussion and Conclusions

The surprisingly wide variety of frequency patterns in $\delta$ Sct stars has
been noted several times (e.g. Bedding et al. 2023). Inspection of many
thousands of periodograms gives the impression that each $\delta$ Sct star has
a unique set of frequencies, unrelated to its location in the H–R diagram.

In this paper, periodograms of $\delta$ Sct stars observed by TESS are
presented. These are shown for stars that lie within small regions of the H–R
diagram. In each region, where stars have essentially the same effective
temperature and luminosity, stars can be found with just a single frequency
peak and also with many frequency peaks. This is in contradiction to pulsation
models which predict that stars with the same global parameters and rotation
rates should pulsate with the same set of frequencies.

In order to obtain an objective measure of similarity of frequency patterns in
a pair of stars, a correlation method was developed. Using non-adiabatic non-
rotating pulsation models, it is found that the frequencies of unstable modes
within a small effective temperature and luminosity range are indeed
correlated. The effect of rotation was simulated by applying rotational
splitting using a 3rd-order approximation formula (Reese et al., 2006). A
spread of rotation periods similar to that expected for mid-A stars was used.
The results show that the effect of rotation is to increase the value of the
correlation. In other words, rotation tends to produce greater overall
similarity between frequency patterns, even though individual frequencies may
differ.

When applied to a large sample of TESS $\delta$ Sct stars with similar
effective temperatures and luminosities, practically no correlation is found.
In other words, the observed frequencies are not related to the location of
the star in the H–R diagram at all.

The possibility that combination frequencies may be responsible for some of
the discrepancy seems unlikely. Although combination frequencies were detected
in about 60 percent of the stars, they constitute only a small fraction of the
total number of frequency peaks observed in a star.

The fact that frequencies in $\delta$ Sct stars with similar effective
temperatures and luminosities vary so widely indicates that there is no
possibility of finding a useful period - luminosity relation. Even for stars
where only one dominant mode is present, the frequency of this mode can vary
widely between stars with essentially the same physical parameters. This means
that the spherical harmonic degree and/or radial order must be different.
However, there is a general trend for the hotter, more luminous stars to have
shorter pulsation periods, in agreement with non-adiabatic pulsation models.

The effect described here suggests that there is a fundamental problem with
current models. The likely presence of starspots and even flares on early-type
stars (Balona, 2021) indicates that the outer stellar layers and atmospheres
of A and B stars have a complex structure. It is possible that local
conditions determine the frequencies that are selected.

Progress might be achieved by including the effect of surface convection, even
in the hottest models of $\delta$ Sct stars. An attempt should be made on
modifying the pulsational parameters in these models to determine if a wide
variety of unstable frequencies can thereby be obtained. The possible effect
of photospheric structures on the number of observed pulsation peaks could be
studied using high-dispersion spectroscopy. Observations of stars with a
single dominant mode and stars with very rich frequency spectra may provide
clues to the difference in mode selection and the extreme sensitivity of
pulsational driving to local conditions.

## Acknowledgments

I thank the National Research Foundation of South Africa for financial support
and Dr. W. Dziembowski for permission to use his code.

This paper includes data collected by the TESS mission. Funding for the TESS
mission is provided by the NASA Explorer Program. Funding for the TESS
Asteroseismic Science Operations Centre is provided by the Danish National
Research Foundation (Grant agreement no.: DNRF106), ESA PRODEX (PEA
4000119301) and Stellar Astrophysics Centre (SAC) at Aarhus University. We
thank the TESS and TASC/TASOC teams for their support of the present work.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium).
Funding for the DPAC has been provided by national institutions, in particular
the institutions participating in the Gaia Multilateral Agreement.

This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France. This research has made use of the VizieR catalogue access
tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original
description of the VizieR service was published in A&AS 143, 23.

The data presented in this paper were obtained from the Mikulski Archive for
Space Telescopes (MAST). STScI is operated by the Association of Universities
for Research in Astronomy, Inc., under NASA contract NAS5-2655.

## Data availability

The data underlying this article can be obtained from
https://sites.google.com/view/tessvariables/home.

## References

  * Asplund et al. [2009] M. Asplund, N. Grevesse, A. J. Sauval, and P. Scott.  The Chemical Composition of the Sun.  _ARA &A_, 47:481–522, September 2009.  10.1146/annurev.astro.46.060407.145222.
  * Balona [2018] L. A. Balona.  Gaia luminosities of pulsating A-F stars in the Kepler field.  _MNRAS_ , 479:183–191, September 2018.  10.1093/mnras/sty1511.
  * Balona [2021] L. A. Balona.  Spots and Flares in Hot Main Sequence Stars Observed by Kepler, K2, and TESS.  _Frontiers in Astronomy and Space Sciences_ , 8:32, 2021\.  10.3389/fspas.2021.580907.  URL https://www.frontiersin.org/article/10.3389/fspas.2021.580907.
  * Balona and Ozuyar [2020] L. A. Balona and D. Ozuyar.  Pulsation among TESS A and B stars and the Maia variables.  _MNRAS_ , 493(4):5871–5879, March 2020.  10.1093/mnras/staa670.
  * Balona [2022] Luis A. Balona.  Identification and classification of TESS variable stars.  _arXiv e-prints_ , art. arXiv:2212.10776, December 2022.
  * Balona [2023] Luis A. Balona.  Maia variables and other anomalies among pulsating stars.  _Frontiers in Astronomy and Space Sciences_ , in press:85, October 2023.  10.3389/fspas.2020.529025.
  * Barac et al. [2022] Natascha Barac, Timothy R. Bedding, Simon J. Murphy, and Daniel R. Hey.  Revisiting bright $\delta$ Scuti stars and their period-luminosity relation with TESS and Gaia DR3.  _MNRAS_ , 516(2):2080–2094, October 2022.  10.1093/mnras/stac2132.
  * Bedding et al. [2023] Timothy R. Bedding, Simon J. Murphy, Courtney Crawford, Daniel R. Hey, Daniel Huber, Hans Kjeldsen, Yaguang Li, Andrew W. Mann, Guillermo Torres, Timothy R. White, and George Zhou.  TESS Observations of the Pleiades Cluster: A Nursery for $\delta$ Scuti Stars.  _ApJ_ , 946(1):L10, March 2023.  10.3847/2041-8213/acc17a.
  * Bowman et al. [2016] D. M. Bowman, D. W. Kurtz, M. Breger, S. J. Murphy, and D. L. Holdsworth.  Amplitude modulation in $\delta$ Sct stars: statistics from an ensemble study of Kepler targets.  _MNRAS_ , 460:1970–1989, August 2016.  10.1093/mnras/stw1153.
  * Bowman et al. [2021] D. M. Bowman, J. Hermans, J. Daszyńska-Daszkiewicz, D. L. Holdsworth, A. Tkachenko, S. J. Murphy, B. Smalley, and D. W. Kurtz.  KIC 5950759: a high-amplitude $\delta$ Sct star with amplitude and frequency modulation near the terminal age main sequence.  _MNRAS_ , 504(3):4039–4053, July 2021.  10.1093/mnras/stab1124.
  * Cohen and Sarajedini [2012] Roger E. Cohen and Ata Sarajedini.  SX Phoenicis period-luminosity relations and the blue straggler connection.  _MNRAS_ , 419(1):342–357, January 2012.  10.1111/j.1365-2966.2011.19697.x.
  * Dziembowski [1977] W. Dziembowski.  Oscillations of giants and supergiants.  _Acta Astron._ , 27:95–126, 1977.
  * Dziembowski [1993] W. A. Dziembowski.  Mode selection and other nonlinear phenomena in stellar oscillations.  In W. W. Weiss and A. Baglin, editors, _IAU Colloq. 137: Inside the Stars_ , volume 40 of _Astronomical Society of the Pacific Conference Series_ , pages 521–534, January 1993.
  * Gaia Collaboration et al. [2016] Gaia Collaboration, T. Prusti, J. H. J. de Bruijne, A. G. A. Brown, A. Vallenari, C. Babusiaux, C. A. L. Bailer-Jones, U. Bastian, M. Biermann, D. W. Evans, and et al.  The Gaia mission.  _A &A_, 595:A1, November 2016.  10.1051/0004-6361/201629272.
  * Gaia Collaboration et al. [2021] Gaia Collaboration, A. G. A. Brown, A. Vallenari, T. Prusti, J. H. J. de Bruijne, C. Babusiaux, M. Biermann, O. L. Creevey, D. W. Evans, and et al.  Gaia Early Data Release 3. Summary of the contents and survey properties.  _A &A_, 649:A1, May 2021.  10.1051/0004-6361/202039657.
  * Garg et al. [2010] A. Garg, K. H. Cook, S. Nikolaev, M. E. Huber, A. Rest, A. C. Becker, P. Challis, A. Clocchiatti, G. Miknaitis, D. Minniti, L. Morelli, K. Olsen, J. L. Prieto, N. B. Suntzeff, D. L. Welch, and W. M. Wood-Vasey.  High-amplitude $\delta$-Scutis in the Large Magellanic Cloud.  _AJ_ , 140(2):328–338, August 2010.  10.1088/0004-6256/140/2/328.
  * Guzik et al. [1998] J. A. Guzik, M. R. Templeton, and P. A. Bradley.  The Effects of Metallicity on Delta Scuti Star Asteroseismology.  In Paul A. Bradley and Joyce A. Guzik, editors, _A Half Century of Stellar Pulsation Interpretation_ , volume 135 of _Astronomical Society of the Pacific Conference Series_ , page 470, January 1998.
  * McNamara [2011] D. H. McNamara.  Delta Scuti, SX Phoenicis, and RR Lyrae Stars in Galaxies and Globular Clusters.  _AJ_ , 142(4):110, October 2011.  10.1088/0004-6256/142/4/110.
  * Mourabit and Weinberg [2023] Mohammed Mourabit and Nevin N. Weinberg.  Resonant Mode Coupling in $\delta$ Scuti Stars.  _ApJ_ , 950(1):6, June 2023.  10.3847/1538-4357/acca16.
  * Paczyński [1970] B. Paczyński.  Evolution of Single Stars. I. Stellar Evolution from Main Sequence to White Dwarf or Carbon Ignition.  _Acta Astron._ , 20:47, 1970.
  * Poleski et al. [2010] R. Poleski, I. Soszyński, A. Udalski, M. K. Szymański, M. Kubiak, G. Pietrzyński, Ł. Wyrzykowski, O. Szewczyk, and K. Ulaczyk.  The Optical Gravitational Lensing Experiment. The OGLE-III Catalog of Variable Stars. VI. Delta Scuti Stars in the Large Magellanic Cloud.  _Acta Astron._ , 60(1):1–16, March 2010.  10.48550/arXiv.1004.0950.
  * Poretti et al. [2008] Ennio Poretti, Gisella Clementini, Enrico V. Held, Claudia Greco, Mario Mateo, Luca Dell’Arciprete, Luca Rizzi, Marco Gullieuszik, and Marcella Maio.  Variable Stars in the Fornax dSph Galaxy. II. Pulsating Stars below the Horizontal Branch.  _ApJ_ , 685(2):947–957, October 2008.  10.1086/591241.
  * Reese et al. [2006] D. Reese, F. Lignières, and M. Rieutord.  Acoustic oscillations of rapidly rotating polytropic stars. II. Effects of the Coriolis and centrifugal accelerations.  _A &A_, 455:621–637, August 2006.  10.1051/0004-6361:20065269.
  * Rogers and Iglesias [1992] F. J. Rogers and C. A. Iglesias.  Radiative atomic Rosseland mean opacity tables.  _ApJS_ , 79:507–568, April 1992.  10.1086/191659.
  * Samus et al. [2017] N. N. Samus, E. V. Kazarovets, O. V. Durlevich, N. N. Kireeva, and E. N. Pastukhova.  General catalogue of variable stars: Version GCVS 5.1.  _Astronomy Reports_ , 61(1):80–88, Jan 2017.  10.1134/S1063772917010085.
  * Smolec [2014] R. Smolec.  Mode selection in pulsating stars.  In J. A. Guzik, W. J. Chaplin, G. Handler, and A. Pigulski, editors, _Precision Asteroseismology_ , volume 301 of _IAU Symposium_ , pages 265–272, February 2014.  10.1017/S1743921313014439.
  * Ziaali et al. [2019] Elham Ziaali, Timothy R. Bedding, Simon J. Murphy, Timothy Van Reeth, and Daniel R. Hey.  The period-luminosity relation for $\delta$ Scuti stars using Gaia DR2 parallaxes.  _MNRAS_ , 486(3):4348–4353, July 2019.  10.1093/mnras/stz1110.


Department of Informatics, University of Bergen,
<EMAIL_ADDRESS>Department of Informatics, University of
Bergen,
NorwayFedor.Fomin@uib.nohttps://orcid.org/0000-0003-1955-4612Department of
Informatics, University of Bergen<EMAIL_ADDRESS>https://orcid.org/0000-0002-2619-2990 Department of Informatics, University of
Bergen<EMAIL_ADDRESS>Department of Informatics, University of
Bergen<EMAIL_ADDRESS>Sayan Bandyapadhyay, Fedor V. Fomin,
Petr A. Golovach, Nidhi Purohit, Kirill Simonov [100]Theory of computation
Design and analysis of algorithms Approximation algorithms analysis Facility
location and clustering The research leading to these results has received
funding from the Research Council of Norway via the project “‘MULTIVAL" (grant
no. 263317) and the European Research Council (ERC) via grant LOPPRE,
reference 819416.John Q. Open and Joan R. Access 2 42nd Conference on Very
Important Topics (CVIT 2016) CVIT 2016 CVIT 2016 December 24–27, 2016 Little
Whinging, United Kingdom 42 23

#  FPT Approximation for Fair Minimum-Load Clustering

Sayan Bandyapadhyay  Fedor V. Fomin  Petr A. Golovach  Nidhi Purohit  Kirill
Simonov

###### Abstract

In this paper, we consider the Minimum-Load $k$-Clustering/Facility Location
(MLkC) problem where we are given a set $P$ of $n$ points in a metric space
that we have to cluster and an integer $k>0$ that denotes the number of
clusters. Additionally, we are given a set $F$ of cluster centers in the same
metric space. The goal is to select a set $C\subseteq F$ of $k$ centers and
assign each point in $P$ to a center in $C$, such that the maximum _load_ over
all centers is minimized. Here the load of a center is the sum of the
distances between it and the points assigned to it.

Although clustering/facility location problems have a rich literature, the
minimum-load objective is not studied substantially, and hence MLkC has
remained a poorly understood problem. More interestingly, the problem is
notoriously hard even in some special cases including the one in line metrics
as shown by Ahmadian et al. [APPROX 2014, ACM Trans. Algorithms 2018]. They
also show APX-hardness of the problem in the plane. On the other hand, the
best known approximation factor for MLkC is $O(k)$, even in the plane.

In this work, we study a fair version of MLkC inspired by the work of
Chierichetti et al. [NeurIPS, 2017], which generalizes MLkC. Here the input
points are colored by one of the $\ell$ colors denoting the group they belong
to. MLkC is the special case with $\ell=1$. Considering this problem, we are
able to obtain a $3$-approximation in $f(k,\ell)\cdot n^{O(1)}$ time. Also,
our scheme leads to an improved $(1+\epsilon)$-approximation in case of
Euclidean norm, and in this case, the running time depends only polynomially
on the dimension $d$. Our results imply the same approximations for MLkC with
running time $f(k)\cdot n^{O(1)}$, achieving the first constant approximations
for this problem in general and Euclidean metric spaces. Our work makes
substantial progress in understanding the structure of the minimum-load
objective, advancing the frontiers of knowledge about the MLkC problem.

###### keywords:

fair clustering, load balancing, parameterized approximation, minimum-load,
$k$-median

###### category:

##  1 Introduction

_Clustering_ is the task of partitioning a set of data items into a number of
groups (or clusters) such that each group contains similar set of items.
Typically, the similarity of the clusters is modeled by a proxy objective
function, which one needs to optimize. Being a fundamental computational
problem in nature, clustering has a host of diverse applications in computer
science and other disciplines. Consequently, the problem has been studied with
several different and possibly independent objectives. Some of these became
notably popular, for example, $k$-means, $k$-median, and $k$-center [21, 18,
5, 25]. In this paper, we consider an objective which is not studied
substantially in the literature. In particular, we consider minimum-load
clustering. Here we are given a set $P$ of points in a metric space that we
have to cluster and an integer $k>0$ that denotes the number of clusters.
Additionally, we are given a set $F$ of cluster centers in the same metric
space. The goal is to select a set $C\subseteq F$ of $k$ centers and assign
each point in $P$ to a center in $C$, such that the maximum _load_ over all
centers is minimized. Here the load of a center is the sum of the distances
between it and the points assigned to it. That is, if $P^{\prime}$ is the set
of points assigned to a center $c$, then its load is $\sum_{p\in
P^{\prime}}d(c,p)$, where $d$ is the given metric. We formally refer to this
problem as Minimum-Load $k$-Clustering (MLkC). MLkC can be used to model
applications where the cost of serving the clients (or points) assigned to a
facility (or center) is incurred by the facility, e.g., assigning jobs to the
$k$ best servers from a pool of servers balancing their loads.

Surprisingly, MLkC is $\operatorname{{\sf NP}}$-hard even if the solution set
of centers $C$ is given, via a reduction from makespan-minimization [2]. In
fact, this assignment version of the problem can be shown to be
$\operatorname{{\sf NP}}$-hard even in line metrics and for $k=2$, via a
simple reduction from the Partition problem [28]. (In Partition, given a set
of integers, the goal is to partition it into two subsets such that the
difference between the sums of the integers in two subsets is minimized.)
Moreover, Ahmadian et al. [2] proved that the problem is strongly
$\operatorname{{\sf NP}}$-hard in line metrics (points on a line) and APX-hard
in the plane. On the positive side, an $O(k)$-approximation follows for this
problem from any existing $O(1)$-approximation for $k$-median [11, 5, 23, 30,
10]. This is true, as $k$-median minimizes the sum of the loads of the
centers. Also, constant-approximations are known for MLkC in some special
cases, e.g., in star metrics and line metrics. Beyond these special cases,
obtaining better than $O(k)$-approximation in polynomial time remained a
notoriously hard question, even in the plane. Indeed, as explicitly pointed
out by Ahmadian et al. [2], MLkC is resilient to attack by the standard
approximation techniques including LP rounding and local search, which has
been fairly successful in obtaining good approximation algorithms for other
clustering problems. Given these difficulties, we investigate whether it is
possible to obtain $O(1)$-approximation for MLkC if we allow time $f(k)\cdot
n^{O(1)}$ instead of only $n^{O(1)}$, for some function $f(.)$ independent of
the input size $n$. Indeed, we study a much more general _fair_ version of the
problem.

Fair clustering was introduced by Chierichetti et al. [13] with the goal of
removing inherent biases from the regular clustering models. In this setting,
we also have a sensitive or protected feature of the data points, e.g., gender
or race. The goal is to obtain a clustering where the fraction of points from
a traditionally underrepresented group (w.r.t. the protected feature) in every
cluster is approximately equal to the fraction of points from this group in
the whole dataset. For simplicity, they assumed that the protected feature can
take only two values and designed fair $k$-center and $k$-median clustering
algorithms in this setting. In particular, here one is given two sets of
points $R$ and $B$ of color red and blue, respectively, and a balance
parameter $t\in[0,1]$. The objective is to find a clustering such that in
every cluster $O$, the ratio between the number of red points and the number
of blue points is at least $t$ and at most $1/t$, i.e., $t\leq\frac{|O\cap
R|}{|O\cap B|}\leq 1/t$. Subsequently, Rösner and Schmidt [32] considered a
general model where the protected feature can take any number of values and
designed fair clustering algorithms for $k$-center objective. Later, Bercea et
al. [8] and Bera et al. [7] independently considered a fair clustering model
that generalizes the models in both [13] and [32]. In this model, we are given
a partition $\\{P_{1}$, $P_{2}$, …, $P_{\ell}\\}$ of the input point set $P$
and balance parameters $0\leq\beta_{i}\leq\alpha_{i}\leq 1$ for each group
$1\leq i\leq\ell$. Then a clustering is called $(\alpha,\beta)$-_fair_ if the
fraction of points from each group $i$ in every cluster is at least
$\beta_{i}$ and at most $\alpha_{i}$. In this paper, we study the
$(\alpha,\beta)$-Fair Minimum-Load $k$-Clustering (FMLkC) problem, where the
goal is to compute an $(\alpha,\beta)$-fair clustering that minimizes the
maximum load. (For a formal definition, please see Section 2.) We note that
the only clustering objectives considered in all the above mentioned works on
fair clustering are $k$-means, $k$-median and $k$-center. To the best of our
knowledge, fair clustering was not studied with the minimum-load objective
before our work.

###  1.1 Our Results and Techniques

Considering the FMLkC problem in general and Euclidean metric spaces we obtain
the following results.

######  Theorem 1.1 (Informal).

There is a $3$-approximation algorithm for $(\alpha,\beta)$-Fair Minimum-Load
$k$-Clustering in general metric spaces that runs in time
$2^{\tilde{O}(k\ell^{2})}n^{O(1)}$. For $d$-dimensional Euclidean spaces,
there is a $(1+\epsilon)$-approximation algorithm for $(\alpha,\beta)$-Fair
Minimum-Load $k$-Clustering with running time
$2^{\tilde{O}(k\ell^{2}/\epsilon^{O(1)})}n^{O(1)}$.

In the above theorem, the $\tilde{O}()$ notation hides logarithmic factors.
Note that all the running times are _fixed-parameter tractable_ (FPT) [15] in
$k$ and $\ell$ for constant $\epsilon$. Moreover, our results imply the same
approximations for Minimum-Load $k$-Clustering with running times FPT in only
$k$, achieving the first constant approximations for this problem in general
and Euclidean metric spaces. Note that in the Euclidean case, the running time
depends only polynomially on the dimension $d$. Recall that no better than
$O(k)$-approximation was known before even in the plane, and this version is
known to be APX-hard. Also, the reduction mentioned before from Partition
eliminates the existence of an exact algorithm in time $f(k)\cdot n^{O(1)}$,
unless P $\neq$ NP, as MLkC in line metrics is already $\operatorname{{\sf
NP}}$-hard when $k=2$. In this sense, our FPT $(1+\epsilon)$-approximation for
Euclidean spaces is tight and the best possible.

Our results are motivated by the recent FPT approximation results for
constrained clustering with popular $k$-median and $k$-means objectives [14,
6]. However, these results are based on coreset construction. A coreset is a
summary of the original dataset from which it is possible to retrieve a near-
optimal clustering. Their main contribution is to show that it is possible to
obtain coresets of size polynomial in $k,\log n$ and $d$. Alternatively, the
input can be compressed to an almost equivalent instance of size $poly(kd\log
n)$. Then one can enumerate all possible $k$-tuples of centers in FPT time
using the coreset and output the $k$-tuple having the minimum clustering cost.
This yields a $(1+\epsilon)$-approximation for Euclidean spaces and a slightly
larger $3$-approximation for general metric spaces due to some technical
reasons. However, such a small-sized coreset is not known for our problems.
Instead, we adapt approaches from [16, 22, 9] used for directly obtaining FPT
approximations for constrained $k$-median and $k$-means clustering. We note
that these schemes were known only in the special Euclidean case until
recently [22]. All these schemes produce in FPT (in $k$) time a list of
$k$-tuples of centers, such that at least one such $k$-tuple is a near-optimal
set of centers. Using the rather tedious similarity of the $k$-median and the
minimum-load objectives, we show these approaches can be adapted for our
problems as well. However, given such a $k$-tuple of centers, assigning the
points to the best centers or finding the optimal clustering, in our case is
still $\operatorname{{\sf NP}}$-hard. Nevertheless, we give a Mixed-Integer
Linear Programming (MILP) based $(1+\epsilon)$-approximation for this
assignment problem that runs in time FPT in $k$ and $\ell$ (in $k$ only for
MLkC). Our MILP is partly motivated by the fair $k$-median MILP [6]. However,
our MILP and its rounding are much more non-trivial compared to that for fair
$k$-median, especially due to the difference in the objectives. For example,
if we forget about the fairness constraints, in that case the assignment
algorithm for $k$-median is trivial. Assign each point to its closest center.
However, even in this case the assignment problem for MLkC is
$\operatorname{{\sf NP}}$-hard. Also, no near-optimal assignment scheme was
known in the literature (a $2$-approximate assignment scheme follows from the
generalized assignment problem (GAP) [33]). Thus, in this case we give a novel
$(1+\epsilon)$-approximate assignment scheme. In this case, we do not need
MILP – rounding of an LP is sufficient to obtain the desired assignment. All
these schemes applied together help us achieve the desired FPT approximations.

###  1.2 Related Work

Even et al. [17] and Arkin et al. [4] studied the MLkC problem under the name
_min-max star cover_ , where $F=P$. In this setting, MLkC can be viewed as a
weighted covering problem where the task is to cover the nodes of a graph by
stars. Both works obtain _bicriteria_ approximation for this problem where the
solution returned has near-optimal load, but uses more than $k$ centers.
Ahmadian et al. [2] studied several special cases of the MLkC problem (under
the name Minimum-Load $k$-Facility Location111MLkC can also be viewed as a
facility location problem with zero facility opening costs where we can still
open only $k$ facilities). They fully resolved the status of the MLkC problem
in line metrics. On the one hand, they designed a PTAS based on dynamic
programming. On the other hand, they proved that this version is strongly
$\operatorname{{\sf NP}}$-hard. They also designed a quasi-PTAS in tree
metrics. Moreover, they studied a variant of the problem with client demands
in star metrics.

The notion of fair clustering introduced by Chierichetti et al. [13] has been
studied extensively in the literature. For $k$-center objective, there are
several polynomial-time true $O(1)$-approximations [32, 8]. For $k$-median and
$k$-means objectives, polynomial-time $O(1)$-approximations are designed by
violating the fairness constraints by an additive factor [8, 7]. On the other
hand, it is possible to obtain true $O(1)$-approximations for these two
objectives if one is allowed to use $f(k,\ell)\cdot n^{O(1)}$ time [6].
Clustering problems have been studied under several other notions of fairness,
e.g., see [3, 12, 26, 27, 20, 31, 1].

#### Organization.

We define some useful notations and our problem formally in Section 2. Then we
describe the assignment algorithm for the FMLkC problem in Section 3. Finally,
in Section 4, we describe the full algorithms for FMLkC in details.

##  2 Preliminaries

We are given a set $P$ of points in a metric space
$(\mathcal{X},d(\cdot,\cdot))$, that we have to cluster. We are also given a
set $F$ of cluster centers in the same metric space. We note that $P$ and $F$
are not-necessarily disjoint, and in fact, $P$ may be equal to $F$. In the
Euclidean version of a clustering problem, $P\subseteq\mathbb{R}^{d}$,
$F=\mathbb{R}^{d}$ and $d(\cdot,\cdot)$ is the Euclidean distance.222Due to
the lack of better notations, we denote the dimension by $d$ and distance
function by $d(\cdot,\cdot)$ In the metric version, we assume that $F$ is
finite. Thus, strictly speaking, the Euclidean version is not a special case
of the metric version. In the metric version, we denote $|P\cup F|$ by $n$ and
in the Euclidean version, $|P|$ by $n$. For any integer $t\geq 1$, we denote
the set $\\{1,2,\ldots,t\\}$ by $[t]$.

For a partition $\mathbb{O}=\\{O_{1},\ldots,O_{k}\\}$ of $P$ and a set of $k$
cluster centers $C=\\{c_{1},\ldots,c_{k}\\}\subset F$, we say that
$\mathbb{O}$ is a _clustering_ of $P$ with the centers $c_{1}$, …, $c_{k}$. We
say that the _minimum-load cost_ of this clustering (also, simply _cost of
clustering_) is $\max_{i\in[k]}\mathrm{cost}_{c_{i}}(O_{i})$. Here
$\mathrm{cost}_{c_{i}}(O_{i})$ denotes the sum-of-distances cost of the
cluster $O_{i}$ with the center $c_{i}$, which is
$\mathrm{cost}_{c_{i}}(O_{i})=\sum_{x\in O_{i}}d(x,c_{i}).$ We use the
following notation to denote the cost of clustering w.r.t. the set of centers
$C$ up to a permutation of the clusters,

$\mathrm{cost}_{C}(\mathbb{O})=\min_{i_{1},\ldots,i_{k}}\max_{j\in[k]}\mathrm{cost}_{c_{j}}(O_{i_{j}}),$

where $i_{1}$, …, $i_{k}$ is a permutation of $[k]$. We also denote by
$\mathrm{cost}(O_{i})$ the cost of a cluster $O_{i}$ with the optimal choice
of a center, that is, $\mathrm{cost}(O_{i})=\min_{c\in
F}\mathrm{cost}_{c}(O_{i})$, and by $\mathrm{cost}(\mathbb{O})$ the optimal
cost of clustering $\mathbb{O}$,

$\mathrm{cost}(\mathbb{O})=\min_{\begin{subarray}{c}C\subset F\\\
|C|=k\end{subarray}}\mathrm{cost}_{C}(\mathbb{O}).$

Alternatively, a clustering with centers in $C\subset F$ can be defined as an
assignment $\varphi:P\to C$. The assignment $\varphi$ then corresponds to a
clustering $\\{\varphi^{-1}(c)\\}_{c\in C}$, and we say that the cost of the
assignment $\varphi$ is $\mathrm{cost}(\varphi)=\max_{c\in C}\sum_{x\in
P:\varphi(x)=c}d(x,c)$.

Now we define the main problem of our interest, where the goal is to find the
minimum-cost clustering that satisfies the fairness constraints.

######  Definition 2.1.

In the $(\alpha,\beta)$-Fair Minimum-Load $k$-Clustering (FMLkC) problem, we
are given a partition $\\{P_{1}$, $P_{2}$, …, $P_{\ell}\\}$ of $P$. We are
also given an integer $k>0$ and two _fairness vectors_
$\alpha,\beta\in[0,1]^{\ell}$, $\alpha=(\alpha_{1},\ldots,\alpha_{\ell})$,
$\beta=(\beta_{1},\ldots,\beta_{\ell})$. The objective is to select a set of
at most $k$ centers $C\subset F$ and an assignment $\varphi:P\to C$ such that
$\varphi$ satisfies the following _fairness constraints_ :

$\displaystyle\left|\\{x\in
P_{i}:\varphi(x)=c\\}\right|\leq\alpha_{i}\cdot\left|\\{x\in
P:\varphi(x)=c\\}\right|,\quad\forall c\in C,\forall i\in[\ell],$
$\displaystyle\left|\\{x\in
P_{i}:\varphi(x)=c\\}\right|\geq\beta_{i}\cdot\left|\\{x\in
P:\varphi(x)=c\\}\right|,\quad\forall c\in C,\forall i\in[\ell],$

and $\mathrm{cost}(\varphi)$ is minimized among all such assignments.

Minimum-Load $k$-Clustering (MLkC) is a restricted case of FMLkC with
$\ell=1$, and hence there is no fairness constraints involved in this case.
The $(\alpha,\beta)$-Fair $k$-median problem is defined identically except
there the cost is $\mathrm{cost}(\varphi)=\sum_{c\in C}\sum_{x\in
P:\varphi(x)=c}d(x,c).$

##  3 Assignment Problem for FMLkC

In the $(\alpha,\beta)$-fair assignment problem, we are additionally given a
set of centers $C\subset F$ and the goal is to find an assignment
$\varphi:P\to C$ such that $\varphi$ satisfies the fairness constraints and
$\mathrm{cost}(\varphi)=\max_{c\in C}\sum_{x\in P:\varphi(x)=c}d(x,c)$ is
minimized.

We refer to an assignment as a fair assignment if it satisfies the fairness
constraints. Also, we denote the optimal cost of $(\alpha,\beta)$-fair
assignment by OPT. In this section, for any $\epsilon>0$, we give a
$(1+\epsilon)$-approximation for this problem in $f(k,\ell,\epsilon)\cdot
n^{O(1)}$ time for some computable function $f$. In particular, we solve a
budgeted version of the problem where we are also given a budget $B$ and the
goal is to decide whether there is a fair assignment of cost at most $B$.

######  Lemma 3.1.

Suppose there is an algorithm $\mathcal{A}$ that given an instance of budgeted
$(\alpha,\beta)$-fair assignment and any $\epsilon>0$, in
$T(n,k,\ell,\epsilon)$ time, either returns a feasible assignment of cost at
most $(1+\epsilon)B$, or correctly detects that there is no feasible
assignment with budget $B$. Then for any $\epsilon>0$, one can obtain a
$(1+\epsilon)$-approximation for $(\alpha,\beta)$-fair assignment in
$(k\ell)^{O(k\ell)}n^{O(1)}+O_{\epsilon}(\log k)\cdot T(n,k,\ell,\epsilon/3)$
time333$O_{\epsilon}()$ notation hides $O(1/\epsilon)$ factor..

######  Proof 3.2.

The idea is to first find a range where OPT belongs and then apply
$\mathcal{A}$ with budget within this range to find a feasible assignment.
Given an instance $I$ of $(\alpha,\beta)$-fair assignment, first we use an
algorithm (Theorem 8.2, [6]) to compute a fair assignment of the points to the
centers of $C$ that minimizes the $(\alpha,\beta)$-fair $k$-median cost. This
algorithm runs in time $(k\ell)^{O(k\ell)}n^{O(1)}$. Let $D$ be the computed
$(\alpha,\beta)$-fair $k$-median cost returned by the algorithm. Then $D\leq
k\cdot$OPT, as the optimal cost of $(\alpha,\beta)$-fair assignment is at
least $1/k$ fraction of the optimal $(\alpha,\beta)$-fair $k$-median cost.
Also, $OPT\leq D$, as optimal $(\alpha,\beta)$-fair assignment cost is at most
the optimal $(\alpha,\beta)$-fair $k$-median cost. Hence $D/k\leq$ OPT $\leq
D$.

Let $\epsilon^{\prime}=\epsilon/3$ and $m$ be the maximum $t$ such that
$(1+\epsilon^{\prime})^{t}\leq D/k$. Also, let $M$ be the minimum $t$ such
that $D\leq(1+\epsilon^{\prime})^{t}$. Thus $(1+\epsilon^{\prime})^{m}\leq$
OPT $\leq(1+\epsilon^{\prime})^{M}$. We run the algorithm $\mathcal{A}$
setting $\epsilon$ to be $\epsilon^{\prime}$ for budget
$B=(1+\epsilon^{\prime})^{i}$ where $m\leq i\leq M$, and terminate it the
first time it returns a feasible assignment for some budget $B$. Let
$B^{\prime}$ be the budget for which this algorithm returns a feasible
assignment. Then $B^{\prime}\leq(1+\epsilon^{\prime})$ OPT, as any instance
with budget $B\geq$ OPT is a yes-instance, and for such a $B$, $\mathcal{A}$
returns a feasible assignment of cost at most $(1+\epsilon)B$. Hence, the cost
of the assignment returned by $\mathcal{A}$ with budget $B^{\prime}$ is at
most $(1+\epsilon^{\prime})^{2}$ OPT $\leq(1+\epsilon)$ OPT. As the algorithm
$\mathcal{A}$ can be used at most
$O_{\epsilon}(\log(M-m+1))=O_{\epsilon}(\log(D/(D/k)))=O_{\epsilon}(\log k)$
times, the whole algorithm runs in time
$(k\ell)^{O(k\ell)}n^{O(1)}+O_{\epsilon}(\log k)\cdot T(n,k,\ell,\epsilon/3)$.

In the following, we design an LP rounding based algorithm for budgeted
$(\alpha,\beta)$-fair assignment with the properties required in the above
lemma. Moreover, this algorithm runs in time
$k^{O(k\ell)}{\ell}^{O((k\ell^{2}/\epsilon)\log(\ell/\epsilon))}n^{O(1)}$.
Hence, we obtain the following theorem.

######  Theorem 3.3.

For any $\epsilon>0$, a $(1+\epsilon)$-approximation for $(\alpha,\beta)$-fair
assignment can be obtained in time
$k^{O(k\ell)}{\ell}^{O((k\ell^{2}/\epsilon)\log(\ell/\epsilon))}n^{O(1)}$.

Next, we design the algorithm for the budgeted version of
$(\alpha,\beta)$-fair assignment. Recall that in the budgeted version, we are
given an instance $I$ containing $\ell$ disjoint groups $\\{P_{i}\\}$ of
$P=\\{p_{1},\ldots,p_{n}\\}$, a set of $k$ centers
$C=\\{c_{1},\ldots,c_{k}\\}$ and the budget $B$. Our algorithm first rounds
each distance to a power of $(1+\epsilon)$. Fix any center $c_{i}\in C$. We
partition the points in $P$ into a number of classes based on their distances
from $c_{i}$. For all $p\in P$, let
$\hat{d}(p,c_{i})=(1+\epsilon)^{t}\epsilon^{2}B$ where
$t=\lceil\log_{1+\epsilon}d(p,c_{i})/(\epsilon^{2}B)\rceil$. Let
$d_{t}=(1+\epsilon)^{t}\epsilon^{2}B$. We refer to the points $p$ with
distance $\hat{d}(p,c_{i})=d_{t}$ as the _distance class_ $t$ with respect to
(w.r.t.) $c_{i}$, which is denoted by $S_{it}$.

For all $p\in P,c\in C$,
${d}(p,c)\leq\hat{d}(p,c)\leq(1+\epsilon)\cdot{d}(p,c)$.

Let $I^{\prime}$ be the new instance of budgeted $(\alpha,\beta)$-fair
assignment with the modified distance $\hat{d}$. As $\hat{d}$ is obtained by
scaling $d$ by at most a factor of $(1+\epsilon)$, we have the following
observation.

If there is a feasible assignment for $I$ with budget $B$, then there is a
feasible assignment for $I^{\prime}$ with budget $(1+\epsilon)B$. Also, if
there is a feasible assignment for $I^{\prime}$ with budget $(1+\epsilon)B$,
then there is a feasible assignment for $I$ with budget $(1+\epsilon)B$.

By the above observation, it is sufficient to consider $\hat{d}$ instead of
$d$ for the purpose of computing an assignment of cost at most
$(1+\epsilon)B$. Henceforth, by distance we mean $\hat{d}$.

Denote by $\varphi^{*}$ a feasible assignment for $I^{\prime}$ of cost at most
$(1+\epsilon)B$ (if any). We define a point $p\in P$ to be _costly_ w.r.t. a
center $c\in C$ if $\hat{d}(p,c)\geq\epsilon^{2}B$. Otherwise, we define the
point to be _cheap_ w.r.t. $c$. Note that the number of costly points that can
be assigned to each center in $\varphi^{*}$ is at most
$(1+\epsilon)/\epsilon^{2}\leq 2/\epsilon^{2}$. The next observation follows
from the fact that $d_{t}=(1+\epsilon)^{t}\epsilon^{2}B$.

For any point $p$ and center $c\in C$ with $\hat{d}(p,c)=d_{t}$, $t<0$ if $p$
is cheap w.r.t. $c$, and $t\geq 0$ if $p$ is costly w.r.t. $c$.

Now, as we are shooting for an assignment of cost at most $(1+\epsilon)B$, we
can discard all the distances $\hat{d}(p,c)$ larger than $(1+\epsilon)B$,
i.e., we can assume that such a $p$ will never be assigned to $c$. Without
loss of generality, we assume that all the distances we have are bounded by
$(1+\epsilon)B$. Let $\Delta$ be the maximum $t$ such that there are $p\in P$
and $c\in C$ with $\hat{d}(p,c)=d_{t}$ for a costly point $p$ w.r.t. $c$. By
our previous assumption, $\hat{d}(p,c)\leq(1+\epsilon)B$. Thus
$\Delta\leq\lceil\log_{1+\epsilon}((1+\epsilon)B/(\epsilon^{2}B))\rceil=O((1/\epsilon)\log(1/\epsilon))$,
which is a constant. For $1\leq i\leq k$, $0\leq t\leq\Delta$ and $1\leq
g\leq\ell$, let $z_{i,t,g}$ be the number of costly points $p\in P_{g}$
assigned to $c_{i}$ in $\varphi^{*}$ such that $\hat{d}(p,c_{i})=d_{t}$. Note
that each $z_{i,t,g}\leq 2/\epsilon^{2}$, as the total number of costly points
assigned to a center is at most $2/\epsilon^{2}$. Thus the total number of
distinct choices for these variables is
$(2/\epsilon^{2})^{k\ell\Delta}=(1/\epsilon)^{O((k\ell/\epsilon)\log(1/\epsilon))}.$

There are $(1/\epsilon)^{O((k\ell/\epsilon)\log(1/\epsilon))}$ distinct
choices for the variables $\\{z_{i,t,g}:1\leq i\leq k,0\leq t\leq\Delta,1\leq
g\leq\ell\\}$.

As we can probe all such possible choices, we assume that we know the exact
values of these variables in $\varphi^{*}$. Next, we describe a Mixed-Integer
Linear Program (MILP) for the budgeted version of the problem, which is partly
motivated by the $(\alpha,\beta)$-fair $k$-median MILP [6]. For every point
$p_{j}$ and center $c_{i}$, we have a fractional variable $x_{ij}$ denoting
the extent up to which $p_{j}$ is assigned to $c_{i}$. For every center
$c_{i}$ and group $g\in\\{1,\ldots,\ell\\}$, we have an integral variable
$y_{gi}$ denoting the “weight” of the points assigned to $c_{i}$ from group
$g$. The constraints of the MILP are described as follows. The Constraint 1
ensures that each point is assigned to the centers up to an extent of 1.
Constraint 2 ensures that the weight assigned from each group $g$ to each
center $c_{i}$ is exactly $y_{gi}$. Constraints 3 and 4 are fairness
constraints. Constraint 5 ensures that the weight of costly points from each
class $t$ and group $g$ assigned to each $c_{i}$ is exactly same as the
guessed value $z_{i,t,g}$. Constraint 6 ensures that the total load assigned
to each $c_{i}$ is bounded by $(1+\epsilon)B$. The first and the second
expressions on the left hand side of this constraint are corresponding to
costly and cheap points, respectively.

$\displaystyle\sum_{1\leq i\leq k}x_{ij}=1$ $\displaystyle\forall j\in[n]$ (1)
$\displaystyle\sum_{j\in[n]:p_{j}\in P_{g}}x_{ij}=y_{gi}$
$\displaystyle\forall i\in[k],\ \forall g\in[\ell],$ (2) $\displaystyle
y_{gi}\geq\beta_{g}\sum_{j\in[n]}x_{ij}$ $\displaystyle\forall i\in[k],\
\forall g\in[\ell]$ (3) $\displaystyle
y_{gi}\leq\alpha_{g}\sum_{j\in[n]}x_{ij}$ $\displaystyle\forall i\in[k],\
\forall g\in[\ell]$ (4) $\displaystyle\sum_{p_{j}\in P_{g}\cap
S_{it}}x_{ij}=z_{i,t,g}$ $\displaystyle\forall i\in[k],\ \forall
t\in\\{0,\ldots,\Delta\\},\ \forall g\in[\ell]$ (5)
$\displaystyle\sum_{g=1}^{\ell}\sum_{0\leq
t\leq\Delta}d_{t}z_{i,t,g}+\sum_{t<0}d_{t}\sum_{p_{j}\in
S_{it}}x_{ij}\leq(1+\epsilon)B$ $\displaystyle\forall i\in[k]$ (6)
$\displaystyle x_{ij}\geq 0$ $\displaystyle\forall j\in[n],\ \forall i\in[k],$
(7) $\displaystyle y_{gi}\in\mathbb{Z}_{\geq 0}$ $\displaystyle\forall
i\in[k],\ \forall g\in[\ell].$ (8)

Let us denote the above MILP by Fair-LP. A solution to Fair-LP is denoted by
$(x,y)$. We note that the assignment $\varphi^{*}$ induces a feasible solution
to Fair-LP. We use the following popular and celebrated result to solve this
MILP.

######  Proposition 3.4 ([29], [24], [19]).

An MILP with $K$ integral variables and encoding size $L$, can be solved in
time $K^{O(K)}L^{O(1)}$.

As Fair-LP has $k\ell$ integral variables $\\{y_{gi}\\}$ and polynomial
encoding size, it can be solved in $(k\ell)^{O(k\ell)}n^{O(1)}$ time. If this
algorithm outputs that there is no feasible solution to Fair-LP, we conclude
that $I$ is a no-instance. Otherwise, let $(x^{*},y^{*})$ denote the feasible
solution returned by this algorithm. Note that although the $y^{*}$ values are
integral, $x^{*}$ values can very well be fractional. Next, we show how to
round these variables to obtain an integral solution to Fair-LP such that the
load of every center is increased by an additive factor of $O(\epsilon\ell B)$
compared to its load in $(x^{*},y^{*})$.

Fix any group $g$. First, we show how to round the variables corresponding to
the points of $P_{g}$. For this purpose, we construct a network
$G_{N}=(V_{N},E_{N})$ with source $S$ and sink $T$ (see Figure 1). For each
point $p_{j}\in P_{g}$, there is a node $v_{j}$ in $V_{N}$. For each distance
class $t$ of every center $c_{i}$, there is a node $w_{it}$. Also, for each
center $c_{i}$, there is a node $u_{i}$. For each $v_{j}\in V_{N}$, there is
an arc $(S,v_{j})$ of capacity 1. For each $p_{j}\in P_{g}$ and center
$c_{i}$, there is an arc $(v_{j},w_{it})$ of capacity 1 where $t$ is the index
such that $p_{j}\in S_{it}$, i.e, $\hat{d}(p_{j},c_{i})=d_{t}$. For center
$c_{i}$ and distance class $t$, let $\lambda_{it}^{g}=\sum_{p_{j}\in P_{g}\cap
S_{it}}x_{ij}^{*}$, i.e, the weight assigned from $P_{g}\cap S_{it}$ to
$c_{i}$. For each node $w_{it}$, there is an arc $(w_{it},u_{i})$ of capacity
$\lceil\lambda_{it}^{g}\rceil$. Lastly, for each center $c_{i}$, there is an
arc $(u_{i},T)$ of capacity $y_{gi}^{*}$.

Note that for each center $c_{i}$, the number of distance classes is at most
the number of points $n$. Hence, the size of $V_{N}$ is a polynomial in $n$.
Also, note that the solution $(x^{*},y^{*})$ projected on the points of
$P_{g}$ induces a feasible fractional solution for the problem of computing a
flow of value $|P_{g}|$ in $G_{N}$.

The network $G_{N}$ has a fractional flow of value $|P_{g}|$.

As all the capacities of the arcs are integral, by integrality of flow, there
exists an integral feasible flow in $G_{N}$ of value $|P_{g}|$. We compute
such a flow $f$ by using any polynomial time flow computation algorithm. This
flow solution $f$ naturally gives us an integral assignment $\varphi_{f}$ of
the points in $P_{g}$ to the centers in $C$.

Figure 1: Figure showing the network $G_{N}$ constructed using the solution
$(x^{*},y^{*})$.

The number of points assigned to each center $c_{i}\in C$ via $\varphi_{f}$ is
exactly $y_{gi}^{*}$.

######  Proof 3.5.

Due to the capacity constraints of the arcs $\\{(u_{i},T)\\}$, the number of
points assigned to $c_{i}$ must be at most $y_{gi}^{*}$. Also, by definition,
$\sum_{i=1}^{k}y_{gi}^{*}=|P_{g}|$. As $f$ has value $|P_{g}|$, the capacity
of the arcs $\\{(u_{i},T)\\}$ must be saturated, which completes the proof.

Next, we analyse the load of $P_{g}$ assigned to each center via
$\varphi_{f}$.

######  Lemma 3.6.

For each center $c_{i}\in C$, $\sum_{p_{j}\in
P_{g}:\varphi_{f}(p_{j})=c_{i}}\hat{d}(p_{j},c_{i})\leq\sum_{t}d_{t}\lambda_{it}^{g}+O(\epsilon)B$.

######  Proof 3.7.

Consider the arcs $\\{(w_{it},u_{i})\\}$. The maximum flow corresponding to
these arcs is bounded by the sum of the capacities
$\sum_{t}\lceil\lambda_{it}^{g}\rceil$. Note that for every $0\leq
t\leq\Delta$, $\lambda_{it}^{g}=z_{i,t,g}$, which is an integer. Now,

$\displaystyle\sum_{t<0}d_{t}=\sum_{t<0}(1+\epsilon)^{t}\epsilon^{2}B<\epsilon^{2}B((1+\epsilon)/\epsilon)=O(\epsilon)B$
(9)

Hence,

$\displaystyle\sum_{t<0}d_{t}\lceil\lambda_{it}^{g}\rceil\leq\sum_{t<0}d_{t}(\lambda_{it}^{g}+1)\leq\sum_{t<0}d_{t}\lambda_{it}^{g}+\sum_{t<0}d_{t}{=}_{[\ref{eq:sumofcheapdistances}]}\sum_{t<0}d_{t}\lambda_{it}^{g}+O(\epsilon)B$

It follows that,

$\displaystyle\sum_{p_{j}\in
P_{g}:\varphi_{f}(p_{j})=c_{i}}\hat{d}(p_{j},c_{i})$
$\displaystyle\leq\sum_{t}d_{t}\lceil\lambda_{it}^{g}\rceil$
$\displaystyle=\sum_{t\geq
0}d_{t}\lambda_{it}^{g}+\sum_{t<0}d_{t}\lceil\lambda_{it}^{g}\rceil$
$\displaystyle=\sum_{t\geq
0}d_{t}\lambda_{it}^{g}+\sum_{t<0}d_{t}\lambda_{it}^{g}+O(\epsilon)B$
$\displaystyle=\sum_{t}d_{t}\lambda_{it}^{g}+O(\epsilon)B$

We repeat the above rounding process for all groups $g$. We combine the
assignment functions $\varphi_{f}$ corresponding to the $\ell$ disjoint groups
to obtain a single assignment for the points in $P$. For simplicity, we also
refer to this combined assignment as $\varphi_{f}$. By Observation 3,
$\varphi_{f}$ is feasible, as for each center $c_{i}$ and each group $g$, the
weight of the points in $P_{g}$ assigned to $c_{i}$ is exactly $y_{gi}^{*}$ as
in $(x^{*},y^{*})$. Next, we analyze the total load of each center.

######  Lemma 3.8.

For each center $c_{i}\in C$,
$\sum_{p_{j}:\varphi_{f}(p_{j})=c_{i}}\hat{d}(p_{j},c_{i})\leq(1+O(\epsilon\ell))B$

######  Proof 3.9.

$\displaystyle\sum_{p_{j}:\varphi_{f}(p_{j})=c_{i}}\hat{d}(p_{j},c_{i})=\sum_{g=1}^{\ell}\sum_{p_{j}\in
P_{g}:\varphi_{f}(p_{j})=c_{i}}\hat{d}(p_{j},c_{i})$
$\displaystyle=\sum_{g=1}^{\ell}\Big{(}\sum_{t}d_{t}\lambda_{it}^{g}+O(\epsilon)B\Big{)}$
$\displaystyle(\text{By Lemma }\ref{lem:singlegrpcost})$
$\displaystyle=\sum_{g=1}^{\ell}\Big{(}\sum_{t\geq
0}d_{t}\lambda_{it}^{g}+\sum_{t<0}d_{t}\lambda_{it}^{g}\Big{)}+O(\epsilon\ell)B$
$\displaystyle=\sum_{g=1}^{\ell}\Big{(}\sum_{t\geq
0}d_{t}z_{i,t,g}+\sum_{t<0}d_{t}(\sum_{p_{j}\in P_{g}\cap
S_{it}}x_{ij}^{*})\Big{)}+O(\epsilon\ell)B$
$\displaystyle=\bigg{(}\sum_{g=1}^{\ell}\sum_{t\geq
0}d_{t}z_{i,t,g}+\sum_{t<0}d_{t}\sum_{p_{j}\in
S_{it}}x_{ij}^{*}\bigg{)}+O(\epsilon\ell)B$
$\displaystyle\leq(1+\epsilon)B+O(\epsilon\ell)B$ $\displaystyle(\text{By
Constraint }\ref{eq:ilp6}\text{ of Fair-LP})$
$\displaystyle=(1+O(\epsilon\ell))B$

By scaling $\epsilon$ down by a factor of $\Omega(\ell)$, we obtain the
desired approximate bound. Thus, by Observation 3, it follows that the number
of distinct possible choices of the $Z_{i,t,g}$ values is
$(\ell/\epsilon)^{O((k\ell^{2}/\epsilon)\log(\ell/\epsilon))}$. For each such
choice, solving the MILP and rounding takes $(k\ell)^{O(k\ell)}n^{O(1)}$ time.
Thus, the algorithm for solving the budgeted version runs in time
$k^{O(k\ell)}{\ell}^{O((k\ell^{2}/\epsilon)\log(\ell/\epsilon))}n^{O(1)}$. The
following lemma completes the proof of Theorem 3.3.

######  Lemma 3.10.

The above MILP based algorithm for budgeted $(\alpha,\beta)$-fair assignment,
in time
$k^{O(k\ell)}{\ell}^{O((k\ell^{2}/\epsilon)\log(\ell/\epsilon))}n^{O(1)}$,
either returns a feasible assignment of budget at most $(1+\epsilon)B$, or
correctly detects that there is no feasible assignment of budget $B$.

##  4 Approximation Algorithms for FMLkC

In this section, we describe the FPT approximation algorithms for the FMLkC
problem, both in the general metric case and in the Euclidean case. In the
general metric case, we aim for a $(3+\epsilon)$-approximation, and in the
Euclidean case for a $(1+\epsilon)$ approximation, for a given $0<\epsilon<1$.
Essentially, we obtain these algorithms as a combination of our assignment
algorithm presented before, and known generic results for constrained
clustering problems that follow the framework of Ding and Xu [16]. Our general
metric algorithm employs the result of Goyal, Jaiswal and Kumar [22], and in
the Euclidean case we use the result of Bhattacharya, Jaiswal and Kumar [9].
Both of these provide algorithms that in FPT time produce a reasonably short
list of candidate sets of $k$ centers, such that for each possible clustering
of the input points one of the sets in the list provides the desired
approximation, with good probability. Note that the results mentioned above
are stated in fact for the $k$-median objective, and not the minimum-load
clustering that we study in this work. However, by tweaking the error
guarantees in the respective proofs we can show that these results hold in the
minimum-load setting as well. Next, we present these in detail.

We start with the Euclidean case and show the following analogue of Theorem 1
in [9] for the minimum-load objective.

######  Theorem 4.1.

Given a set of $n$ points $P\subset\mathbb{R}^{d}$, parameters $k$ and
$0<\epsilon<1$, there is a randomized algorithm that in time
$2^{\tilde{O}(k/\epsilon^{O(1)})}nd$ outputs a list $L$ of
$2^{\tilde{O}(k/\epsilon^{O(1)})}$ sets of centers of size $k$ such that for
any partition $\mathbb{P}^{*}=\\{P_{1}^{*}$, …, $P_{k}^{*}\\}$ of $P$ the
following event occurs with probability at least $1/2$: there is a set $C$ in
$L$ such that

$cost_{C}(\mathbb{P}^{*})\leq(1+\epsilon)\max_{i\in[k]}\mathrm{cost}(P_{i}^{*}).$

######  Proof 4.2.

The algorithm proceeds exactly as Algorithm 5.1 in [9]. For the analysis, we
observe that Bhattacharya et al. prove the following statement (follows
immediately from invariant $P(i)$ in [9]): With constant probability, there is
a set of centers $C=\\{c_{1},\ldots,c_{k}\\}$ in the output of the algorithm
and the permutation $i_{1}$, …, $i_{k}$ of the clusters in $\mathbb{P}^{*}$
such that for each $j\in[k]$,

$\mathrm{cost}_{c_{j}}(P_{i_{j}}^{*})\leq(1+\frac{\epsilon}{2})\cdot\mathrm{cost}(P_{i_{j}}^{*})+\frac{\epsilon}{2k}\cdot\sum_{i=1}^{k}\mathrm{cost}(P_{i}^{*}).$

From here it easily follows that the set of centers $C$ achieves
$(1+\epsilon)$-approximation of the cost of $\mathbb{P}^{*}$ with respect to
the minimum-load objective:

$\mathrm{cost}_{C}(\mathbb{P}^{*})\leq\max_{j\in[k]}\mathrm{cost}_{c_{j}}(P_{i_{j}}^{*})\leq(1+\frac{\epsilon}{2})\cdot\max_{j\in[k]}\mathrm{cost}(P_{i_{j}}^{*})+\frac{\epsilon}{2k}\cdot\sum_{i=1}^{k}\mathrm{cost}(P_{i}^{*})\leq(1+\epsilon)\cdot\max_{j\in[k]}\mathrm{cost}(P_{i_{j}}^{*}),$

since $\sum_{i=1}^{k}\mathrm{cost}(P_{i}^{*})\leq
k\max_{j\in[k]}\mathrm{cost}(P_{i_{j}}^{*})$.

In the general metric case, a similar result can be shown, however with the
approximation factor of $(3+\epsilon)$. Specifically, we show an analogue of
Theorem 5 in [22] for the minimum-load objective. Similarly to Theorem 4.1,
the algorithm and the analysis is identical to what is presented in [22], up
to a different view on the cost upper bound.

######  Theorem 4.3.

Given a set of $n$ points $P$ in a metric space, parameters $k$ and
$0<\epsilon<1$, there is a randomized algorithm that in time
$(k/\epsilon)^{O(k)}n$ outputs a list $L$ of $(k/\epsilon)^{O(k)}$ sets of
centers of size $k$ such that for any partition $\mathbb{P}^{*}=\\{P_{1}^{*}$,
…, $P_{k}^{*}\\}$ of $P$ the following event occurs with probability at least
$1/2$: there is a set $C$ in $L$ such that

$\mathrm{cost}_{C}(\mathbb{P}^{*})\leq(3+\epsilon)\max_{i\in[k]}\mathrm{cost}(P_{i}^{*}).$

######  Proof 4.4.

The algorithm proceeds exactly as Algorithm 1 in [22]. For the analysis, we
observe that Goyal et al. prove the following statement (encapsulated by
Property-I in [22]): With constant probability, there is a set of centers
$C=\\{c_{1},\ldots,c_{k}\\}$ in the output of the algorithm and the
permutation $i_{1}$, …, $i_{k}$ of the clusters in $\mathbb{P}^{*}$ such that
for each $j\in[k]$,

$\mathrm{cost}_{c_{j}}(P_{i_{j}}^{*})\leq(3+\frac{\epsilon}{2})\cdot\mathrm{cost}(P_{i_{j}}^{*})+\frac{\epsilon}{2k}\cdot\sum_{i=1}^{k}\mathrm{cost}(P_{i}^{*}).$

Now, analogously to the proof of Theorem 4.1, it follows that the set of
centers $C$ achieves $(3+\epsilon)$-approximation of the cost of
$\mathbb{P}^{*}$ with respect to the minimum-load objective:

$\mathrm{cost}_{C}(\mathbb{P}^{*})\leq\max_{j\in[k]}\mathrm{cost}_{c_{j}}(P_{i_{j}}^{*})\leq(3+\frac{\epsilon}{2})\cdot\max_{j\in[k]}\mathrm{cost}(P_{i_{j}}^{*})+\frac{\epsilon}{2k}\cdot\sum_{i=1}^{k}\mathrm{cost}(P_{i}^{*})\leq(3+\epsilon)\cdot\max_{j\in[k]}\mathrm{cost}(P_{i_{j}}^{*}),$

since $\sum_{i=1}^{k}\mathrm{cost}(P_{i}^{*})\leq
k\max_{j\in[k]}\mathrm{cost}(P_{i_{j}}^{*})$.

Now, Theorem 4.1 and Theorem 4.3 imply that for the Minimum-Load
$k$-Clustering problem with any given set of constraints on the desired
clustering, there exists a $(1+\epsilon)$-approximation algorithm in the
Euclidean case, and a $(3+\epsilon)$-approximation algorithm in the general
metric case. The running time is $2^{\tilde{O}(k/\epsilon^{O(1)})}(nd+T)$ for
both algorithms, where $T$ is the running time of an algorithm solving the
respective assignment problem, either exact or $(1+\epsilon)$-approximate. In
particular, combining the theorems with our approximation algorithm for
$(\alpha,\beta)$-fair assignment (Theorem 3.3), for the FMLkC problem we
obtain a $(1+\epsilon)$-approximation in $\mathbb{R}^{d}$ and a
$(3+\epsilon)$-approximation in general metric in FPT time when parameterized
by the number of clusters $k$ and the number of protected groups $\ell$.

######  Theorem 4.5.

For any $0<\epsilon<1$, there exists a randomized $(1+\epsilon)$-approximation
algorithm for $(\alpha,\beta)$-Fair Minimum-Load $k$-Clustering in
$\mathbb{R}^{d}$ with running time

$2^{\tilde{O}(k\ell^{2}/\epsilon^{O(1)})}n^{O(1)}d.$

The same holds in general metric with the approximation factor of
$(3+\epsilon)$, where the running time becomes

$2^{\tilde{O}(k\ell^{2}/\epsilon)}n^{O(1)}.$

######  Proof 4.6.

First, we deal with the Euclidean case. Fix an optimal fair min-load
$k$-clustering $\mathbb{P}^{*}=\\{P_{1}^{*}$, …, $P_{k}^{*}\\}$ of $P$. Run
the algorithm of Theorem 4.1 on $P$ with error parameter $\epsilon_{0}$ to
obtain the list $L$ of candidate sets of centers, here $\epsilon_{0}$ is such
that $(1+\epsilon_{0})^{2}\leq(1+\epsilon)$. In the following, assume that the
event described in the statement of Theorem 4.1 occurs for the clustering
$\mathbb{P}^{*}$, by constant number of repetitions the probability of this
can be lifted arbitrarily close to one. That is, there exists a set of $k$
centers $C^{\prime}$ in $L$ such that

$cost_{C^{\prime}}(\mathbb{P}^{*})\leq(1+\epsilon_{0})\max_{i\in[k]}\mathrm{cost}(P_{i}^{*}).$
(10)

Now, for each set of centers $C$ in $L$, run the
$(1+\epsilon_{0})$-approximate assignment algorithm given by Theorem 3.3 on
$(P,C)$, and choose the set of centers $C^{\prime\prime}$ that gives the best
assignment cost among the considered sets, denote the computed assignment from
$P$ to $C^{\prime\prime}$ by $\varphi$. In what follows, we show that the set
of centers $C^{\prime\prime}$ and the assignment $\varphi:P\to
C^{\prime\prime}$ provide $(1+\epsilon)$-approximate solution to the given
FMLkC instance. Denote by $\psi$ the assignment from $P$ to $C^{\prime}$ that
the algorithm outputs,

$\mathrm{cost}(\varphi)\leq\mathrm{cost}(\psi)\leq(1+\epsilon_{0})\mathrm{cost}_{C^{\prime}}(\mathbb{P}^{*})\leq(1+\epsilon_{0})^{2}\max_{i\in[k]}\mathrm{cost}(P_{i}^{*})\leq(1+\epsilon)\max_{i\in[k]}\mathrm{cost}(P_{i}^{*}),$

where the first inequality is by the choice of $C^{\prime\prime}$ and
$\varphi$, the second is by Theorem 3.3, and the third inequality is by (10).

Finally, we show that the running time bound holds. Invoking Theorem 4.1 takes
time $2^{\tilde{O}(k/\epsilon_{0}^{O(1)})}nd$, and produces a list of
$2^{\tilde{O}(k/\epsilon_{0}^{O(1)})}$ sets of centers. On each of them,
running the algorithm of Theorem 3.3 takes time
$k^{O(k\ell)}{\ell}^{O((k\ell^{2}/\epsilon_{0})\log(\ell/\epsilon_{0}))}n^{O(1)}d$.
Since $\epsilon_{0}=O(\epsilon)$, the total running time can be bounded as

$2^{\tilde{O}(k/\epsilon^{O(1)})}\left(nd+k^{O(k\ell)}{\ell}^{O((k\ell^{2}/\epsilon)\log(\ell/\epsilon))}n^{O(1)}d\right)=2^{\tilde{O}(k\ell^{2}/\epsilon^{O(1)})}n^{O(1)}d.$

The general metric case is identical, but to obtain the list of candidate sets
of centers we use Theorem 4.3 instead of Theorem 4.1. The final cost bound
changes to

$\mathrm{cost}(\varphi)\leq\mathrm{cost}(\psi)\leq(1+\epsilon_{0})\mathrm{cost}_{C^{\prime}}(\mathbb{P}^{*})\leq(3+\epsilon_{0})\cdot(1+\epsilon_{0})\max_{i\in[k]}\mathrm{cost}(P_{i}^{*})\leq(3+\epsilon)\max_{i\in[k]}\mathrm{cost}(P_{i}^{*}),$

where $\epsilon_{0}$ is chosen so that
$(3+\epsilon_{0})\cdot(1+\epsilon_{0})\leq(3+\epsilon)$.

## References

  * [1] Mohsen Abbasi, Aditya Bhaskara, and Suresh Venkatasubramanian.  Fair clustering via equitable group representations.  In Madeleine Clare Elish, William Isaac, and Richard S. Zemel, editors, FAccT ’21: 2021 ACM Conference on Fairness, Accountability, and Transparency, Virtual Event / Toronto, Canada, March 3-10, 2021, pages 504–514. ACM, 2021.
  * [2] Sara Ahmadian, Babak Behsaz, Zachary Friggstad, Amin Jorati, Mohammad R. Salavatipour, and Chaitanya Swamy.  Approximation algorithms for minimum-load _k_ -facility location.  ACM Trans. Algorithms, 14(2):16:1–16:29, 2018.
  * [3] Sara Ahmadian, Alessandro Epasto, Ravi Kumar, and Mohammad Mahdian.  Clustering without over-representation.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 267–275, 2019.
  * [4] Esther M. Arkin, Refael Hassin, and Asaf Levin.  Approximations for minimum and min-max vehicle routing problems.  J. Algorithms, 59(1):1–18, 2006.
  * [5] Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and Vinayaka Pandit.  Local search heuristics for k-median and facility location problems.  SIAM J. Comput., 33(3):544–562, 2004.
  * [6] Sayan Bandyapadhyay, Fedor V. Fomin, and Kirill Simonov.  On coresets for fair clustering in metric and euclidean spaces and their applications.  CoRR, abs/2007.10137, 2020.
  * [7] Suman Bera, Deeparnab Chakrabarty, Nicolas Flores, and Maryam Negahbani.  Fair algorithms for clustering.  In Advances in Neural Information Processing Systems, pages 4954–4965, 2019.
  * [8] Ioana O Bercea, Martin Groß, Samir Khuller, Aounon Kumar, Clemens Rösner, Daniel R Schmidt, and Melanie Schmidt.  On the cost of essentially fair clusterings.  In Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques (APPROX/RANDOM 2019). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2019.
  * [9] Anup Bhattacharya, Ragesh Jaiswal, and Amit Kumar.  Faster Algorithms for the Constrained k-means Problem.  Theory of Computing Systems, 62(1):93–115, 2018.  doi:10.1007/s00224-017-9820-7.
  * [10] Jaroslaw Byrka, Thomas W. Pensyl, Bartosz Rybicki, Aravind Srinivasan, and Khoa Trinh.  An improved approximation for _k_ -median and positive correlation in budgeted optimization.  ACM Trans. Algorithms, 13(2):23:1–23:31, 2017.
  * [11] Moses Charikar, Sudipto Guha, Éva Tardos, and David B. Shmoys.  A constant-factor approximation algorithm for the _k_ -median problem (extended abstract).  In Proceedings of the 31st Annual ACM Symposium on Theory of Computing, pages 1–10, 1999.
  * [12] Xingyu Chen, Brandon Fain, Liang Lyu, and Kamesh Munagala.  Proportionally fair clustering.  In International Conference on Machine Learning, pages 1032–1041, 2019.
  * [13] Flavio Chierichetti, Ravi Kumar, Silvio Lattanzi, and Sergei Vassilvitskii.  Fair clustering through fairlets.  In Advances in Neural Information Processing Systems, pages 5029–5037, 2017.
  * [14] Vincent Cohen-Addad and Jason Li.  On the fixed-parameter tractability of capacitated clustering.  In Christel Baier, Ioannis Chatzigiannakis, Paola Flocchini, and Stefano Leonardi, editors, 46th International Colloquium on Automata, Languages, and Programming, ICALP 2019, July 9-12, 2019, Patras, Greece, volume 132 of LIPIcs, pages 41:1–41:14. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2019.
  * [15] Marek Cygan, Fedor V Fomin, Łukasz Kowalik, Daniel Lokshtanov, Dániel Marx, Marcin Pilipczuk, Michał Pilipczuk, and Saket Saurabh.  Parameterized algorithms, volume 5.  Springer, 2015.
  * [16] Hu Ding and Jinhui Xu.  A unified framework for clustering constrained data without locality property.  Algorithmica, 82(4):808–852, 2020.
  * [17] Guy Even, Naveen Garg, Jochen Könemann, R. Ravi, and Amitabh Sinha.  Covering graphs using trees and stars.  In Sanjeev Arora, Klaus Jansen, José D. P. Rolim, and Amit Sahai, editors, Approximation, Randomization, and Combinatorial Optimization: Algorithms and Techniques, 6th International Workshop on Approximation Algorithms for Combinatorial Optimization Problems, APPROX 2003 and 7th International Workshop on Randomization and Approximation Techniques in Computer Science, RANDOM 2003, Princeton, NJ, USA, August 24-26, 2003, Proceedings, volume 2764 of Lecture Notes in Computer Science, pages 24–35. Springer, 2003.
  * [18] Tomás Feder and Daniel Greene.  Optimal algorithms for approximate clustering.  In Proceedings of the twentieth annual ACM symposium on Theory of computing, pages 434–444, 1988.
  * [19] András Frank and Éva Tardos.  An application of simultaneous diophantine approximation in combinatorial optimization.  Combinatorica, 7(1):49–65, March 1987.  doi:10.1007/bf02579200.
  * [20] Mehrdad Ghadiri, Samira Samadi, and Santosh S. Vempala.  Socially fair k-means clustering.  In Madeleine Clare Elish, William Isaac, and Richard S. Zemel, editors, FAccT ’21: 2021 ACM Conference on Fairness, Accountability, and Transparency, Virtual Event / Toronto, Canada, March 3-10, 2021, pages 438–448. ACM, 2021.
  * [21] Teofilo F Gonzalez.  Clustering to minimize the maximum intercluster distance.  Theoretical Computer Science, 38:293–306, 1985.
  * [22] Dishant Goyal, Ragesh Jaiswal, and Amit Kumar.  FPT Approximation for Constrained Metric k-Median/Means.  In Yixin Cao and Marcin Pilipczuk, editors, 15th International Symposium on Parameterized and Exact Computation (IPEC 2020), volume 180 of Leibniz International Proceedings in Informatics (LIPIcs), pages 14:1–14:19, Dagstuhl, Germany, 2020. Schloss Dagstuhl–Leibniz-Zentrum für Informatik.  URL: https://drops.dagstuhl.de/opus/volltexte/2020/13317, doi:10.4230/LIPIcs.IPEC.2020.14.
  * [23] Kamal Jain and Vijay V. Vazirani.  Approximation algorithms for metric facility location and _k_ -median problems using the primal-dual schema and lagrangian relaxation.  J. ACM, 48(2):274–296, 2001.
  * [24] Ravi Kannan.  Minkowski’s convex body theorem and integer programming.  Mathematics of Operations Research, 12(3):415–440, August 1987\.  doi:10.1287/moor.12.3.415.
  * [25] Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman, and Angela Y. Wu.  A local search approximation algorithm for k-means clustering.  Comput. Geom., 28(2-3):89–112, 2004.
  * [26] Matthäus Kleindessner, Pranjal Awasthi, and Jamie Morgenstern.  Fair k-center clustering for data summarization.  In 36th International Conference on Machine Learning, ICML 2019, pages 5984–6003. International Machine Learning Society (IMLS), 2019.
  * [27] Matthäus Kleindessner, Samira Samadi, Pranjal Awasthi, and Jamie Morgenstern.  Guarantees for spectral clustering with fairness constraints.  arXiv preprint arXiv:1901.08668, 2019.
  * [28] Richard E. Korf.  A complete anytime algorithm for number partitioning.  Artificial Intelligence, 106(2):181–203, 1998.
  * [29] H. W. Lenstra.  Integer programming with a fixed number of variables.  Mathematics of Operations Research, 8(4):538–548, November 1983\.  doi:10.1287/moor.8.4.538.
  * [30] Shi Li and Ola Svensson.  Approximating k-median via pseudo-approximation.  SIAM J. Comput., 45(2):530–547, 2016.
  * [31] Yury Makarychev and Ali Vakilian.  Approximation algorithms for socially fair clustering.  CoRR, abs/2103.02512, 2021.
  * [32] Clemens Rösner and Melanie Schmidt.  Privacy preserving clustering with constraints.  In 45th International Colloquium on Automata, Languages, and Programming (ICALP 2018). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2018\.
  * [33] David B. Shmoys and Éva Tardos.  An approximation algorithm for the generalized assignment problem.  Math. Program., 62:461–474, 1993.  doi:10.1007/BF01585178.