# 2D-subband spectra variations under persistent tunnelling photoconductivity
condition in tunnel delta-GaAs/Al structures

S. E. Dizhur<EMAIL_ADDRESS>Institute of Radioengineering and Electronics
RAS, 101999 Moscow, Russia  I. N. Kotel’nikov  Institute of Radioengineering
and Electronics RAS, 101999 Moscow, Russia  V. A. Kokin  Institute of
Radioengineering and Electronics RAS, 101999 Moscow, Russia  F. V. Shtrom
Institute of Radioengineering and Electronics RAS, 101999 Moscow, Russia

###### Abstract

2DEG states of Al/$\delta$-GaAs structures were investigated in the persistent
tunnelling photoconductivity (PTPC) regime at low temperatures. ”Thickening”
of the unoccupied subbands to the ground state of 2DEG was observed at $T=4.2$
$K$. It was found that there is a uniform shift of the subband levels at
$T=77$ $K$. The behaviour of the persistent 2D states in tunnelling spectra
after various illumination and at different temperatures allows us to assume
that there are of two PTPC mechanisms, namely, the broadening of the potential
well profile of the $\delta$-layer and photoionization of the DX-centers.

##  I Introduction

We began to take interest in persistent tunnelling photoconductivity effect
during experiments concerning the intersubband resonant polaron in the
Al/$\delta$-GaAs structures Kotel’nikov  _et al._ (2000), in which the
diamagnetic shift of the 2D subbands in the $\delta$-layer was used. The
positions of the 2D subbands can be changed by persistent photoconductivity
(PPC) too. This effect well known for the case of occupied subbands of the
$\delta$-layers in the bulk of GaAs (see, for instance, A.Zrenner (1987);
Koenraad  _et al._ (1990)). The concentration $n_{2D}$ of the two dimensional
electron gas (2DEG) increases after the pulse illumination of such structures
at liquid helium temperatures. Occupations of the 2D subbands in the PPC
regime are determined from magnetotransport measurements. The origin of the
increase of $n_{2D}$ of the $\delta$-doped structures in the PPC regime is not
clear yet. This effect is explained by the influence of the DX-centers,
interface states, deep levels of the GaAs gap, etc Arscott, Missous, and
Dobaczewski (1992); de Oliveira  _et al._ (1995); Li and Jagadish (1997);
Valyaev  _et al._ (1998). In recent years communications about new features of
the PPC in $\delta$-GaAs at low-doping densities Chen  _et al._ (1998) and
about the effects of charge correlations in the samples of $\delta$-GaAlAs
Piotrzkowski  _et al._ (1999) appear. As far as we know, the investigations of
PPC in the near-surface $\delta$-layers ($\gtrsim 30$ $nm$ from GaAs surface)
were conducted by the capacitance-voltage profiles technique and by the deep
level transient spectroscopy only (see, e.g. Arscott, Missous, and Dobaczewski
(1992)). But such methods do not allow to analyze the energy spectrum of 2DEG
in the $\delta$-layers. Though the near-surface $\delta$-layer can be
characterized by tunnelling spectroscopy Kotel’nikov  _et al._ (2000), these
experiments in PPC regime were not carried out. The tunnelling spectroscopy,
as distinct from the traditional magnetotransport measurements, allows to
obtain spectra of both occupied and unoccupied 2D subbands directly. Therefore
the modification of the self-consistent potential profile near $\delta$-doped
layer can be studied under the PPC conditions. The first persistent tunnelling
photoconductivity experiments in the Al/$\delta$-GaAs structures demonstrated
that a significant variation of the 2D subband spectrum occurs in the
unoccupied subbands 7 (2001).

##  II Samples and experiment

###  II.1 Samples

The tunnel structures Al/$\delta$-GaAs were prepared in IRE RAS on semi-
insulating (100) GaAs substrates by the molecular beam epitaxy (MBE) method.
The $\delta$-doped layer with the density of Si atoms $5.2\times 10^{12}$
$cm^{-2}$ was formed at $570~{}{{}^{\circ}}C$ at the distance of $L\simeq 20$
$nm$ from Al/GaAs interface. Deposition of Al from the Knudsen cell took place
directly in the MBE chamber Kotel’nikov  _et al._ (2000). The width of a
semitransparent Al film was 80–100 $nm$. Al/$\delta$-GaAs tunnel junctions
were photolitographically defined and etched into samples. The ohmic contacts
to the $\delta$-layer were made by Au-Ge-Ni alloy. The sketchy sectional view
of this structure is shown at Figure 1. Three different geometric of samples
(’a’, ’b’, ’c’) were used in our measurements as shown in the same figure.

Figure 1:  Microphotographs of the different samples ’a’, ’b’, ’c’ and sketchy
sectional view of Al/$\delta$-GaAs structure. The dimensions of the octagonal
Al electrodes of the samples ’a’ and ’b’ are 0.7 and 0.2 $mm$ respectively,
the width of Al strips of the sample ’c’ is 10 $\mu m$.

The charge of the impurity of the $\delta$-layer forms a strong asymmetric
(due to its closeness to the interface) V-shaped potential well with 2DEG. 2D
electron sheet concentration in our samples was in the range of $0.5\times
10^{12}$ to $1.5\times 10^{12}$ $cm^{-2}$, as it was found from tunnel spectra
(TS); 2DEG mobility was equal approx to 1500 $cm^{2}/(Vs)$. The band diagram
of such QW is shown on Figure 2, where the potential profile was calculated
self-consistently. Here $\Phi_{S}$ is a barrier height at Al/GaAs interface;
$E_{F}$ – Fermi energy in GaAs; $E_{0},E_{1}\ldots$ – subbands levels in
$\delta$-layer QW.

Figure 2:  Band diagram of $\delta$-doped structures. Self-consistent
calculation of the potential well profile was made with the following
parameters (best approximation for the samples d9): donor (Si) concentration
in $\delta$-layer $N^{\delta}_{Si}=4.7\times 10^{12}$ $cm^{-2}$; distance
between $\delta$-layer and Al/GaAs interface $z_{\delta}=20$ $nm$;
$\delta$-layer width $\Delta z_{\delta}=2.5$ $nm$; $\Phi_{S}=0.87$ $eV$; bias
voltage at tunnel junction $U=0$; acceptor concentrations in GaAs
$N_{a}=5.5\times 10^{15}$ and $1\times 10^{15}$ $cm^{-3}$ (dashed line). The
insert shows potential profile in the bulk of GaAs.

###  II.2 Experiment

The experiments were made at temperatures from 4.2 to 300 $K$. TS was obtained
by standard lock-in technique using the 4-contact method of measurement. Two
lock-in amplifiers simultaneously detected first and second harmonics of the
signal, proportional to $dU/dI$ and $d^{2}U/dI^{2}$ of the U-I curve
respectively. Tunnel spectrum $S(U)=d\ln\sigma/dU$, where $\sigma=dI/dU$, was
calculated from experimental data. Tunnel conductance $\sigma(U)$ and tunnel
spectrum $S(U)$ at temperatures $T=4.2$ and 77 $K$ are shown in Figure 3.

Figure 3:  Tunnel conductance $\sigma$ (bottom curves and left axis) and its
logarithmic derivative $S=d\ln\sigma/dU$ (top curves and right axis) at
temperatures $T=4.2$ and 77 $K$ for sample d9a11 ($n_{2D}=1.1\times 10^{12}$
$cm^{-2}$).

Similar curves were obtained for different geometric of the samples (Figure 1)
at the same $n_{2D}$. We make use of tunnel spectra in $S(U)$ form, because of
convenience from several points of view. First, this value is independent of
absolute magnitude of conductance, that allow us to compare TS of the
different metal electrodes diameters. Second, $S(U)$ is measured in universal
units of $V^{-1}$. When conductance is the exponential function of the bias
voltage $U$ then $S(U)=const$. For our junctions Al/$\delta$-GaAs potential
barrier is extremely close to linear function of the coordinate (see Figure
2). It is easy to show Kotel’nikov  _et al._ (2000) (for quasi-classical
approximation and one occupied subband energy $E_{0}$) that at bias voltages
$U\simeq+100$ $mV\ll(\Phi_{S}+E_{F})/e\simeq 1$ $V$ the tunnel conductance

$\sigma\sim
D(E_{0},eU)\sim\exp\left(\frac{eU-2(\Phi_{S}+E_{F})}{\mathcal{E}_{00}}\right),$

where

$\mathcal{E}_{00}=\frac{3}{2}\frac{\hbar}{\sqrt{2m}}\frac{\sqrt{\Phi_{S}+E_{F}}}{z_{\delta}}$

.

Figure 3 shows that such a simplified model describes well TS behaviour at
large positive $U$. The valuation $e/\mathcal{E}_{00}\simeq 14$ $V^{-1}$ is
easy determine from the shown spectra in the region $U\simeq 100$ $mV$, and
gives quantity of the characteristic energy $\mathcal{E}_{00}\simeq 71$ $meV$.
The calculation by formula above gives good agreement with the experiment for
parameters of junction $\Phi_{S}+E_{F}=1.1$ $eV$ and $z_{\delta}=17$ $nm$.
Each strong dip in TS is related to the bottom $E_{i}$ of the two-dimensional
subbands in the ${\delta}$-layer. Tunnel conductance in 2DEG tunnel structure
is combined of the conductance of 2D subbands:
$\sigma(U)=\sum_{i}\sigma_{i}\propto\sum_{i}\rho_{\|i}(E_{i}-eU)\times
D(E_{i},E_{F}-eU)$, where $i$ is the subband’s number,
$\rho_{\|i}=(m/\pi\hbar^{2})\times\Theta(E_{F}-E_{i}-eU)$ is the two-
dimensional density of states (2DOS), $D$ is the barrier transmission, $E_{F}$
is Fermi energy and $U$ is bias voltage on the tunnel junction; function
$\Theta$, by definition, is $\Theta\equiv 1$ for $x\geq 0$ and zero otherwise.
The derivation of the stepped 2DOS gives the dips seen in TS. Thus the
positions $U_{i}$ of the dips in $d\ln\sigma/dU$ give the subband energies
$E_{i}$ in 2DEG (with the accuracy of weak dependance $E_{i}(U)$ Kotel’nikov
_et al._ (2000)). Figure 3 shows that the amplitudes of the dips decrease at
$T=77$ $K$. It is associated with degradation of the tunnel spectroscopy
resolution with increasing temperature. But the positions $U_{i}$ of the dips
are still rather well measured at $T=77$ $K$. Note that the dips at $U>0$ and
$U<0$ correspond to occupied and unoccupied subbands, respectively. The many-
body features in TS at 4.2 $K$ (near $U=0$ and $U=\hbar\omega_{LO}/e=36.5$
$mV$) will not be discussed in this paper.

In PTPC experiments, the samples were cooled to the temperature 4.2 $K$ under
dark conditions, then they were illuminated by the pulses of the red LED
mounted in cryostat ($h\nu>E_{g}$, where $E_{g}$ is GaAs energy gap).

##  III Results and discussion

Figure 4 demonstrates TS of the junction Al/$\delta$-GaAs with concentration
$n_{2D}=1.1\times 10^{12}$ $cm^{-2}$ at $T=4.2$ $K$ after cooling under dark
conditions (upper curve) and different illuminations (bottom curves).

Figure 4:  a) tunnel spectra at the dark (top curve) and different
illuminations (bottom curves). Signs $E_{i}$ show subbands positions. Spectra
at illumination time $t=0.2$ and 2 $s$ are shifted from dark spectrum at 50
and 90 $V^{-1}$. Dashed line on the bottom curve shows polynomial
interpolation of one of the subbands. b) tunnel conductance at the dark (solid
line) and after illumination time $t=2$ $s$ (dash line). Temperature is 4.2
$K$. Sample d9b91.

According to Figure 4, only one subband $E_{0}$ is occupied in this sample and
the value of Fermi energy is $eU_{0}=-E_{0}\approx 40$ $meV$. In all tested
samples only $E_{0}$ subband was partially filled under dark conditions and
the position $U_{0}$ in TS at different 2D electron sheet concentrations was
changed from 20 to 50 $mV$. Illumination of the samples by LED at liquid
helium temperatures in time $t$ led to variation of the tunnel conductance and
of the dips positions in TS (bottom curves on Figure 4) 7 (2001). It is well
seen from Figure 4 that the main consequence of illumination is the change of
the 2D subband spectra, instead of the conductance change. Observed PTPC and
subband spectra variation depend upon ”dark” concentration of 2D electrons and
reaches saturation that corresponds to the illumination time. The persistent
state of the illuminated tunnel structure was not actually changed in several
hours: for example, the tunnel conductance decreased for less than 0.1% from
the initial value in one hour. But nevertheless, the dip positions were not
changed within the accuracy of our measurements.

Heating of the sample to room temperature with following cooling in the dark
returned the structures parameters to the initial (”dark”) values. The ground
state position weakly depends on $t$ that corresponds to magnetotransport
measurements for the occupied subbands of the ”bulk” $\delta$-layers at PPC
regime. The difference $E_{0}(0)-E_{0}(t)$ at saturation corresponds to
increase of 2D electron concentration at 10–15%. Figure 5 demonstrates
subbands energies (measured from $E_{F}$) as functions of $t$, calculated from
dip positions in TS.

Figure 5:  Subband energies $E_{i}$ of the samples d7c52 and d5b41 as
functions of the illumination time $t$ at $T=4.2$ $K$. The thicknesses of Al
films of these junctions were 80 and 100 $nm$ respectively. Different signs
corresponds to different series of experiments.
$\varepsilon_{LO}=\hbar\omega_{LO}$ is the LO-phonon energy.

For more precise determination of minima of the dips, the approximation with
polynomial of power 3–5 was used (see, for instance, the dashed curve for
$E_{0}$ on Figure 4). It is seen that for empty subbands the shifts are 10–50
$meV$ and for filled level $E_{0}$ shift is 2–5 $meV$. Note that at some $t$
on Figure 5 the resonant conditions for intersubband polaron interaction are
realized too. We will discuss theses effects in future papers.

