\subsection{Energy gap, critical point and scaling with system size}

\begin{figure}[h]
\includegraphics[trim=0cm 3cm 0cm 3cm,clip=true,width=0.5\textwidth]{supp1.pdf}
\caption{\label{suppfig1} (a) Energy gap $\Delta$ between ground and first excited state for a 10-spin chain. The energy gap decreases as $\Delta=1-g/g_c(p)+O(g^2)$, where $g_c(p)$ is the critical interaction strength separating the paramagnetic and ferromagnetic phases. In the finite-sized system, we define $g_c(p)$ to be the point where $\partial^2\Delta/\partial g^2$, shown in (b), is a maximum. Inset in (a) shows the dependence of $\Delta$ on $N$ for $g=0.4\omega$.}
\end{figure}
The energy gap $\Delta$ between the ground and first excited state is shown in Fig.~\ref{suppfig1}(a). The energy gap decreases approximately linearly with increasing $g/g_c(p)$. In an infinite system with nearest neighbour interactions, the exact energy gap is $\Delta=|\omega-g|$~\cite{pfeuty1970} and hence $\partial^2\Delta/\partial g^2=\delta(\omega-g)$. Finite-size effects regularize the divergence of $\partial^2\Delta/\partial g^2$, however we still observe a clear peak at a critical value of $g$, see Fig.~\ref{suppfig1}(b). We define $g_c(p)$ to be the value of $g$ corresponding to this peak. The critical value is close to $\omega$ for $p=\infty$ and decreases for decreasing $p$. The critical value will depend on $\omega$. In the main text, we fix $g_c(p)$ by defining this to be the critical value at $\omega=\omega_0$.

The dependence of work and efficiency on chain size $N$ is shown in Fig.~\ref{suppfig2}(a) for $g=0.4\omega_0$. For $p=\infty,3$ the performance scales extensively for $N\gtrsim 6$. Finite size effects are slightly more pronounced for $p=2$, however the scaling is close to extensive. For $p=1$ the performance scales super-extensively with $g$. This is due to a dependence of $\Delta$ on $N$, see inset to Fig.~\ref{suppfig1}(a). The critical point $g_c$ also depends on $N$ for $p=1$. In the cold, perturbative limit (see next section), $\Delta\sim \omega-g(\gamma+\ln N)$ with $\gamma$ the Euler-Mascheroni constant. The two-level approximation [Eq.~4 in the main text] then gives $W\sim N^{1+\beta_H g}$. This approximation captures the dependence of work on particle number for small $g/g_c$, see Fig.~\ref{suppfig2}(a,ii),(b,ii). Increasing either $g$ or $N$ increases $g/g_c$ and the perturbative approximation $W\sim N^{1+\beta_H g}$ breaks down. Scaling $g$ by $g_c$ effectively Kac renormalizes the interactions, rendering the long-range system extensive~\cite{kac1963,kac1969,cannas1996,gonzalezlazo2021,zhang2017}.

\begin{figure}[H]
\includegraphics[trim=0cm 8cm 0cm 8cm,clip=true,width=\textwidth]{supp2.pdf}
\caption{\label{suppfig2} The (a,i) work output and (a,ii) efficiency scale approximately extensively for $p=2,3,\infty$ whereas they scale non-extensively for $p=1$. The non-extensive scaling for $p=1$ agrees with the perturbative prediction $W\sim N^{1+\beta_H g}$ (black dashed lines). Results for $g=0.4\omega_0$. (b) The non-extensive scaling of $W$ and $\eta$ for $p=1$ and $g/\omega_0=0.2,0.3,0.4,0.5$ along with the perturbative prediction $W\sim N^{1+\beta_H g}$ (matching dashed black lines). The gray lines in (a,ii),(b,ii) are the Carnot efficiency $\eta_\mathrm{C}$. All results are for $\beta_H=10\omega_0^{-1}$, $\beta_C=2\beta_H$ and $r=r_\mathrm{NI}^\mathrm{max}$.}
\end{figure}

\subsection{Performance of the transverse Ising model}
The dimensionless free energy (free energy divided by $k_BT$) of the transverse Ising model confined to the first $N$ excited states is
\begin{equation}
\begin{split}
\ln Z&=\sum_{k=0}^{N-1} e^{-\beta\sqrt{\omega^2+g^2-2\omega g\cos(2\pi k/N)}},\\
&\approx\frac{N}{\pi}\int_0^\pi e^{-\beta\sqrt{\omega^2+g^2-2\omega g\cos\theta}}\,d\theta.
\end{split}
\end{equation}
Here and below we work in a frame where the ground state energy is zero. In the adiabatic limit, the system energy at the four points in the cycle in Fig.~1(a) of the main text are
\begin{equation}
\begin{split}
(1)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}e^{-\beta_H\sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}}\,d\theta,\\
(2)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}e^{-\beta_H\sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}}\,d\theta,\\
(3)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}e^{-\beta_C\sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}}\,d\theta,\\
(4)&\hspace{0.5cm}\langle E\rangle=\frac{N}{\pi}\int_0^\pi \sqrt{r^2\omega_0^2+g^2-2r\omega_0 g\cos\theta}e^{-\beta_C\sqrt{\omega_0^2+g^2-2\omega_0 g\cos\theta}}\,d\theta.
\end{split}
\end{equation}
The heat and work outputs can be computed from these, see Fig.~1(b) in the main text.

\subsection{Perturbative treatment}

We assume a large chain so that boundary effects can be ignored and hence we approximate the system as being translationally invariant (equivalently we can impose periodic boundary conditions). The spin operators can be converted to expressions in terms of bosonic operators via a Holstein-Primakoff transformation~\cite{holstein1940},
\begin{equation}
\hat{\sigma}_x^{(i)}\rightarrow \frac{\left(\sqrt{1-\hat{a}_i^\dagger\hat{a}_i}\right)\hat{a}_i+\hat{a}_i^\dagger\left(\sqrt{1-\hat{a}_i^\dagger\hat{a}_i}\right)}{2},\hspace{1cm}\hat{\sigma}_z^{(i)}\rightarrow\frac{1}{2}-\hat{a}_i^\dagger\hat{a}_i.
\end{equation}
In the low excitation regime we expand the Hamiltonian [Eq.~1 in the main text] to quadratic order in the bosonic operators $\hat{a}_i$, $\hat{a}_i^\dagger$. This gives~\cite{holstein1940,vogl2020,konig2021}
\begin{equation}
\hat{H}=-\frac{\omega N}{2}+\omega\sum_{i=1}^N\hat{a}_i^\dagger\hat{a}_i-\frac{g}{4}\sum_{\substack{i,j=1\\(j\ne i)}}^NJ_{ij}\left(\hat{a}_i+\hat{a}_i^\dagger\right)\left(\hat{a}_j+\hat{a}_j^\dagger\right).
\end{equation}
Fourier transforming and carrying out a Bogoliubov transformation gives~\cite{cevolani2016}
\begin{equation}
\hat{H}=\epsilon_0+\sum_{k=0}^{N-1} \left(\omega\operatorname{sign}\left(\omega-g\tilde{J}_k\right)\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right)b_k^\dagger b_k,
\end{equation}
with $\epsilon_0$ the ground-state energy, $b_k$ bosonic operators for the quasiparticle modes and
\begin{equation}\label{Jk}
\tilde{J}_k=\sum_{m=1}^N \frac{\cos\left(2\pi k m/N\right)}{m^p}=C_p^N\left(\frac{2\pi k}{N}\right),
\end{equation}
where $C_p^N(\theta)=\sum_{m=1}^N\cos(m\theta)/m^p$ is the finite-$N$ generalized Clausen function~\cite{wu2010}. For nearest neighbour interactions, we have $\tilde{J}_k=\cos (2\pi k/N)$.

From hereon we assume $\operatorname{sign}(\omega-g\tilde{J}_k)>0$. The thermodynamic properties of the spin chain can be determined from the partition function,
\begin{equation}
Z=\prod_{k=0}^{N-1}\frac{1}{1-\exp\left(-\beta\omega\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right)},
\end{equation}
with logarithm (dimensionless free energy),
\begin{equation}
\begin{split}
\ln Z&=\sum_{k=0}^{N-1} \ln \left[\frac{1}{1-\exp\left(-\beta\omega\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right)}\right],\\
&\approx\sum_{k=0}^{N-1}\exp\left(-\beta\omega\sqrt{1-\frac{2g}{\omega}\tilde{J}_k}\right).
\end{split}
\end{equation}
The latter approximation assumes low temperature. We have ignored the ground state energy, which adds an inconsequential constant to $\ln Z$. To proceed analytically, we presume small $g\tilde{J}_k/\omega$ and approximate $\sqrt{1-2g\tilde{J}_k/\omega}\approx 1-g\tilde{J}_k\omega$. Hence
\begin{equation}\label{Zapprox}
\ln Z\approx e^{-\beta\omega}\sum_{k=0}^{N-1}\exp\left(m\beta g\tilde{J}_k\right).
\end{equation}

