# Simultaneous Trapping of Two Optical Pulses in an Atomic Ensemble as
Stationary Light Pulses

U-Shin Kim<EMAIL_ADDRESS>Yoon-Ho Kim<EMAIL_ADDRESS>Department of Physics, Pohang University of Science and Technology (POSTECH),
Pohang 37673, Korea

###### Abstract

The stationary light pulse (SLP) refers to a zero-group-velocity optical pulse
in an atomic ensemble prepared by two counter-propagating driving fields.
Despite the uniqueness of an optical pulse trapped within an atomic medium
without a cavity, observations of SLP so far have been limited to trapping a
single optical pulse due to the stringent SLP phase-matching condition, and
this has severely hindered the development of SLP-based applications. In this
paper, we first show theoretically that the SLP process in fact supports two
phase-matching conditions and we then utilize the result to experimentally
demonstrate simultaneous SLP trapping of two optical pulses for the duration
from 0.8 $\mu$s to 2.0 $\mu$s. The characteristic dissipation time, obtained
by the release efficiency measurement from the SLP trapping state, is 1.22
$\mu$s, which corresponds to an effective Q-factor of $2.9\times 10^{9}$. Our
work is expected to bring forth interesting SLP-based applications, such as,
efficient photon-photon interaction, spatially multi-mode coherent quantum
memory, creation of exotic photonic gas states, etc.

A stationary light pulse (SLP) is a zero-group-velocity optical pulse trapped
in an atomic ensemble prepared by two counter-propagating classical driving
fields Hansen07 ; Otterbach09 ; Nikoghosyan09 . In contrast to the light
storage based on electromagnetically induced transparency (EIT) in which the
photonic excitation is fully mapped to the atomic coherence Fleishhauer02 ;
Liu01 ; Phillips01 ; Cho10a ; Cho10b ; Cho16 , the SLP process effectively
traps the optical pulse in the atomic ensemble while retaining its
electromagnetic field nature Hansen07 ; Otterbach09 ; Nikoghosyan09 ; Zimmer08
; Yan12 ; Iakoupov16 . In the experiment, SLP for a classical light pulse has
been reported in a warm atomic vapor system Bajcsy03  and in a cold atomic
ensemble Lin09 ; Everett17 . Recently, based on a cold atomic ensemble,
quantum SLP has been demonstrated Park18  and the formation of quantum SLP
from a free-propagating single-photon has also been reported Kim22 .

The fact that a zero-group-velocity optical pulse can be trapped in an atomic
ensemble without the need for a cavity has brought forth to proposals on novel
SLP-based applications, such as, efficient cavity-free nonlinear optics
Andre05 ; Chen12 ; Murray17 , photonic quantum gates Iakoupov18 ,
fermionization of polaritons Chang08 , mimicking interacting relativistic
theories with SLPs Angelakis13 , Bose-Einstein condensation of stationary-
light polaritons Fleischhauer08 , etc. Although the majority of such proposed
SLP-based applications require interactions of multiple SLPs occupying
different spatial modes, observations of SLP so far have been limited to
trapping a single optical pulse due to the stringent SLP phase-matching
condition Zhou11 , which has been one of the bottlenecks for realizing various
theoretical proposals on SLP-based physics. Extending the SLP physics to
support multiple spatial modes would therefore promote the development of, in
addition to the above-mentioned potential applications, multiplexed quantum
memory, optical image buffer based on SLP, and efficient photon-photon
nonlinear optics.

In this work, we demonstrate the simultaneous SLP trapping of two optical
pulses in an atomic ensemble. We first show in theory that the SLP process in
fact supports two phase-matching conditions and we then utilize the result to
experimentally demonstrate simultaneous SLP trapping of two optical pulses for
the duration from 0.8 $\mu$s to 2.0 $\mu$s. The formation of SLPs is verified
by comparing the probe waveforms of SLP with those of EIT. We have also
carried out the release efficiency measurement from the SLP trapping state,
and the resulting characteristic dissipation time 1.22 $\mu$s, which
corresponds to an effective Q-factor of $2.9\times 10^{9}$.

Propagation of the photon-atom excitation through the atomic ensemble is
described by a quasi-particle $\hat{\psi}$ described as Moiseev14 ,

$\hat{\psi}=[\mathcal{\hat{E}}^{+}\mathrm{cos}\phi+\mathcal{\hat{E}}^{-}\mathrm{sin}\phi]\mathrm{cos}\theta-\hat{\mathcal{S}}\mathrm{sin}\theta,$
(1)

where $\mathcal{\hat{E}}^{+}$ and $\mathcal{\hat{E}}^{-}$ are the slowly
varying field operators for the forward-propagating and backward-propagating
optical fields, respectively, and $\hat{\mathcal{S}}$ is the collective atomic
spin operator. Here,
$\mathrm{tan}\phi=\Omega_{\mathrm{BWC}}/\Omega_{\mathrm{FWC}}$, where
$\Omega_{\mathrm{FWC}}$ and $\Omega_{\mathrm{BWC}}$ respectively are the Rabi
frequencies due to the forward coupling (FWC) and backward coupling (BWC)
classical driving fields for the atomic medium, and
$\mathrm{tan}^{2}\theta=g^{2}n/\Omega^{2}$, where $g$ is the atom-photon
coupling constant, $n$ is the number of atoms, and
$\Omega^{2}=\Omega_{\mathrm{FWC}}^{2}+\Omega_{\mathrm{BWC}}^{2}$. The
polariton propagates through the atomic medium with the group velocity
$v_{g}=c_{0}\mathrm{cos}^{2}\theta\mathrm{cos}2\phi$, where $c_{0}$ is the
speed of light in vacuum Otterbach14 . Note that, when the Rabi frequencies
due to FWC and BWC driving fields are balanced, i.e., $\phi=45^{\circ}$, the
SLP state with $v_{g}=0$ is formed within the atomic medium.

Figure 1: (a) Typical phase-matching condition of the SLP trapping process.
The momentum conservation conditions,
$\vec{k}^{+}_{1}=\vec{k}_{s}+\vec{k}_{\mathrm{FWC}}$ and
$\vec{k}^{-}_{1}=\vec{k}_{s}+\vec{k}_{\mathrm{BWC}}$, and energy conservation
condition, $|\vec{k}^{+}_{1}|=|\vec{k}^{-}_{1}|$, must be simultaneously
satisfied, leading to the angle restriction between $\vec{k}_{1}$ and
$\vec{k}_{\mathrm{FWC}}$. (b) The mirror symmetry with respect to
$\vec{k}_{\mathrm{FWC}}$ helps to identify another equally-likely SLP phase-
matching condition, enabling simultaneous trapping of two optical pulses.
Note, however, that the three vectors $\vec{k}_{1}$, $\vec{k}_{2}$, and
$\vec{k}_{\mathrm{FWC}}$ need not lie on the same plane.

In addition to balancing FWC and BWC Rabi frequencies, the formation of SLP
from an optical pulse injected into the atomic medium requires a stringent
phase-matching condition, see Fig. 1(a). First, a probe pulse, $\vec{k}_{1}$,
injected into the atomic medium dressed by FWC, $\vec{k}_{\mathrm{FWC}}$, is
fully mapped into the collective atomic spin state $\hat{\mathcal{S}}$ with
the wave vector $\vec{k}_{s}=\vec{k}_{1}-\vec{k}_{\mathrm{FWC}}$ by
adiabatically turning off FWC, which corresponds to the EIT storage process
Fleischhauer00 . If FWC is turned back on after a specified storage time, the
propagating optical pulse is retrieved from the collective atomic spin state.
For the SLP trapping of the retrieved pulse, BWC is simultaneously turned on
with FWC. At this step, the wave vectors of trapped forward-propagating field,
$\vec{k}^{+}_{1}$, and backward-propagating field, $\vec{k}^{-}_{1}$, must
simultaneously satisfy the following conditions:
$\vec{k}^{+}_{1}=\vec{k}_{s}+\vec{k}_{\mathrm{FWC}}$ and
$\vec{k}^{-}_{1}=\vec{k}_{s}+\vec{k}_{\mathrm{BWC}}$ for the momentum
conservation, and $|\vec{k}^{+}_{1}|=|\vec{k}^{-}_{1}|$ for the energy
conservation. This phase-matching condition determines the injection angle
between probe 1 and FWC to be 0.345∘ for our experimental conditions Park18 ;
Kim22 , in stark contrast to EIT storage which does not require any phase-
matching for the incoming probe beams Cho10b .

Simultaneous trapping of two SLPs requires another input angle that satisfies
the phase-matching condition and it can be found by considering the mirror
symmetry with respect to $\vec{k}_{\mathrm{FWC}}$. As shown in Fig. 1(b), a
second probe pulse $\vec{k}_{2}$ injected into the atomic medium symmetrically
with respect to $\vec{k}_{\mathrm{FWC}}$ can create the phase-matching
condition which is a mirrored image of the typical SLP phase-matching
condition shown in Fig. 1(a). Note that the wave vector for the spin-wave in
the case of Fig. 1(b) is opposite to the case of Fig. 1(a). Therefore, for
simultaneous SLP trapping of two optical pulses, the second probe pulse,
$\vec{k}_{2}$, must be injected into the atomic medium at the angle of
$-0.345^{\circ}$ with respect to $\vec{k}_{\mathrm{FWC}}$, while the first
probe pulse, $\vec{k}_{1}$, is injected at the angle of $+0.345^{\circ}$ with
respect to $\vec{k}_{\mathrm{FWC}}$.

We note that the phase-matching condition shown in Fig. 1 can support the
formation of multi-mode SLPs for more than two optical pulses because the
conditions Fig. 1(a) and Fig. 1(b) can be satisfied separately. As the only
restriction is the angles between the probe pulses $\vec{k}_{1}$ and
$\vec{k}_{2}$ and the FWC beam $\vec{k}_{\mathrm{FWC}}$, the wave vectors
$\vec{k}_{1}$ and $\vec{k}_{2}$ form a cone around $\vec{k}_{\mathrm{FWC}}$.
The three vectors need not line on the same plane. This feature can be
utilized to enable multi-mode photon-photon interaction for interesting
quantum information applications, somewhat analogous to multi-mode EIT quantum
memory pu17 . In our experiment, we demonstrate the pair-wise formation of
SLPs on the plane parallel to the optical table.

Figure 2: (a) Experimental schematic for the simultaneous SLP trapping of two
optical pulses. (b) The level diagram of 87Rb used in the experiment. The
detuning $\Delta$=$2\pi\times 4$ MHz is introduced to suppress higher-order
Raman excitations, which hinders the formation of SLP. APD: avalanched photo-
diode, FWC: forward coupling beam, BWC: backward coupling beam, MOT: magneto-
optical trap.