The obtained results reveal two effects in our structures at liquid helium
temperatures due to illumination: the ”thickening” of 2D levels in QW of the
$\delta$-layer and the increase of the 2DEG concentration. Let us discuss
these effects. Our experiments show that the change $\Delta\sigma_{S}$ of the
tunnelling conductivity at the PTPC saturation regime does not exceed 10% (see
Figure 4). Since $\sigma$ depends on barrier height $\Phi_{S}$ at the
interface Al/GaAs as $\sigma\propto\exp(-2(\Phi_{S}+E_{F})/\mathcal{E}_{00})$,
the observed behaviour of $\sigma(t)$ shows negligible variation of $\Phi_{S}$
due to illumination. Therefore, the ”thickening” of the 2D subbands of the
$\delta$-layer at PTPC regime can be explained only by the broadening of the
potential well profile to the bulk GaAs direction (dashed potential profile in
insert to Figure 2). It can be associated with the decrease of the charged
acceptors concentration $N_{a}$ in the bulk GaAs caused by the photoexiting
electrons captured on acceptor levels beyond the $\delta$-doped region. A
self-consistent calculation of the dependence $E_{i}(N_{a})$ qualitatively
proves this assumption (Figure 6).

Figure 6:  Theoretical calculation of the 2D-subband positions at different
acceptor concentrations at $T=4.2$ $K$ and $U=0$ (other calculation parameters
are the same as at the caption to the Figure 2).

The calculations were made at constant concentration of Si in $\delta$-layer
and constant bias $U=0$ at tunnel junction. However, the behaviour of the
theoretical dependances $E_{i}(N_{a})$ reflects main features of the
experimental curves $E_{i}(t)$: subbands are ”thickening” to $E_{0}$, and
value $|E_{0}|$ is increasing by several $meV$, i.e. 2DEG concentration of the
$\delta$-layer is growing due to reduction of the compensative charge on the
acceptor states in bulk GaAs.

Note that Figure 3 shows that resolution of tunnelling spectroscopy at $T=77$
$K$ is still high enough to observe 2D subbands in TS. Accurate measurements
of the TS in PTPC regime at 77 $K$ demonstrate the absence of the effect of
”thickening” of the 2D subbands: 2D levels shift to the conduction band
direction approximately by the same value $\sim 1$ $meV$ (Figure 7). This
behaviour can be associated with increase of active impurities concentration
in $\delta$-layer (because of DX-center photoionization in $\delta$-layer
area, for instance). The model self-consistent calculations of the behaviour
of $E_{i}$ as the function of the impurity concentration $N_{Si}$ in
$\delta$-layer ($N_{a}=const$) demonstrate almost uniform diminution of the
subband energies. As a result the growth of $N_{Si}$ on $\simeq 2\%$ leads to
level shifts by $\simeq 1$ $meV$. And the intersubband energies practically do
not change as well as experimental curves $E_{i}(t)$ at $T=77$ $K$ (Figure 7).

Figure 7:  Persistent tunnelling photoconductivity effect at $T=77$ $K$. Error
bars show uncertainty of measurement, solid lines are polynomial
interpolations of the data. Sample d9a11.

Calculations demonstrate that this behaviour can be explained by variation of
donor concentrations in $\delta$-layer. Thus, according to tunnel spectroscopy
data at 77 $K$ at PTPC regime the observed growth of 2DEG concentration is due
to increase of active impurities number in $\delta$-layer only. At liquid
helium temperatures, on the contrary, the main effect of the PTPC is the
change of occupied acceptors states beyond $\delta$-layer.

The data on the temperature dependence of $S_{0}=d\ln\sigma/dU$ at $U=0$ shows
(Figure 8) the presence of two annealing temperatures of the PTPC.

Figure 8:  $S_{0}=d\ln\sigma/dU$ at $U=0$ as a function of temperature:
cooling in darkness (thick curve), reverse warming after illumination right up
to saturation (thin curve). The insert shows the behaviour of $S_{0}$ around
120 $K$. Sample d5b42.

At $T\simeq 50$ $K$ the main changes in $S_{0}$ vanish. This apparently
corresponds to thermally activation of charged centers in bulk GaAs. The
further heating of the sample suppresses change of the active donors in
$\delta$-layer, caused by illumination (at $T\geq 120$ $K$ the curve
$S_{0}(T)$ coincides with the ”dark” curve). As $T\simeq 120$ $K$ corresponds
to the DX-center’s activation temperature in GaAs Mooney (1990), we suppose
that the later effect relates to the DX-centers in the $\delta$-layer region.
The annealing temperature $T=50$ $K$ of the PPC effect in a ”bulk”
$\delta$-layer was reported in Chen  _et al._ (1998), where it was related to
deep centers in GaAs, which are DX-like.

##  IV Conclusions

Thus, our results show that the tunnelling spectroscopy of the near-surface
delta-layer permits to obtain qualitative new information about persistent
photoconductivity effect in $\delta$-layer structures on the base of GaAs. The
self-consistent and strong asymmetric potential profile (such as delta-layer
QW) depends on both number of charged impurity in bulk GaAs and active donor
concentration in $\delta$-layer. Each of these factors have different effects
on 2D subband’s spectra in $\delta$-layer. The decrease of acceptors
concentration beyond $\delta$-layer led to ”thickening” of the 2D subbands,
whereas the increase of the active donor number in $\delta$-layer gives rise
to almost uniform shift of the dimensional quantization levels. Note that in
both above-mentioned cases the electron concentration in $\delta$-layer is
augmented. We can separate these cases through character of 2DEG spectra
transformation.

According to TS data obtained, persistent effects depend upon temperature, at
which the samples are illuminated. In 50–120 $K$ temperature range the
reorganization of the TS at PTPC regime evidence of 2D electron concentration
increase in QW due to donor (DX-center) photoionization in $\delta$-layer. At
temperatures less than 50 $K$ 2DEG concentration increases in addition and,
besides, the spreading of the QW’s potential profile of the $\delta$-layer is
observed. Additional persistent effect is associated with photoionization of
the acceptor states in bulk GaAs in this case.

Besides, the results obtained show that PTPC permits to provide the resonance
conditions for detecting polaron effects in the quantum well of the near-
surface $\delta$-layer.

##  V Acknowledgements

We are grateful to Yu.V. Fedorov and A.V. Hook for preparation of the samples,
V.A. Volkov for helpful discussions and V.G. Mokerov for attention. This work
was supported by RFBR and Russian Federal Programs PAS, FTNS.

## References

  * Kotel’nikov  _et al._ (2000) I. N. Kotel’nikov, V. A. Kokin, Y. V. Fedorov, A. V. Hook, and D. T. Talbaev, “Intersubband resonance polarons in al/$\delta$-gaas tunneling junctions,” Journal of Experimental and Theoretical Physics Letters 71, 387–390 (2000).
  * A.Zrenner (1987) A. A.Zrenner, Ph.D. thesis, Technische Universität München (1987).
  * Koenraad  _et al._ (1990) P. Koenraad, A. Voncken, J. Singleton, F. Blom, C. Langerak, M. Leys, J. Perenboom, S. Spermon, W. van Vleuten, and J. Wolter, “Characterization of silicon $\delta$-doped gaas grown by mbe at various temperatures,” Surface Science 228, 538 – 541 (1990).
  * Arscott, Missous, and Dobaczewski (1992) S. Arscott, M. Missous, and L. Dobaczewski, “Observation of persistent photoconductivity in delta-doped GaAs,” Semiconductor Science and Technology 7, 620–623 (1992).
  * de Oliveira  _et al._ (1995) A. G. de Oliveira, G. M. Ribeiro, D. A. W. Soares, J. A. Corrêa F., M. I. N. da Silva, and H. Chacham, “Photoconductivity and n- to p-type transition in silicon planar-doped gaas structures with a doped cap layer,” Journal of Applied Physics 78, 2659–2665 (1995), https://doi.org/10.1063/1.360127 .
  * Li and Jagadish (1997) G. Li and C. Jagadish, “Effect of illumination on the subband electronic structure of si $\delta$-doped gaas,” Applied Physics Letters 70, 90–92 (1997), https://doi.org/10.1063/1.119317 .
  * Valyaev  _et al._ (1998) V. V. Valyaev, V. L. Gurtovoi, D. Y. Ivanov, S. V. Morozov, V. V. Sirotkin, Y. V. Dubrovskii, S. Y. Shapoval, Y. N. Khanin, E. E. Vdovin, and A. N. Pustovit, “Conductivity and persistent photoconductivity in gaas epitaxial films containing single and double delta-doped layers,” Journal of Experimental and Theoretical Physics 86, 383–387 (1998).
  * Chen  _et al._ (1998) C. Y. Chen, T. Thio, K. L. Wang, K. W. Alt, and P. C. Sharma, “Persistent photoconductivity in si delta-doped gaas at low doping concentration,” Applied Physics Letters 73, 3235–3237 (1998), https://doi.org/10.1063/1.122729 .
  * Piotrzkowski  _et al._ (1999) R. Piotrzkowski, E. Litwin-Staszewska, F. Bosc, J. Sicart, and J. Robert, “New type of persistent photoconductivity related to dx-center: the study of interband ppc in si-doped algaas,” Physica B: Condensed Matter 273-274, 792 – 795 (1999).
  * 7 (2001) _Proc. of All-Russian Workshop ”Nanophotonics”_ (IPM RAS, N.Novgorod, 2001) in Russian.
  * Mooney (1990) P. M. Mooney, “Deep donor levels (dx centers) in iii-v semiconductors,” Journal of Applied Physics 67, R1–R26 (1990), https://doi.org/10.1063/1.345628 .


# Exact and approximate solutions to the minimum of $1+x+\cdots+x^{2n}$

Aaron Hendrickson
<EMAIL_ADDRESS>

Claude F. Leibovici
<EMAIL_ADDRESS>
###### Abstract

The polynomial $f_{2n}(x)=1+x+\cdots+x^{2n}$ and its minimizer on the real
line $x_{2n}=\operatorname*{arg\,inf}f_{2n}(x)$ for $n\in\mathbb{N}$ are
studied. Results show that $x_{2n}$ exists, is unique, corresponds to
$\partial_{x}f_{2n}(x)=0$, and resides on the interval $[-1,-1/2]$ for all
$n$. It is further shown that $\inf f_{2n}(x)=(1+2n)/(1+2n(1-x_{2n}))$ and
$\inf f_{2n}(x)\in[1/2,3/4]$ for all $n$ with an exact solution for $x_{2n}$
given in the form of a finite sum of hypergeometric functions of unity
argument. Perturbation theory is applied to generate rapidly converging and
asymptotically exact approximations to $x_{2n}$. Numerical studies are carried
out to show how many terms of the perturbation expansion for $x_{2n}$ are
needed to obtain suitably accurate approximations to the exact value.

Keywords: Lagrange inversion, Tschirnhaus transformation, perturbation series,
transcendental roots, hypergeometric function.
2010 Mathematics Subject Classification: Primary 65H04, Secondary 65Q30.

##  1 Introduction

The inspiration for this work came from a question posted by Wang [11] on the
_Mathematics Stack Exchange_ discussion board March 13, 2021, which sought a
solution to the minimum of the polynomial $1+x+\cdots+x^{2n}$ for
$n\in\mathbb{N}$. In the question is was noted that the minimum appeared to
correspond to a vanishing derivative and thus could be found by solving for
the real roots of $\partial_{x}(1+x+\cdots x^{2n})$. When $n=1,2$ these roots
are algebraic with their exact forms being recovered using the standard
formulae for linear and cubic equations. However for $n\geq 3$, the work of
Abel and Galois shows no general algebraic solution exists; hence, motivating
the need for more powerful methods [1]. Given the broad and pervasive
applications of geometric sums in the literature, further study of this
polynomial and its minimum is a worthwhile venture.

##  2 Preliminaries

Throughout this work we define $\mathbb{N}=\\{1,2,\dots\\}$,
$\mathbb{N}_{0}=\mathbb{N}\cup\\{0\\}$, and $\mathbb{E}=\\{2,4,\dots\\}$ to be
the sets of positive integers, nonnegative integers, and positive even
integers, respectively. For the sake of brevity we shall denote $m=2n$ so that
the the polynomial of interest and its minimizer becomes

$f_{m}(x)\coloneqq 1+x+\cdots+x^{m},\quad m\in\mathbb{E}$

and

$x_{m}\coloneqq\operatorname*{arg\,inf}_{x\in\mathbb{R}}f_{m}(x).$

The following definitions and relations will also be used.

######  Definition 1 (Gamma function).

$\Gamma(s)\coloneqq\int_{0}^{\infty}t^{s-1}e^{-t}\,\mathrm{d}t,\quad\Re s>0$

######  Definition 2 (Factorial power (falling factorial)).

$(s)^{(n)}\coloneqq\frac{\Gamma(s+1)}{\Gamma(s-n+1)}$

######  Definition 3 (Pochhammer symbol (rising factorial)).

$(s)_{n}\coloneqq\frac{\Gamma(s+n)}{\Gamma(s)}$

######  Definition 4 (Generalized hypergeometric series).

${{}_{p}F_{q}}\left({a_{1},\dots,a_{p}\atop
b_{1},\dots,b_{q}};z\right)\coloneqq\sum_{k=0}^{\infty}\frac{(a_{1})_{k}\cdots(a_{p})_{k}}{(b_{1})_{k}\cdots(b_{q})_{k}}\frac{z^{k}}{k!}$

######  Definition 5 ($k$-gamma function and Pochhammer $k$-symbol [2]).

The $k$-gamma function and Pochhammer $k$-symbol are given by

$\Gamma_{k}(x)\coloneqq k^{\frac{x}{k}-1}\Gamma\left(\frac{x}{k}\right)$

and

$(x)_{n,k}\coloneqq\frac{\Gamma_{k}(x+nk)}{\Gamma_{k}(x)},$

respectively.

######  Relation 1.

If $n\in\mathbb{N}_{0}$ then $(\alpha)_{n,k}=k^{n}(\alpha/k)_{n}$ [7, Prop.
3.1].

With these definitions at hand we are ready to begin studying the properties
of $f_{m}$ and $x_{m}$.

##  3 Properties of $f_{m}$ and $x_{m}$

Our first goal is to establish the existence and uniqueness of $x_{m}$. To
accomplish this it will be helpful to use the closed-form for geometric sums
and write $f_{m}$ in the form

$f_{m}(x)=\frac{1-x^{m+1}}{1-x}.$ (1)

######  Lemma 1.

The polynomial $f_{m}(x)$ is strictly convex on $x\in\mathbb{R}$ for all
$m\in\mathbb{E}$.