We reserve a discussion of $p=1$ for later and for now assume $p>1$. For large $N$, $C_p^N(\theta)\rightarrow C_p(\theta)=\sum_{m=1}^\infty \cos(m\theta)/m^p$ and we can replace the sum over $k$ in Eq.~\eqref{Zapprox} by an integral
\begin{equation}
\sum_{k=0}^{N-1}\exp\left(\beta g\tilde{J}_k\right)\rightarrow \frac{N}{\pi}\int_0^\pi \exp\left(\beta g C_p(\theta)\right)\,d\theta.
\end{equation}
For even integers $p$,
\begin{equation}
C_p(\theta)=-\frac{(-1)^{p/2}(2\pi)^p}{2p!}B_p\left(\frac{\theta}{2\pi}\right),
\end{equation}
with $B_p(x)$ the Bernoulli polynomials, which are polynomials of order $p$. We obtain the following results for $p=\infty$ and $p=2$,
\begin{equation}\label{p12}
\begin{split}
\frac{1}{\pi}\int_0^\pi\exp\left(\beta g C_\infty(\theta)\right)\,d\theta&=\frac{1}{\pi}\int_0^\pi\exp\left(\beta g \cos(\theta)\right)\,d\theta=I_0(\beta g)\overset{\beta\rightarrow\infty}\sim \frac{e^{\beta g}}{\sqrt{2\pi \beta g}},\\
\frac{1}{\pi}\int_0^1\exp\left(\beta g C_2(\theta)\right)\,d\theta&=\frac{2 e^{\beta g\pi^2/6}}{\sqrt{\beta g\pi^2}}D\left(\frac{\sqrt{\beta g\pi^2}}{2}\right)\overset{\beta\rightarrow\infty}\sim \frac{e^{\beta g\pi^2/6}}{3\beta g\pi^2/6},
\end{split}
\end{equation}
where $I_0(x)$ is the modified Bessel function and $D(x)$ is Dawson's function (expressible in terms of the error function $\operatorname{erf}(x)$ via $D(x)=-(\sqrt{\pi}/2)e^{-x^2}i\operatorname{erf}(ix)$).

In general, for $p>3$ the integral $\int_0^1 \exp((m\beta gC_p(2\pi x)))\,dx$ can be approximated at low temperatures using the method of steepest descent,
\begin{equation}
\frac{1}{\pi}\int_0^\pi \exp(\beta gC_p(\theta))\,d\theta=\sqrt{\frac{1}{2\pi \beta gC_{p-2}(0)}}\exp(\beta gC_p(0)),
\end{equation}
where $C_p(0)=\sum_{m=1}^\infty m^{-p}=\zeta(p)$. This does not work for $p=2$ since $C_2'(\theta)|_{\theta=0}\ne 0$. This reflects that the spectrum is linear rather than quadratic around the lowest energy state. Nor does it work for $p=3$, since $C_3^{''}(\theta)|_{\theta=0}$ diverges. We observe numerically that $\exp(\beta gC_3(\theta))$ is dominated by its small $\theta$ behaviour for large $\beta g$. Hence for $p=3$ we expand the Clausen function in a power series around the maximum $\theta=0$,
\begin{equation}\label{p3}
\begin{split}
\frac{1}{\pi}\int_0^\pi \exp(\beta gC_3(\theta))\,d\theta &\approx \frac{\exp(\beta g\zeta(3))}{\pi}\int_0^\pi\exp\left(-\frac{\beta g \theta^2}{4}\left(3-\ln\theta^2\right)\right)\,d\theta,\\
&=\frac{\exp(\beta g\zeta(3))}{\pi}\int_0^\pi\exp\left(-\frac{\beta g \theta^2}{4}\left(3+\ln(\beta g)-\ln(\beta g\theta^2)\right)\right)\,d\theta,\\
&=\frac{\exp(\beta g\zeta(3))}{\pi\sqrt{\beta g}}\int_0^{\pi\sqrt{\beta g}}u^{u^2/2}\exp\left(-\frac{1}{4}u^2\left(\ln(\beta g)+3\right)\right)\,du,\\
&\overset{\beta\rightarrow\infty}\sim \frac{\exp(\beta g\zeta(3))}{\pi\sqrt{\beta g}}\int_0^{\pi\sqrt{\beta g}}\exp\left(-\frac{1}{4}u^2(\ln(\beta g)+3)\right)\,du,\\
&=\frac{\exp(\beta g\zeta(3))}{\sqrt{\pi \beta g(3+\ln\beta g)}}\operatorname{erf}\left(\frac{\pi}{2}\sqrt{3\beta g+\beta g\ln\beta g}\right),\\
&\overset{\beta\rightarrow\infty}\sim \frac{\exp(\beta g\zeta(3))}{\sqrt{\pi \beta g(3+\ln\beta g})}.
\end{split}
\end{equation}
Hence we find a logarithmic correction $\ln \beta g$ to the partition function for $p=3$, which interpolates between $p=2$ and $p>3$. Note that $3>2\ln \pi\approx 2.29$ and so the expansion around $\theta=0$ creates no convergence problems for large $\theta>1$. The asymptotic behaviour in the fourth line follows by replacing $u^{u^2/2}$ by $\lim_{u\rightarrow 0^+}u^{u^2/2}=1$ due to $\exp(-\frac{1}{4}u^2(\ln(\beta g)+3))$ being sharply peaked around the origin for large $\beta g$. The validity of this approximation was confirmed numerically. Equations~\eqref{p12} and \eqref{p3} give $\ln Z=N\mathcal{G}_p(\beta g)e^{-\beta \Delta}$ with $\Delta=\omega-g\zeta(p)$ and $\mathcal{G}_p(\beta g)$ arising from thermal fluctuations and dependent on $p$.

For $p=1$, Eq.~\eqref{Jk} at $k=0$ diverges with $N$ as $\tilde{J}_0\approx \ln N+\gamma$, with $\gamma$ the Euler-Mascheroni constant. We separate out this term in Eq.~\eqref{Zapprox}. For $k\ne 0$, Eq.~\eqref{Jk} converges and we can use $\sum_{m=0}^\infty \cos(m\theta)/m=C_1(\theta)=-\ln |2\sin(\theta/2)|$. Hence
\begin{equation}
\ln Z\approx Ne^{-\beta\omega}\left(e^{\beta g (\ln N+\gamma)}+2^{1-\beta g}N^{-1}\sum_{k=1}^{(N-1)/2}\left(\sin \frac{\pi k}{N}\right)^{-\beta g}\right),
\end{equation}
where we have assumed $N$ is odd for simplicity (assuming even $N$ will give the same final result below). Noting that the sum is dominated by small $k$ terms, we can use the small-angle approximation $\sin(\theta)\approx\theta$,
\begin{equation}
\begin{split}
\ln Z&\approx Ne^{-\beta\omega}\left(e^{\beta g(\ln N+\gamma)}+\frac{1}{\pi}\left(\frac{N}{2\pi}\right)^{\beta g-1}\sum_{k=1}^{(N-1)/2}k^{-\beta g}\right),\\
&\approx Ne^{-\beta(\omega-g(\ln N+\gamma))}\left(1+\frac{2}{N}\left(\frac{1}{2\pi e^\gamma}\right)^{\beta g}\zeta\left(\beta g\right)\right),\\
&\approx Ne^{-\beta(\omega-g(\ln N+\gamma))},\\
&=e^{-\beta \omega}e^{\beta g\gamma}N^{1+\beta g}.
\end{split}
\end{equation}
Hence for $p=1$ fluctuations are suppressed for low temperatures and large $N$, and the dimensionless free energy is $\ln Z\approx Ne^{-\beta\Delta}$ with $\Delta\approx \omega-g(\ln N+\gamma)$. The low-temperature work output within this approximation is then (see Fig.~\ref{suppfig2}),
\begin{equation}\label{Wapprox}
W\approx N\omega_0(r-1)\left(e^{-\beta_H\Delta}-e^{-\beta_C\Delta}\right)\approx \omega_0(r-1)e^{-\beta_H\omega_0}e^{\beta_Hg\gamma}N^{1+\beta_Hg}.
\end{equation}