The schematic of the experiment to demonstrate simultaneous SLP trapping of
two optical pulses is shown in Fig. 2(a). We exploit a cigar-shaped 87Rb cold
atomic ensemble consisting of a pair of rectangular anti-Helmholtz coils and
1-inch-diameter trapping beams of 130 mW Cho14 ; Lee16 ; Park17 ; Zhao19 . The
duration of each experimental sequence is 200 ms, of which 190 ms is used for
preparing the atomic medium starting with the plain magneto-optical trap (MOT)
loading phase of 173.5 ms. It is then followed by the temporal dark MOT phase
of 15 ms to trap more atoms by gradually ramping down the trapping beam
powers. The decreased beam powers help to increase the optical depth (OD) of
the atomic ensemble by preventing the trapped atoms from reheating by the
intense trapping beams Ketterle92 ; Hsiao14 . At the start of the temporal
dark MOT phase, the magnetic field compression phase is simultaneously
initiated by ramping up the current of the anti-Helmholtz coils from 5.8 A to
13.6 A (SM100-AR-75, DELTA ELEKTRONIKA). After 13 ms of the magnetic field
compression phase, the anti-Helmholtz coils are then switched off using an
insulated-gate bipolar transistor (IGBT) (1MBI400V-120-50, Fuji Electronic) to
remove residual magnetic field stray. Finally, the atoms are prepared in the
ground state $|5^{2}S_{1/2},F=1\rangle$ by shining a 1.5 ms optical pulse
which is resonant to $|5^{2}S_{1/2},F=2\rangle$ to $|5^{2}P_{3/2},F=2\rangle$
transition. The SLP trapping experiment is carried out during the experimental
window of 10 ms and the relevant energy level diagram is shown in Fig. 2(b).
The probes and the coupling beams are prepared from a 795 nm-centered
external-cavity diode laser (ECDL) combined with a tapered amplifier system
(TA Pro, Toptica). The ECDL is locked at the 120 MHZ red-detuned frequency
from $|5^{2}S_{1/2},F=2\rangle$ to $|5^{2}P_{1/2},F=2\rangle$ transition. FWC
and BWC are then modulated using acousto-optic modulators (AOM) (ATM-1201A2,
IntraAction) for the experimental sequence by 120 MHz and 116 MHz,
respectively. Probe 1 and probe 2 are prepared by up-shifting the optical
frequency of the ECDL using an electro-optic modulator (EOM) (Visible Phase
Modulator 4851, New Focus) by 6.79 GHz. Finally, the probes are modulated
using AOMs (ATM-801A2, IntraAction) by 80 MHz to have a 2 $\mu$s temporal
width and the frequency resonant to $|5^{2}S_{1/2},F=1\rangle$ to
$|5^{2}P_{1/2},F=2\rangle$ transition as shown in Fig. 2(b). As a result, the
atomic ensemble prepared for our experiment has an OD of 60, the dephasing
rate between the two ground states of $2\pi\times 60$ kHz, and the phase
mismatch for SLP process, $\Delta_{k}L$, of 0.05.

For the available 10 ms experimental window, for optimal results, we use the
first 10 $\mu$s for the SLP experiment. First, probe 1 and probe 2 with 2
$\mu$s temporal widths are simultaneously injected to the atomic ensemble
dressed by FWC with the phase-matching angle $\pm 0.345^{\circ}$ with respect
to $\vec{k}_{\mathrm{FWC}}$ as shown in Fig. 2(a). Probe 1 and probe 2 are
then mapped into the collective spin states by adiabatically turning off FWC,
which corresponds to the EIT storage process. After the 2 $\mu$s EIT storage
process, FWC and BWC are turned on simultaneously to form SLPs of probe 1 and
probe 2 inside the atomic ensemble. The detuning between FWC and BWC,
$\Delta=2\pi\times 4$ MHz, is introduced to avoid higher-order Raman
excitations, which hinder the formation of SLP Wu10 . After a specified SLP
trapping duration, the probe pulses are released from the SLP states by
turning off BWC and detected by two avalanched photo-diodes (APDs) (APD120A/M,
Thorlabs). All traces from the APDs are acquired by averaging 50 experimental
sequences using an 1 GHz-bandwidth oscilloscope (DPO7104, Tektronix).

Figure 3 shows, for probe 1 and probe 2, the experimental results of the EIT
storage of 2 $\mu$s only and those of the EIT storage of 2 $\mu$s followed by
the SLP trapping of 1 $\mu$s. The prepared cold atom ensemble (for the
duration of 190 ms prior to the experimental window) is dressed with FWC and,
at the time of 1.6 $\mu$s, probe 1 and probe 2 are injected into the atomic
medium as shown in Fig. 2. Due to EIT formed by FWC, when FWC is turned off
between 3.6 $\mu$s and 5.6 $\mu$s, probe 1 and probe 2 are simultaneously
mapped to the collective atomic spin states. After the EIT storage of 2
$\mu$s, FWC is turned back on at the time of 5.6 $\mu$s to retrieve the two
probe pulses as propagating optical pulses from the collective atomic spin
states, as evidenced in the EIT traces in Fig. 3. If both FWC and BWC are
turned on simultaneously at 5.6 $\mu$s to 6.6 $\mu$s, the retrieved optical
pulses are immediately trapped within the atomic medium as SLPs, and, when BWC
is later turned off, SLPs are converted into propagating optical pulses. The
suppressed emission during the time both FWC and BWC are simultaneously turned
on and the subsequent emission of the optical pulses when BWC is turned off
clearly indicate the formation of simultaneous SLPs for probe 1 and probe 2
for the duration of 1 $\mu$s.

Figure 3: EIT storage versus SLP trapping for (a) probe 1 and (b) probe 2. EIT
storage occurs by turning off FWC for 2 $\mu$s. When only FWC is turned back
on, probe pulses are immediately retrieved, but if both FWC and BWC are turned
on simultaneously, photon retrieval is suppressed and delayed until BWC is
turned off, signaling SLP formation. (c) Temporal sequence of FWC and BWC
laser beams.

The EIT/SLP efficiency differences observed in Fig. 3(a) and Fig. 3(b) for
probe 1 and probe 2 come from the fact that, in our experimental setup, each
probe beam experiences a slightly different OD and has a different overlap
with FWC. This can be confirmed by comparing their efficiencies and group
delays during the EIT storage processes shown in Fig. 3. Let us first consider
the retrieval efficiencies of the two probe beams after the 2 $\mu$s EIT
storage process. The retrieval efficiency of probe 1 is measured to be 80.8%,
and that of probe 2 is measured to be 59.1%. Since higher OD gives higher
retrieval efficiency under the same storage duration, we conclude that probe 1
experiences higher OD than probe 2 does. Let us now compare the group delays
(due to the slow-light effect) experienced by the two probe beams after the 2
$\mu$s EIT storage process. The group delay is calculated by subtracting the
storage time of 2 $\mu$s from the total time delay measured at the peak of the
observed probe waveform. The group delay of probe 1 is measured to be 1.9
$\mu$s, and that of probe 2 is measured to be 2.4 $\mu$s from Fig. 3. In
general, higher OD and lower $\Omega_{\mathrm{FWC}}$ result in a larger group
delay. Since probe 1 experiences higher OD than probe 2, we can conclude that
the spatial overlap between probe 2 and FWC is slightly worse than the spatial
overlap between probe 1 and FWC.

Figure 4: Simultaneous SLP trapping of (a) probe 1 and (b) probe 2 from 0.8
$\mu$s to 2.0 $\mu$s. (c) Temporal sequence of FWC and BWC laser beams.

In Fig. 4, we show the experimental results of simultaneous SLP trapping of
probe 1 and probe 2 for different SLP trapping times by varying the duration
of the BWC turn-on time. As before, SLPs are formed for probe 1 and probe 2
after the EIT storage process of 2 $\mu$s. The SLP trapping time is varied
from 0.8 $\mu$s to 2.0 $\mu$s in 0.2 $\mu$s steps. The experimental results in
Fig. 4 clearly display the suppressed emission of the probes while FWC and BWC
are both turned on, and the delayed release of the probe pulses when BWC is
turned off after a specified SLP trapping time.

As the SLP trapping times become longer, the peaks of the released optical
pulses become smaller as shown in Fig. 4. In Fig. 5, we summarize the release
efficiency measurements of probe 1 and probe 2, as a function of the SLP
trapping time, from the simultaneous SLP trapping shown in Fig. 4. The release
efficiencies of probe 1 and probe 2 exhibit the same exponential energy
dissipation with the characteristic dissipation time of $\tau=1.22$ $\mu$s.
Since probe 1 and probe 2 are trapped in the atomic ensemble while retaining
their electromagnetic field nature, the simultaneous SLP trapping process can
be understood with the analogy of an optical cavity system that confines the
electromagnetic field within a restricted volume. The Q-factor of an optical
cavity system is typically given by $Q={2\pi f_{0}E}/{P}=2\pi f_{0}\tau$,
where $f_{0}$ is the resonant frequency, $E$ is the stored energy in the
cavity, $P=-{dE}/{dt}$ is the energy dissipation rate, and $\tau$ is the
characteristic energy dissipation time Saleh07 ; anderson84 . The
characteristic dissipation time of $\tau=1.22$ $\mu$s corresponds to the
optical cavity system with an effective Q-factor of $2.9\times 10^{9}$ Chen12
. Note that the effective Q-factor of $2.9\times 10^{9}$ is comparable to
other ultrahigh Q-factor optical cavity systems such as the silicon or
crystalline resonator system.

Figure 5: The observed SLP trapping efficiency as the trapping times are
varied. Both probe 1 and probe 2, as they undergo simultaneous SLP processes,
exhibit the same exponential decay with the characteristic dissipation time of
1.22 $\mu$s represented as the solid line. The characteristic dissipation time
of 1.22 $\mu$s corresponds to the effective Q-factor of $2.9\times 10^{9}$.

The high effective Q-factor of the SLP trapping process demonstrated in this
work implies that there is a strong atom-photon coupling. The atom-photon
coupling can be estimated by $N$-atom cooperativity, $C_{N}$, as the photons
are effectively trapped within the atomic ensemble, similarly to a cavity
quantum electrodynamics (QED) system. The $N$-atom cooperativity is given as
$C_{N}=4g^{2}N/\kappa\Gamma$, where $g$ is the atom-field coupling constant,
$\kappa$ is the effective linewidth of the cavity, $\Gamma$ is the decay rate
of the atom, and $N$ is the number of atoms beck16 ; beck14 . In our SLP
experimental setup, the relevant parameters are estimated as
$\\{g,\kappa,\Gamma\\}=2\pi\times\\{0.24,0.13,5.8\\}\ \mathrm{MHz}$ from the
experimental parameters (the optical depth $\mathrm{OD}=60$, the group delay
$\tau_{g}\approx 2\ \mu\mathrm{s}$, the Q-factor $Q=2.9\times 10^{9}$, and the
length of atomic ensemble $L=10\ \mathrm{mm}$) using the relations
$\tau_{g}=\mathrm{OD}\times\Gamma/|\Omega_{\mathrm{FWC}}|^{2}$ fleischhauer05
and
$v_{g}={L}/{\tau_{g}}={c_{0}|\Omega_{\mathrm{FWC}}|^{2}}/({|\Omega_{\mathrm{FWC}}|^{2}+g^{2}N})$
Fleischhauer00 . This leads to the $N$-atom cooperativity of our SLP system as
$C_{N}\sim 8\times 10^{6}$, which is larger by an order of magnitude than
those reported in recent cavity QED experiments beck16 ; beck14 ; chen13 ;
hosseini16 .

In summary, we have for the first time demonstrated simultaneous SLP trapping
of two optical pulses in a cold atomic ensemble. We have shown theoretically
that the SLP process in fact supports two phase-matching geometries, and
experimentally demonstrated simultaneous SLP trapping of two optical pulses
for the duration from 0.8 $\mu$s to 2.0 $\mu$s. Although the SLP trapping
processes for the two optical pulses showed slightly different efficiencies
due to the experimental geometry, both SLP trapping processes exhibited the
same characteristic dissipation time of 1.22 $\mu$s, which corresponds to an
effective Q-factor of $2.9\times 10^{9}$. As our SLP system reports a very
large effective $N$-atom cooperativity, $C_{N}\sim 8\times 10^{6}$, our work
is expected to bring forth interesting SLP-based applications, such as,
cavity-free nonlinear optics Andre05 ; Chen12 ; Murray17 , photonic quantum
gates Iakoupov18 , fermionization of polaritons Chang08 , mimicking
relativistic theories with SLPs Angelakis13 , and Bose-Einstein condensates of
stationary-light polaritons Fleischhauer08 .