###### Proof.

To establish strict convexity it is sufficient to show
$f_{m}^{\prime\prime}(x)>0$ everywhere on $x\in\mathbb{R}$. It is trivial to
show $f_{m}^{\prime\prime}(x)>0$ holds for $x\geq 0$ so all that is left to do
is to consider the complementary case $x<0$. Equating the second derivative
with zero we find $f_{m}^{\prime\prime}(x)=0\iff h_{m}(x)=0$, where

$h_{m}(x)=(m-1)mx^{m+1}-2(m^{2}-1)x^{m}+m(m+1)x^{m-1}-2.$

The signs of the coefficients of $h_{m}(-x)$ in order of descending variable
exponent yields the sequence $(-1,-1,-1,-1)$, which are all negative. It
follows from Descartes’ rule of signs that $f_{m}^{\prime\prime}(x)$ has zero
roots on the interval $x\in(-\infty,0)$. But,
$f_{m}^{\prime\prime}(-1)=\frac{1}{2}m^{2}>0$; thus, we conclude
$f_{m}^{\prime\prime}(x)>0$ also holds for all $x<0$. The proof is now
complete. ∎

######  Theorem 1.

The minimizer $x_{m}$ exists, is unique, and resides on the interval
$[-1,-1/2]$ for all $m\in\mathbb{E}$, .

###### Proof.

We begin by establishing that $f_{m}^{\prime}(x)$ has exactly one real root.
It is immediately obvious that $f_{m}^{\prime}(x)>0$ for all $x\geq 0$. Now
assuming $x<0$, we deduce $f_{m}^{\prime}(x)=0\iff g_{m}(x)=0$, where
$g_{m}(x)=mx^{m+1}-(m+1)x^{m}+1$. The signs of the coefficients of $g_{m}(-x)$
in order of descending variable exponent gives the sequence $(-1,-1,+1)$;
revealing a single variation in sign. Again appealing to Descartes’ rule of
signs we conclude $f_{m}^{\prime}(x)$ must have exactly one real root on the
interval $x\in(-\infty,0)$. However, $f_{m}^{\prime}(-1)=-\frac{1}{2}m$ and
$f_{m}^{\prime}(-1/2)=\frac{1}{9}2^{1-m}(2^{m+1}-3m-2)\geq 0$ with the latter
inequality following from induction on $m\in\mathbb{E}$. Consequently,
$f_{m}^{\prime}(x)$ has a single root on the real line contained in the
interval $[-1,-1/2]$ for all $m\in\mathbb{E}$. Furthermore, the strict
convexity of $f_{m}$ proven in Lemma 1 implies that the solution to
$f_{m}^{\prime}(x)=0$ also corresponds to the unique global minimum of
$f_{m}$, which completes the proof. ∎

With the existence and uniqueness of $x_{m}$ established, we turn to finding a
simple formula for the minimum of $f_{m}$ as a function of $x_{m}$.

######  Lemma 2.

Let $x_{m}\in[-1,-1/2]$ denote the unique minimizer of $f_{m}$ such that
$f_{m}(x_{m})=\inf_{x\in\mathbb{R}}f_{m}(x)$. Then,

$f_{m}(x_{m})=\frac{1+m}{1+m(1-x_{m})}$

and $f_{m}(x_{m})\in[1/2,3/4]$ for all $m\in\mathbb{E}$ with
$f_{2}(x_{2})=3/4$ and $\lim_{m\to\infty}f_{m}(x_{m})=1/2$.

###### Proof.

From Theorem 1 we know that $x_{m}$ satisfies
$mx_{m}^{m+1}-(m+1)x_{m}^{m}+1=0$, which can be rewritten as
$x_{m}^{m+1}=x_{m}/(1+m(1-x_{m}))$. Substituting this expression for
$x_{m}^{m+1}$ into (1) yields the desired form for $f_{m}(x_{m})$. The bounds
on $f_{m}(x_{m})$ are then found by minimizing and maximizing
$f(m,x)=(1+m)/(1+m(1-x))$ on $(m,x)\in\mathbb{E}\times[-1,-1/2]$. We find
$\inf f(m,x)=\lim_{m\to\infty}f(m,-1)=1/2$ and $\sup f(m,x)=f(2,-1/2)=3/4$,
which are equivalent to $\lim_{m\to\infty}f_{m}(x_{m})$ and $f_{2}(x_{2})$,
respectively. The proof is now complete. ∎

##  4 Explicit expression for $x_{m}$

In the previous section we showed that the minimizer $x_{m}$ exists, is
unique, and resides on the interval $[-1,-1/2]$ for all $m\in\mathbb{E}$.
Furthermore, we were able to establish a very simple expression for $\inf
f_{m}$ as a function of this minimizer so that the problem of evaluating $\inf
f_{m}$ is equivalent to finding $x_{m}$. For $m=2,4$ we may apply the standard
equations for roots of linear and cubic equations to derive exact algebraic
expressions for $x_{m}$. Furthermore, as $m\to\infty$ we find for $|x|<1$:
$f_{m}(x)\to(1-x)^{-1}$, which is strictly increasing on $x\in(-1,-1/2]$.
Bringing these observations together we have

$\displaystyle x_{2}$ $\displaystyle=-\frac{1}{2}$ $\displaystyle x_{4}$
$\displaystyle=-\frac{1}{4}\left(1+\sqrt[3]{5/9}\left(\sqrt[3]{9+4\sqrt{6}}-\sqrt[3]{4\sqrt{6}-9}\right)\right)$
$\displaystyle\ \,\vdots$ $\displaystyle x_{\infty}$ $\displaystyle=-1.$

While a general algebraic solution for $x_{m}$ with $m\geq 6$ does not exist,
methods for expressing exact solutions to higher-order polynomial roots have
been thoroughly studied [10]. For example, the work of Hermite shows that
$x_{6}$ can be solved exactly in terms of nonelementary functions [5]. One way
this is accomplished is by reducing the quintic equation
$\partial_{x}f_{6}(x)=0$ to its Bring–Jerrard normal form and then using
series reversion to express $x_{6}$ in terms of hypergeometric functions.
Using this approach as a clue, Theorem 2 presents an exact and general
solution for $x_{m}$ based on an adaptation of the method used for solving
trinomial equations [4, 3, 8].

######  Theorem 2.

For all $m\in\mathbb{E}$

$x_{m}=\sum_{k=1}^{m}\frac{(-m)^{k-2}}{(1+m)^{\frac{mk+k}{m}-1}}\frac{\Gamma\left(\frac{mk+k}{m}-1\right)}{\Gamma\left(\frac{m+k}{m}\right)\Gamma(k)}{{}_{m+2}F_{m+1}}\left({1,\\{\frac{k}{m}+\frac{\ell-1}{m+1}\\}_{\ell=0}^{m}\atop\frac{m+k}{m},\\{\frac{k+\ell}{m}\\}_{\ell=0}^{m-1}};1\right).$

###### Proof.

From Theorem 1 we know $x_{m}$ satisfies

$mx_{m}^{m+1}-(m+1)x_{m}^{m}+1=0,\quad m\in\mathbb{E}.$

Performing the substitution $x_{m}\mapsto-\zeta^{-\frac{1}{m}}$ we obtain the
transformed expression

$\zeta=1+m+m\phi(\zeta),$ (2)

with $\phi(\zeta)=\zeta^{-\frac{1}{m}}$. By Lagrange’s inversion theorem it
follows for a function $F$ analytic in a neighborhood of the root of (2) that

$F(\zeta)=F(1+m)+\sum_{n=1}^{\infty}\frac{m^{n}}{n!}\left[\partial_{w}^{n-1}F^{\prime}(w)\phi^{n}(w)\right]_{w=1+m}.$

Choosing $F(\zeta)=-\zeta^{-\frac{1}{m}}$ we subsequently obtain

$x_{m}=-(1+m)^{-\frac{1}{m}}+\sum_{n=1}^{\infty}\frac{m^{n-1}}{n!}\left[\partial_{w}^{n-1}w^{-\frac{m+n+1}{m}}\right]_{w=1+m},$

which upon further noting that
$\partial_{w}^{n-1}w^{-s}=(-1)^{n-1}(s)_{n-1}w^{-s-n+1}$ yields after some
algebraic manipulation

$x_{m}=-\frac{(1+m)^{-\frac{1}{m}}}{m}\sum_{n=0}^{\infty}\frac{\Gamma\left(\frac{mn+n+1}{m}\right)}{\Gamma\left(\frac{m+n+1}{m}\right)}\frac{\bigl{(}-m(1+m)^{-\frac{m+1}{m}}\bigr{)}^{n}}{n!}.$
(3)

To evaluate the series (3) we write
$x_{m}=\sum_{n=0}^{\infty}c_{n}=\sum_{k=1}^{m}\sum_{n=0}^{\infty}c_{mn+k-1}$;
resulting in $m$ new series containing Pochhammer symbols of the form
$(\cdot)_{(m+1)n}$ and $(\cdot)_{mn}$. Then using the identity [7, Eq. 2.13]

$(\alpha)_{rn}=r^{rn}\prod_{j=0}^{r-1}\left(\frac{\alpha+j}{r}\right)_{n},\quad
r\in\mathbb{N}$

we arrive at

$x_{m}=\sum_{k=1}^{m}\frac{(-m)^{k-2}}{(1+m)^{\frac{mk+k}{m}-1}}\frac{\Gamma\left(\frac{mk+k}{m}-1\right)}{\Gamma\left(\frac{m+k}{m}\right)\Gamma(k)}\sum_{n=0}^{\infty}\frac{(1)_{n}\prod_{\ell=0}^{m}\bigl{(}\frac{k}{m}+\frac{\ell-1}{m+1}\bigr{)}_{n}}{\left(\frac{m+k}{m}\right)_{n}\prod_{\ell=0}^{m-1}\left(\frac{k+\ell}{m}\right)_{n}}\frac{1}{n!},$

which is the desired result. ∎

To demonstrate the validity of the closed-form for $x_{m}$ given by Theorem 2
we substitute $m=2$ and find

$x_{2}=\frac{1}{9}{{}_{3}F_{2}}\left({1,\frac{2}{3},\frac{4}{3}\atop
2,\frac{3}{2}};1\right)-\frac{1}{\sqrt{3}}{{}_{2}F_{1}}\left({\frac{1}{6},\frac{5}{6}\atop\frac{3}{2}};1\right).$

The ${}_{3}F_{2}(\cdot)$ term is reduced to ${}_{2}F_{1}(\cdot)$ by way of [6,
Eq. 07.27.03.0120.01]

${{}_{3}F_{2}}\left({1,\beta,\gamma\atop
2,\epsilon};z\right)=\frac{\epsilon-1}{(\beta-1)(\gamma-1)z}\left({{}_{2}F_{1}}\left({\beta-1,\gamma-1\atop\epsilon-1};z\right)-1\right).$

Gauss’s hypergeometric summation theorem

${{}_{2}F_{1}}\left({\alpha,\beta\atop\gamma};1\right)=\frac{\Gamma(\gamma)\Gamma(\gamma-\alpha-\beta)}{\Gamma(\gamma-\alpha)\Gamma(\gamma-\beta)},\quad\Re(\gamma-\alpha-\beta)>0$

then permits us to write the remaining ${}_{2}F_{1}(\cdot)$ terms as ratios of
gamma functions. After some simplification we find

$x_{2}=-\frac{1}{2},$

which is the exact value of $x_{2}$.

For $m\geq 4$, reducing the closed-form for $x_{m}$ to more elementary
functions in this manner becomes very cumbersome if not impossible. Without
the ability to reduce the hypergeometric functions present in $x_{m}$, this
expression also becomes difficult to implement numerically, especially as $m$
becomes large. To obtain approximations we could turn to the series given by
(3); however, the slow convergence of this series renders it impractical. For
example, substituting $m=2$ and adding up the first one-hundred terms of (3)
we obtain $x_{2}\approx-0.499885$, which corresponds to an absolute relative
error of $2.3\times 10^{-4}$. Given that numerical root finding methods can
achieve more accurate approximations in just a few iterations we find this
means of approximation to be less than satisfactory.

##  5 Perturbation series expansion of $x_{m}$

In the previous section we were able to find an exact expression for $x_{m}$
but this expression was not useful for the purpose of computing numerical
approximations. Here, we apply the methods of perturbation theory to obtain a
faster converging series expansion for this purpose.

We begin by recalling from Theorem 1 that $x_{m}$ satisfies

$g_{m}(x_{m})=0,\quad\text{with}\quad
g_{m}(x)=x^{m}\left(1-x+\tfrac{1}{m}\right)-\tfrac{1}{m}$

and $x_{m}\to-1$ as $m\to\infty$. The fact that $x_{m}+1$ vanishes as $m$
becomes large suggests we instead study the perturbed problem

$g_{m,\epsilon}(x_{m,\epsilon})=0,\quad\text{with}\quad
g_{m,\epsilon}(x)=x^{m}\left(2-(1+x)\epsilon+\tfrac{1}{m}\right)-\tfrac{1}{m},$

where

$x_{m,\epsilon}=\sum_{k=0}^{\infty}a_{k}\epsilon^{k}.$ (4)

Upon inspection we observe $g_{m,1}(x)=g_{m}(x)$ and so it follows that
$x_{m}$ can be recovered by evaluating the perturbation series (4) at
$\epsilon=1$. To determine the coefficients $a_{k}$ we first consider the
well-known result for integer powers of series to express powers of
$x_{m,\epsilon}$ as

$x_{m,\epsilon}^{p}=\sum_{k=0}^{\infty}c_{k,p}\epsilon^{k},\quad
p\in\mathbb{N}$

with

$\displaystyle c_{0,p}$ $\displaystyle=a_{0}^{p}$ $\displaystyle c_{k,p}$
$\displaystyle=\frac{1}{a_{0}k}\sum_{\ell=1}^{k}((p+1)\ell-k)a_{\ell}c_{k-\ell,p}.$

Using Faá di Bruno’s formula we may also obtain a closed-form for the
coefficients $c_{k,p}$ as

$c_{k,p}=\frac{1}{k!}\sum_{\ell=1}^{k}(p)^{(\ell)}a_{0}^{p-l}B_{k,\ell}(1!\,a_{1},\dots,(k-\ell+1)!\,a_{k-\ell+1}),$