Summarising, to lowest order in $g/\omega$ and for large $N$ and $\beta$, we obtain the dimensionless free energies,
\begin{equation}
\ln Z\approx \left\{\begin{array}{ll}\sqrt{\frac{1}{2\pi\beta g\zeta(p-2)}}e^{-\beta\Delta},&\hspace{0.5cm}p>3,\\
\sqrt{\frac{1}{\pi\beta g(3+\ln\beta g)}}e^{-\beta\Delta},&\hspace{0.5cm}p=3,\\
\frac{1}{3\beta g\zeta(2)}e^{-\beta\Delta},&\hspace{0.5cm}p=2,\\
e^{-\beta\Delta},&\hspace{0.5cm}p=1,\end{array}\right.
\end{equation}
with $\zeta(s)=\sum_{m=1}^\infty 1/m^s$ the Riemann-zeta function and $\Delta(\omega)=\omega(1-\omega_0 g/\omega g_c)$.

\subsection{Calculation of entanglement}
Within the quadratic approximation in the previous section, the first excited state is a state with one excitation uniformly spread across all spins,
\begin{equation}
\ket{1}=\frac{1}{\sqrt{N}}\sum_{i=1}^N\hat{\sigma}_+^{(i)}\ket{0}.
\end{equation}
We consider a partition dividing the chain in half. The entanglement entropy with respect to this partition and state $\ket{1}$ is
\begin{equation}
\mathcal{S}=-\operatorname{Tr}_R(\rho_L\ln\rho_L),
\end{equation}
with $\rho_L=\operatorname{Tr}_L(\ket{1}\bra{1})$ and $\operatorname{Tr}_{R(L)}$ a partial trace over states in the right(left) half of the chain. It is straightforward to show that this gives $\mathcal{S}=\ln (N/2)$.

For low temperatures, we can qualitatively describe the system by the thermal state
\begin{equation}
\rho_T^{01}=\frac{1}{1+e^{-\beta\Delta}}\left(\ket{0}\bra{0}+e^{-\beta\Delta}\ket{1}\bra{1}\right),
\end{equation}
where $\Delta=\omega(1-g/g_c)$ within the quadratic approximation. The Peres-Horodecki criterion states that a separable density matrix has positive partial trace~\cite{peres1996,horodecki1996}. Hence, if the density matrix does not have positive partial trace, the state is entangled (``non-PPT entanglement''). Due to the symmetry of the partition and the state, a sufficient condition for $\rho_T^{01}$ to exhibit non-PPT entanglement is $\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})<0$, with $\hat{M}$ any Hermitian operator acting on either the left ($\hat{M}\otimes$) or right ($\otimes \hat{M}$) partition of the chain~\cite{toth2009}. We choose $\hat{M}\otimes\hat{M}=\prod_{i=1}^N\hat{m}_i$ with
\begin{equation}
\hat{m}_i=\frac{1}{(1+\alpha)^{1/N}}\left[\frac{\hat{\mathcal{I}}}{2}-\hat{\sigma}_z^{(i)}-\alpha\left(\frac{\hat{\mathcal{I}}}{2}+\hat{\sigma}_z^{(i)}\right)\right],
\end{equation}
with $\hat{\mathcal{I}}$ the identity operator and $\alpha>0$ a parameter to be chosen. This gives
\begin{equation}
\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})=\frac{1-\alpha e^{-\beta\Delta}}{(1+\alpha)(1+e^{-\beta\Delta})}.
\end{equation}
We can make $\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})$ arbitrarily close to $-(1+e^{\beta\Delta})^{-1}$ by choosing $\alpha$ to be large ($\alpha\gg e^{\beta\Delta}$). This choice of $\hat{M}$ gives negative $\operatorname{Tr}(\hat{M}\otimes\hat{M}\rho_T^{01})$ and hence the state exhibits non-PPT entanglement.

\subsection{Calculation of counterdiabatic drive}
Here we derive Eq.~6 in the main text. We find the coefficients $C_{ij}$ in $\hat{H}_\mathrm{cd}=\sum_{ij(j\ne i)}C_{ij}\hat{\sigma}_x^{(i)}\hat{\sigma}_y^{(j)}$ by minimising the action
\begin{equation}
S=\operatorname{Tr}(G(\hat{H}_\mathrm{cd})^2),
\end{equation}
with $G=\frac{\partial \hat{H}}{\partial t}+i[\hat{H}_\mathrm{cd},\hat{H}]$~\cite{sels2017,kolodrubetz2017}. Hence we need to solve
\begin{equation}
\frac{\partial S}{\partial C_{ij}}=2\operatorname{Tr}\left(G\frac{\partial G}{\partial C_{ij}}\right)=0.
\end{equation}
For notational simplicity, we set $J_{ii}=C_{ii}=0$. We have
\begin{equation}
\begin{split}
G&=-\omega'(t)\sum_p\sigma_z^p-i\omega(t)\sum_{i,j,p}C_{ij}[\sigma_x^i\sigma_y^j,\sigma_z^p]-ig\sum_{i,j,p,q}C_{ij}J_{pq}\sigma_x^i[\sigma_y^j,\sigma_x^p\sigma_x^q],\\
&=-\omega'(t)\sum_p\sigma_z^p+\omega(t)\sum_{i,j,p}C_{ij}(\sigma_x^i\sigma_x^j\delta_{pj}-\sigma_y^i\sigma_y^j\delta_{pi})-g\sum_{i,j,p,q}C_{ij}J_{pq}\sigma_x^i(\sigma_z^p\sigma_x^q\delta_{jp}+\sigma_x^p\sigma_z^q\delta_{jq}),\\
&=-\omega'(t)\sum_p\sigma_z^p+\omega(t)\sum_{i,j}C_{ij}(\sigma_x^i\sigma_x^j-\sigma_y^i\sigma_y^j)-2g\sum_{i,p,q}J_{pq}C_{ip}\sigma_x^i\sigma_x^q\sigma_z^p.
\end{split}
\end{equation}
Hence
\begin{equation}
\frac{\partial G}{\partial C_{mn}}=\omega(t)(\sigma_x^m\sigma_x^n-\sigma_y^m\sigma_y^n)-2g\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n.
\end{equation}
We now want to calculate the trace of $G\partial G/\partial C_{mn}$ ($m\ne n$). Considering the three terms in $G$ and the two terms in $\partial G/\partial C_{mn}$ gives a total of six terms. However three of these are zero, since the trace of terms with an odd number of spin-1/2 operators is zero. The three remaining terms are as follows.

Firstly,
\begin{equation}
\operatorname{Tr}\left[g\omega'(t)\sum_p\sigma_z^p\sum_{\ell}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right]=2^N\frac{g\omega'(t)J_{mn}}{16}.
\end{equation}

Next,
\begin{equation}
\operatorname{Tr}\left[\omega(t)^2\sum_{p,q}C_{pq}(\sigma_x^p\sigma_x^q-\sigma_y^p\sigma_y^q)(\sigma_x^m\sigma_x^n-\sigma_y^m\sigma_y^n)\right]=2^N\frac{\omega(t)^2C_{mn}}{4}.
\end{equation}