This work was supported in part by the National Research Foundation of Korea
(2019R1A2C3004812) and the ITRC support program (IITP-2022-2020-0-01606).

## References

  * (1) K. R. Hansen and K. Mølmer, Phys. Rev. A 75, 053802 (2007).
  * (2) J. Otterbach, R. G. Unanyan, and M. Fleischhauer, Phys. Rev. Lett. 102, 063602 (2009).
  * (3) G. Nikoghosyan and M. Fleischhauer, Phys. Rev. A 80, 013818 (2009).
  * (4) M. Fleischhauer and M. D. Lukin, Phys. Rev. A 65, 022314 (2002).
  * (5) C. Liu, Z. Dutton, C. H. Behroozi, and L. V. Hau, Nature 409, 490-493 (2001).
  * (6) D. F. Phillips, A. Fleischhauer, A. Mair, R. L. Walsworth, and M. D. Lukin, Phys. Rev. Lett. 86, 783 (2001).
  * (7) Y.-W. Cho and Y.-H. Kim, Opt. Express 18, 25786 (2010).
  * (8) Y.-W. Cho and Y.-H. Kim, Phys. Rev. A82, 033830 (2010).
  * (9) Y.-W. Cho, G. T. Campbell, J. L. Everett, J. Bernu, D. B. Higginbottom, M. T. Cao, J. Geng, N. P. Robins, P. K. Lam, and B. C. Buchler, Optica 3, 100-107 (2016).
  * (10) F. E. Zimmer, J. Otterbach, R. G. Unanyan, B. W. Shore, and M. Fleischhauer, Phys. Rev. A 77, 063823 (2008).
  * (11) Y. Zhang, Y. Zhang, X.-H. Zhang, M. Yu, C.-L. Cui, J.-H. Wu, Phys. Lett. A 376, 656-661 (2012).
  * (12) I. Iakoupov, J. R. Ott, D. E. Chang, and A. S. Sørensen, Phys. Rev. A 94, 053824 (2016).
  * (13) M. Bajcsy, A. S. Zibrov, and M. D. Lukin, Nature 426, 638-641 (2003).
  * (14) Y.-W. Lin, W.-T. Liao, T. Peters, H.-C. Chou, J.-S. Wang, H.-W. Cho, P.-C. Kuan, and I. A. Yu, Phys. Rev. Lett. 102, 213601 (2009).
  * (15) J. L. Everett, G. T. Campbell, Y.-W. Cho, P. V.-Gris, D. B. Higginbottom, O. Pinel, N. P. Robins, P. K. Lam and B. C. Buchler, Nat. Phys. 13, 68-73 (2017).
  * (16) K.-K. Park, Y.-W. Cho, Y.-T. Chough, and Y.-H. Kim, Phys. Rev. X 8, 021016 (2018).
  * (17) U.-S. Kim, Y. S. Ihn, C.-H. Lee, and Y.-H. Kim, AVS Quantum Sci. 4, 021403 (2022).
  * (18) A. André, M. Bajcsy, A. S. Zibrov, and M. D. Lukin, Phys. Rev. Lett. 94, 063902 (2005).
  * (19) Y.-H. Chen, M.-J. Lee, W. Hung, Y.-C. Chen, Y.-F. Chen, and I. A. Yu, Phys. Rev. Lett. 108, 173603 (2012).
  * (20) C. R. Murray and T. Pohl, Phys. Rev. X 7, 031007 (2017).
  * (21) I. Iakoupov, J. Borregaard, and A. S. Sørenson, Phys. Rev. Lett. 120, 010502 (2018).
  * (22) D. E. Chang , V. Gritsev, G. Morigi, V. Vuletić, M. D. Lukin, and E. A. Demler, Nat. Phys. 4, 884-889 (2008).
  * (23) D. G. Angelakis, M.-X. Huo, D. Chang, L. C. Kwek, and V. Korepin, Phys. Rev. Lett. 110, 100502 (2013).
  * (24) M. Fleischhauer, J. Otterbach, and R. G. Unanyan, Phys. Rev. Lett. 101, 163601 (2008).
  * (25) H.-T. Zhou, D.-W. Wang, D. Wang, J.-X. Zhang, and S.-Y. Zhu, Phys. Rev. A 84, 053835 (2011).
  * (26) S. A. Moiseev, A. I. Sidorova, and B. S. Ham, Phys. Rev. A 89, 043802 (2014).
  * (27) J. Otterbach, J. Ruseckas, R.G. Unanyan, G. Juzeli$\mathrm{\bar{u}}$nas, and M. Fleischhauer, Phys. Rev. Lett. 104, 033903 (2010).
  * (28) M. Fleischhauer and M. D. Lukin, Phys. Rev. Lett. 84, 5094 (2000).
  * (29) Y.-F. Pu, N. Jiang, W. Chang, C. Li, and L.-M. Duan, Nat. Commun. 8, 15359 (2017).
  * (30) Y.-W. Cho, K.-K. Park, J.-C. Lee, and Y.-H. Kim, Phys. Rev. Lett. 113, 063602 (2014).
  * (31) J.-C. Lee, K.-K. Park, T.-M. Zhao, and Y.-H. Kim, Phys. Rev. Lett. 117, 250501 (2016).
  * (32) K.-K. Park, Y.-W. Cho, J.-H. Kim, T.-M. Zhao, and Y.-H. Kim, Optica 4, 1293 (2017).
  * (33) T.-M. Zhao, Y. S. Ihn, and Y.-H. Kim, Phys. Rev. Lett. 122, 123607 (2019).
  * (34) W. Ketterle, K. B. Davis, M. A. Joffe, A. Martin, and D. E. Pritchard, Phys. Rev. Lett. 70, 2253 (1993)
  * (35) Y.-F. Hsiao, H.-S. Chen, P.-J. Tsai, and Y.-C. Chen, Phys. Rev. A 90, 055401 (2014).
  * (36) J.-H. Wu, M. Artoni, and G. C. La Rocca, Phys. Rev. A 82, 013807 (2010).
  * (37) B. E. A. Saleh and M. C. Teich, Fundamentals of photonics (Wiley-Interscience, New Jersey, 2007), 2nd ed., p. 376.
  * (38) D. Z. Anderson, J. C. Frisch, and C. S. Masser, Appl. Opt. 23, 1238-1245 (1984).
  * (39) K. M. Beck, M. Hosseini, Y. Duan, and V. Vuletić, PNAS 113, 35 (2016).
  * (40) K. M. Beck, W. Chen, Q. Lin, M. Gullans, M. D. Lukin, and V. Vuletić, Phys. Rev. Lett. 113, 113603 (2014).
  * (41) M. Fleischhauer, A. Imamoglu, and J. P. Marangos, Rev. Mod. Phys. 77, 2 (2005).
  * (42) W. Chen, K. M. Beck, R. Bücker, M. Gullans, M. D. Lukin, H. T.-Suzuki, and V. Vuletić, Science 341, 768-770 (2013).
  * (43) M. Hosseini, K. M. Beck, Y. Duan, W. Chen, and V. Vuletić, Phys. Rev. Lett. 116, 033602 (2016).


#  Global Algorithms for Mean-Variance Optimization in Markov Decision
Processes

Li Xia, Shuai Ma L. Xia and S. Ma are both with the School of Business, Sun
Yat-Sen University, Guangzhou 510275, China. (email: xiali5@sysu.edu.cn).

###### Abstract

Dynamic optimization of mean and variance in Markov decision processes (MDPs)
is a long-standing challenge caused by the failure of dynamic programming. In
this paper, we propose a new approach to find the globally optimal policy for
combined metrics of steady-state mean and variance in an infinite-horizon
undiscounted MDP. By introducing the concepts of pseudo mean and pseudo
variance, we convert the original problem to a bilevel MDP problem, where the
inner one is a standard MDP optimizing pseudo mean-variance and the outer one
is a single parameter selection problem optimizing pseudo mean. We use the
sensitivity analysis of MDPs to derive the properties of this bilevel problem.
By solving inner standard MDPs for pseudo mean-variance optimization, we can
identify worse policy spaces dominated by optimal policies of the pseudo
problems. We propose an optimization algorithm which can find the globally
optimal policy by repeatedly removing worse policy spaces. The convergence and
complexity of the algorithm are studied. Another policy dominance property is
also proposed to further improve the algorithm efficiency. Numerical
experiments demonstrate the performance and efficiency of our algorithms. To
the best of our knowledge, our algorithm is the first that efficiently finds
the globally optimal policy of mean-variance optimization in MDPs. These
results are also valid for solely minimizing the variance metrics in MDPs.

Keywords: Markov decision process, mean-variance optimization, bilevel MDP,
pseudo mean, pseudo variance, global optimum

##  1 Introduction

Mean-variance optimization is an important model for the risk control in
finance engineering, which was first proposed by Markowitz (1952) for single-
period portfolio management. Extending to multi-period scenarios is a natural
but challenging research topic. This is because the variance criterion in
multi-period is not additive, which induces the time inconsistency and the
failure of dynamic programming. This important topic attracts research
attention over past decades (Dai et al., 2021; Gao and Li, 2013; Hernandez-
Lerma et al., 1999; Sobel, 1994, 1982), while it is not completely solved yet.

Since Markov models are widely used to study multi-period stochastic systems,
there is rich literature on Markov decision processes (MDPs) with variance
related criteria, either for discounted or undiscounted, discrete-time or
continuous-time, discrete-state or continuous-state, finite-horizon or
infinite-horizon MDPs. Excellent works can be referred to Chung (1994); Filar
and Lee (1985); Haskell and Jain (2013); Hernandez-Lerma et al. (1999); Sobel
(1982, 1994); Guo and Song (2009), just to name a few. Many of these works
study the variance minimization of accumulated rewards in a policy set, in
which the mean performance has already been optimized. In such scenarios, the
variance minimization problem can be equivalently converted to another
standard MDP with a new cost function (Guo et al., 2012; Huang, 2018; Sobel,
1982; Xia, 2018). These approaches are not applicable to directly optimize
variance or mean-variance combined metrics in MDPs when the mean performance
is not optimized. Another method to study the mean-variance optimization of
MDPs is to reformulate these problems as mathematical programming models and
to do further analytical investigations (Chung, 1994; Haskell and Jain, 2013;
Sobel, 1994). How to efficiently solve these mathematical programs is
challenging.

Another research stream on multi-period mean-variance optimization is from the
perspective of stochastic control. The seminal work by Li and Ng (2000); Zhou
and Li (2000) formulated the mean-variance portfolio selection problem as a
linear quadratic (LQ) control problem and used an embedding method to develop
an iterative procedure to analytically solve this problem. There are numerous
works following this research line (Gao and Li, 2013; Zhou and Yin, 2004; Zhu
et al., 2004) and interested audience can refer to a recent survey paper (Cui
et al., 2022). However, these works use an LQ model with linear state
transitions, which may properly characterize the portfolio selection problem
but lack much generalization compared with Markov models.

Recently there are also some works that study mean-variance optimization in
the regime of reinforcement learning. Although the principle of dynamic
programming fails, gradient-based algorithms for parameterized policies
(represented by neural networks) still work. Most of these studies focus on
improving the sampling efficiency for learning the gradient estimators for
variance related metrics (Borkar, 2010; Prashanth and Ghavamzadeh, 2013; Tamar
et al., 2012). A recent progress is to reformulate mean-variance optimization
with Fenchel duality (Xie et al., 2018), and to adopt gradient-based
algorithms to find local optima (Bisi et al., 2020; Zhang et al., 2021).
However, all these gradient-based learning algorithms suffer from slow
convergence speed and trap into local optima. Globally solving the mean-
variance optimization problem in MDPs is still an unanswered question.