where $B_{n,k}(x_{1},\dots,x_{n-k+1})$ is the partial Bell-polynomial. Using
these results we substitute $x_{m,\epsilon}$ into $g_{m,\epsilon}$ and collect
terms by powers of $\epsilon$ yielding

$g_{m,\epsilon}(x_{m,\epsilon})=\left(2+\tfrac{1}{m}\right)a_{0}^{m}-\tfrac{1}{m}+\sum_{k=1}^{\infty}\left[\left(2+\tfrac{1}{m}\right)c_{k,m}-c_{k-1,m}-c_{k-1,m+1}\right]\epsilon^{k}.$

Since $g_{m,\epsilon}(x_{m,\epsilon})=0$ we equate the coefficients of
$\epsilon^{k}$ with zero to yield an infinite system of equations that recover
the coefficients $a_{k}$. Setting the constant term equal to zero gives
$a_{0}^{m}=(1+2m)^{-1}$. Knowing that $x_{m}\in[-1,-1/2]$ and $m\in\mathbb{E}$
we take the negative solution to this equation and set the higher-order
coefficients of $\epsilon^{k}$ to zero yielding

$a_{0}=-(1+2m)^{-\frac{1}{m}},\quad(1+2m)c_{k,m}-m(c_{k-1,m}+c_{k-1,m+1})=0.$
(5)

Evaluating the first several coefficients we are able to conjecture a closed-
form for $a_{k}$, which leads to the following result.

######  Theorem 3.

For all $m\in\mathbb{E}$

$x_{m}=\sum_{k=0}^{\infty}a_{k},$

where

$\displaystyle a_{0}$ $\displaystyle=-(1+2m)^{-\frac{1}{m}}$ $\displaystyle
a_{k}$
$\displaystyle=\sum_{\ell=0}^{k}\frac{(\ell+m+1)_{k-1,m}}{\ell!(k-\ell)!}a_{0}^{mk+\ell+1}.$

###### Proof.

We begin by considering the closed-form for $x_{m}$ claimed in the statement
of Theorem 2, which consists of a sum of $m$ hypergeometric functions.
Denoting $\\{a_{j,k}\\}_{j=1}^{m+2}$ as the top parameters and
$\\{b_{j,k}\\}_{j=1}^{m+1}$ as the bottom parameters of the hypergeometric
function in the $k$th term we find
$\gamma_{k}=(b_{1}+\cdots+b_{m+1})-(a_{1}+\dots+a_{m+2})=\frac{1}{2}$ for all
$k=1,\dots,m$. Since $\gamma_{k}>0$, each of the $m$-terms of $x_{m}$ can be
written as an absolutely convergent series; hence, the entire expression
representing $x_{m}$ must also be absolutely convergent. Now using the
conjectured closed-form for $a_{k}$ we write

$x_{m}=\sum_{k=0}^{\infty}\sum_{\ell=0}^{k}\frac{m^{k-1}}{\ell!(k-\ell)!}\frac{\Gamma\left(k+\frac{\ell+1}{m}\right)}{\Gamma\left(1+\frac{\ell+1}{m}\right)}a_{0}^{mk+\ell+1}.$

If this expression is equal to that given in the statement of Theorem 2, then
it is also absolutely convergent and permits rearrangement of its terms.
Interchanging the order of summation we find after some simplification

$x_{m}=\frac{a_{0}}{m}\sum_{\ell=0}^{\infty}\frac{\Gamma\left(\frac{m\ell+\ell+1}{m}\right)}{\Gamma\left(1+\frac{\ell+1}{m}\right)}\frac{(ma_{0}^{m+1})^{\ell}}{\ell!}\sum_{k=0}^{\infty}\left(\tfrac{m\ell+\ell+1}{m}\right)_{k}\frac{(ma_{0}^{m})^{k}}{k!}.$

The interior sum over $k$ can now be evaluated in terms of
${{}_{1}F}_{0}(\alpha;-;z)=(1-z)^{-\alpha}$. Reintroducing $a_{0}$ yields

$x_{m}=-\frac{(1+m)^{-\frac{1}{m}}}{m}\sum_{\ell=0}^{\infty}\frac{\Gamma\left(\frac{m\ell+\ell+1}{m}\right)}{\Gamma\left(\frac{m+\ell+1}{m}\right)}\frac{\bigl{(}-m(1+m)^{-\frac{m+1}{m}}\bigr{)}^{\ell}}{\ell!},$

which is the series expansion for $x_{m}$ given in (3). By the uniqueness of
Taylor series it follows that the conjectured form for $a_{k}$ must be
correct. ∎

So does the perturbation series for $x_{m}$ converge faster than that given by
(3)? Substituting $m=2$ and adding the first one-hundred terms we find for the
absolute relative error $5.6\times 10^{-64}$, which is a significant
improvement on the absolute relative error of $2.3\times 10^{-4}$ obtained
from the first one hundred terms of (3).

We conclude this section with an important property of approximations for
$x_{m}$ obtained via the perturbation series of Theorem 3.

######  Corollary 1.

If $\tilde{x}_{m,n}=\sum_{k=0}^{n}a_{k}$, then $x_{m}\sim\tilde{x}_{m,n}$ as
$m\to\infty$.

###### Proof.

Using the expression for $a_{k}$ given in Theorem 3 we have
$\lim_{m\to\infty}a_{0}=-1$ and $\lim_{m\to\infty}a_{k}=0$ for all $k\geq 1$;
thus, $\lim_{m\to\infty}\tilde{x}_{m,n}=-1$ for all $n\in\mathbb{N}_{0}$.
Since $\lim_{m\to\infty}x_{m}=-1$ the result follows. ∎

##  6 Numerical results

From Corollary 1 we know $x_{m}\sim\tilde{x}_{m,n}$ as $m\to\infty$ and so we
expect the number $n$ needed to guarantee $|x_{m}-\tilde{x}_{m,n}|<\epsilon$
should decrease as $m$ increases. Since we have closed-forms for $x_{2}$ and
$x_{4}$, which can be computed to arbitrary precision, our first task will be
to study the convergence of $\tilde{x}_{m,n}\to x_{m}$ as a function of $n$
for $m=2,4$. Given that we expect less terms will be needed for larger values
of $m$, the results of this exercise should give us a worst case scenario for
how large $n$ must be to obtain the desired accuracy in our approximation.

Using Mathematica software, we evaluated $\tilde{x}_{m,n}$ for $m=2,4$ and
$n=0,1,\dots,100$. To compare the approximation to the exact values we used
the absolute relative error

$R_{m}(n)=\left\lvert\frac{\tilde{x}_{m,n}}{x_{m}}-1\right\rvert,$

the results of which are plotted in Figure 1. From the figure we see the error
decreases exponentially with $n$ and that $R_{4}(n)<R_{2}(n)$ for each value
of $n$. Working with the data for $R_{4}(n)$ we further determined

$R_{4}(n)<5\times 10^{-(2+0.759n)},$

which suggests setting $n$ equal to

$n^{\ast}=\max\left\\{0,\left\lceil\tfrac{q-2}{0.759}\right\rceil\right\\}$

is sufficient to guarantee $\tilde{x}_{m,n^{\ast}}$ agrees with $x_{m}$ to at
least $q$ significant digits for all $m\geq 4$.

To test this hypothesis we first note that $x_{m}\in(-1,-1/2]$ for all finite
$m\in\mathbb{E}$. Since the leading exponent in the decimal expansion of
$x_{m}$ is always negative one it follows that $\tilde{x}_{m,n}$ has $p$
significant digits of $x_{m}$ if $|x_{m}-\tilde{x}_{m,n}|\leq 5\times
10^{-(p+1)}$. Furthermore, we know $x_{m}$ is the unique real root of
$g_{m}(x)=x^{m}(1-x+\frac{1}{m})-\frac{1}{m}$ with $g_{m}(x_{m}-\epsilon)$ and
$g_{m}(x_{m}+\epsilon)$ differing in sign; hence a lower bound on the number
of significant digits obtained by $\tilde{x}_{m,n}$ is found by determining
the largest nonnegative integer $p$ such that

$g_{m}(\tilde{x}_{m,n}-5\times 10^{-(p+1)})g_{m}(\tilde{x}_{m,n}+5\times
10^{-(p+1)})\leq 0.$

For the sake of example we chose $q=10$ for the number of desired significant
digits resulting in $n^{\ast}=11$. Using the above mentioned procedure, the
value $p$ was computed for $m=4,6,\dots,100$ with the results presented in
Figure 2. From the figure we observe $p>q$ for each value of $m$ as is
expected. Finally, Table 1 presents numerical values for $x_{m}$ and
$f_{m}(x_{m})$ computed using $\tilde{x}_{m,n}$ and the formula in Lemma 2.

Figure 1: Absolute relative error incurred from the approximation
$\tilde{x}_{m,n}$ versus $n$ for $m=2,4$. Plot produced with matlab2tikz [9].
Figure 2: Number of significant figures obtained by the approximation
$\tilde{x}_{m,n^{\ast}}(q)$ for $q=10$ versus $m$. Table 1: Numerical values
for $x_{m}$ and $f_{m}(x_{m})$.

$m$ | $x_{m}$ | $f_{m}(x_{m})$
---|---|---
2 | -0.5000000000 | 0.7500000000
4 | -0.6058295862 | 0.6735532235
6 | -0.6703320476 | 0.6350938940
8 | -0.7145377272 | 0.6115666906
10 | -0.7470540749 | 0.5955429324
12 | -0.7721416355 | 0.5838576922
14 | -0.7921778546 | 0.5749221276
16 | -0.8086048979 | 0.5678463037
18 | -0.8223534102 | 0.5620909079
20 | -0.8340533676 | 0.5573090540
22 | -0.8441478047 | 0.5532669587
24 | -0.8529581644 | 0.5498010211
26 | -0.8607238146 | 0.5467931483
28 | -0.8676269763 | 0.5441558518
30 | -0.8738090154 | 0.5418228660
32 | -0.8793814184 | 0.5397430347
34 | -0.8844333818 | 0.5378762052
36 | -0.8890371830 | 0.5361903986
38 | -0.8932520563 | 0.5346598151
40 | -0.8971270425 | 0.5332633990
42 | -0.9007031162 | 0.5319837878
44 | -0.9040147981 | 0.5308065300
46 | -0.9070913919 | 0.5297194951
48 | -0.9099579456 | 0.5287124219
50 | -0.9126360054 | 0.5277765690
52 | -0.9151442141 | 0.5269044410
54 | -0.9174987898 | 0.5260895727
56 | -0.9197139122 | 0.5253263565
58 | -0.9218020367 | 0.5246099035
60 | -0.9237741513 | 0.5239359311
62 | -0.9256399895 | 0.5233006711
64 | -0.9274082062 | 0.5227007942
66 | -0.9290865244 | 0.5221333471
68 | -0.9306818591 | 0.5215957008
70 | -0.9322004214 | 0.5210855067
72 | -0.9336478067 | 0.5206006599
74 | -0.9350290699 | 0.5201392683
76 | -0.9363487901 | 0.5196996259

$m$ | $x_{m}$ | $f_{m}(x_{m})$
---|---|---
78 | -0.9376111258 | 0.5192801905
80 | -0.9388198625 | 0.5188795643
82 | -0.9399784542 | 0.5184964771
84 | -0.9410900592 | 0.5181297723
86 | -0.9421575717 | 0.5177783938
88 | -0.9431836485 | 0.5174413759
90 | -0.9441707340 | 0.5171178332
92 | -0.9451210804 | 0.5168069528
94 | -0.9460367670 | 0.5165079864
96 | -0.9469197164 | 0.5162202447
98 | -0.9477717091 | 0.5159430910
100 | -0.9485943966 | 0.5156759367
102 | -0.9493893132 | 0.5154182363
104 | -0.9501578860 | 0.5151694840
106 | -0.9509014444 | 0.5149292100
108 | -0.9516212282 | 0.5146969770
110 | -0.9523183955 | 0.5144723780
112 | -0.9529940289 | 0.5142550329
114 | -0.9536491420 | 0.5140445872
116 | -0.9542846846 | 0.5138407092
118 | -0.9549015479 | 0.5136430885
120 | -0.9555005690 | 0.5134514340
122 | -0.9560825347 | 0.5132654729
124 | -0.9566481855 | 0.5130849489
126 | -0.9571982191 | 0.5129096209
128 | -0.9577332933 | 0.5127392623
130 | -0.9582540286 | 0.5125736594
132 | -0.9587610115 | 0.5124126108
134 | -0.9592547961 | 0.5122559267
136 | -0.9597359069 | 0.5121034274
138 | -0.9602048403 | 0.5119549436
140 | -0.9606620669 | 0.5118103146
142 | -0.9611080328 | 0.5116693885
144 | -0.9615431615 | 0.5115320215
146 | -0.9619678551 | 0.5113980772
148 | -0.9623824957 | 0.5112674259
150 | -0.9627874469 | 0.5111399447
$\infty$ | -1.0000000000 | 0.5000000000

##  7 Conclusions

In this note, we were able to establish many useful facts about the polynomial
$f_{m}(x)=1+x+\cdots+x^{m}$ and its minimum value on the real line. In
particular, we were able to show
$\operatorname*{arg\,inf}f_{m}(x)\in[-1,-1/2]$ and $\inf f_{m}(x)\in[1/2,3/4]$
for all $m\in\mathbb{E}$ as well as provide a very simple formula for the
minimum as a function of the minimizer $x_{m}$. Lagrange inversion and
perturbation theory were applied to derive two different series expansions for
$x_{m}$, which lead to a closed-form in terms of hypergeometric functions.
Furthermore, numerical studies were conducted which gave a rule of thumb for
how large $n$ must be to achieve a desired accuracy in approximating $x_{m}$
with $\tilde{x}_{m,n}$.