Finally,
\begin{equation}
\begin{split}
&\operatorname{Tr}\left[g^2\sum_{p,q,r(q,r\ne p)}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=\operatorname{Tr}\left[g^2\sum_{p,q,r(q,r\ne p)}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=\operatorname{Tr}\left[g^2\sum_{p,q}J_{pq}C_{pq}(\sigma_x^q)^2\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right]+\operatorname{Tr}\left[g^2\sum_{p,q,r(\mathrm{distinct})}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=2^N \frac{g^2\sum_q J_{nq}C_{nq}J_{nm}}{64}+\operatorname{Tr}\left[g^2\sum_{p,q,r(\mathrm{distinct})}J_{pq}C_{pr}\sigma_x^r\sigma_x^q\sigma_z^p\sum_{\ell\ne m,n}J_{n\ell}\sigma_x^m\sigma_x^\ell\sigma_z^n\right],\\
&=2^N \frac{g^2\sum_q J_{mn}J_{nq}C_{nq}}{64}+2^N\frac{g^2\sum_{\ell\ne m,n}\left(J_{n\ell}^2C_{mn}+J_{mn}J_{n\ell}C_{n\ell}\right)}{64},\\
&=2^N \frac{g^2(J_{mn}(JC)_{nn}+(J^2)_{nn}C_{mn}/2-J_{mn}^2C_{mn})}{32}.
\end{split}
\end{equation}

Combining terms gives the coupled linear equations that determine $C_{mn}$,
\begin{equation}
g\omega'(t)J_{mn}+2\omega(t)^2C_{mn}+g^2(J_{mn}(JC)_{nn}+(J^2)_{nn}C_{mn}/2-J_{mn}^2C_{mn})=0.
\end{equation}
This gives
\begin{equation}\label{Cmn}
C_{mn}=-\frac{g\omega'(t)J_{mn}+g^2J_{mn}(JC)_{nn}}{2\omega(t)^2+g^2(J^2)_{nn}/2-g^2J_{mn}^2}.
\end{equation}
Multiplying both sides by $g^2J_{nm}$ and summing over $m$ gives
\begin{equation}\label{Cnn}
g^2(JC)_{nn}=-(g\omega'(t)+g^2(JC)_{nn})f_n,
\end{equation}
with
\begin{equation}
f_n=\sum_m\frac{g^2J_{mn}^2}{2\omega(t)^2+g^2\left[(J^2)_{nn}/2-J_{mn}^2\right]}.
\end{equation}
Equation~\eqref{Cnn} can be rearranged to give
\begin{equation}
g^2(JC)_{nn}=-\frac{g\omega'(t)f_n}{1+f_n}.
\end{equation}
Substituting this into Eq.~\eqref{Cmn} gives
\begin{equation}
C_{mn}=-\frac{g\omega'(t)\left(\frac{1}{1+f_n}\right)J_{mn}}{2\omega(t)^2+g^2\left[(J^2)_{nn}/2-J_{mn}^2\right]}.
\end{equation}
For an infinite chain with $J_{mn}=1/|m-n|^p$, $(J^2)_{nn}/2=\zeta(2p)$ with $\zeta(s)$ the Riemann zeta function and $f_n=f$ is independent of $n$. Furthermore, we find that for integer values of $p$ we can evaluate $f$ analytically using Mathematica (for example, for $p=1$, $f=1-x\cot x$ with $x=g\pi/\sqrt{g^2\zeta(2)+64\omega(t)^2}$). For nearest neighbour interactions, $f=g^2/\omega(t)^2$, and
\begin{equation}
C_{mn}=-\frac{g\omega'(t)}{2(\omega(t)^2+g^2)}\delta_{n,m+1}.
\end{equation}


\section{Introduction}
Submodular function refers to a broad class of functions which satisfy the natural diminishing returns property: adding an additional item to a larger existing subset is less beneficial. A wide range of machine learning and AI problems, including exemplar-based clustering \citep{dueck2007non}, feature selection \citep{das2008algorithms},  active learning \citep{golovin2011adaptive}, influence maximization in social networks \citep{tang2020influence}, recommender system \citep{el2011beyond}, and diverse data summarization \citep{sipos2012temporal}, can be formulated as a submodular maximization problem. This problem, whose goal is to select a set of items to maximize a submodular function, has been extensively studied in the literature subject to various constraints, including cardinality, matroid, or knapsack-type restrictions.

 We notice that in practise, items or individuals are often associated with different groups based on various attributes, such as gender, race, age, religion, or other factors. Existing algorithms might exhibit bias if left unchecked, for example, some of the groups might be over- or under-represented in the final selected subset. Therefore, it becomes increasingly important to design \emph{fairness-aware} algorithms to mitigate such issues. Towards this end, we propose and study the classic non-monotone submodular maximization problem subject to novel group fairness constraints.
 Our goal is to select a \emph{balanced} set of items that maximizes a non-monotone submodular function, such that the ratio of selected items from each group to its size is within a desired range, as determined by the decision maker. Formally, we consider a set $V$ of items (e.g., datapoints) which are partitioned into $m$ groups: $V_1, V_2, \cdots, V_m$ such that items from the same group share same attributes (e.g., gender). We say that a set $S\subseteq V$ of items is $(\alpha,\beta)$-fair if for all groups $i \in[m]$, it holds that $\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i| \leq \lfloor\beta |V_i|\rfloor$.  Using our model, it allows for the decision maker to specify the desired level of fairness by setting appropriate values of $\alpha$ and $\beta$. Specifically, setting $\alpha=\beta$ leads to the highest level of fairness in that the number of selected items is \emph{strictly} proportional to its group size; if we set $\alpha=0$ and $\beta=1$, there are no fairness constraints. Our goal is to find such a $(\alpha,\beta)$-fair subset of items that maximizes a submodular objective function. Our definition of fairness, which balances solutions with respect to sensitive attributes, has gained widespread acceptance in the academic community, as demonstrated by its frequent use in previous studies \citep{celis2018multiwinner,el2020fairness,chierichetti2019matroids}. There are several other notations of fairness that can be captured by our formulation such as the $80\%$-rule \citep{biddle2017adverse}, statistical parity \citep{10.1145/2090236.2090255} and proportional representation \citep{monroe1995fully}.

\subsection{Our Contributions}

\begin{itemize}
\item Our study breaks new ground by examining the classic (non-monotone) submodular maximization problem under $(\alpha,\beta)$-fairness constraints. Our model offers flexibility in capturing varying degrees of fairness as desired by the decision maker, by adjusting the values of $\alpha$ and $\beta$.
\item We develop the first constant-factor approximation algorithm for this problem. We observe that the parameter $\alpha$ is closely linked to the complexity of solving the $(\alpha,\beta)$-fair non-monotone submodular maximization problem. In particular, when $\alpha\leq 1/2$, we design a $\frac{\gamma}{2}$-approximation algorithm and when $\alpha> 1/2$, we develop a $\frac{\gamma}{3}$-approximation algorithm, where $\gamma$ is the approximation ratio of the current best algorithm for matroid-constrained submodular maximization. We also extend the basic model to incorporate an additional global size constraint on the total number of selected items. We provide approximation algorithms that have a constant-factor approximation ratio for this extended model.
\end{itemize}

\subsection{Additional Related Works}
In recent years, there has been a growing awareness of the importance of  fair and unbiased decision-making systems. This has led to an increased interest in the development of fair algorithms in a wide range of applications, including  influence maximization \citep{tsang2019group}, classification \citep{zafar2017fairness}, voting \citep{celis2018multiwinner}, bandit learning \citep{joseph2016fairness},  and data summarization \citep{celis2018fair}. Depending on the specific context and the type of bias that one is trying to mitigate, existing studies  adopt different metrics of fairness. This can lead to different optimization problems and different fair algorithms that are tailored to the specific requirements of the application. Our notation of fairness is general enough to capture many existing notations such as the $80\%$-rule \citep{biddle2017adverse}, statistical parity \citep{10.1145/2090236.2090255} and proportional representation \citep{monroe1995fully}.  \cite{el2020fairness} develop fair algorithms in the context of streaming submodular maximization. Their proposed notation of fairness is more general than ours, leading to a more challenging optimization problem which does not admit any constant-factor approximation algorithms. Very recently, \cite{tang2022group} extend studies of fair algorithms to a more complicated adaptive setting and they propose a new metric of fairness called group equality.

\section{Preliminaries and Problem Statement}

We consider a set $V$ of $n$ items. There is
a non-negative submodular utility function $f: 2^V\rightarrow \mathbb{R}_+$. Denote by $f(e\mid S)$ the marginal utility of $e\in V$ on top of $S\subseteq V$. We say a function  $f: 2^V\rightarrow \mathbb{R}_+$ is submodular if for any two sets $X, Y\subseteq V$ such that $X\subseteq Y$ and any item $e \in V\setminus Y$,
\[f(e\mid Y) \leq f(e\mid X).\]

Submodular maximization with group fairness. Assume $V$ is partitioned into $m$ groups: $V_1, V_2,\cdots, V_m$. We assume that there is a given
lower and upper bound on the fraction of items of each group that must be contained in a feasible solution. These two bounds, namely $\alpha$ and $\beta$, represent group fairness constraints. The problem of $(\alpha, \beta)$-fair submodular maximization problem  (labelled as $\textbf{P.0}$) can be written as follows.

 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.0}$
$\max f(S)$ \\
\textbf{subject to:}  \\
$\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i| \leq \lfloor\beta |V_i|\rfloor, \forall i \in[m]$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

One can adjust the degree of group fairness in a feasible solution through choosing appropriate values of $\alpha$ and $\beta$. I.e., strict group fairness is achieved at $\alpha=\beta$ in which case every feasible solution must contain the same $\alpha$ fraction of items from each group; if we set $\alpha=0$ and $\beta=1$, then there is no group fairness constraints. We next present the hardness result of this problem.

\begin{lemma}
Problem $\textbf{P.0}$ is NP-hard.
\end{lemma}
\emph{Proof:} We prove this through reduction to the classic \emph{cardinality constrained submodular maximization problem} which we define below.
\begin{definition} The input of  cardinality constrained submodular maximization problem is a group of items $U$, a submodular function $h: 2^U\rightarrow \mathbb{R}_+$, and a cardinality constraint $b$; we aim to select a group of items $S\subseteq U$ such that $h(S)$ is maximized and $|S|\leq b$.
\end{definition}

We next show a reduction from cardinality constrained submodular maximization problem to $\textbf{P.0}$.
Consider any given instance of cardinality constrained submodular maximization problem, we construct a corresponding instance of $\textbf{P.0}$ as follows: Let $V=U$, $f=h$, assume there is only one group, i.e., $V=V_1$, and let $\alpha=0$, $\beta=b/|U|$. It is easy to verify that these two instances are equivalent. This finishes the proof of the reduction. $\Box$

\section{Non-monotone Submodular Maximization with Group Fairness}
\label{sec:adaptive-general}
\paragraph{Warm-up: Monotone Utility Function} For the special case of our problem when $f$ is monotone and submodular, it is easy to verify that $\textbf{P.0}$ simplifies to $\textbf{P.1}$. This is because in this case, increasing the size of a solution by adding more items will not decrease its utility. As a result, the lower bound constraints in $\textbf{P.0}$, which state that $\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i| \leq \lfloor\beta |V_i|\rfloor$ for all $i \in[m]$, can always be met by adding sufficient items to the solution.

 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.1}$
$\max f(S)$ \\
\textbf{subject to:}  \\
$|S \cap V_i| \leq \lfloor\beta |V_i|\rfloor, \forall i \in[m]$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

Since $f$ is a monotone submodular function, $\textbf{P.1}$ is a well-known problem of maximizing a monotone submodular function subject to matroid constraints\footnote{A matroid is a pair $\mathcal{M} = (V, \mathcal{I})$ where $\mathcal{I} \subseteq 2^V$ and
1. $\forall Y \in \mathcal{I}, X \subseteq Y \rightarrow X \in \mathcal{I}$.
2. $\forall X, Y \in \mathcal{I}; |X| < |Y| \rightarrow \exists e\in Y \setminus X; X \cup\{ e\} \in \mathcal{I}$.}. This problem has a $(1-1/e)$-approximation algorithm.

We then proceed to develop approximation algorithms for non-monotone functions. We will examine two scenarios, specifically when $\alpha\leq 1/2$ and when $\alpha> 1/2$.
\subsection{The case when $\alpha\leq 1/2$}
In the scenario where $\alpha\leq 1/2$, we use the solution of $\textbf{P.1}$ as a building block to construct our algorithm. First, it is easy to verify that $\textbf{P.1}$ is relaxed version of  $\textbf{P.0}$ with lower bound constraints $\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i|$ in $\textbf{P.0}$ being removed. Because $f$ is a submodular function, $\textbf{P.1}$ is a classic problem of maximizing a submodular function subject to matroid constraints. There exist effective solutions for $\textbf{P.1}$. Now we are ready to present the design of our algorithm as below.
\begin{enumerate}
\item Apply the state-of-the-art algorithm $\mathcal{A}$ for matroid constrained submodular maximization to solve $\textbf{P.1}$ and obtain a solution $A^{P.1}$.

\item Note that $A^{P.1}$ is not necessarily a feasible solution to $\textbf{P.0}$ because it might violate the lower bound constraints $\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i|$ for some groups. To make it feasible, we add additional items to $A^{P.1}$. Specifically, for each group $i\in[m]$ such that $|A^{P.1} \cap V_i|< \lfloor\alpha |V_i|\rfloor$, our algorithm selects a backup set $B_i$ of size $\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|$, by randomly sampling  $\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|$ items from $V_i\setminus A^{P.1}$. Define $B_i=\emptyset$ if $|A^{P.1} \cap V_i|\geq \lfloor\alpha |V_i|\rfloor$.
\item At the end, add $\cup_{i\in [m]}B_i$ to $A^{P.1}$ to build the final solution $A^{approx}$, i.e., $A^{approx}= A^{P.1}\cup(\cup_{i\in [m]}B_i)$.
\end{enumerate}

\begin{algorithm}[hptb]
\caption{Approximation Algorithm for $\textbf{P.0}$ when $\alpha\leq1/2$}
\label{alg:2}
\begin{algorithmic}[1]
\STATE Apply  $\mathcal{A}$  to solve $\textbf{P.1}$ and obtain a solution $A^{P.1}$
\FOR {every group $i\in[m]$}
\IF {$|A^{P.1} \cap V_i|< \lfloor\alpha |V_i|\rfloor$}
\STATE  select a random backup set $B_i$ of size $\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|$ from $V_i\setminus A^{P.1}$
\ELSE
\STATE $B_i\leftarrow \emptyset$
\ENDIF
\ENDFOR
\STATE $A^{approx}\leftarrow A^{P.1}\cup(\cup_{i\in [m]}B_i)$
\RETURN $A^{approx}$
\end{algorithmic}
\end{algorithm}

The pseudocode of this approximation algorithm is given as Algorithm \ref{alg:2}. Observe that   $A^{P.1}$ is a feasible solution to $\textbf{P.1}$, hence, $A^{P.1}$  satisfies upper bound constraints of $\textbf{P.1}$ and hence $\textbf{P.0}$, i.e., $|S \cap V_i| \leq \lfloor\beta |V_i|\rfloor, \forall i \in[m]$. According to the construction of $B_i$, it is easy to verify that  adding $\cup_{i\in [m]}B_i$ to $A^{P.1}$ does not violate the upper bound constraints because $\cup_{i\in [m]}B_i$ are only supplemented to those groups which do not satisfy the lower bound constraints of  $\textbf{P.0}$, i.e., $\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i|$. Moreover, adding $\cup_{i\in [m]}B_i$ to $A^{P.1}$ makes it satisfy lower bound constraints of   $\textbf{P.0}$. Hence, $A^{approx}$ is a feasible solution to $\textbf{P.0}$.
 \begin{lemma}
 $A^{approx}$ is a feasible solution to $\textbf{P.0}$.
 \end{lemma}
\subsubsection{Performance Analysis}
\label{sec:performance1}

We next analyze the performance of Algorithm \ref{alg:2}. We first introduce a useful lemma from  \citep{buchbinder2014submodular}.

\begin{lemma}If $f$ is submodular and $S$ is a random subset of $V$, such that each item in $V$ is contained in
$S$ with probability at most $p$, then
$\mathbb{E}_A[f(A)] \geq (1-p)f(\emptyset)$.
\label{lem:ama}
\end{lemma}

The next lemma states that if $A^{P.1}$ is a $\gamma$-approximate solution of $\textbf{P.1}$, then $f(A^{P.1})$ is at least $\gamma$ fraction of the optimal solution of $\textbf{P.0}$.
 \begin{lemma}
 \label{lem:seq}
Suppose $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint. Let $OPT$ denote the optimal solution of $\textbf{P.0}$, we have $f(A^{P.1})\geq \gamma f(OPT)$.
 \end{lemma}
\emph{Proof:} Because $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint, we have $f(A^{P.1})\geq \gamma f(O^{P.1})$ where $O^{P.1}$ denotes the optimal solution of $\textbf{P.1}$. Moreover, because $\textbf{P.1}$ is a relaxed version of $\textbf{P.0}$, we have  $ f(O^{P.1})\geq f(OPT)$. Hence, $f(A^{P.1})\geq \gamma f(OPT)$. $\Box$

 We next show that augmenting $A^{P.1}$ with  items from the random set $\cup_{i\in [m]}B_i$ reduces its utility  by a factor of at most $1/2$ in expectation. Here the expectation is taken over the distribution of $\cup_{i\in [m]}B_i$.

 \begin{lemma}
 \label{lem:seq2}
Suppose  $\alpha\leq 1/2$, we have $\mathbb{E}_{A^{approx}}[f(A^{approx})]\geq \frac{1}{2}f(A^{P.1})$ where $A^{approx}= A^{P.1}\cup(\cup_{i\in [m]}B_i)$.
 \end{lemma}
\emph{Proof:} Recall that $B_i= \emptyset$ for all $i\in[m]$ such that $|A^{P.1} \cap V_i|\geq \lfloor\alpha |V_i|\rfloor$, hence, adding those $B_i$ to $A^{P.1}$ does not affect its utility. In the rest of the proof we focus on those $B_i$ with
\begin{eqnarray}
\label{eq:aamas}
|A^{P.1} \cap V_i|< \lfloor\alpha |V_i|\rfloor.
 \end{eqnarray}
 Recall that for every $i\in[m]$ such that $|A^{P.1} \cap V_i|< \lfloor\alpha |V_i|\rfloor$, $B_i$ is a random set of size $\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|$ that is sampled from $V_i\setminus A^{P.1}$. It follows that each item in  $V_i\setminus A^{P.1}$ is contained in  $B_i$  with probability at most
\begin{eqnarray}
\label{eq:haha}
\frac{\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|}{|V_i\setminus A^{P.1}|}.
\end{eqnarray}

We next give an upper bound of (\ref{eq:haha}). First,
\begin{eqnarray}
\label{eq:haha1}
\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|\leq \lfloor\alpha |V_i|\rfloor\leq  |V_i|/2,
\end{eqnarray} where the second inequality is by the assumption that $\alpha\leq 1/2$.
Moreover,
\begin{eqnarray}
&&|V_i\setminus A^{P.1}|=|V_i|- |A^{P.1} \cap V_i|=(\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|)+ (|V_i|-\lfloor\alpha |V_i|\rfloor)\\
&&\geq (\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|)+|V_i|/2,\label{eq:haha2}
\end{eqnarray}
where the inequality is by the assumption that $\alpha\leq 1/2$.

Hence,
 \begin{eqnarray}
 (\ref{eq:haha}) \leq \frac{\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|}{ (\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|)+|V_i|/2}\leq \frac{|V_i|/2}{|V_i|/2+|V_i|/2}=1/2,
  \end{eqnarray}
 where the first inequality is by  (\ref{eq:haha2}); the second inequality is by (\ref{eq:haha1}) and the assumption that $\lfloor\alpha |V_i|\rfloor-|A^{P.1} \cap V_i|>0$ (listed in (\ref{eq:aamas})).
 That is, the probability that each item  in $V_i\setminus A^{P.1}$ is contained in $B_i$ is at most $1/2$. It follows that the probability that each item in $V \setminus A^{P.1}$ is contained in $\cup_{i\in [m]}B_i$ is at most $1/2$. Moreover,
Lemma \ref{lem:ama} states that if $f$ is submodular and $S$ is a random subset of $V$, such that each item in $V$ appears in
$S$ with probability at most $p$, then
$\mathbb{E}_A[f(A)] \geq (1-p)f(\emptyset)$. With the above discussion and the observation that $f(A^{P.1}\cup\cdot)$ is submodular, it holds that $\mathbb{E}_{A^{approx}}[f(A^{approx})]=\mathbb{E}_{\cup_{i\in [m]}B_i}[f(A^{P.1}\cup(\cup_{i\in [m]}B_i))]\geq  (1-\frac{1}{2})f(A^{P.1}\cup\emptyset)=\frac{1}{2}f(A^{P.1})$.
$\Box$

Our main theorem as below follows from Lemma \ref{lem:seq} and Lemma \ref{lem:seq2}.

 \begin{theorem}
 \label{thm:1}
Suppose $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint and $\alpha\leq 1/2$, we have $\mathbb{E}_{A^{approx}}[f(A^{approx})]\geq \frac{\gamma}{2}f(OPT)$.
 \end{theorem}

One option of $\mathcal{A}$  is the continuous double greedy algorithm proposed in \citep{feldman2011unified} which gives
 a $1/e-o(1)$-approximation solution, that is, $\gamma\geq 1/e-o(1)$. This, together with Theorem \ref{thm:1}, implies that $\mathbb{E}_{A^{approx}}[f(A^{approx})]\geq \frac{1/e-o(1)}{2}f(OPT)$.

\subsection{The case when $\alpha> 1/2$}
We next consider the case when  $\alpha> 1/2$. We first introduce a new utility function $g: 2^V\rightarrow \mathbb{R}_+$ as below:
\begin{eqnarray}
g(\cdot) = f(V\setminus \cdot).
\end{eqnarray} We first present a well-known result, which states that   submodular functions maintain their submodularity property when taking their complement.

\begin{lemma}
\label{lem:submodular}
If $f$ is submodular, then $g$ must be submodular.
\end{lemma}

With utility function $g$, we present a new optimization problem $\textbf{P.2}$ as below:

 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.2}$
$\max g(S)$ \\
\textbf{subject to:}  \\
$|V_i|-\lfloor\beta |V_i|\rfloor \leq |S \cap V_i| \leq |V_i|-\lfloor\alpha |V_i|\rfloor, \forall i \in[m]$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

$\textbf{P.2}$ is a flipped version of the original problem $\textbf{P.0}$ in the sense that if there is a $\gamma$-approximate solution $A^{P.2}$ to $\textbf{P.2}$, it can be easily verified that $V\setminus A^{P.2}$ is a $\gamma$-approximate solution to $\textbf{P.0}$. As a result, we will focus on solving $\textbf{P.2}$ for the rest of this section.

To solve $\textbf{P.2}$, we introduce another problem (labeled as $\textbf{P.3}$) as follows:

 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.3}$
$\max g(S)$ \\
\textbf{subject to:}  \\
$|S \cap V_i| \leq |V_i|-\lfloor\alpha |V_i|\rfloor, \forall i \in[m]$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

$\textbf{P.3}$ is relaxed version of  $\textbf{P.2}$ with lower bound constraints $|V_i|-\lfloor\beta |V_i|\rfloor \leq |S \cap V_i|$ in $\textbf{P.2}$ being removed. Because $g$ is a submodular function, $\textbf{P.3}$ is a classic problem of maximizing a submodular function subject to matroid constraints. Now we are ready to present the design of our algorithm.
\begin{enumerate}
\item Apply the state-of-the-art algorithm $\mathcal{A}$ for matroid constrained submodular maximization to solve $\textbf{P.3}$ and obtain a solution $A^{P.3}$.

\item Note that $A^{P.3}$ is not necessarily a feasible solution to $\textbf{P.2}$ because it might violate the lower bound constraints $|V_i|-\lfloor\beta |V_i|\rfloor \leq |S \cap V_i|$  for some groups. We add additional items to $A^{P.3}$ to make it feasible. Specifically, for each group $i\in[m]$ such that $|A^{P.3} \cap V_i|< |V_i|-\lfloor\beta |V_i|\rfloor$, our algorithm selects a backup set $B_i$ of size $|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|$, by randomly sampling  $|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|$ items from $V_i\setminus A^{P.3}$. Define $B_i=\emptyset$ if $|A^{P.1} \cap V_i|\geq |V_i|-\lfloor\beta |V_i|\rfloor$.
\item Add $\cup_{i\in [m]}B_i$ to $A^{P.3}$ to build $A^{approx}$, i.e., $A^{approx}= A^{P.3}\cup(\cup_{i\in [m]}B_i)$. Return $V\setminus A^{approx}$ as the final solution.
\end{enumerate}

\begin{algorithm}[hptb]
\caption{Approximation Algorithm for $\textbf{P.0}$ when $\alpha>1/2$}
\label{alg:3}
\begin{algorithmic}[1]
\STATE Apply  $\mathcal{A}$  to solve $\textbf{P.3}$ and obtain a solution $A^{P.3}$
\FOR {every group $i\in[m]$}
\IF {$|A^{P.3} \cap V_i|< |V_i|-\lfloor\beta |V_i|\rfloor$}
\STATE  select a random backup set $B_i$ of size $|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|$ from $V_i\setminus A^{P.3}$
\ELSE
\STATE $B_i\leftarrow \emptyset$
\ENDIF
\ENDFOR
\STATE $A^{approx}\leftarrow A^{P.3}\cup(\cup_{i\in [m]}B_i)$
\RETURN $V\setminus A^{approx}$
\end{algorithmic}
\end{algorithm}

The pseudocode of this approximation algorithm is given as Algorithm \ref{alg:3}. Observe that $A^{P.3}$ satisfies upper bound constraints of  $\textbf{P.3}$ and hence $\textbf{P.2}$ because  $A^{P.3}$ is a feasible solution to $\textbf{P.3}$.  According to the construction of $B_i$, adding $\cup_{i\in [m]}B_i$ to $A^{P.1}$ does not violate the upper bound constraints because $\cup_{i\in [m]}B_i$ are added to meet the lower bound constraints of  $\textbf{P.2}$ if necessary.  Moreover, adding $\cup_{i\in [m]}B_i$ to $A^{P.3}$ makes it satisfy lower bound constraints of   $\textbf{P.2}$. Hence, $A^{approx}$ is a feasible solution to $\textbf{P.2}$.
 \begin{lemma}
 $A^{approx}$ is a feasible solution to $\textbf{P.2}$.
 \end{lemma}

\subsubsection{Performance Analysis}
\label{sec:performance2}

 We first introduce a technical lemma which states that if $A^{P.3}$ is a $\gamma$-approximate solution of $\textbf{P.3}$, then $f(A^{P.3})$ is at least $\gamma$ fraction of the optimal solution of $\textbf{P.2}$. This lemma follows from the observation that $\textbf{P.3}$ is a relaxation of  $\textbf{P.2}$ .
 \begin{lemma}
 \label{lem:jian}
Suppose $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint. Let $O^{P.2}$ denote the optimal solution of $\textbf{P.2}$, it holds that $g(A^{P.3})\geq \gamma g(O^{P.2})$.
 \end{lemma}

 We next show that augmenting $A^{P.3}$ with  items from $\cup_{i\in [m]}B_i$ reduces its utility  by a factor of at most $2/3$ in expectation.

 \begin{lemma}
 \label{lem:ti}
Suppose  $\alpha> 1/2$, $\mathbb{E}_{A^{approx}}[g(A^{approx})]\geq \frac{1}{3}g(A^{P.3})$ where $A^{approx}= A^{P.3}\cup(\cup_{i\in [m]}B_i)$.
 \end{lemma}
\emph{Proof:} Recall that $B_i= \emptyset$ for all $i\in[m]$ such that $|A^{P.3} \cap V_i|\geq |V_i|-\lfloor\beta |V_i|\rfloor$, hence, adding those $B_i$ to $A^{P.3}$ does not affect its utility. Therefore, we focus on those groups $i\in[m]$ with $|A^{P.3} \cap V_i|< |V_i|-\lfloor\beta |V_i|\rfloor$ in the rest of the proof. Let $M=\{i\mid |A^{P.3} \cap V_i|< |V_i|-\lfloor\beta |V_i|\rfloor\}$ denote the set containing the indexes of all such groups and we assume $M\neq \emptyset$ to avoid trivial cases. We next show that it is safe to assume $\min_{i\in M}|V_i|>1$ without loss of generality, i.e., the smallest group in $M$ contains at least two items. To prove this, we consider two cases, depending on the value of $\beta$. If $\beta=1$, then $|A^{P.3} \cap V_i|< |V_i|-\lfloor\beta |V_i|\rfloor$ does not hold for any group $i$ such that $|V_i|=1$, that is, $\min_{i\in M}|V_i|>1$. If $\beta<1$, then according to the group fairness constraints listed in  $\textbf{P.0}$, we are not allowed to select any items from those groups with $|V_i|=1$. Hence, removing all groups with size one from consideration does not affect  the quality of the optimal solution.

With the assumption that $\min_{i\in M}|V_i|>1$, we are now in position to prove this lemma.  Recall that for every $i\in M$, $B_i$ is a random set of size $|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|$ that is sampled from $V_i\setminus A^{P.3}$. It follows that each item in  $V_i\setminus A^{P.3}$ appears in $B_i$ with probability at most
\begin{eqnarray}
\label{eq:haha3}
\frac{|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|}{|V_i\setminus A^{P.3}|}.
\end{eqnarray}

We next give an upper bound of (\ref{eq:haha3}). Because we assume $\alpha> 1/2$, we have $\beta\geq\alpha>1/2$. This, together with the assumption that $\min_{i\in M}|V_i|>1$, implies that for all $i\in M$,
\begin{eqnarray}
\label{eq:haha4}
|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i| \leq  |V_i|-\lfloor\beta |V_i|\rfloor \leq 2|V_i|/3.
\end{eqnarray}
Moreover,
\begin{eqnarray}
&&|V_i\setminus A^{P.3}|=|V_i|-|A^{P.3} \cap V_i|\\
&&=(|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|)+(|V_i|-(|V_i|-\lfloor\beta |V_i|\rfloor)) \\
&&= (|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|)+\lfloor\beta |V_i|\rfloor \\
&&\geq (|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|)+|V_i|/3, \label{eq:haha5}
\end{eqnarray}
where the inequality is by the observation that $\beta>1/2$. It follows that
\begin{eqnarray}
(\ref{eq:haha3})\leq \frac{|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|}{(|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|)+|V_i|/3}\leq\frac{2|V_i|/3}{2|V_i|/3+|V_i|/3}=2/3,
\end{eqnarray}
where the first inequality is by (\ref{eq:haha5}) and the second inequality is by (\ref{eq:haha4}) and the assumption that $|V_i|-\lfloor\beta |V_i|\rfloor -|A^{P.3} \cap V_i|>0$. That is, each item in $V_i\setminus A^{P.3}$ appears in  $B_i$ with probability at most $2/3$. Lemma \ref{lem:ama} and the observation that $g(A^{P.3}\cup\cdot)$ is submodular imply that $\mathbb{E}_{A^{approx}}[g(A^{approx})]=\mathbb{E}_{\cup_{i\in [m]}B_i}[g(A^{P.3}\cup(\cup_{i\in [m]}B_i))]\geq (1-\frac{2}{3})g(A^{P.3}\cup\emptyset)=  \frac{1}{3}g(A^{P.3})$. $\Box$

Lemma \ref{lem:jian} and Lemma \ref{lem:ti} together imply that $\mathbb{E}_{A^{approx}}[g(A^{approx})]\geq \frac{1}{3}g(A^{P.3})\geq  \frac{\gamma}{3}g(O^{P.2})$. By the definition of function $g$, we have  $\mathbb{E}_{A^{approx}}[f(V\setminus A^{approx})]=\mathbb{E}_{A^{approx}}[g(A^{approx})] \geq  \frac{\gamma}{3}g(O^{P.2})=\frac{\gamma}{3}f(OPT)$ where the last equality is by the observation that $\textbf{P.2}$ and $\textbf{P.0}$ share the same value of the optimal solution. Hence, the following main theorem holds.

 \begin{theorem}
 \label{thm:2}
Suppose $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint and $\alpha> 1/2$, we have $\mathbb{E}_{A^{approx}}[f(V\setminus A^{approx})]\geq \frac{\gamma}{3}f(OPT)$.
 \end{theorem}

If we adopt the continuous double greedy algorithm  \citep{feldman2011unified} as $\mathcal{A}$ to compute  $A^{P.3}$,  it gives
 a $1/e-o(1)$-approximation solution, that is, $\gamma\geq 1/e-o(1)$. This, together with Theorem \ref{thm:2}, implies that $\mathbb{E}_{A^{approx}}[f(V\setminus A^{approx})]\geq \frac{1/e-o(1)}{3}f(OPT)$.
\section{Extension: Incorporating Global Cardinality Constraint}
\label{sec:non-adaptive-extended}

In this section, we extend $\textbf{P.0}$ to incorporate a global cardinality constraint. A formal definition of this problem is listed in $\textbf{P.A}$. Our objective is to find a best $S$ subject to a group fairness constraint $(\alpha,\beta)$ and an additional cardinality constraint $c$.

 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.A}$
$\max f(S)$ \\
\textbf{subject to:}  \\
$\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i| \leq \lfloor\beta |V_i|\rfloor, \forall i \in[m]$.\\
$|S|\leq c$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

\subsection{The case when $\alpha\leq 1/2$}
We first consider the case when $\alpha\leq 1/2$. We introduce a new optimization problem $\textbf{P.B}$ as follows:
 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.B}$
$\max f(S)$ \\
\textbf{subject to:}  \\
$|S \cap V_i| \leq \lfloor\beta |V_i|\rfloor, \forall i \in[m]$.\\
$\sum_{i\in[m]}\max\{\lfloor\alpha |V_i|\rfloor , |S\cap V_i|\} \leq c$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

It is easy to verify that $\textbf{P.B}$ is a relaxation of $\textbf{P.A}$ in the sense that every feasible solution to $\textbf{P.A}$ is also a feasible solution to  $\textbf{P.B}$. Hence, we have the following lemma.

 \begin{lemma}
 \label{lem:ti1}
Let $OPT$ denote the optimal solution of $\textbf{P.A}$ and $O^{P.B}$ denote the optimal solution of $\textbf{P.B}$, we have $f(O^{P.B})\geq f(OPT)$.
 \end{lemma}

It has been shown that the constraints in $\textbf{P.B}$ gives rise to a matroid \citep{el2020fairness}. This, together with the assumption that $f$ is a submodular function, implies that $\textbf{P.B}$ is a classic problem of maximizing a submodular function subject to  matroid constraints. Now we are ready to present the design of our algorithm.
\begin{enumerate}
\item Apply the state-of-the-art algorithm $\mathcal{A}$ for matroid constrained submodular maximization to solve $\textbf{P.B}$ and obtain a solution $A^{P.B}$.

\item Note that $A^{P.B}$ is not necessarily a feasible solution to $\textbf{P.A}$ because it might violate the lower bound constraints $\lfloor\alpha |V_i|\rfloor \leq |S \cap V_i|$ for some groups. To make it feasible, we add additional items to $A^{P.B}$. Specifically, for each group $i\in[m]$ such that $|A^{P.B} \cap V_i|< \lfloor\alpha |V_i|\rfloor$, our algorithm selects a backup set $B_i$ of size $\lfloor\alpha |V_i|\rfloor-|A^{P.B} \cap V_i|$, by randomly sampling  $\lfloor\alpha |V_i|\rfloor-|A^{P.B} \cap V_i|$ items from $V_i\setminus A^{P.B}$. Define $B_i=\emptyset$ if $|A^{P.1} \cap V_i|\geq \lfloor\alpha |V_i|\rfloor$.
\item  At the end, add $\cup_{i\in [m]}B_i$ to $A^{P.B}$ to build the final solution $A^{approx}$, i.e., $A^{approx}= A^{P.B}\cup(\cup_{i\in [m]}B_i)$.
\end{enumerate}

\begin{algorithm}[hptb]
\caption{Approximation Algorithm for $\textbf{P.A}$ when $\alpha\leq 1/2$}
\label{alg:21}
\begin{algorithmic}[1]
\STATE Apply  $\mathcal{A}$  to solve $\textbf{P.B}$ and obtain a solution $A^{P.B}$
\FOR {every group $i\in[m]$}
\IF {$|A^{P.B} \cap V_i|< \lfloor\alpha |V_i|\rfloor$}
\STATE  select a random backup set $B_i$ of size $\lfloor\alpha |V_i|\rfloor-|A^{P.B} \cap V_i|$ from $V_i\setminus A^{P.B}$
\ELSE
\STATE $B_i\leftarrow \emptyset$
\ENDIF
\ENDFOR
\STATE $A^{approx}\leftarrow A^{P.B}\cup(\cup_{i\in [m]}B_i)$
\RETURN $A^{approx}$
\end{algorithmic}
\end{algorithm}

The pseudocode of this approximation algorithm is given as Algorithm \ref{alg:21}. Observe that $A^{P.B}$ satisfies the group-wise upper bound constraints of  $\textbf{P.A}$ because  $A^{P.B}$ meets the first set of constraints in $\textbf{P.B}$. According to the construction of $B_i$, adding $\cup_{i\in [m]}B_i$ to $A^{P.B}$ does not violate the group-wise upper bound constraints of  $\textbf{P.A}$ because $\cup_{i\in [m]}B_i$ are added to meet the lower bound constraints of  $\textbf{P.A}$ if necessary.   Moreover, adding $\cup_{i\in [m]}B_i$ to $A^{P.B}$ does not violate the global cardinality constraint of  $\textbf{P.A}$ because $A^{P.B}$ meets the second set of constraints in $\textbf{P.B}$. At last, it is easy to verify that adding $\cup_{i\in [m]}B_i$ to $A^{P.B}$ makes it satisfy the lower bound constraints of   $\textbf{P.A}$. Hence, $A^{approx}$ is a feasible solution to $\textbf{P.A}$.
 \begin{lemma}
 $A^{approx}$ is a feasible solution to $\textbf{P.A}$.
 \end{lemma}

Following the same proof of Theorem \ref{thm:1}, we have the following theorem.
 \begin{theorem}
 \label{thm:11}
Suppose $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint and $\alpha\leq 1/2$, we have $\mathbb{E}_{A^{approx}}[f(A^{approx})]\geq \frac{\gamma}{2}f(OPT)$.
 \end{theorem}

 \subsection{The case when $\alpha> 1/2$}
We next consider the case when  $\alpha> 1/2$. Recall that $g(\cdot) = f(V\setminus \cdot)$.  We first present a flipped formation of $\textbf{P.A}$ as below:

 \begin{center}
\framebox[0.5\textwidth][c]{
\enspace
\begin{minipage}[t]{0.5\textwidth}
\small
$\textbf{P.C}$
$\max g(S)$ \\
\textbf{subject to:}  \\
$|V_i|-\lfloor\beta |V_i|\rfloor \leq |S \cap V_i| \leq |V_i|-\lfloor\alpha |V_i|\rfloor, \forall i \in[m]$.\\
$|S|\geq n-c$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

Suppose there is a $\gamma$-approximate solution $A^{P.C}$ to $\textbf{P.C}$, it is easy to verify that $V\setminus A^{P.C}$ is a $\gamma$-approximate solution to $\textbf{P.A}$. We focus on solving $\textbf{P.C}$ in the rest of this section. We first introduce a new optimization problem (labeled as $\textbf{P.D}$) as follows:

 \begin{center}
\framebox[0.4\textwidth][c]{
\enspace
\begin{minipage}[t]{0.4\textwidth}
\small
$\textbf{P.D}$
$\max g(S)$ \\
\textbf{subject to:}  \\
$|S \cap V_i| \leq |V_i|-\lfloor\alpha |V_i|\rfloor, \forall i \in[m]$.
\end{minipage}
}
\end{center}
\vspace{0.1in}

$\textbf{P.D}$ is relaxed version of  $\textbf{P.C}$ with both group-wise lower bound constraints $|V_i|-\lfloor\beta |V_i|\rfloor \leq |S \cap V_i|$ and global  lower bound constraints $|S|\geq n-c$ in $\textbf{P.C}$ being removed. Hence, we have the following lemma.

 \begin{lemma}
 \label{lem:ti2}
Let $O^{P.C}$ denote the optimal solution of $\textbf{P.C}$  and $O^{P.D}$ denote the optimal solution of $\textbf{P.D}$, we have $g(O^{P.D})\geq g(O^{P.C})$.
 \end{lemma}

Recall that if $f$ is submodular, $g$ must be submodular (by Lemma \ref{lem:submodular}). Hence, $\textbf{P.D}$ is a classic problem of maximizing a submodular function subject to matroid constraints. We next present the design of our algorithm.
\begin{enumerate}
\item Apply the state-of-the-art algorithm $\mathcal{A}$ for matroid constrained submodular maximization to solve $\textbf{P.D}$ and obtain a solution $A^{P.D}$.

\item Note that $A^{P.D}$ is not necessarily a feasible solution to $\textbf{P.C}$ because it might violate the group-wise or the global lower bound constraints of  $\textbf{P.C}$. We add additional items to $A^{P.D}$ to make it feasible. Specifically, for each group $i\in[m]$, our algorithm selects a backup set $B_i$ of size $|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|$, by randomly sampling  $|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|$ items from $V_i\setminus A^{P.D}$. Define $B_i=\emptyset$ if $|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|=0$.
\item Add $\cup_{i\in [m]}B_i$ to $A^{P.D}$ to build  $A^{approx}$, i.e., $A^{approx}= A^{P.D}\cup(\cup_{i\in [m]}B_i)$. Return $V\setminus A^{approx}$ as the final solution.
\end{enumerate}

\begin{algorithm}[hptb]
\caption{Approximation Algorithm for $\textbf{P.A}$ when $\alpha>1/2$}
\label{alg:31}
\begin{algorithmic}[1]
\STATE Apply  $\mathcal{A}$  to solve $\textbf{P.D}$ and obtain a solution $A^{P.D}$
\FOR {every group $i\in[m]$}
\IF {$|A^{P.D} \cap V_i|< |V_i|-\lfloor\beta |V_i|\rfloor$}
\STATE  select a random backup set $B_i$ of size  $|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|$  from $V_i\setminus A^{P.D}$
\ELSE
\STATE $B_i\leftarrow \emptyset$
\ENDIF
\ENDFOR
\STATE $A^{approx}\leftarrow A^{P.D}\cup(\cup_{i\in [m]}B_i)$
\RETURN $V\setminus A^{approx}$
\end{algorithmic}
\end{algorithm}

The pseudocode of this approximation algorithm is given as Algorithm \ref{alg:31}. Observe that adding $\cup_{i\in [m]}B_i$ to $A^{P.D}$ ensures that each group contributes exactly $|V_i|-\lfloor\alpha |V_i|\rfloor$ number of items to the solution. Because $c\leq \sum_{i\in[m]}(|V_i|-\lfloor\alpha |V_i|\rfloor)$, $A^{P.D}\cup(\cup_{i\in [m]}B_i)$ must satisfy all constraints in  $\textbf{P.C}$. Hence, we have the following lemma.
 \begin{lemma}
 $A^{approx}$ is a feasible solution to $\textbf{P.C}$.
 \end{lemma}

We next analyze the performance of $A^{approx}$. The following lemma states that adding $\cup_{i\in [m]}B_i$ to  $A^{P.D}$ reduces its utility  by a factor of at most $2/3$ in expectation.

 \begin{lemma}
 \label{lem:ti3}
Suppose  $\alpha> 1/2$, we have $\mathbb{E}_{A^{approx}}[g(A^{approx})]\geq \frac{1}{3}g(A^{P.D})$.
 \end{lemma}
\emph{Proof:} Observe that $B_i= \emptyset$ for all $i\in[m]$ such that $|A^{P.D} \cap V_i|= |V_i|-\lfloor\alpha |V_i|\rfloor$, hence, adding those $B_i$ to $A^{P.3}$ does not affect its utility. Therefore, we focus on those groups $i\in[m]$ with $|A^{P.D} \cap V_i|< |V_i|-\lfloor\alpha |V_i|\rfloor$ in the rest of the proof. Let $Z=\{i\in[m]\mid |A^{P.D} \cap V_i|< |V_i|-\lfloor\alpha |V_i|\rfloor\}$ denote the set containing the indexes all such groups. We assume $Z\neq \emptyset$ to avoid trivial cases.  We next show that it is safe to assume $\min_{i\in Z}|V_i|>1$ without loss of generality, i.e., the smallest group in $Z$ contains at least two items. To prove this, we consider two cases, depending on the value of $\alpha$. If $\alpha=1$, then $|A^{P.D} \cap V_i|< |V_i|-\lfloor\alpha |V_i|\rfloor$ does not hold for any group $i$ such that $|V_i|=1$. Hence, $\min_{i\in Z}|V_i|>1$. If $\alpha<1$, then according to the group fairness constraints listed in  $\textbf{P.A}$, we are not allowed to select any items from those groups with $|V_i|=1$. Hence, removing all groups with size one from consideration does not affect the quality of the optimal solution.

With the assumption that $\min_{i\in Z}|V_i|>1$, we are now ready to prove this lemma.  Recall that for every $i\in Z$, $B_i$ is a random set of size $|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|$ that is sampled from $V_i\setminus A^{P.D}$. It follows each item in  $V_i\setminus A^{P.D}$ appears in  $B_i$  with probability at most
\begin{eqnarray}
\label{eq:haha31}
\frac{|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|}{|V_i\setminus A^{P.D}|}.
\end{eqnarray}

We next give an upper bound of (\ref{eq:haha31}). Because we assume $\alpha> 1/2$ and $\min_{i\in Z}|V_i|>1$, it holds that for all $i\in M$,
\begin{eqnarray}
\label{eq:haha41}
|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i| \leq  |V_i|-\lfloor\alpha |V_i|\rfloor \leq 2|V_i|/3.
\end{eqnarray}
Moreover,
\begin{eqnarray}
&&|V_i\setminus A^{P.D}|= |V_i|-|A^{P.D} \cap V_i|\\
&&=(|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|)+(|V_i|-(|V_i|-\lfloor\alpha |V_i|\rfloor)) \\
&&= (|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|)+\lfloor\alpha |V_i|\rfloor \\
&&\geq (|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|)+|V_i|/3, \label{eq:haha51}
\end{eqnarray}
where the inequality is by the assumption that $\alpha>1/2$. It follows that
\begin{eqnarray}
(\ref{eq:haha31})\leq \frac{|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.3} \cap V_i|}{(|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.3} \cap V_i|)+|V_i|/3}\leq\frac{2|V_i|/3}{2|V_i|/3+|V_i|/3}=2/3,
\end{eqnarray}
where the first inequality is by (\ref{eq:haha51}) and the second inequality is by (\ref{eq:haha41}) and the assumption that $|V_i|-\lfloor\alpha |V_i|\rfloor -|A^{P.D} \cap V_i|>0$. That is, each item in  $V_i\setminus A^{P.D}$ appears in $B_i$ with probability at most $2/3$. Lemma \ref{lem:ama} and the observation that $g(A^{P.D}\cup\cdot)$ is submodular imply that $\mathbb{E}_{A^{approx}}[g(A^{approx})]=\mathbb{E}_{\cup_{i\in [m]}B_i}[g(A^{P.D}\cup(\cup_{i\in [m]}B_i))]\geq (1-\frac{2}{3})g(A^{P.D}\cup\emptyset)=\frac{1}{3}g(A^{P.D})$. $\Box$

Suppose  $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint, we have  $\mathbb{E}_{A^{approx}}[g(A^{approx})]\geq \frac{1}{3}g(A^{P.D})\geq \frac{\gamma}{3}g(O^{P.D})$ where the first inequality is by Lemma \ref{lem:ti3}. This, together with $g(O^{P.D})\geq g(O^{P.C})$ (as proved in Lemma \ref{lem:ti2}), implies that  $\mathbb{E}_{A^{approx}}[g(A^{approx})]\geq \frac{\gamma}{3}g(O^{P.C})$. By the definition of function $g$, we have $\mathbb{E}_{A^{approx}}[f(V\setminus A^{approx})] = \mathbb{E}_{A^{approx}}[g(A^{approx})]\geq \frac{\gamma}{3}g(O^{P.C})=\frac{\gamma}{3}f(OPT)$ where the last equality is by the observation that $\textbf{P.A}$ and $\textbf{P.C}$ share the same value of the optimal solution. Hence, the following main theorem holds.

 \begin{theorem}
Suppose $\mathcal{A}$ is a $\gamma$-approximate algorithm for non-monotone submodular
maximization subject to a matroid constraint and $\alpha> 1/2$, we have $\mathbb{E}_{A^{approx}}[f(V\setminus A^{approx})]\geq \frac{\gamma}{3}f(OPT)$.
 \end{theorem}

\section{Conclusion}
This paper presents a comprehensive investigation of the non-monotone submodular maximization problem under group fairness constraints. Our main contribution is the development of several constant-factor approximation algorithms for this problem. In the future, we plan to expand our research to explore alternative fairness metrics.
\bibliographystyle{ijocv081}