In this paper, we study global algorithms for the mean-variance optimization
problem in an infinite-horizon discrete-time undiscounted MDP. The mean and
variance of rewards are measured in a steady-state environment, similar to
those in the works by Bisi et al. (2020); Chung (1994); Sobel (1994); Xia
(2016). By introducing an auxiliary variable called pseudo mean
$y\in\mathbb{R}$, we convert the steady-state mean-variance optimization
problem to a bilevel MDP problem, where the inner level is a standard MDP
$\mathcal{M}(y)$ optimizing the so-called pseudo mean-variance and the outer
level is a single parameter selection problem optimizing the pseudo mean $y$.
With the sensitivity analysis of MDPs, we show that the optimal value of the
pseudo mean-variance optimization problem $\mathcal{M}(y)$ is a convex
piecewise quadratic function with respect to $y$ and its global optimum equals
the optimum of the mean-variance optimization problem. We further discover
policy dominance properties which help us discard the worse policies dominated
by the optimal policy of $\mathcal{M}(y)$. Thus, the optimization complexity
can be significantly reduced. Based on these properties, we develop an
iterative algorithm which is shown to find the global optimum of the mean-
variance optimization problem after a finite number of iterations. The
computation complexity and some variants of the algorithm are also studied.
Compared with the literature work only capable of finding a local optimum of
mean-variance optimization in MDPs (Xia, 2020), our algorithms guarantee a
global convergence. The performance and efficiency of our algorithms are also
demonstrated by numerical experiments. To the best of our knowledge, our work
is the first to compute the globally optimal policies of mean-variance
optimization in MDPs.

The rest of the paper is organized as follows. In Section 2, we give the MDP
formulation for the mean-variance optimization problem. Section 3 presents the
main results of this paper, including the policy dominance property and the
algorithmic analysis. Numerical experiments are conducted in Section 4 to
demonstrate the performance of our algorithms. Finally, we conclude this paper
in Section 5.

##  2 Problem Formulation

Consider an infinite-horizon discrete-time MDP denoted by a tuple
$\mathcal{M}=\langle\mathcal{S},\mathcal{A},\mathcal{P},\bm{r}\rangle$, where
$\mathcal{S}=\\{1,2,\dots,S\\}$ is the state space,
$\mathcal{A}=\\{a_{1},a_{2},\dots,a_{A}\\}$ is the action space,
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\overset{D}{\mapsto}\mathcal{S}$ is
the state transition probability kernel with element $p(j|i,a)$ where
$\overset{D}{\mapsto}$ represents a mapping to the distribution on
$\mathcal{S}$, and $\bm{r}:\mathcal{S}\times\mathcal{A}\mapsto\mathbb{R}$ is
the reward function with element $r(i,a)$, $i,j\in\mathcal{S}$,
$a\in\mathcal{A}$. When the system is in state $i$ and action $a$ is adopted,
it will transit to the next state $j$ with probability $p(j|i,a)$ and a reward
$r(i,a)$ is incurred. Since deterministic policies can attain optimal mean and
variance in MDPs (Haskell and Jain, 2013; Xia, 2020), we only consider
stationary deterministic policies $d:\mathcal{S}\mapsto\mathcal{A}$, where
$d(i)\in\mathcal{A}$ indicates the action adopted in state $i$. The
corresponding policy space is denoted by $\mathcal{D}$ and we assume that the
MDP with any policy $d\in\mathcal{D}$ is a unichain. When a policy $d$ is
adopted, the state transition probability matrix is denoted by $\bm{P}^{d}$
and its $(i,j)$-th element is $p(j|i,d(i))$, $i,j\in\mathcal{S}$. The
associated steady-state distribution is denoted by an $S$-dimensional row
vector $\bm{\pi}^{d}:=(\pi^{d}(i))_{i\in\mathcal{S}}$. Obviously, we have

$\bm{\pi}^{d}\bm{P}^{d}=\bm{\pi}^{d},\quad\bm{\pi}^{d}\bm{e}=1,\quad\bm{P}^{d}\bm{e}=\bm{e},$

where $\bm{e}$ is a column vector of 1’s with a proper dimension size. We
consider long-run performance metrics of this MDP, which are independent of
the initial state at time 0. The _long-run average (mean)_ reward of the MDP
under policy $d$ is defined as

$\mu^{d}:=\lim\limits_{T\rightarrow\infty}\mathbb{E}^{d}\left\\{\frac{1}{T}\sum_{t=0}^{T-1}r(X_{t},A_{t})\right\\}=\bm{\pi}^{d}\bm{r}^{d},$
(1)

where $\mathbb{E}^{d}$ indicates the expectation under policy $d$, $X_{t}$ is
the system state at time $t$, $A_{t}=d(X_{t})$ is the action adopted at time
$t$, $\bm{r}^{d}$ is an $S$-dimensional column vector whose element is
$r(i,d(i))$, $i\in\mathcal{S}$. Similarly, the _long-run variance_ (or steady-
state variance) of the MDP under policy $d$ is defined as (Xia, 2016, 2020)

$\sigma^{d}:=\lim\limits_{T\rightarrow\infty}\mathbb{E}^{d}\left\\{\frac{1}{T}\sum_{t=0}^{T-1}[r(X_{t},A_{t})-\mu^{d}]^{2}\right\\}=\bm{\pi}^{d}(\bm{r}^{d}-\mu^{d}\bm{e})_{\odot}^{2},$
(2)

where $(\bm{r}^{d}-\mu^{d}\bm{e})_{\odot}^{2}$ is the component-wise square of
vector $(\bm{r}^{d}-\mu^{d}\bm{e})$, i.e.,

$(\bm{r}^{d}-\mu^{d}\bm{e})_{\odot}^{2}:=((r(1,d(1))-\mu^{d})^{2},\
(r(2,d(2))-\mu^{d})^{2},\ \dots,\ (r(S,d(S))-\mu^{d})^{2})^{\mathrm{T}}.$

When the finite Markov chain is a unichain, we can view $r(X_{t},A_{t})$ as a
random variable whose value realization set is
$\\{r(i,d(i)):i\in\mathcal{S}\\}$ and distribution is
$(\bm{P}^{d})^{t}\bm{\nu}$, where $\bm{\nu}$ is the vector of initial state
distribution. We can verify that

$\displaystyle\mu^{d}$ $\displaystyle=$
$\displaystyle\lim\limits_{t\rightarrow\infty}\mathbb{E}[r(X_{t},A_{t})],$
$\displaystyle\sigma^{d}$ $\displaystyle=$
$\displaystyle\lim\limits_{t\rightarrow\infty}\mathbb{VAR}[r(X_{t},A_{t})].$

Mean-variance optimization was originally proposed by Markowitz (1952) for
portfolio selection, where decision makers usually aim at maximizing the mean
return while minimizing the variance risk, which is a multi-objective
optimization problem. Usually, the Pareto frontier composed of Pareto
efficient solutions is the optimization goal, which is illustrated by Fig. 1.
A common way of obtaining Pareto optima is to optimize the combined objective

$\eta^{d}:=\beta\sigma^{d}-\mu^{d},$ (3)

where $\beta\geq 0$ is the tradeoff weight between mean and variance.
Therefore, our goal is to solve the _mean-variance optimization_ problem:

$\mbox{(P0): \hskip 113.81102pt
}\begin{array}[]{cc}&\eta^{*}=\min\limits_{d\in\mathcal{D}}\\{\beta\sigma^{d}-\mu^{d}\\},\\\
&d^{*}\in\operatorname*{argmin}\limits_{d\in\mathcal{D}}\\{\beta\sigma^{d}-\mu^{d}\\}.\end{array}$
(4)

Note that $\eta^{*}$ and $d^{*}$ depend on $\beta$, and we may also use
$\eta^{*}(\beta)$ and $d^{*}(\beta)$ if necessary. In Fig. 1, the red star
points are Pareto efficient solutions which dominate the black dot solutions.
The dashed curve is the Pareto frontier which can be obtained by solving (4)
with different $\beta\geq 0$. We can also observe that the dashed line is
tangent to the Pareto frontier, where the slope is $\beta$ and the tangent
point is $(\sigma^{d^{*}(\beta)},\mu^{d^{*}(\beta)})$.

Figure 1: Illustration of the Pareto frontier in mean-variance optimization.

How to efficiently solve (4) is the key of the mean-variance optimization in
MDPs. Since the variance function $(r(i,d(i))-\mu^{d})^{2}$ depends on history
and future behaviors through $\mu^{d}$, it is not either additive or
Markovian. The mean-variance optimization problem (4) does not fit a standard
model of MDPs and the principle of dynamic programming fails (Chung, 1994;
Sobel, 1994; Xia, 2016). Although there is a recent progress on this problem
by using the technique of sensitivity-based optimization instead of the
traditional dynamic programming (Xia, 2020), it can only find a local optimum
of this mean-variance optimization problem. A local optimum is not guaranteed
as a Pareto efficient solution. Thus, finding the global optimum of (4) is
still an unsolved problem in the mean-variance optimization of MDPs and we
accomplish this challenge in the rest of this paper.

##  3 Main Results

First, we introduce the concept of _pseudo mean_ and _pseudo variance_ of an
MDP under policy $d\in\mathcal{D}$ (Xia, 2016, 2020):

$\tilde{\sigma}^{d}(y)=\bm{\pi}^{d}(\bm{r}^{d}-y\bm{e})^{2}_{\odot},\qquad
y\in\mathbb{R},$ (5)

where $\tilde{\sigma}^{d}(y)$ is called the pseudo variance of the MDP with
the pseudo mean $y$. We can derive that the difference between the pseudo
variance and the real variance is

$\displaystyle\Delta^{d}(y)$ $\displaystyle:=$
$\displaystyle\tilde{\sigma}^{d}(y)-\sigma^{d}=\bm{\pi}^{d}(\bm{r}^{d}-y\bm{e})^{2}_{\odot}-\bm{\pi}^{d}(\bm{r}^{d}-\mu^{d}\bm{e})^{2}_{\odot}$
(6) $\displaystyle=$
$\displaystyle\sum_{i\in\mathcal{S}}\pi^{d}(i)[(r(i,d(i))-y)^{2}-(r(i,d(i))-\mu^{d})^{2}]$
$\displaystyle=$ $\displaystyle(y-\mu^{d})^{2}\geq 0.$

We call $\Delta^{d}(y)$ the _variance distortion_ caused by the pseudo mean
$y$. Interestingly, we observe that the pseudo variance
$\tilde{\sigma}^{d}(y)$ is a _convex quadratic_ function of $y$, since
$\tilde{\sigma}^{d}(y)=\sigma^{d}+(y-\mu^{d})^{2}$. When the pseudo mean $y$
equals the real mean $\mu^{d}$, the variance distortion is zero and the pseudo
variance attains its minimum which is exactly the real variance, i.e.,

$\sigma^{d}=\min_{y\in\mathbb{R}}\tilde{\sigma}^{d}(y)=\tilde{\sigma}^{d}(y^{*})\Big{|}_{y^{*}=\mu^{d}}.$
(7)

Remark 1. The above property of variance is analogous to CVaR (Conditional
Value at Risk) discovered by Rockafellar and Uryasev (2000): the CVaR of
random variable $X$ at probability level $\alpha$ is defined as
$\mathrm{CVaR}_{\alpha}(X):=\mathbb{E}[X|X\geq F^{-1}_{X}(\alpha)]$, and
equals
$\min\limits_{y\in\mathbb{R}}\mathbb{E}[y+\frac{1}{1-\alpha}[X-y]^{+}]$, where
$F^{-1}_{X}(\cdot)$ is the inverse distribution function of $X$,
$[X-y]^{+}:=\max\\{0,X-y\\}$, $\mathbb{E}[y+\frac{1}{1-\alpha}[X-y]^{+}]$ is a
convex function of $y$ and its minimum attains at $y^{*}=F^{-1}_{X}(\alpha)$.