## References

  * [1] N. H. Abel.  Beweis der Unmöglichkeit, algebraische Gleichungen von höheren Graden als dem vierten allgemein aufzulösen.  J. Reine Angew. Math., 1:65–84, 1826.
  * [2] Rafael Díaz and Eddy Pariguan.  On hypergeometric functions and Pochhammer $k$-symbol.  Divulgaciones Matemáticas, 15(2), 2007.
  * [3] Albert Eagle.  Series for all the roots of a trinomial equation.  The American Mathematical Monthly, 46(7):422–425, 1939.
  * [4] M.L. Glasser.  Hypergeometric functions and the trinomial equation.  Journal of Computational and Applied Mathematics, 118(1):169–173, 2000.
  * [5] Charles Hermite.  Sur la résolution de l’équation du cinquieme degré.  Comptes rendus de l’Académie des Sciences, 46(1858):508–515, 1858.
  * [6] Wolfram Research Inc.  The Wolfram functions site.  Visited on 2021-04-28.
  * [7] Shahid Mubeen and Abdur Rehman.  A note on $k$-gamma function and pochhammer $k$-symbol.  Journal of Informatics and Mathematical Sciences, 6(2), 2014.
  * [8] Daniele Ritelli and Giulia Spaletta.  Trinomial equation: the hypergeometric way.  Open Journal of Mathematical Sciences, 5(1):236–247, 2021.
  * [9] Nico Schlömer.  matlab2tikz: A script to convert MATLAB/Octave into TikZ figures for easy and consistent inclusion into LaTeX.  GitHub.  URL: https://github.com/matlab2tikz/matlab2tikz (retrieved May 8, 2021).
  * [10] Hiroshi Umemura.  Resolution of algebraic equations by theta constants, pages 261–270.  Birkhäuser Boston, Boston, MA, 2007.
  * [11] Yiwei Wang.  The minimum of $f(x)=1+x+\cdots+x^{2n}$.  Mathematics Stack Exchange.  URL: https://math.stackexchange.com/q/4060608 (version: 2021-04-30).


1

# Algorithms for constrained optimal transport

Martin Corless and Anthony Quinn and Sarah Boufelja Y. and Robert Shorten.
Affiliations: Martin Corless (corless@purdue.edu), is with the School of
Astronautics and Aeronautics, Purdue University, West Lafayette. The other
authors are with the Dyson School of Design Engineering, Imperial College
London. Anthony Quinn is also with the Dept. of Electronic & Electrical
Engineering, Trinity College Dublin.

###### Abstract

We derive iterative scaling algorithms of the Sinkhorn-Knopp (SK) type for
constrained optimal transport. The constraints are in the form of prior-
imposed zeroes in the transport plan. Based on classical Bregman arguments, we
prove asymptotic convergence of our algorithms to a unique optimal solution.
New insights obtained from the convergence proof are highlighted. An example
from electrical vehicle charging in a smart city context is outlined, in which
the prior zero-constraints prevent energy from being transported from some
providers to some vehicles.

##  I Introduction: Optimal Transport (OT)

Optimal Transport (OT) has emerged as a key enabling mathematical technology
which is driving a growing number of contemporary engineering applications in
fields such as machine learning, image processing, and in optimization and
control [1], [2], [3]. The field is an old one, dating back to Monge in the
18th century. Nevertheless, OT is now attracting accelerated attention from
both theoreticians and practitioners alike, with widespread application in
engineering practice. OT refers to a class of problems in which we seek to
transport a finite resource (e.g. mass) from a distributed source to a
distributed target, in some optimal manner. For example, the classical
problems of Monge and Kantorovitch relate to the transportation of a
distributed pile of sand [4] from one location to another. The contemporary
analogues of these formulations (for example, in economics) arise in resource
(re)allocation problems, and in supply-demand matching problems [5]. Closely
related problems arise in key machine learning (ML) settings, such as domain
adaptation. Here, source data are optimally adapted (i.e. transported) for use
in a related, but distinct, target domain [6]. An example is “data repair”, in
which OT is used to de-bias ML tasks in order to improve their fairness
properties. Classifiers and recommender systems—trained on data with fairness
deficiencies such as representation bias, dependence on protected attributes,
etc.—can be fairness-repaired using methods of OT [7].

In this paper, we confine ourselves to the discrete case, in which the
resource (mass) is distributed across a finite (though typically large) number
of states in both the source and target. For every pair of source-target
states, the designer specifies a cost per unit mass of transporting the
resource from this source state to this target state. The OT problem then
involves the design of a transport plan that specifies the amount of resource
to be transported for every one of these source-target state pairs, in a way
that minimizes the total cost of transport. A feature of our setting of the OT
problem in this paper is that we impose constraints in which some source
states are not allowed to transport to some target states. These prior zero-
transport constraints may arise for a variety of reasons. For example, in
economics, certain producers may not be allowed to supply a product to certain
consumers. In fairness-aware domain adaptation for ML (above), data for a
particular demographic group may not be mapped to a different demographic
group. These OT problems give rise to linear programs of very large scale [8].
Linear programming solvers are known to scale poorly, and much of the current
OT work focuses on regularizing the classical OT problem to ensure
tractability at scale. The best known of these approaches is the entropic
regularisation of the OT problem [8]. In essence, this amounts to the addition
of an extra term in the (total cost) objective, which penalizes a specified
divergence of the transport plan from a prior-defined ideal plan [9]. The
resulting optimization problem is then strongly convex on a compact convex
set, and can be solved efficiently using a Sinkhorn-Knopp-type (SK) iterative
scaling algorithm, with its favourable large-scale computational properties
[10], [11]. These computational OT results and guarantees—based on the SK
algorithm— have been central to the widespread adoption of the OT paradigm in
applications, particularly in ML [4], [6], [12], [13].

In contrast to previous results in the literature, our purpose in this paper
is to consider OT problems in which some elements of the transport plan are
forced to be zero a priori, and to derive iterative scaling algorithms of the
SK-type to solve these problems. We use Bregman-type arguments [14] to prove
convergence of our algorithms to the (unique) solution of the OT problems
under consideration. We outline an application of our results in a smart
cities context, where prior zeroes in the transport plan are an important
constraint.

##  II Prior work on OT and the SK Algorithm

OT problems are often cast in the language of probability, where the resource
(mass) to be transported is a probability measure, and so the net resource
(mass) is normalized to one. The transport plan is thus a bivariate
probability mass function (pmf), and the source and target distributions are
its two marginals. We speak of a transport plan that moves the source
distribution (marginal) to the target distribution (marginal). In this way, OT
can be viewed as a mathematical framework for measuring the distance between
(source and target marginal) probability distributions, this distance being
the minimal expected cost associated with the optimal transport plan [4]. This
Kantorovich-Rubinstein (also called Wasserstein or earth-mover’s) norm [4] for
classical entropy-regularized OT gives rise to a unique solution of a strongly
convex optimization problem on a compact convex support [15], [16], with
linear convergence of the SK iterative scaling algorithm to this unique
solution [17]. If, as in this paper, we impose prior zero-constraints on some
elements of the transport plan, there may not exist a plan that satisfies
exact marginal constraints at the source and target, as already noted. Another
situation in which a transport plan satisfying exact marginal constraints does
not exist is unbalanced OT (UOT) [16], in which the prescribed net masses of
the source and target are not the same. Earlier, in [18], an SK-type iterative
scaling algorithm was derived for unconstrained UOT (i.e. without prior zero-
constraints on the transport plan), by relaxing the source and target marginal
constraints via extra terms in the cost that penalize noncompliance with the
marginal constraints. This was shown to converge linearly to a unique
solution. In [16], the authors study the UOT case with prior zero-constraints
on the transport plan. They prove that the SK algorithm has two convergent
subsequences in this case and they propose a transport plan which is the
component-wise geometric mean of the limits plans of the two convergent
subsequences

The paper is organized as follows: in Section III, our OT problem is specified
with prior constraints on the pattern of the transport plan. After developing
background material, two iterative scaling algorithms are presented in Section
IV for OT problems with prior zero-constraints. Section V presents some
properties of the resulting OT plans. Section VI provides the proof of
convergence or our algorithms. Finally, a use-case from the domain of smart
mobility is outlined in Section VII, illustrating the efficacy of our new
algorithms.

##  III The Optimal transport (OT) problem with prior zero-constraints

Consider an OT problem in which we have a finite number, $m$, of source agents
(i.e. source states), $i=1,2,\dots,m$, and a finite number, $n$, of target
agents (i.e. target states), $j=1,2,\dots,n$. Suppose each source agent, $i$,
has capacity or mass, $\tilde{u}_{i}>0$, and each target agent, $j$, desires a
capacity or mass, $\tilde{v}_{j}>0$. We wish to transport the masses of the
source agents to the target agents. For the time being, we impose mass balance
(conservation), i.e.

$\displaystyle\sum_{i=1}^{m}\tilde{u}_{i}=\sum_{j=1}^{n}\tilde{v}_{j}$ (1)

that is, ${\bm{1}}^{\prime}\tilde{u}={\bm{1}}^{\prime}\tilde{v}$ where
${\bm{1}}$ is a vector of ones, of appropriate length. Let $t_{ij}\geq 0$ be
the mass or capacity that source agent, $i$, transports to target agent, $j$.
We will refer to the $m\times n$ matrix, $T=\left\\{t_{ij}\right\\}$, as the
transport plan. As a consequence of the capacity constraints of the source and
target agents, we must have

$\displaystyle\sum_{j=1}^{n}t_{ij}$ $\displaystyle=\tilde{u}_{i},\qquad
i=1,\dots,m$ (2)

or, $T{\bm{1}}=\tilde{u}$, and

$\displaystyle\sum_{i=1}^{m}t_{ij}=\tilde{v}_{j},\qquad j=1,\dots,n$ (3)

or, $T^{\prime}{\bm{1}}=\tilde{v}$. Note that constraints (1), (2) and (3)
imply that

$\sum_{i=1}^{m}\sum_{j=1}^{n}t_{ij}=\sum_{i=1}^{m}\tilde{u}_{i}=\sum_{j=1}^{n}\tilde{v}_{j}\
$ (4)

Next, suppose that a specific source agent, $i$, never transports to a
specific target agent, $j$, so that $t_{ij}=0$. We let $\mathcal{Z}$ be the
set of these prior non-transporting index pairs:

$t_{ij}=0\text{ if }(i,j)\in\mathcal{Z}$ (5)

The set of allowable transport plans is given by

$\displaystyle\mathcal{T}\equiv\left\\{T\in\mathbb{R}^{m\times
n}:\begin{array}[]{ll}t_{ij}\geq 0&\forall(i,j)\notin\mathcal{Z}\\\
t_{ij}=0&\forall(i,j)\in\mathcal{Z}\end{array}\right\\}$ (8)

We assume that for every $i$, that there is at least one
$(i,j)\notin\mathcal{Z}$ and for every $j$, there is also at least one
$(i,j)\notin\mathcal{Z}$; i.e. if $t_{ij}>0$ for all $(i,j)\notin\mathcal{Z}$
then, $T$ has no zero-rows or zero-columns. If $\mathcal{Z}$ is empty, i.e. if
every source agent can transport to every target agent, then the marginal
constraints (2) and (3) are feasible (i.e. there is at least one solution in
$\mathcal{T}$, for example,
$T\equiv\tilde{u}\tilde{v}^{\prime}/{\bm{1}}^{\prime}\tilde{v}$). However, if
$\mathcal{Z}$ is non-empty, then (2) and (3) may not be feasible (i.e. there
may not be any solution in $\mathcal{T}$). An example is given by $m=n=2$ and
$\mathcal{Z}\equiv\\{(2,2)\\}$:

$T=\left(\begin{array}[]{cc}t_{11}&t_{12}\\\
t_{21}&0\end{array}\right)\qquad\mbox{with}\qquad\tilde{u}=\tilde{v}=\left(\begin{array}[]{c}1\\\
a\end{array}\right),\;\;a>1$

#### Transport cost.

Suppose that there is a cost, $c_{ij}\in\mathbb{R}$, of moving a unit mass
from source, $i$, to target, $j$. Then, the total cost of transport is
$\sum_{(i,j)\notin\mathcal{Z}}c_{ij}t_{ij}$ This gives rise to the constrained
OT problem:

$\displaystyle\min_{T\in\mathcal{U}(\tilde{u},\tilde{v})}\sum_{(i,j)\notin\mathcal{Z}}c_{ij}t_{ij},$
(9) $\displaystyle\mathcal{U}(\tilde{u},\tilde{v})$ $\displaystyle=$
$\displaystyle\left\\{T\in\mathcal{T}:\small{\text{constraints
\eqref{eq:rowconstraints} and \eqref{eq:colconstraints} hold}}\right\\}$ (10)

#### An ideal transport plan and regularization

Next, suppose that it is also desirable for $T$ to be close to some ideal
(desired) transport plan, $\tilde{T}\in\mathcal{T}$ (8) with
$\tilde{t}_{ij}>0$ for $(i,j)\notin\mathcal{Z}$. We therefore modify the cost
to

$\sum_{(i,j)\notin\mathcal{Z}}c_{ij}t_{ij}+\gamma_{0}KL(T|\tilde{T})$ (11)

where $\gamma_{0}>0$ is a pre-assigned regularization constant, and the
Kullback-Leibler ($KL$) divergence [9] of $T$ from $\tilde{T}$ is defined as

$\displaystyle
KL(T|\tilde{T})\equiv\sum_{(i,j)\notin\mathcal{Z}}kl(t_{ij}|\tilde{t}_{ij})$
(12)

where, for any two scalar, $t\geq 0$ and $\tilde{t}>0$, we define