With this property (7), we can convert the original mean-variance optimization
problem to a bilevel MDP problem and directly derive the following lemma.

######  Lemma 1 (Bilevel MDP).

The mean-variance optimization problem (4) is equivalent to a bilevel MDP
problem where the inner one is a standard MDP with cost function
$\beta(\bm{r}-y\bm{e})^{2}_{\odot}-\bm{r}$:

$\eta^{*}=\min_{d\in\mathcal{D}}\\{\beta\sigma^{d}-\mu^{d}\\}=\min_{y\in\mathbb{R}}\min_{d\in\mathcal{D}}\\{\beta\tilde{\sigma}^{d}(y)-\mu^{d}\\}.$
(8)

The above bilevel MDP formulation is similar to the Fenchel duality
formulation (Xie et al., 2018), while our formulation (8) naturally comes from
(6) of pseudo variance which was originally discovered by Xia (2016). The
inner problem
$\min\limits_{d\in\mathcal{D}}\\{\beta\tilde{\sigma}^{d}(y)-\mu^{d}\\}$ aims
to optimize the pseudo mean-variance, which is a standard MDP denoted by tuple
$\mathcal{M}(y):=\langle\mathcal{S},\mathcal{A},\mathcal{P},\beta(\bm{r}-y\bm{e})^{2}_{\odot}-\bm{r}\rangle$.
We can use traditional dynamic programming to solve this MDP. For different
outer variable of pseudo mean $y$, we have to solve different MDP
$\mathcal{M}(y)$. The number of solving inner MDPs is equal to the number of
$y\in\mathbb{R}$, which is computationally intractable. Therefore, efficiently
solving this bilevel MDP problem (8) is challenging.

With (7), we see that the optimal $y^{*}$ in (8) satisfies $y^{*}=\mu^{d^{*}}$
for the optimal policy $d^{*}\in\mathcal{D}$. Therefore, we can restrict $y$’s
value domain from $y\in\mathbb{R}$ to a much smaller set $y\in\mathcal{Y}$,
where $\mathcal{Y}:=\\{\mu^{d}:\forall d\in\mathcal{D}\\}$. Although
$\mathcal{Y}$ is still computationally intractable, we know that
$\mathcal{Y}\subset[\underline{r},\overline{r}]$, where
$\underline{r}:=\min\limits_{i\in\mathcal{S},a\in\mathcal{A}}\\{r(i,a)\\}$ and
$\overline{r}:=\max\limits_{i\in\mathcal{S},a\in\mathcal{A}}\\{r(i,a)\\}$.
Therefore, the bilevel MDP problem (8) can be rewritten as

$\eta^{*}=\min_{d\in\mathcal{D}}\\{\beta\sigma^{d}-\mu^{d}\\}=\min_{y\in[\underline{r},\overline{r}]}\min_{d\in\mathcal{D}}\\{\beta\tilde{\sigma}^{d}(y)-\mu^{d}\\}.$
(9)

Since the inner problem $\mathcal{M}(y^{*})$ is a standard MDP, we can derive
a concise proof about the optimality of deterministic policies (detailed
proofs can also be referred to Haskell and Jain (2013); Xia (2020)): Suppose
$(y^{*},d^{*})$ is an optimal solution to (9). It is well known that there
exists a deterministic policy $d_{0}$ which attains the minimum of standard
MDP $\min_{d}\\{\beta\tilde{\sigma}^{d}(y^{*})-\mu^{d}\\}$. It is obvious that
$y^{*}$ must be the real mean of the MDP under policy $d_{0}$. Thus,
$\beta\tilde{\sigma}^{d_{0}}(y^{*})-\mu^{d_{0}}=\beta\sigma^{d_{0}}-\mu^{d_{0}}$,
which indicates that the deterministic policy $d_{0}$ attains the minimum of
mean-variance performance.

When the pseudo mean $y$ is fixed, the inner standard MDP $\mathcal{M}(y)$ is
an auxiliary problem, and its long-run average performance under policy $d$ is
a combined performance
$\tilde{\eta}^{d}(y)=\beta\tilde{\sigma}^{d}(y)-\mu^{d}$. We also call
$\mathcal{M}(y)$ a _pseudo mean-variance optimization_ problem:

$(\mathcal{M}(y)):\hskip
85.35826pt\begin{array}[]{cc}&\tilde{\eta}^{*}(y)=\min\limits_{d\in\mathcal{D}}\\{\beta\tilde{\sigma}^{d}(y)-\mu^{d}\\}.\\\
&\tilde{d}^{*}(y)\in\operatorname*{argmin}\limits_{d\in\mathcal{D}}\\{\beta\tilde{\sigma}^{d}(y)-\mu^{d}\\}.\end{array}$
(10)

For notation simplicity, sometimes we may omit $y$, and use $\tilde{\eta}^{*}$
and $\tilde{d}^{*}$ if no confusion caused. Therefore, the bilevel MDP (9) for
mean-variance optimization can be rewritten as below.

$\eta^{*}=\min_{y\in[\underline{r},\overline{r}]}\\{\tilde{\eta}^{*}(y)\\}.$

If we plot a curve of $\tilde{\eta}^{*}(y)$ with respect to $y$, we can
observe that $\eta^{*}$ is the global minimum of this curve at point $y^{*}$
and the corresponding $\tilde{d}^{*}(y^{*})$ is the optimal policy of the
original problem (4). From the sensitivity analysis of MDPs, we can derive the
following lemma.

######  Lemma 2 (Critical points).

There exists a series of intervals $[y^{c}_{k-1},y^{c}_{k}]$ with
$\bigcup\limits_{k=1,\dots,K}[y^{c}_{k-1},y^{c}_{k}]=[\underline{r},\overline{r}]$,
in which the optimal policy of $\mathcal{M}(y)$ remains unvaried as
$\tilde{d}^{*}_{k}:=\tilde{d}^{*}(y)$ when $y\in[y^{c}_{k-1},y^{c}_{k}]$.

###### Proof.

We rewrite the standard MDP problem $\mathcal{M}(y)$ as a linear programming
(LP) model:

$\begin{array}[]{rcl}\tilde{\eta}^{*}(y)&=&\min\limits_{\bm{x}}\left\\{\sum\limits_{i\in
S}\sum\limits_{a\in\mathcal{A}}[\beta(r(i,a)-y)^{2}-r(i,a)]x(i,a)\right\\}\\\
&\mbox{s.t., }&\sum\limits_{a\in\mathcal{A}}x(i,a)=\sum\limits_{j\in
S}\sum\limits_{a\in\mathcal{A}}p(i|j,a)x(j,a),\qquad\forall
i\in\mathcal{S},\\\ &&\sum\limits_{i\in
S}\sum\limits_{a\in\mathcal{A}}x(i,a)=1,\\\ &&x(i,a)\geq 0,\qquad\forall
i\in\mathcal{S},a\in\mathcal{A}.\end{array}$ (11)

The above problem can be represented as a standard LP model:

$\tilde{\eta}^{*}(y)-\beta
y^{2}=\min_{\bm{x}}\\{(\bm{c}+y\bm{c}^{\prime})^{\mathrm{T}}\bm{x}|\bm{A}\bm{x}=\bm{b},\bm{x}\geq\bm{0}\\},$
(12)

where we utilize the fact $\sum\limits_{i\in\mathcal{S},a\in\mathcal{A}}\beta
y^{2}x(i,a)=\beta y^{2}$, the $S$-by-$SA$ matrix $\bm{A}$ and the
$S$-dimensional column vector $\bm{b}$ are determined by the constraint
equations in (11), $\bm{c}=\beta\bm{r}^{2}_{\odot}-\bm{r}$,
$\bm{c}^{\prime}=-2\beta\bm{r}$, $\bm{r}$ and $\bm{x}$ are $SA$-dimensional
column vector with element $r(i,a)$ and $x(i,a)$, respectively. We observe
that the right-hand-side of (12) is a _parametric linear programming_ (PLP)
(Gal and Greenberg, 1997; Tan and Hartman, 2011) with a linear parameter $y$.
Below we do sensitivity analysis for this PLP problem. For a given $y$,
suppose $\bm{x}^{*}_{k}$ is the optimal solution of (12) and its associated
basis matrix is $\bm{B}$. We can verify that $\bm{x}_{k}^{*}$ in this LP is
equivalent to the optimal policy $\tilde{d}^{*}_{k}$ of the MDP
$\mathcal{M}(y)$, where the optimal action in state $i$ is
$\tilde{d}^{*}_{k}(i)\in\operatorname*{argmax}\limits_{a\in\mathcal{A}}\\{x^{*}_{k}(i,a)\\}$,
$i\in\mathcal{S}$. With the terminology of LP, we denote
$\bm{A}=[\bm{B},\bm{N}]$, $\bm{x}=[\bm{x}_{B};\bm{x}_{N}]$,
$\bm{c}=[\bm{c}_{B};\bm{c}_{N}]$, and
$\bm{c}^{\prime}=[\bm{c}^{\prime}_{B};\bm{c}^{\prime}_{N}]$. The optimality
test of the simplex method requires that all the test coefficients should be
nonpositive, i.e.,

$(\bm{c}_{B}+y\bm{c}^{\prime}_{B})^{\mathrm{T}}\bm{B}^{-1}\bm{A}-(\bm{c}+y\bm{c}^{\prime})^{\mathrm{T}}=(\bm{c}_{B}^{\mathrm{T}}\bm{B}^{-1}\bm{A}-\bm{c}^{\mathrm{T}})+y({\bm{c}^{\prime}_{B}}^{\mathrm{T}}\bm{B}^{-1}\bm{A}-{\bm{c}^{\prime}}^{\mathrm{T}})\leq\bm{0}.$

For notation simplicity, we denote the $SA$-dimensional test coefficients
vector as

$\displaystyle\bm{\zeta}^{\mathrm{T}}$ $\displaystyle:=$
$\displaystyle\bm{c}_{B}^{\mathrm{T}}\bm{B}^{-1}\bm{A}-\bm{c}^{\mathrm{T}},$
(13) $\displaystyle{\bm{\zeta}^{\prime}}^{\mathrm{T}}$ $\displaystyle:=$
$\displaystyle{\bm{c}^{\prime}_{B}}^{\mathrm{T}}\bm{B}^{-1}\bm{A}-{\bm{c}^{\prime}}^{\mathrm{T}}.$
(14)

In order to find the interval $[y^{c}_{k-1},y^{c}_{k}]$ that any $y$ therein
makes $\bm{x}^{*}_{k}$ remain optimal, we only need to solve $y$ satisfying
$\bm{\zeta}+y\bm{\zeta}^{\prime}\leq\bm{0}$. It is easy to verify that the
solution is

$\displaystyle
y^{c}_{k-1}=\max_{i\in\mathcal{S},a\in\mathcal{A}}\left\\{-\frac{\zeta(i,a)}{\zeta^{\prime}(i,a)}\Bigg{|}\zeta^{\prime}(i,a)<0\right\\}\qquad(\max\\{\varnothing\\}=-\infty),$
(15) $\displaystyle
y^{c}_{k}=\min_{i\in\mathcal{S},a\in\mathcal{A}}\left\\{-\frac{\zeta(i,a)}{\zeta^{\prime}(i,a)}\Bigg{|}\zeta^{\prime}(i,a)>0\right\\}\qquad(\min\\{\varnothing\\}=+\infty).\
$ (16)

Obviously, we can first let $y=\underline{r}$, and use (15) and (16) to obtain
$y^{c}_{0}=-\infty$ and $y^{c}_{1}$, respectively. Other $y^{c}_{k}$’s can be
computed sequentially, and $y^{c}_{K}=+\infty$. The lemma is proved. ∎

We call such $y^{c}_{k}$’s in Lemma 2 _critical points_ for the sensitivity
analysis of MDP $\mathcal{M}(y)$, and $K+1$ is the number of critical points.
Actually, by using the specific structures of $\bm{A},\bm{b},\bm{r}$ of the LP
for $\mathcal{M}(y)$, we can verify that

$[\bm{B}^{-1}]^{\mathrm{T}}=(\bm{I}-\bm{P}^{d}+\bm{e}\bm{e}^{\mathrm{T}})^{-1},$
(17)

which is a generalized _fundamental matrix_ in MDPs (Xia and Glynn, 2016),
where the policy $d$ corresponds to the vector of feasible basic variables
$\bm{x}_{B}$ of the basis matrix $\bm{B}$. The associated vector $\bm{b}$
equals $\bm{e}$. The matrix $\bm{A}$ has a similar structure of $\bm{B}$,
i.e., $\bm{A}=\bm{I}_{e}-\bm{P}^{\mathrm{T}}_{e}+\bm{e}\bm{e}^{\mathrm{T}}$,
where $\bm{I}_{e}$ is an $S$-by-$SA$ matrix whose element $I_{e}(i,(j,a))=1/0$
if $i=j$/otherwise, $\bm{P}_{e}$ is an $SA$-by-$S$ matrix whose element
$P_{e}((i,a),j)=p(j|i,a)$, $\bm{e}\bm{e}^{\mathrm{T}}$ is an $S$-by-$SA$
matrix of 1’s. We also observe that $\bm{c}_{B}$ is the column vector of the
cost function of the MDP under policy $d$ (associated with $\bm{x}_{B}$). We
can derive that $\bm{c}^{\mathrm{T}}_{B}\bm{B}^{-1}$ is equal to the
_performance potential_ or _relative value function_ in MDPs (Cao, 2007;
Puterman, 1994). Thus, we can verify that $\bm{c}_{B}^{\mathrm{T}}\bm{B}^{-1}$
in (13) coincides with the _Poisson equation_ in MDPs

$\displaystyle\bm{g}^{d}$ $\displaystyle=$
$\displaystyle[\bm{B}^{-1}]^{\mathrm{T}}\bm{c}_{B}=(\bm{I}-\bm{P}^{d}+\bm{e}\bm{e}^{\mathrm{T}})^{-1}(\beta{\bm{r}^{d}}^{2}_{\odot}-\bm{r}^{d}),$
$\displaystyle{\bm{g}^{\prime}}^{d}$ $\displaystyle=$
$\displaystyle[\bm{B}^{-1}]^{\mathrm{T}}\bm{c}^{\prime}_{B}=-2\beta(\bm{I}-\bm{P}^{d}+\bm{e}\bm{e}^{\mathrm{T}})^{-1}\bm{r}^{d},$

where $\bm{g}^{d}$ and ${\bm{g}^{\prime}}^{d}$ are $S$-dimensional column
vector of performance potentials for the MDP under policy $d$ with cost
function $\beta{\bm{r}^{d}}^{2}_{\odot}-\bm{r}^{d}$ and $-2\beta\bm{r}^{d}$,
respectively. Thus, we can rewrite the element of (13) and (14) as

$\displaystyle\zeta(i,a)$ $\displaystyle:=$ $\displaystyle
g^{d}(i)-\sum_{j\in\mathcal{S}}p(j|i,a)g^{d}(j)+\sum_{j\in\mathcal{S}}g^{d}(j)+r(i,a)-\beta{r(i,a)}^{2}=:-\tilde{A}^{d}(i,a),$
(18) $\displaystyle\zeta^{\prime}(i,a)$ $\displaystyle:=$
$\displaystyle{g^{\prime}}^{d}(i)-\sum_{j\in\mathcal{S}}p(j|i,a){g^{\prime}}^{d}(j)+\sum_{j\in\mathcal{S}}{g^{\prime}}^{d}(j)+2\beta
r(i,a)=:-\tilde{A^{\prime}}^{d}(i,a),$ (19)

where we utilize a fact that $\bm{e}^{\mathrm{T}}\bm{g}^{d}$ equals the long-
run average performance of the MDP under policy $d$, which can be verified
from the Poisson equation. Therefore, we can substitute (18) and (19) into
(15) and (16) to compute all the critical points $y^{c}_{k}$’s.

Remark 2. Equations (18) and (19) interestingly indicate that the test
coefficient $\zeta(i,a)$ in LP coincides with the _advantage function_
$\tilde{A}^{d}(i,a)$ which is a key quantity widely used in reinforcement
learning (Sutton and Barto, 2018). $\tilde{A}^{d}(i,a)<0$ means that action
$a$ at state $i$ induces a smaller long-run average cost than the current
policy $d$ in MDPs, which hints that $\zeta(i,a)>0$ and the corresponding
variable $x(i,a)$ should be an entering basic variable in LP.

By using (6), we can obtain the relation between the pseudo and real mean-
variance combined performances as

$\tilde{\eta}^{d}(y)=\beta[\sigma^{d}+(y-\mu^{d})^{2}]-\mu^{d}=\eta^{d}+\beta(y-\mu^{d})^{2},\quad\forall
d\in\mathcal{D},y\in\mathbb{R}.$

Since the optimal policy of $\mathcal{M}(y)$ remains the same as
$\tilde{d}^{*}_{k}$ for any $y\in[y^{c}_{k-1},y^{c}_{k}]$, we have

$\begin{array}[]{rcll}\tilde{\eta}^{*}(y)&=&\eta^{\tilde{d}^{*}_{k}}+\beta(y-\mu^{\tilde{d}^{*}_{k}})^{2},&\forall
y\in[y^{c}_{k-1},y^{c}_{k}],\\\
\tilde{\eta}^{*}(y)&=&\min\limits_{k=1,\dots,K}\left\\{\eta^{\tilde{d}^{*}_{k}}+\beta(y-\mu^{\tilde{d}^{*}_{k}})^{2}\right\\},&\forall
y\in\mathbb{R}.\end{array}$ (20)

That is, each piece of curves in Fig. 2 is a convex quadratic function of $y$,
and the whole curve is the minimum of all these quadratic functions. With
(20), it is interesting to note that all the piecewise curves have the same
shape (the same term of $\beta y^{2}$) but different positions. At each
critical point $y^{c}_{k}$, we can validate that both $\tilde{d}^{*}_{k}$ and
$\tilde{d}^{*}_{k+1}$ are optimal policies of MDP $\mathcal{M}(y^{c}_{k})$, so
$\tilde{\eta}^{*}(y)$ is continuous in $y$. Therefore, we directly derive the
following lemma.

######  Lemma 3.

The pseudo mean-variance performance $\tilde{\eta}^{*}(y)$ is a convex
piecewise quadratic function and continuous in $y$, and its global minimum is
the optimal solution of (4).

Figure 2: Illustration of the convex piecewise quadratic curve of the pseudo
mean-variance performance $\tilde{\eta}^{*}(y)$ with respect to the pseudo
mean $y$.

From Fig. 2, we can observe that $\min\limits_{y}\\{\tilde{\eta}^{*}(y)\\}$ is
difficult to solve, because $\tilde{\eta}^{*}(y)$ may have multiple local
optima $\hat{y}^{*}$ which has a zero derivative, i.e.,

$\frac{\partial\tilde{\eta}^{*}(y)}{\partial
y}\Big{|}_{\hat{y}^{*}}=2\beta(\hat{y}^{*}-\mu^{\tilde{d}^{*}_{k}})=0,\qquad\mbox{if
}\hat{y}^{*}\in[y^{c}_{k-1},y^{c}_{k}].$

This indicates that a local optimum $\hat{y}^{*}$ must satisfy the following
_fixed point equation_

$y=\mu^{\tilde{d}^{*}(y)}.$ (21)

Note that the fixed point solutions of (21), as indicated by “inverted”
triangles in Fig. 2, are not necessarily local optima of
$\tilde{\eta}^{*}(y)$. The reason is when a fixed point is also a critical
point $y^{c}_{k}$, we can verify that its left-derivative is 0 (or positive)
and its right-derivative is negative (or 0), and such point is not a local
optimum $\hat{y}^{*}$ of $\tilde{\eta}^{*}(y)$, as illustrated by Fig. 2. We
can verify that the policies indicated by all these fixed point solutions of
(21) are exactly the so-called local optimal policies in mixed or randomized
policy space of MDPs, as discovered by Xia (2020). These two kinds of local
optima are different: local optima $\hat{y}^{*}$ in this paper are included by
local optimal policies (fixed point solutions) defined in Xia (2020), as
illustrated by triangles and “inverted” triangles in Fig. 2, respectively.

Note that $\tilde{d}^{*}$ is optimal only for the pseudo problem (10), not for
the original problem (4). Fortunately, we discover that $\tilde{d}^{*}$ has a
better mean-variance performance than some other policies, which is described
by the following lemma.

######  Lemma 4 (Policy dominance).

For any $y\in\mathbb{R}$, $\tilde{d}^{*}$ is an optimal policy of the MDP
$\mathcal{M}(y)$ in (10). If a policy $d\in\mathcal{D}$ satisfying
$\mu^{d}\in[y-|y-\mu^{\tilde{d}^{*}}|,\ y+|y-\mu^{\tilde{d}^{*}}|]$, then

$\beta\sigma^{\tilde{d}^{*}}-\mu^{\tilde{d}^{*}}\leq\beta\sigma^{d}-\mu^{d}.$
(22)

###### Proof.

Since $\tilde{d}^{*}$ is an optimal policy of the standard MDP problem
$\mathcal{M}(y)$, we have

$\beta\tilde{\sigma}^{\tilde{d}^{*}}(y)-\mu^{\tilde{d}^{*}}\leq\beta\tilde{\sigma}^{d}(y)-\mu^{d},\quad\forall
d\in\mathcal{D}.$ (23)

With (6), we derive

$\tilde{\sigma}^{d}(y)=\sigma^{d}+(y-\mu^{d})^{2},\quad\forall
d\in\mathcal{D}.$

Substituting the above equation into (23), we have

$\beta\sigma^{\tilde{d}^{*}}-\mu^{\tilde{d}^{*}}+\beta(y-\mu^{\tilde{d}^{*}})^{2}\leq\beta\sigma^{d}-\mu^{d}+\beta(y-\mu^{d})^{2},\quad\forall
d\in\mathcal{D}.$ (24)

For any policy $d$ satisfying $\mu^{d}\in[y-|y-\mu^{\tilde{d}^{*}}|,\
y+|y-\mu^{\tilde{d}^{*}}|]$, we have

$(y-\mu^{\tilde{d}^{*}})^{2}\geq(y-\mu^{d})^{2}.$ (25)

Substituting (25) into (24), we directly obtain (22), and the lemma is proved.
∎

Moreover, if $d$ satisfies $\mu^{d}\in(y-|y-\mu^{\tilde{d}^{*}}|,\
y+|y-\mu^{\tilde{d}^{*}}|)$, we have
$(y-\mu^{\tilde{d}^{*}})^{2}>(y-\mu^{d})^{2}$, and the inequality in (22)
strictly holds. Therefore, Lemma 4 indicates that $\tilde{d}^{*}$ dominates
all the policies whose means lie in the interval $[y-|y-\mu^{\tilde{d}^{*}}|,\
y+|y-\mu^{\tilde{d}^{*}}|]$, and these dominated policies can be removed from
the policy space $\mathcal{D}$ to save computation. We illustrate this
property by Fig. 3, where we can see that the shadow area can be discarded
since the policies therein are always dominated by $\tilde{d}^{*}$. Thus, we
can utilize Lemma 4 to significantly reduce the complexity of the mean-
variance problem (4).