$kl(t|\tilde{t})\equiv\left\\{\begin{array}[]{ccc}\tilde{t}&\text{if}&t=0\\\
\displaystyle{t\log\left(\frac{t}{\tilde{t}}\right)}-t+\tilde{t}&\text{if}&t>0\end{array}\right.$

One may readily show that, for all $t\geq 0$ and $\tilde{t}>0$,
$kl(t|\tilde{t})\geq 0$; also $kl(t|\tilde{t})=0$ iff $t=\tilde{t}$. Hence,
for all $T\in\mathcal{T}$, $KL(T|\tilde{T})\geq 0$ and $KL(T|\tilde{T})=0$ iff
$T=\tilde{T}$.

###### Remark 1

Note again that original optimization problem (9) is a linear programming
problem and may be computationally burdensome to solve. So, a common solution
approach is approximately to solve the original problem (9) by minimizing the
regularized cost in (11) with $\gamma_{0}$ small but postive. This problem can
be solved more efficiently than the original problem.

Next, note that
$c_{ij}t_{ij}=\gamma_{0}t_{ij}\log\left(\frac{1}{\exp(-c_{ij}/\gamma_{0})}\right)$
and define

$\displaystyle k_{ij}$
$\displaystyle\equiv\left\\{\begin{array}[]{ccl}\tilde{t}_{ij}\exp(-c_{ij}/\gamma_{0})&\text{if}&(i,j)\notin\mathcal{Z}\\\
0&\text{if}&(i,j)\in\mathcal{Z}\end{array}\right.$ (15)

Then, whenever $t_{ij}\neq 0$,

$\displaystyle c_{ij}t_{ij}+\gamma_{0}kl(t_{ij}|\tilde{t}_{ij})$
$\displaystyle=\gamma_{0}t_{ij}\log\left(\frac{\tilde{t}_{ij}}{k_{ij}}\right)+\gamma_{0}t_{ij}\log\left(\frac{t_{ij}}{\tilde{t}_{ij}}\right)-\gamma_{0}t_{ij}+\gamma_{0}\tilde{t}_{ij}$
$\displaystyle=\gamma_{0}t_{ij}\log\left(\frac{t_{ij}}{k_{ij}}\right)-\gamma_{0}t_{ij}+\gamma_{0}k_{ij}+\gamma_{0}\tilde{t}_{ij}-\gamma_{0}k_{ij}$
$\displaystyle=\gamma_{0}kl(t_{ij}|k_{ij})+\gamma_{0}\tilde{t}_{ij}-\gamma_{0}k_{ij}$
(16)

(16) also holds for $t_{ij}=0$ , $(i,j)\notin\mathcal{Z}$. Hence, cost (11)
can be expressed as
$\gamma_{0}KL(T|K)+\sum_{(i,j)\notin\mathcal{Z}}\gamma_{0}(\tilde{t}_{ij}-k_{ij})$.
Therefore, we consider the new optimization problem:

$\min_{T\in{\mathcal{U}}(\tilde{u},\tilde{v})}KL(T|K)$ (17)

where $\mathcal{U}(\tilde{u},\tilde{v})$ is given by (10). When the
optimization problem (17) is feasible (i.e. $\mathcal{U}(\tilde{u},\tilde{v})$
is non-empty), its solution can be obtained using the Sinkhorn-Knopp (SK)
Algorithm [17].

Sinkhorn-Knopp (SK) Algorithm. Initialize $d_{2j}(0)=1$. Iterate for
$l=0,1,\dots$

$\displaystyle d_{1i}(l+1)$
$\displaystyle=\frac{\tilde{u}_{i}}{\Sigma_{j=1}^{n}k_{ij}d_{2j}(l)}$ (18)
$\displaystyle d_{2j}(l+1)$
$\displaystyle=\frac{\tilde{v}_{j}}{\Sigma_{i=1}^{m}d_{1i}(l+1)k_{ij}}$ (19)
$\displaystyle t_{ij}(l+1)$ $\displaystyle=d_{1i}(l+1)k_{ij}d_{2j}(l+1)$ (20)

${\hskip 220.50885pt}\square$

If (17) is feasible, the sequence, $\\{T(l)\\}$ converges to a limit, $T^{*}$,
which is the unique minimizer for (17) [16]. If $t^{*}_{ij}=0$ for some
$(i,j)\notin\mathcal{Z}$, then either $\\{d_{ii}(l)\\}$ or $\\{d_{2j}(l)\\}$
do not converge. If (17) is not feasible, the sequence, $\\{T(l)\\}$ has two
convergent subsequences with different limits [16].

#### When desired marginals constraints are not achievable

If $\mathcal{Z}$ is non-empty, it may not be possible to simultaneously
satisfy the marginal constraints, (2) and (3). If—as already imposed in
(8)—$T$ has no enforced rows of zeroes, the row constraints can always be
achieved, as explained in Remark 2, below. Therefore, consider the situation
in which row constraints (2) are satisfied but column constraints (3) are not
necessarily satisfied. Instead of requiring satisfaction of the column
constraints, we add a KL-based penalty term to the cost function, which is a
measure of non-compliance with the column constraints, $\tilde{v}$. We
therefore consider a new optimization problem given by

$\min_{T\in\mathcal{V}(\tilde{u})}KL(T|K)+\gamma KL(v_{T}|\tilde{v})$ (21)

where $KL(v_{T}|\tilde{v})=\sum_{j=1}^{n}kl(v_{T_{j}}|\tilde{v}_{j})$,
$v_{T_{j}}{\equiv}\sum_{i=1}^{m}t_{ij}$, $\gamma>0$,

$\displaystyle\mathcal{V}(\tilde{u})$
$\displaystyle\equiv\left\\{T\in\mathcal{T}:\sum_{j=1}^{n}t_{ij}=\tilde{u}_{i},\;i=1,\dots,m\right\\}.$
(22)

###### Remark 2

For any admissible $\tilde{u}$, one can always obtain a $T$ in
$\mathcal{V}(\tilde{u})$ with $t_{ij}>0$ for all $(i,j)\notin\mathcal{Z}$.
Simply choose any $T\in\mathcal{T}$ with $t_{ij}>0$ for all
$(i,j)\notin\mathcal{Z}$, and scale each of its rows to yield a new matrix in
$\mathcal{V}(\tilde{u})$. Therefore, $\mathcal{V}(\tilde{u})$ is non-empty,
and since we are minimizing a continuous strictly convex function (21) over a
compact convex set (22), a unique minimizer exists.

${\hskip 220.50885pt}\square$

###### Remark 3

[19] considers a generalization of the problem considered here. To highlight
the novelty of our work we make the following comments.

  * (a)

It is true that our problem is a special case of the initial general problem
considered in [19] by, among other things, specializing to the bivariate (so,
bi-marginal) case, and taking $c_{ij}=\infty$ for all entries in the transport
plan that are to be zero.

  * (b)

However, our algorithms—which have been designed, among other things, to be
relevant in important applications of OT, such as the sharing economy
application of our Section VII—are not presented in [19].

  * (c)

More substantially, we note the following difference. In order to prove
convergence of their algorithm to solve their general problem, it is assumed
that all the elements of their cost tensor $\mathbf{C}$ are finite (Assumption
C on page 5, left column, of [19]). Under this assumption, our problem is not
a special case of the problem solved by the algorithm in [19], because for our
problem to be a special case, some of the elements of their cost tensor
$\mathbf{C}$ must be infinite. This assumption is used in the proof of Lemma
III.8 in [19]. Note also that prior zeros are not imposed on their transport
tensor, $M$ (corresponding to our transport matrix, $T$).

##  IV Main result: a modified SK algorithm

We now present the main results of the paper. These results yield iterative
scaling algorithms that produce a sequence, $\\{T(l\\}$, converging to the
optimal transport plan, $T^{*}$, for (21). In particular, two algorithms are
presented. Algorithm 1 is obtained using results in [14]. Algorithm 2 is
equivalent to Algorithm I and is presented for comparison to other related
algorithms such as the SK Algorithm. A third algorithm, the Chizat Algorithm
is a related algorithm from the literature that is included to provide context
for our contributions.

Algorithm 1. Initialize $T(0)=K,v(0)=\tilde{v}$. Iterate for $l=0,1,\dots$

$\displaystyle c_{1i}(l+1)$
$\displaystyle=\frac{\tilde{u}_{i}}{\Sigma_{j=1}^{n}t_{ij}(l)}$ (23)
$\displaystyle c_{2j}(l+1)$
$\displaystyle=\left(\frac{v_{j}(l)}{\Sigma_{i=1}^{m}c_{1i}(l+1)t_{ij}(l)}\right)^{\frac{\gamma}{1+\gamma}}$
(24) $\displaystyle t_{ij}(l+1)$
$\displaystyle=c_{1i}(l+1)t_{ij}(l)c_{2j}(l+1)$ (25) $\displaystyle
v_{j}(l+1)$ $\displaystyle=c_{2j}(l+1)^{-1/\gamma}v_{j}(l)$ (26)

${\hskip 220.50885pt}\square$

The following theorem provides the main result of this paper. A proof is
provided in Section VI.

###### Theorem 1

Consider the sequence, $\\{T(l),v(l)\\}$, generated by Algorithm 1. This
sequence converges to the (unique) limit, $(T^{*},v^{*})$, and $T^{*}$ is the
minimizer for the optimization problem given by (21).

###### Remark 4

Since

$\displaystyle\Sigma_{j=1}^{n}t^{*}_{ij}=\tilde{u}_{i}$ (27)

from (23), $\lim_{l\rightarrow\infty}c_{1i}(l)=1$ for all $i$. In Lemma 2 (see
Section V), it will be shown that $t^{*}_{ij}>0$ for all
$(i,j)\notin\mathcal{Z}$. Hence, $T^{*}$ is guaranteed not to contain any row
of zeroes, and it follows from (25) that
$\lim_{l\rightarrow\infty}c_{2j}(l)=1$ for all $j$. From (24), we see that

$\displaystyle\Sigma_{i=1}^{m}t^{*}_{ij}=v^{*}_{j}$ (28)

and, since $T^{*}$ does not contain any column of zeroes, we must have
$v^{*}_{j}>0$, $\forall j$. It now follows, from (27) and (28), that
$\sum_{j=1}^{n}v_{j}^{*}=\sum_{i=1}^{m}\tilde{u}_{i}$ and so $v^{*}$ and
$\tilde{u}$ are balanced  (i.e. mass-conserving). The above result holds even
in the unbalanced (i.e. non-mass-conserving) problem,
$\sum_{i=1}^{m}\tilde{u}_{i}\neq\sum_{j=1}^{n}\tilde{v}_{j}$. In this case,
application of the SK algorithm produces a sequence, $\\{T(l)\\}$, which has
two convergent subsequences with different limits [16].

${\hskip 220.50885pt}\square$

In order to obtain an algorithm for comparison with the SK algorithm, we
introduce scaling parameters, $d_{1i}(l)$ and $d_{2j}(l)$, defined by

$\displaystyle d_{1i}(l+1)$ $\displaystyle\equiv c_{1i}(l+1)d_{1i}(l),\qquad
d_{1i}(1)\equiv c_{i1}(1)$ (29) $\displaystyle d_{2j}(l+1)$
$\displaystyle\equiv c_{2j}(l+1)d_{2j}(l),\qquad d_{2j}(0)\equiv 1$ (30)

Then

$\displaystyle t_{ij}(l+1)$ $\displaystyle=d_{1i}(l+1)k_{ij}d_{2j}(l+1)$ (31)
$\displaystyle v_{j}(l+1)$
$\displaystyle=d_{2j}(l+1)^{-1/\gamma}\tilde{v}_{j}$ (32)

The sequence, $\\{T(l)\\}$, obtained from Algorithm 1, can then also be
obtained from the following algorithm.

Algorithm 2. Initialize $d_{2j}(0)=1.$ Iterate for $l=0,1,\dots$

$\displaystyle d_{1i}(l+1)$
$\displaystyle=\frac{\tilde{u}_{i}}{\Sigma_{j=1}^{n}k_{ij}d_{2j}(l)}$ (33)
$\displaystyle d_{2j}(l+1)$
$\displaystyle=\left(\frac{\tilde{v}_{j}}{\Sigma_{i=1}^{m}d_{1i}(l+1)k_{ij}}\right)^{\frac{\gamma}{1+\gamma}}$
(34) $\displaystyle t_{ij}(l+1)$ $\displaystyle=d_{1i}(l+1)k_{ij}d_{2j}(l+1)$
(35)

${\hskip 220.50885pt}\square$

###### Remark 5

Since $\tilde{v}_{j},v^{*}_{j}>0$ for all $j$, it follows from (32) that the
sequence $\\{d_{2j}(l)\\}$ has a limit $d^{*}_{2j}$ which is non-zero for all
$j$. Now (31) implies that that the sequence $\\{d_{1i}(l)\\}$ has a limit
$d^{*}_{1i}$ for all $i$. Moreover (33) and (34) show that these limits
satisfy

$d^{*}_{1i}=\frac{\tilde{u}_{i}}{\sum_{j=1}^{n}k_{ij}d^{*}_{2j}}\qquad
d^{*}_{2j}=\left(\frac{\tilde{v}_{j}}{\sum_{i=1}^{m}d^{*}_{1i}k_{ij}}\right)^{\frac{\gamma}{1+\gamma}}$
(36)

and (35) results in

$t^{*}_{ij}=d^{*}_{1i}k_{ij}d^{*}_{2j}$ (37)

${\hskip 220.50885pt}\square$

###### Remark 6

If one considers the limit, as $\gamma\rightarrow\infty$, in Algorithm 2, one
obtains the SK algorithm [17], whose sequence, $\\{T(l)\\}$, converges to the
minimizer for optimization problem (17), provided this problem is feasible.

${\hskip 220.50885pt}\square$

###### Remark 7

Before proceeding we note that [18] considers a problem which is related (but
different from) to the problem considered here. Namely, to solve

$\min_{T\in\mathcal{T}}KL(T|K)+\gamma_{1}KL(u_{T}|\tilde{u})+\gamma_{2}KL(v_{T}|\tilde{v})$
(38)

where $u_{T_{i}}=\sum_{j=1}^{n}t_{ij}$ for all $i$,
$k_{ij}=\exp(-c_{ij}/\gamma_{0})\text{ for all }(i,j)$ and
$\gamma_{1},\gamma_{2}>0$. In [18], neither marginal constraint is enforced
and none of the elements of $T$ are constrained to be zero. This is in
contrast to the current paper, in which the source marginal is constrained to
be $\tilde{u}$, and some of the elements of $T$ are enforced to be zero,
specifically $t_{ij}=0$ for $(i,j)\in\mathcal{Z}$. [18] shows that the
sequence, $\\{T(l)\\}$, generated by the following algorithm converges to a
limit which is the minimizer for optimization problem (38).

Chizat Algorithm. Initialize $d_{2j}(0)=1.$ Iterate for $l=0,1,\dots$

$\displaystyle d_{1i}(l+1)$
$\displaystyle=\left(\frac{\tilde{u}_{i}}{\Sigma_{j=1}^{n}k_{ij}d_{2j}(l)}\right)^{\frac{\gamma_{1}}{1+\gamma_{1}}}$
(39) $\displaystyle d_{2j}(l+1)$
$\displaystyle=\left(\frac{\tilde{v}_{j}}{\Sigma_{i=1}^{m}d_{1i}(l+1)k_{ij}}\right)^{\frac{\gamma_{2}}{1+\gamma_{2}}}$
(40) $\displaystyle t_{ij}(l+1)$ $\displaystyle=d_{1i}(l+1)k_{ij}d_{2j}(l+1)$
(41)

${\hskip 220.50885pt}\square$

###### Remark 8

Remarkably, if one considers the limit as $\gamma_{1}\rightarrow\infty$, then
the Chizat algorithm reduces to our Algorithm 2 in the case where there are no
constraints on $T$ and all the elements of $\tilde{T}$ equal one.

##  V Properties of the minimizer, $T^{*}$

###### Lemma 1

If $T^{*}$ is a minimizer for (21), then $t^{*}_{ij}>0$ for all
$(i,j)\notin\mathcal{Z}$.

###### Proof:

Suppose that $T^{*}$ is a minimizer for (21). From Remark 2, there is a
$\hat{T}\in\mathcal{V}(\tilde{u})$ with $\hat{t}_{ij}>0$ for all
$(i,j)\notin\mathcal{Z}$. Since $\mathcal{V}(\tilde{u})$ is convex,
$(1-\lambda)T^{*}+\lambda\hat{T}$ is in $\mathcal{V}(\tilde{u})$ for all
$\lambda\in[0,1]$. Also, there are bounds, $\beta_{1}$ and $\beta_{2}$, such
that, if $t^{*}_{ij}>0$ then for any $(i,j)\notin\mathcal{Z}$,
$0<\beta_{1}\leq t_{ij}\leq\beta_{2}$ for
$T=(1-\lambda)T^{*}+\lambda\hat{T}\text{ and }\lambda\in[0,1]$. The function
to be minimized in (21) can be expressed as

$f(T)=\sum_{(i,j)\notin{\mathcal{Z}}}kl(t_{ij}|k_{ij})+\gamma\sum_{j=1}^{n}kl\left(\sum_{l=1}^{m}t_{lj}\big{|}\tilde{v}_{j}\right)$

Consider any $\lambda\in(0,1]$ and $T=(1-\lambda)T^{*}+\lambda\hat{T}$. By the
mean value theorem, $\exists$ $\underline{\lambda}\in(0,\lambda)$ such that

$\displaystyle
f(T)=f(T^{*})+\lambda\\!\sum_{(i,j)\notin{\mathcal{Z}}}\frac{\partial
f}{\partial t_{ij}}(\underline{T})(\hat{t}_{ij}-t^{*}_{ij})$ (42)

where $\underline{T}=(1-\underline{\lambda})T^{*}+\underline{\lambda}\hat{T}$
and

$\frac{\partial f}{\partial
t_{ij}}(\underline{T})=\log\left(\frac{\underline{t}_{ij}}{k_{ij}}\right)+\gamma\log\left(\frac{\sum_{l=1}^{m}\underline{t}_{lj}}{\tilde{v}_{j}}\right)$
(43)

If $t^{*}_{ij}>0$ then, for all $\lambda\in(0,1]$, we have
$0<\beta_{1}\leq\underline{t}_{ij}\leq\beta_{2}$, and
$0<\beta_{1}\leq\sum_{l=1}^{m}\underline{t}_{lj}\leq\beta_{3}$ for some
$\beta_{3}$. Hence $\frac{\partial f}{\partial
t_{ij}}(\underline{T})(\hat{t}_{ij}-t^{*}_{ij})\leq\gamma_{ij}$ for some
$\gamma_{ij}$. Suppose that $t^{*}_{ij}=0$ for some $(i,j)\notin\mathcal{Z}$.
Then $\lim_{\lambda\rightarrow 0}\underline{t}_{ij}=t^{*}_{ij}=0$; and
$\lim_{\lambda\rightarrow 0}\frac{\partial f}{\partial
t_{ij}}(\underline{T})(\hat{t}_{ij}-t^{*}_{ij})=-\infty$. This implies that,
for $\lambda>0$ sufficiently small,
$\sum_{(i,j)\notin\mathcal{Z}}\frac{\partial f}{\partial
t_{ij}}(\underline{T})(\hat{t}_{ij}-t^{*}_{ij})<0$ which along with (42)
yields the contradiction, $f(T)<f(T^{*}).$

The following result can be obtained from the discussion of Algorithm 2; see
Remark 5. However we wish to provide a proof which is independent of any
algorithm.

###### Lemma 2

A matrix, $T^{*}$, solves the OT problem given by (21) iff $\exists$ positive
scalars, $d_{11},\dots,d_{1m}$ and $d_{21},\dots,d_{2n}$, such that, for all
$(i,j)$ (37) and (36) hold.

###### Proof:

If $(i,j)\in\mathcal{Z}$, we have $t^{*}_{ij}=k_{ij}=0$ and the expression for
$t^{*}_{ij}$ in (37) holds. The Lagrangian associated with this optimization
problem is

$\displaystyle L(T,\alpha)$
$\displaystyle=\sum_{(i,j)\notin{\mathcal{Z}}}kl(t_{ij}|k_{ij})+\gamma\sum_{j=1}^{n}kl\left(\sum_{l=1}^{m}t_{lj}|\tilde{v}_{j}\right)$
$\displaystyle+\sum_{i=1}^{m}\alpha_{i}\left(\sum_{l=1}^{n}t_{il}-\tilde{u}_{i}\right)$

For $(i,j)\notin\mathcal{Z}$ and $t_{ij}\neq 0$,

$\displaystyle\frac{\partial L}{\partial t_{ij}}(T,\alpha)$
$\displaystyle=\log\left(\frac{t_{ij}}{k_{ij}}\right)+\gamma\log\left(\frac{\sum_{l=1}^{m}t_{lj}}{\tilde{v}_{j}}\right)+\alpha_{i}$

From Lemma 1, we know that $t^{*}_{ij}>0$ for all $(i,j)\notin\mathcal{Z}$.
Hence $T^{*}$ is a minimizer iff $\exists$ scalars,
$\alpha_{1},\dots,\alpha_{m}$, such that $\frac{\partial L}{\partial
t_{ij}}(T^{*},\alpha)=0$ for all $(i,j)\notin\mathcal{Z}$, that is,
$\log\left(\frac{t^{*}_{ij}}{k_{ij}}\right)=-\alpha_{i}-\beta_{j}$, where

$\displaystyle\beta_{j}\equiv\gamma\log\left(\frac{\sum_{l=1}^{m}t^{*}_{lj}}{\tilde{v}_{j}}\right)$
(44)

Hence

$\displaystyle t^{*}_{ij}$
$\displaystyle=k_{ij}\exp\left(-\alpha_{i}-\beta_{j}\right)=\exp(-\alpha_{i})k_{ij}\exp(-\beta_{j})$
$\displaystyle=d_{1i}k_{ij}d_{2j}$ (45)

where $d_{1i}\equiv\exp(-\alpha_{i})$ and $d_{2j}\equiv\exp(-\beta_{j})$.
Using the row constraints, we have
$\tilde{u}_{i}=\sum_{j=1}^{m}t^{*}_{ij}=\sum_{j=1}^{m}d_{1i}k_{ij}d_{2j}$.
Hence, $d_{1i}=\frac{\tilde{u}_{i}}{\sum_{j=1}^{m}k_{ij}d_{2j}}$. It follows,
from (44) and (45), that

$\displaystyle
d_{2j}\sum_{i=1}^{m}d_{1i}k_{ij}=\sum_{i=1}^{m}t^{*}_{ij}=\tilde{v}_{j}e^{\beta_{j}/\gamma}=\tilde{v}_{j}d_{2j}^{-1/\gamma}$

that is,
$d_{2j}^{\frac{1+\gamma}{\gamma}}\sum_{i=1}^{m}d_{1i}k_{ij}=\tilde{v}_{j}$, or

$d_{2j}=\left(\frac{\tilde{v}_{j}}{\sum_{i=1}^{m}d_{1i}k_{ij}}\right)^{\frac{\gamma}{1+\gamma}}$

##  VI Convergence of the Modified SK algorithms

To prove Theorem 1, we need a result from [14]. Let

$\mathcal{X}=\left\\{x\in\mathbb{R}^{q}:x_{i}>0\text{ for
}i=1,2,\dots,q\right\\}$

and, for a fixed $\tilde{x}\in\mathcal{X}$, consider the strictly convex
function, $f:\bar{\mathcal{X}}\rightarrow\mathbb{R}$, given by

$f(x)=KL(x|\tilde{x})=\sum_{i=1}^{q}kl(x_{i}|\tilde{x}_{i})$ (46)

where $\overline{\mathcal{X}}$ is the closure of $\mathcal{X}$, that is

$\overline{\mathcal{X}}=\mathbb{R}^{q}_{+}=\left\\{x\in\mathbb{R}^{q}:x_{i}\geq
0\text{ for }i=1,2,\dots,q\right\\}$

Also, let $A_{i}\in\mathbb{R}^{m_{i}\times q}$ and
$b_{i}\in\mathbb{R}^{m_{i}}$, for $i=1,2,\dots,N$, for some positive integers,
$N$ and $m_{i}$, and consider the following optimization problem:

$\displaystyle\min_{x\in\overline{\mathcal{X}}}f(x)\qquad s.t.\quad
A_{i}x=b_{i}\quad i=1,\dots,N$ (47)

Let $\mathcal{C}_{i}$ be the closed convex set,
$\\{x\in\mathbb{R}^{q}:A_{i}x=b_{i}\\}$, and assume that
$\mathcal{C}\equiv\bigcap_{i=1}^{N}\mathcal{C}_{i}$ is non-empty. We also
require the following assumption.

###### Assumption 1

For each $i=1,2,\dots,N$ and $x\in\mathcal{X}$, $\exists$
$y^{*}\in\mathcal{X}\cap\mathcal{C}_{i}$ such that

$\displaystyle KL(y^{*}|x)=\min_{y\in\mathcal{X}\cap\mathcal{C}_{i}}KL(y|x)$
(48)

Note that, in this assumption, optimization is over
$\mathcal{X}\cap\mathcal{C}_{i}$ and not over
$\overline{\mathcal{X}}\cap\mathcal{C}_{i}$. We will denote the point $y^{*}$
above by $P_{i}(x)$ and refer to it as the $KL$-projection of $x$ onto
$\mathcal{C}_{i}$. Let $p$ be a permutation on $\\{1,2,\dots,N\\}$, that is

$p^{k}(i)\neq i\text{ for }k=1,\dots,N-1\qquad p^{N}(i)=i$

where $p^{k}$ is the application of $p$, $k$ times.

The following result may be gleaned from [14].

###### Theorem 2

Suppose $\mathcal{C}\cap\overline{\mathcal{X}}$ is non-empty, that Assumption
1 holds, and

$\displaystyle x(l+1)$ $\displaystyle=P_{i_{l}}(x(l))\quad\text{with}\quad
x(0)\in{\mathcal{X}}$ (49)

where $i_{l+1}=p(i_{l})$. Then,
$\lim_{l\rightarrow\infty}x(l)=x^{*}\in{\mathcal{C}}\cap\overline{\mathcal{X}}$.
Moreover, if

$\nabla f(x(0))\text{ is in the range of
}\left(A^{\prime}_{1}\;A^{\prime}_{2}\,\cdots\,A^{\prime}_{m}\right)$ (50)

where ′ denotes transpose, then

$\displaystyle f(x^{*})=\min_{x\in\overline{\mathcal{X}}}f(x)\quad s.t.\quad
A_{i}x=b_{i}\quad i=1,\dots,N$ (51)

The algorithm in Theorem 2 initially chooses some point, $x(0)\in\mathcal{X}$.
It then cycles indefinitely through each index, $i$, and projects onto
$\mathcal{C}_{i}$. It can be viewed as an alternating projection algorithm.
The resulting sequence, $\\{x(l)\\}_{l=0}^{\infty}$, converges to a point
which is common to all of the sets, $\mathcal{C}_{i}$ and
$\overline{\mathcal{X}}$. In addition, if $\nabla f(x(0))$ satisfies (50),
then $x^{*}$ is a minimizer for optimization problem (47). This algorithm is
very useful when one can readily solve the optimization problems in (48).

The optimization problem of this paper (21) can be rewritten as
$\min_{(T,v)\in\overline{\mathcal{TV}}}J(T,v)$, subject to

$\displaystyle\sum_{j=1}^{n}t_{ij}=\tilde{u}_{i},\,i=1,\dots,m,\qquad\sum_{i=1}^{m}t_{ij}=v_{j}\,,j=1,\dots
n,$ (52)

where

$\displaystyle J(T,v)$ $\displaystyle=KL(T|K)+\gamma
KL(v|\tilde{v})=KL(T|K)+KL(\gamma v|\tilde{\gamma v})$ (53)

$\gamma>0$ and

$\displaystyle\mathcal{S}=$
$\displaystyle\left\\{(T,v):T\in\mathbb{R}^{m\times
n},v\in\mathbb{R}^{n}\text{ and
}t_{ij}=0\;\forall\;(i,j)\in\mathcal{Z}\right\\}$ $\displaystyle\mathcal{TV}=$
$\displaystyle\\{(T,v)\in\mathcal{S}:t_{ij}>0\,\forall\,(i,j)\notin\mathcal{Z}\text{
and }v_{j}>0\,\forall\,i\\}$

Let $q=n_{T}+n$, where $n_{T}$ is the number of index pairs $(i,j)$ not in
$\mathcal{Z}$. Then, by appropriate definition of $x$ in $\mathbb{R}^{q}$, one
can associate each element of $x$ to an element of $T$ or $\gamma v$. We
denote this by $x=\text{vec}(T,\gamma v)$ and the objective function in (53)
can be written as $f(x)=KL(x|\tilde{x})$, where
$\tilde{x}=\text{vec}(K,\gamma\tilde{v})$. Also, the constraints in (52) can
be expressed as $A_{1}x=b_{1},\;\;A_{2}x=b_{2}$, where

$\displaystyle A_{1}x$
$\displaystyle=b_{1}\quad\text{iff}\quad(T,v)\in\tilde{\mathcal{C}}_{1}\equiv\left\\{(T,v)\in\mathcal{S}:T{\bm{1}}=\tilde{u}\right\\}$
$\displaystyle A_{2}x$
$\displaystyle=b_{2}\quad\text{iff}\quad(T,v)\in\tilde{\mathcal{C}}_{2}\equiv\left\\{(T,v)\in\mathcal{S}:T^{\prime}{\bm{1}}=v\right\\}$

If $(T,v),(S,w)\in\mathcal{TV}$, $y=\text{vec}(T,\gamma v)$,
$x=\text{vec}(S,\gamma w)$:

$KL(y|x)=KL(T|S)+\gamma KL(v|w)$

###### Lemma 3

If $(S,w)\in\mathcal{TV}$ then

$(T^{*},v^{*})=\text{argmin}_{(T,v)\in\mathcal{TV}\cap\tilde{\mathcal{C}}_{1}}KL(T|S)+\gamma
KL(v|w)$

iff $t^{*}_{ij}=c_{1i}s_{ij},\;\;v^{*}_{j}=w_{j}$, where

$c_{1i}=\frac{\tilde{u}_{i}}{\Sigma_{j=1}^{n}s_{ij}}$ (54)

###### Proof:

Since $t_{ij}=0$ for $(i,j)\in\mathcal{Z}$, the Lagrangian associated with
this optimization problem is

$\displaystyle L(T,v,\alpha)$ $\displaystyle=KL(T|S)+\gamma
KL(v|w)+\sum_{i=1}^{m}\alpha_{i}\left(\sum_{j=1}^{n}t_{ij}-\tilde{u}_{i}\right)$

When $(i,j)\notin\mathcal{Z}$,

$\displaystyle\frac{\partial L}{\partial
t_{ij}}=\log\left(\frac{t_{ij}}{s_{ij}}\right)+\alpha_{i}\qquad\frac{\partial
L}{\partial v_{j}}=\gamma\log\left(\frac{v_{j}}{w_{j}}\right)$

Setting these to zero yields

$\displaystyle\log\left(\frac{t^{*}_{ij}}{s_{ij}}\right)=-\alpha_{i}\qquad\log\left(\frac{v^{*}_{j}}{w_{j}}\right)=0$

Hence, $t^{*}_{ij}=s_{ij}\exp\left(-\alpha_{i}\right)=c_{1i}s_{ij}$, where
$c_{1i}=\exp(-\alpha_{i})$ and $v^{*}_{j}=w_{j}$. Also,
$\sum_{j=1}^{n}t^{*}_{ij}=\tilde{u}_{i}$ implies (54).

###### Lemma 4

If $(S,w)\in\mathcal{TV}$ then

$(T^{*},v^{*})=\text{argmin}_{(T,v)\in\mathcal{TV}\cap\tilde{\mathcal{C}}_{2}}KL(T|S)+\gamma
KL(v|w)$

iff $t^{*}_{ij}=c_{2j}s_{ij},\;\;v^{*}_{j}=c_{2j}^{-\frac{1}{\gamma}}w_{j}$,
where

$\displaystyle
c_{2j}=\left(\frac{w_{j}}{\Sigma_{i=1}^{m}s_{ij}}\right)^{\frac{\gamma}{1+\gamma}}$
(55)

###### Proof:

Since $t_{ij}=0$ for $(i,j)\in\mathcal{Z}$, the Lagrangian associated with
this optimization problem is given by

$\displaystyle L(T,v,\beta)$ $\displaystyle=KL(T|S)+\gamma
KL(v|w)+\sum_{j=1}^{n}\beta_{j}\left(\sum_{i=1}^{m}t_{ij}-v_{j}\right)$

When $(i,j)\notin{\mathcal{Z}}$,

$\displaystyle\frac{\partial L}{\partial
t_{ij}}=\log\left(\frac{t_{ij}}{s_{ij}}\right)+\beta_{j}\qquad\frac{\partial
L}{\partial v_{j}}=\gamma\log\left(\frac{v_{j}}{w_{j}}\right)-\beta_{j}$

Setting these to zero results in

$\displaystyle\log\left(\frac{t^{*}_{ij}}{s_{ij}}\right)=-\beta_{j}\qquad\quad\log\left(\frac{v^{*}_{j}}{w_{j}}\right)=\beta_{j}/\gamma$

Hence, $t^{*}_{ij}=s_{ij}c_{2j}$ where $c_{2j}=\exp(-\beta_{j})$. Also,
$v^{*}_{j}=w_{j}\exp(\beta_{j}/\gamma)=c_{2j}^{-1/\gamma}w_{j}$. Due to the
column sum constraints, we must have
$c_{2j}\sum_{i=1}^{m}s_{ij}=\sum_{i=1}^{m}t^{*}_{ij}=v^{*}_{j}=c_{2j}^{-1/\gamma}w_{j}$,
that is, $c_{2j}^{\frac{1+\gamma}{\gamma}}\sum_{i=1}^{m}s_{ij}=w_{j}$, which
implies (55).

###  VI-A Proof of Theorem 1

It follows, from Theorem 2 and Lemmas 3 and 4, that the sequence,
$\\{(T(l),v(l)\\}$, converges to a limit, $(T^{*},v^{*})$, with
$T^{*}\in{\mathcal{T}}$, $v^{*}\in\mathbb{R}^{n}_{+}$, and this limit
satisfies the constraints in (52). We now only need to prove that $\nabla
f(x(0))$ satisfies (50), where $x(0)=\text{vec}(K,\gamma\tilde{v})$. With $J$
given by (53), we have, for $t_{ij}>0$ and $v_{j}>0$,

$\frac{\partial J}{\partial
t_{ij}}(T,v)=\log\left(\frac{t_{ij}}{k_{ij}}\right),\qquad\frac{\partial
J}{\partial\gamma v_{j}}(T,v)=\log\left(\frac{\gamma
v_{j}}{\gamma\tilde{v}_{j}}\right)$

Hence, $\frac{\partial J}{\partial t_{ij}}(K,\tilde{v})=0,\;\;\frac{\partial
J}{\partial\gamma v_{j}}(K,\tilde{v})=0$, and $\nabla f(x(0))_{i}=0$ for
$i=1,2,\dots,q$. Therefore, $\nabla f(x(0))$ satisfies (50).

${\hskip 220.50885pt}\blacksquare$

##  VII Example: resource allocation in the sharing economy

SK-type iterations make sense for large-scale problems. Smart cities are a
natural place to look for such problems, and, in particular, sharing economy
[20] applications, since these are precisely the problem domains where
allocation of resources at scale emerge, and where the scale of the problem is
subject to temporal variations. One such problem arises in the context of
charging $m$ electric vehicles (EVs) overnight in a city such as London. With
the advent of vehicle-to-grid (V2G), vehicle-to-vehicle (V2V) and widespread
availability of solar, it is likely that many entities that currently consume
energy will become prosumers in the near future; i.e. most homes, and even
cars, will consume energy and also make energy available, depending on the
circumstance. Such a scenario is clearly very large scale, with a potentially
very large number, $m$, of EVs (there are currently more than 2.5M cars
registered in London), and a large number, $n$, of energy providers
(consisting of conventional utilities, energy brokerages, households and even
other cars, and constituting the target agents of our OT setup in Section
III). In the setting of this paper (Section III), each EV specifies the energy
it requires in KWhrs, and then this demand for energy is communicated (i.e.
‘transported’) to a set of providers. Each provider may set a cost based on
their type of energy generation, their proximity to the car being charged, and
the type of vehicle being charged. In addition, some providers may prohibit
certain types of vehicles, for example plug-in hybrid vehicles (PHEVs), or
vehicles that are very large in size. Our OT formulation (21) captures the
realistic scenario in which the $n$ energy providers make available flexible
amounts of energy, nominally $\tilde{v}$, to $m$ EVs whose demands are exactly
$\tilde{u}$. Furthermore, energy transfer between specified provider-EV pairs,
$(i,j)\in{\mathcal{Z}}$ (5), are not allowed a priori.

To simulate this scenario, consider $m=10,000$ EVs requiring charging and
$n=10$ energy providers. We simulate the specified energy charging
requirement, $\tilde{u}_{i}$, $i=1,\ldots,m$, of each EV via independent and
uniform draws in the range $(0,1)$ KWhrs. We simulate the nominal available
energy of each supplier, $\tilde{v}_{j}$, $j=1,\ldots,n$, in the same way.
Therefore, (1) may not hold, i.e. the transport problem may be unbalanced, in
that the total energy required by the EVs may differ from the total nominal
energy made available by the providers. This UOT problem is also solved by our
algorithms. To effect prior zero constraints, we assume that even-indexed cars
are PHEVs, and that even-indexed providers will not supply these PHEVs. This
defines a transport plan with $\frac{mn}{4}$ pre-specified zeroes. The non-
zero elements of $K$ in (15) are obtained with $\tilde{t}_{ij}=1$ and
$\gamma_{0}=1.99$, and the transport cost, $c_{ij}$,
$\forall(i,j)\notin\mathcal{Z}$ are again iid uniformly drawn from $(0,1)$. We
use Algorithm 1, with $\gamma=1.005$ (21), to obtain the OT plan, $T^{*}$.
Figure 1 illustrates the convergence of the algorithm for 10 simulation runs.
Let $T(l-1)$ and $T(l)$ be the transport plans obtained from two consecutive
iterations of the algorithm in any one simulation, and let $\Delta(l)$ be the
matrix with $(i,j)$th entry defined as
$\delta_{ij}(l)=|t_{ij}(l)-t_{ij}(l-1)|$. Figure 1 plots the log of the sum of
the $\delta_{ij}$s, normalised by the first log-sum, for each of the 10
simulations, i.e.
$\frac{\log\left(\sum_{(i,j)\notin{\mathcal{Z}}}|t_{ij}(l)-t_{ij}(l-1)|\right)}{\log\left(\sum_{(i,j)\notin{\mathcal{Z}}}|t_{ij}(1)-t_{ij}(0)|\right)}$.
We have normalized in order to facilitate comparison of the different
simulations.

Figure 1: Convergence of the modified algorithm (Algorithm 1)

##  VIII Conclusions

In this paper, we have presented SK-type algorithms for constrained optimal
transport. Specifically, our algorithms allow for transport plans that force
some entries to be zero a priori. The convergence proof is based on Bregman-
type ideas. An example in resource allocation for the sharing economy is
provided, pointing to situations in which our algorithm is relevant. Future
work will investigate extension of this work to problems that incorporate
other forms of regularisation terms [21], and other potential use-cases.

## References

  * [1] Y. Chen, G. Tryphon, and M. Pavon, “Optimal transport in systems and control,” Annual Review of Control, Robotics, and Autonomous Systems, vol. 4, no. 1, pp. 89–113, 2021.
  * [2] M. Kusner, Y. Sun, N. Kolkin, and K. Weinberger, “From word embeddings to document distances,” in Proceedings of the 32nd International Conference on Machine Learning (F. Bach and D. Blei, eds.), vol. 37 of Proceedings of Machine Learning Research, (Lille, France), pp. 957–966, PMLR, 07–09 Jul 2015.
  * [3] K. Saenko, B. Kulis, M. Fritz, and T. Darrell, “Adapting visual category models to new domains,” in Computer Vision – ECCV 2010 (K. Daniilidis, P. Maragos, and N. Paragios, eds.), (Berlin, Heidelberg), pp. 213–226, Springer Berlin Heidelberg, 2010.
  * [4] C. Villani, Optimal Transport: Old and New.  Springer, 2008.
  * [5] A. Galichon, Optimal transport methods in economics.  Princeton University Press, 2016.
  * [6] N. Courty, R. Flamary, D. Tuia, and A. Rakotomamonjy, “Optimal transport for domain adaptation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 39, no. 9, pp. 1853–1865, 2017.
  * [7] P. Gordaliza, E. Barrio, G. Fabrice, and J. Loubes, “Obtaining fairness using optimal transport theory,” in ICML, pp. 214–223, 2019.
  * [8] G. Peyré and M. Cuturi, Computational Optimal Transport: With Applications to Data Science.  Foundations and Trends in Machine Learning, 2019.
  * [9] S. Boufelja Y., A. Quinn, M. Corless, and R. Shorten, “Fully probabilistic design for optimal transport,” arXiv, no. 2212.13912, 2023.
  * [10] M. Cuturi, “Sinkhorn distances: Lightspeed computation of optimal transport.,” in NIPS, pp. 2292–2300, 2013.
  * [11] S. Marino and A. Gerolin, “An optimal transport approach for the Schrödinger bridge problem and convergence of Sinkhorn algorithm,” Journal of Scientific Compututing, vol. 85, no. 27, 2020.
  * [12] M. Arjovsky, S. Chintala, and L. Bottou, “Wasserstein generative adversarial networks,” in ICML, pp. 214–223, 2017.
  * [13] J. Solomon, F. Goes, G. Peyré, M. Cuturi, A. Butscher, A. Nguyen, T. Du, and L. Guibas, “Convolutional Wasserstein distances: Efficient optimal transportation on geometric domains,” ACM Transactions on Graphics, vol. 34, 07 2015.
  * [14] L. M. Bregman, “The relaxation method of finding the common point of convex sets and its application to the solution of problems in convex programming,” USSR Computational Mathematics and Mathematical Physics, vol. 7, pp. 200–217, 1967.
  * [15] R. A. Brualdi, “Convex sets of non-negative matrices,” Canadian Journal of Mathematics, vol. 20, p. 144–157, 1968.
  * [16] A. Baradat and E. Ventre, “Convergence of the Sinkhorn algorithm when the Schrodinger problem has no solution,” arXiv, no. 2207.02977, 2022.
  * [17] P. A. Knight, “The Sinkhorn–Knopp algorithm: Convergence and applications,” SIAM Journal on Matrix Analysis and Applications, vol. 30, no. 1, pp. 261–275, 2008.
  * [18] L. Chizat, G. Peyré, B. Schmitzer, and F.-X. Vialard, “Scaling algorithms for unbalanced optimal transport problems,” Mathematics of Computation, vol. 87, pp. 2563–2609, November 2018.
  * [19] A. Ringh and J. K. Haasler I. Y. Chen, “Graph-structured tensor optimization for nonlinear density control and mean field games,” arXiv, no. 2112.05645v2, 2022.
  * [20] E. Crisostomi, B. Ghaddar, F. Haeusler, J. Naoum-Sawaya, G. Russo, and R. Shorten, Analytics for the Sharing Economy: Mathematics, Engineering and Business Perspectives.  Springer, 2020.
  * [21] S. Ferradans, N. Papadakis, G. Peyré, and J.-F. Aujol, “Regularized discrete optimal transport,” SIAM Journal of Imaging Sciences, vol. 7, no. 3, pp. 1853–1882, 2014.