Figure 3: Illustration of the dominated policy area indicated by Lemma 4.

With Lemma 4, we can develop an algorithm to iteratively solve the bilevel MDP
problem (8), which is described by Algorithm 1. The key idea is to solve a
series of auxiliary problems $\mathcal{M}(y)$’s, and repeatedly reduce the
auxiliary variable $y$’s value domain $\mathbb{Y}$ by using Lemma 4. When
$\mathbb{Y}$ is shrunk to be empty, the best-so-far solution of
$\mathcal{M}(y)$’s is the global optimum of the bilevel MDP problem (8) or
(9). The global convergence of Algorithm 1 is guaranteed by the following
theorem.

Algorithm 1 An iterative algorithm to find globally optimal policies of mean-
variance MDPs.

1:Input: MDP parameters
$\langle\mathcal{S},\mathcal{A},\mathcal{P},\bm{r}\rangle$ and coefficient
$\beta$.

2:Output: The globally optimal policy $d^{*}$ and its optimum $\eta^{*}$.

3: _initialization_ : set $\mathbb{Y}=[\underline{r},\overline{r}]$, $l=0$,
and $\eta^{*}=+\infty$.

4:while $\mathbb{Y}\neq\varnothing$ do

5: _set pseudo mean_ : $\mathbb{Y}$ is composed of multiple intervals,
descendingly sort these intervals as
$\mathbb{Y}=\\{\mathbb{Y}_{1},\mathbb{Y}_{2},\dots\\}$, set $y_{l}$ as the
median value of the first interval of $\mathbb{Y}$, i.e.,
$y_{l}=(\max\\{\mathbb{Y}_{1}\\}+\min\\{\mathbb{Y}_{1}\\})/2$.

6: _solve auxiliary problem_ : solve the standard MDP $\mathcal{M}(y_{l})$ by
using traditional dynamic programming methods, such as policy iteration or
value iteration, and obtain an optimal policy $\tilde{d}^{*}(y_{l})$ of
$\mathcal{M}(y_{l})$.

7: _domain shrinking_ : use Lemma 4 to remove the interval
$[y_{l}-|y_{l}-\mu^{\tilde{d}^{*}(y_{l})}|,\
y_{l}+|y_{l}-\mu^{\tilde{d}^{*}(y_{l})}|]$ from the value domain $\mathbb{Y}$,
i.e., $\mathbb{Y}=\mathbb{Y}-[y_{l}-|y_{l}-\mu^{\tilde{d}^{*}(y_{l})}|,\
y_{l}+|y_{l}-\mu^{\tilde{d}^{*}(y_{l})}|]$.

8: _update parameters_ :

9: if $\eta^{*}>\beta\sigma^{\tilde{d}^{*}(y_{l})}-\mu^{\tilde{d}^{*}(y_{l})}$
then

10: $\eta^{*}=\beta\sigma^{\tilde{d}^{*}(y_{l})}-\mu^{\tilde{d}^{*}(y_{l})}$.

11: $d^{*}=\tilde{d}^{*}(y_{l})$.

12: end if

13: set $l:=l+1$.

14:end whilereturn $d^{*}$ and $\eta^{*}$.

######  Theorem 1.

Algorithm 1 converges to the global optimum of the mean-variance optimization
problem after a finite number of iterations.

###### Proof.

To prove the convergence of Algorithm 1, we only need to prove that the value
domain $\mathbb{Y}$ is reduced to an empty set after a finite number of
iterations. From the algorithm procedure, we can observe that for each
iteration of solving an auxiliary problem $\mathcal{M}(y_{l})$, we will derive
a policy $\tilde{d}^{*}(y_{l})$ and remove a square area with y-axis interval
$[y_{l}-|y_{l}-\mu^{\tilde{d}^{*}(y_{l})}|,\
y_{l}+|y_{l}-\mu^{\tilde{d}^{*}(y_{l})}|]$, as stated by Lemma 4. This removed
area at least contains the policy $\tilde{d}^{*}(y_{l})$, as illustrated by
the 1st and 2nd iterations in Fig. 4. Usually, it contains multiple policies
dominated by the policy $\tilde{d}^{*}(y_{l})$, as illustrated by Fig. 3. If
the current policy $\tilde{d}^{*}(y_{l})$ has already been removed by previous
domain shrinking operations, the current domain shrinking operation will
remove at least the interval $\mathbb{Y}_{1}$, as illustrated by the 3rd, 4th,
and 5th iterations in Fig. 4, which can be verified by the fact of $y_{l}$
being the median value of $\mathbb{Y}_{1}$ and Lemma 4. In summary, each
domain shrinking operation will either delete at least one policy (not deleted
previously) or delete at least one interval $\mathbb{Y}_{1}$. It is easy to
verify that the largest number of intervals is $|\mathcal{D}|+1$, where each
iteration only deletes a very small area around $\tilde{d}^{*}(y_{l})$.
Therefore, in the worst case, we need $|\mathcal{D}|$ iterations to delete
every policy and $|\mathcal{D}|+1$ iterations to delete every interval. The
algorithm will stop after at most $2|\mathcal{D}|+1$ iterations.

Since each $\tilde{d}^{*}(y_{l})$ dominates all the other policies located in
the shrunk area of the $l$-th iteration, the best-so-far solution among all
$\tilde{d}^{*}(y_{l})$’s is the global optimum of the original mean-variance
optimization problem. This completes the proof. ∎

Figure 4: A worst-case illustration of domain shrinking procedure of Algorithm
1 for an example of a policy space with only two solutions.

Fig. 4 gives an illustration of the worst case for an example of a policy
space with only 2 solutions, it requires $2\times 2+1=5$ iterations to cover
the whole interval $[\underline{r},\overline{r}]$. From the proof of Theorem
1, we directly derive the following corollary about the computational
complexity of Algorithm 1.

######  Corollary 1.

The computational complexity of Algorithm 1 is $2|\mathcal{D}|+1$ times of
solving $\mathcal{M}(y)$.

Although the above computational complexity is not attractive, it accounts for
the worst case. Numerical experiments in Section 4 demonstrate that the
convergence speed of Algorithm 1 is very fast in most cases.

Remark 3. By changing the update rule of $y_{l}$ in Algorithm 1, we can obtain
different version of algorithms. One example is to let
$y_{l+1}=\mu^{\tilde{d}^{*}(y_{l})}$, i.e., the pseudo mean $y_{l+1}$ is set
as the real mean of the optimal policy $\tilde{d}^{*}(y_{l})$ of
$\mathcal{M}(y_{l})$. Such revised algorithm is very similar to the policy
iteration algorithm for solving local optimality equation in Xia (2020), both
converge to a fixed point solution to (21).

Besides Lemma 4, we may further improve the shrinking efficiency of dominated
areas by using other properties. From the viewpoint of bi-objective
optimization in Fig. 1, we observe that the minimization of objective
$\beta\sigma^{d}-\mu^{d}$ is interpreted to find the last solution $d^{*}$
tangent with the line of slope $\beta$ when the line is moving toward top-
left. All the solutions located at the down-right side of this line have a
worse objective $\beta\sigma^{d}-\mu^{d}$ than that of the solution $d^{*}$.
This fact is illustrated by Fig. 5.

Figure 5: Illustration of the dominated policy area indicated by Lemma 5.

Therefore, based on an optimal policy $\tilde{d}^{*}$ by solving the pseudo
mean-variance MDP $\mathcal{M}(y)$, we directly derive the following lemma
about the shrinkage of dominated areas.

######  Lemma 5.

For any policy $\tilde{d}^{*}\in\mathcal{D}$, the policies in the following
areas are dominated by $\tilde{d}^{*}$ and can be discarded:
① any policy $d\in\mathcal{D}$ satisfying $\mu^{d}\in(-\infty,\
\mu^{\tilde{d}^{*}}-\beta\sigma^{\tilde{d}^{*}}]$;
② any policy $d\in\mathcal{D}$ satisfying
$\mu^{d}\in[\mu^{\tilde{d}^{*}}-\beta\sigma^{\tilde{d}^{*}},\ +\infty)$ and
$\sigma^{d}\in[\sigma^{\tilde{d}^{*}}+\frac{1}{\beta}(\mu^{d}-\mu^{\tilde{d}^{*}}),\
+\infty)$.

The area ① in Lemma 5 is similar to the area discarded by Lemma 4, both are
square areas and have no constraints on variances. Therefore, we can utilize
the rule ① in Lemma 5 to speed up the domain shrinking of $\mathbb{Y}$ in
Algorithm 1. That is, at line 7 of Algorithm 1, we can add an extra operation
to discard the area ① indicated by Lemma 5:

$\mathbb{Y}=\mathbb{Y}-(-\infty,\
\mu^{\tilde{d}^{*}(y_{l})}-\beta\sigma^{\tilde{d}^{*}(y_{l})}].$ (26)

We call such algorithm revision _Algorithm 1-Plus_, whose performance is
compared in our numerical experiments in Section 4. For example, in Fig. 4,
when $\beta$ is relatively small, the Iteration 5 will be saved if we apply
(26) for policy $\tilde{d}^{*}(y_{1})$ at the Iteration 2. This demonstrates
that Algorithm 1-Plus is computationally saving compared with Algorithm 1.

Remark 4. It is easy to verify that all the results in this paper can be
extended to solely minimizing the steady-state variance of MDPs. One trivial
method is to let the coefficient $\beta$ in (4) large enough to approximate
the variance minimization of MDPs. Actually, if we replace the mean-variance
objective $\beta\sigma^{d}-\mu^{d}$ in (4) with the variance $\sigma^{d}$, we
can rigorously prove that all the previous results hold for this variance
minimization problem. Algorithm 1 also works to find the optimal policy that
attains the global minimum of the variance in MDPs.

##  4 Numerical Experiments

In this section, we validate the proposed algorithms with a multi-period
inventory control problem, where we consider both the steady-state mean and
variance of rewards. This problem is modeled as an infinite-horizon discrete-
time undiscounted MDP. The inventory capacity is $C\in\mathbb{N}^{+}$. At each
epoch $t=0,1,\dots$, the inventory level
$s_{t}\in\mathcal{S}=\\{0,1,\dots,C\\}$ is reviewed and an order
$a_{t}\in\mathcal{A}(s_{t})=\\{0,\dots,C-s_{t}\\}$ is made. The demands
$\xi_{t}\sim B(C,p)$ are independent and identically distributed, where
$B(C,p)$ is a binomial distribution and $p$ is the probability of success.
There is no lead time and the next inventory level is determined as
$s_{t+1}=[s_{t}+a_{t}-\xi_{t}]^{+}$. The reward function is
$r(s_{t},a_{t})=\mathbb{E}[r(s_{t+1}|s_{t},a_{t})]=-\mathbb{E}\\{ba_{t}+hs_{t+1}+l[\xi_{t}-s_{t}-a_{t}]^{+}\\}$,
where $b,h$ and $l$ are ordering, holding and shortage costs per unit,
respectively. By default, we set $C=4,p=0.6,b=1,h=0.7,l=2.9,$ and $\beta=10$.
We run algorithms 50 replications for statistical analysis.

Fig. 6 illustrates an example of the convergence process of Algorithm 1, where
the interval $[\underline{r},\overline{r}]$ is covered iteratively and the
global optimum is found after only 6 iterations. This demonstrates the
efficiency of Algorithm 1, although the policy space is large as
$|\mathcal{D}|=(C+1)!$.

Figure 6: A convergence process of Algorithm 1, where the blue cross signs
represent the optimal policies of auxiliary problems $\mathcal{M}(y)$ and the
red stars represent the best-so-far policies.

As a comparison, we also implement the local optimization algorithm proposed
by Xia (2020). Considering that the mean-variance optimization of this problem
usually has multiple local optima, we illustrate the performance comparison of
these two algorithms in Fig. 7, where different problem sizes
$C\in\\{4,7,10,20,30,50\\}$ are used. We can see that our global algorithm has
much better performance and the local algorithm by Xia (2020) may converge to
different local optima shown by the whiskers of standard deviations.

Figure 7: Performance comparisons between Algorithm 1 and the local
optimization algorithm by Xia (2020), with different inventory capacities.

Fig. 8 shows the curves of optimal pseudo mean-variance $\tilde{\eta}^{*}(y)$
with respect to the pseudo mean $y$. For capacity $C=4$, the global optimum is
$\eta^{*}=4.500$ and the other two local optima are 5.376 and 6.382, which
coincide with the left pair of bars in Fig. 7. The pseudo mean corresponding
to $\eta^{*}$ is $y^{*}=-3.891$, which also equals the mean of the star point
in the last subfigure of Fig. 6. All these demonstrate that our Algorithm 1
truly finds the global optimum and the local algorithm by Xia (2020) randomly
converges to different local optima. Moreover, when the capacity increases,
the curve of $\tilde{\eta}^{*}(y)$ has more local optima and the local
algorithm is more possibly trapped in a worse local optimum. This also
explains the big performance gaps in Fig. 7 when the capacity is large.

Figure 8: Curves of the optimal pseudo mean-variance $\tilde{\eta}^{*}(y)$
under different capacities.

Furthermore, we study the effect of risk coefficient $\beta$ on the curve
$\tilde{\eta}^{*}(y)$, as illustrated in Fig. 9. We observe that the problem
complexity is increasing with respect to $\beta$. When $\beta$ is small, the
curve has only a single local optimum, which indicates that the problem is
easy to solve. This is because a mean-variance problem with a small $\beta$ is
approximately equivalent to only optimizing the mean performance, which is a
standard MDP easy to solve. Oppositely, when $\beta$ is large, the curve has
multiple local optima and the associated optimization problem is difficult to
solve.

Figure 9: Curves of the optimal pseudo mean-variance $\tilde{\eta}^{*}(y)$
under different $\beta$’s.

Finally, we study the effect of Lemma 5 on the algorithm efficiency. We
compare the performance difference between Algorithm 1 and Algorithm 1-Plus
under different capacities and $\beta$’s. We observe that Algorithm 1-Plus can
achieve a significant efficiency improvement when the problem size (capacity)
is large, as shown in Fig. 10(a). When $\beta$ is changed, there are three
cases as shown in Fig. 10(b):

  1. 1.

When $\beta$ is relatively small ($\leq 0.1$), the variance is trivial, and
the mean-variance optimization is approximately equivalent to a mean
optimization problem which is a standard MDP. The problem is relatively easy,
and these two algorithms have similar efficiency;

  2. 2.

When $\beta$ is relatively large ($\geq 100$), Lemma 5 may rarely remove areas
with means smaller than $\underline{r}$, which can be illustrated by the
intercept in Fig. 5 when the line slope is large. Thus, these two algorithms
also have similar efficiency in this case;

  3. 3.

In other cases, Lemma 5 significantly improves the algorithm convergence
speed, and Algorithm 1-Plus is quite more efficient than Algorithm 1.

(a) Effect of Lemma 5 under different capacities.

(b) Effect of Lemma 5 under different $\beta$’s.

Figure 10: Efficiency comparisons between Algorithm 1 and Algorithm 1-Plus
under different capacities and $\beta$’s.

##  5 Discussion and Conclusion

This paper proposes the global algorithms for solving multi-period mean-
variance optimization in the framework of MDPs, which is a long-standing
challenge caused by the failure of dynamic programming. We convert this
problem to a bilevel MDP formulation, where the inner optimization is a
standard MDP $\mathcal{M}(y)$ for pseudo mean-variance optimization and the
outer one is a single parameter selection problem optimizing pseudo mean $y$.
Interestingly, the optimal value of $\mathcal{M}(y)$ is a convex piecewise
quadratic function of $y$. By the square form difference between the real
variance and the pseudo variance, we discover policy dominance properties to
help remove worse policy spaces iteratively. The global optimum can be found
by repeatedly removing these dominated policy spaces. The convergence and
efficiency of our algorithms are studied both theoretically and
experimentally.

Our work demonstrates a promising approach to globally optimize the steady-
state mean-variance metrics in undiscounted MDPs. It is meaningful to further
extend our approach to mean-variance optimization of discounted MDPs. Another
interesting topic is to develop reinforcement learning algorithms based on our
global optimization approach, which can make our approach implementable in a
data-driven environment.

## References

  * Bisi et al. (2020) Bisi, L., Sabbioni, L., Vittori, E., Papini, M., and Restelli, M. (2020). Risk-averse trust region optimization for reward-volatility reduction. _Proceedings of the 29th International Joint Conference on Artificial Intelligence (IJCAI’2020)_ , Special Track on AI in FinTech, 4583-4589.
  * Borkar (2010) Borkar, V. (2010). Learning algorithms for risk-sensitive control. _Proceedings of the 19th International Symposium on Mathematical Theory of Networks and Systems (MTNS’2010)_ , July 5-9, 2010, Budapest, Hungary, 1327-1332.
  * Cao (2007) Cao, X. R. (2007). _Stochastic Learning and Optimization – A Sensitivity-Based Approach_. New York: Springer.
  * Chung (1994) Chung, K. J. (1994). Mean-variance tradeoffs in an undiscounted MDP: The unichain case. _Operations Research_ 42, 184-188.
  * Cui et al. (2022) Cui, X. Y., Gao, J. J., Li, X., and Shi, Y. (2022). Survey on multi-period mean-variance portfolio selection model. _Journal of the Operations Reserach Society of China_ 10, 599-622.
  * Dai et al. (2021) Dai, M., H. Jin, S. Kou, Y. Xu. (2021). A dynamic mean-variance analysis for log returns. _Management Science_ 67(2), 1093-1108.
  * Filar and Lee (1985) Filar, J. A. and Lee, H. M. (1985). Gain/variability tradeoffs in undiscounted Markov decision processes. _Proceedings of the 24th IEEE Conference on Decision and Control (CDC’1985)_ , 1106-1112.
  * Gal and Greenberg (1997) Gal, T. and Greenberg, H. J. (1997). _Advances in Sensitivity Analysis and Parametric Programming_. Kluwer, Dordrecht.
  * Gao and Li (2013) Gao, J. J. and Li, D. (2013). Optimal cardinality constrained portfolio selection. _Operations Research_ 61(3), 745-761.
  * Guo and Song (2009) Guo, X. and Song, X. Y. (2009). Mean-variance criteria for finite continuous-time Markov decision processes. _IEEE Transactions on Automatic Control_ 54, 2151-2157.
  * Guo et al. (2012) Guo, X., Ye, L., and Yin, G. (2012). A mean-variance optimization problem for discounted Markov decision processes. _European Journal of Operational Research_ 220, 423-429.
  * Haskell and Jain (2013) Haskell, W. B. and Jain, R. (2013). Stochastic dominance-constrained Markov decision processes. _SIAM Journal on Control and Optimization_ 51(1), 273-303.
  * Hernandez-Lerma et al. (1999) Hernandez-Lerma, O., Vega-Amaya, O., and Carrasco, G. (1999). Sample-path optimality and variance-minimization of average cost Markov control processes. _SIAM Journal on Control and Optimization_ 38, 79-93.
  * Huang (2018) Huang, Y. (2018). Finite horizon continuous-time Markov decision processes with mean and variance criteria. _Discrete Event Dynamic Systems: Theory and Applications_ 28(4), 539-564.
  * Li and Ng (2000) Li, D. and Ng, W.-L. (2000). Optimal dynamic portfolio selection: Multiperiod mean-variance formulation. _Mathematical Finance_ 10(3), 387-406.
  * Markowitz (1952) Markowitz, H. (1952). Portfolio selection. _The Journal of Finance_ 7, 77-91.
  * Prashanth and Ghavamzadeh (2013) Prashanth, L. A. and Ghavamzadeh, M. (2013). Actor-critic algorithms for risk-sensitive MDPs. _Advances in Neural Information Processing Systems (NIPS’2013)_ , 252-260.
  * Puterman (1994) Puterman, M. L. (1994). _Markov Decision Processes: Discrete Stochastic Dynamic Programming_. New York: John Wiley & Sons.
  * Rockafellar and Uryasev (2000) Rockafellar, R. T. and Uryasev, S. (2000). Optimization of conditional value-at-risk. _Journal of Risk_ 2(3), 21-42.
  * Sobel (1982) Sobel, M. J. (1982). The variance of discounted Markov decision processes. _Journal of Applied Probability_ 19(4), 794-802.
  * Sobel (1994) Sobel, M. J. (1994). Mean-variance tradeoffs in an undiscounted MDP. _Operations Research_ 42, 175-183.
  * Sutton and Barto (2018) Sutton, R. S. and Barto, A. G. (2018). _Reinforcement Learning: An Introduction, 2nd Edition_. MIT Press, Cambridge, MA.
  * Tamar et al. (2012) Tamar, A., Castro, D. D., and Mannor, S. (2012). Policy gradients with variance related risk criteria. _Proceedings of the 29th International Conference on Machine Learning (ICML’2012)_ , Edinburgh, Scotland, 387-396.
  * Tan and Hartman (2011) Tan, C. H. and Hartman, J. C. (2011). Sensitivity analysis in Markov decision processes with uncertain reward parameters. _Journal of Applied Probabilities_ 48, 954-967.
  * Xia (2016) Xia, L. (2016). Optimization of Markov decision processes under the variance criterion. _Automatica_ 73, 269-278.
  * Xia (2018) Xia, L. (2018). Mean-variance optimization of discrete time discounted Markov decision processes. _Automatica_ 88, 76-82.
  * Xia (2020) Xia, L. (2020). Risk-sensitive Markov decision processes with combined metrics of mean and variance. _Production and Operations Management_ 29(12), 2808-2827.
  * Xia and Glynn (2016) Xia, L. and Glynn, P. W. (2016). A generalized fundamental matrix for computing fundamental quantities of Markov systems. arXiv:1604.04343.
  * Xie et al. (2018) Xie, T., Liu, B., Xu, Y., Ghavamzadeh, M., Chow, Y., Lyu, D., and Yoon, D. (2018). A block coordinate ascent algorithm for mean-variance optimization. _In Advances in Neural Information Processing Systems (NIPS’2018)_ , 1065-1075.
  * Zhang et al. (2021) Zhang, S., Liu, B., and Whiteson, S. (2021). Mean-variance policy iteration for risk-averse reinforcement learning. _Proceedings of the AAAI Conference on Artificial Intelligence (AAAI’2021)_ 35(12), 10905-10913.
  * Zhou and Li (2000) Zhou, X. Y. and Li, D. (2000). Continuous-time mean-variance portfolio selection: A stochastic LQ framework. _Applied Mathematics and Optimization_ 42, 19-33.
  * Zhou and Yin (2004) Zhou, X. Y. and Yin, G. (2004). Markowitz’s mean-variance portfolio selection with regime switching: A continuous-time model. _SIAM Journal on Control and Optimization_ 42, 1466-1482.
  * Zhu et al. (2004) Zhu, S. S., Li, D., and Wang, S. Y. (2004). Risk control over bankruptcy in dynamic portfolio selection: A generalized mean-variance formulation. _IEEE Transactions on Automatic Control_ 49(3), 447-457.