##  0.1 Graph Shift Operator

By focusing on analogies with the classical signal processing time shift or
delay operator as a fundamental operation that shifts a signal along the time
axis without altering its shape, the graph shift operator (GSO) was initially
defined in [sandryhailaDiscreteSignalProcessing2013] as a matrix
$S\in\mathbb{R}^{N\times N}$ with $S_{ij}=0$ for $i\neq j$ and
$i,i\notin\mathcal{E}$. This means that $S$ can take nonzero values at the
edges of $\mathcal{G}$ or its diagonal, although this is not required. $A$
[shumanEmergingFieldSignal2013] and $L$
[sandryhailaDiscreteSignalProcessing2013] were used in the development of GSP,
and [ortegaGraphSignalProcessing2018] later showed them to be equivilent.
[heimowitzUnifiedViewDiffusion2017] showed how $P$ can be used as a GSO,
allowing one to define a graph signal as a function of diffusion embedding
vectors instead of a function over the nodes of a graph. Another approach in
[shumanEmergingFieldSignal2013] observes that the translation of a signal x is
the convolution of x with the Dirac mass $\delta_{\tau}$. In recent years, the
GSO has been reframed as an algebraic descriptor of the graph topology in
[mateosConnectingDotsIdentifying2019] and extended to consider energy
preservation shifts [gaviliShiftOperatorGraph2017] and isometric shifts
[giraultTranslationGraphsIsometric2015].

Expanding on the original definition, the current definition of a GSO $S$,
associated with a graph structure $\mathcal{G}$, is a linear mapping between
graph signals $S:\mathcal{X}(\mathcal{G})\rightarrow\mathcal{X}(\mathcal{G})$
such that the output graph signal $y=Sx$ is calculated as
$y_{v}=\sum_{v\in\mathcal{N}(v)}S_{uv}X_{u}$. The GSO is determined by the
adjacency rule that assigns the values $S_{uv}$ for every
$(u,v)\in\mathcal{E}$. [giraultSignalProcessingGraphs] summarizes three
properties111The paper also defines Localization preserving: $\forall i\exists
k_{i}:S(\delta_{i})=\delta_{k_{i}}$ (delta signal is mapped to another delta
signal) as another property. But, the paper then shows that all four
properties of the time shift can be preserved if and only if the graph Fourier
transform is a DFT, which is unlikely, so there is no usable operator on graph
signals verifying all four properties. Thus, they choose to enforce only the
first three: linearity, convolution, and isometric concerning the chosen
energy definition. that a graph shift operator should share with the time
shift:

  * •

Linear: $S(\alpha x+y)=\alpha S(x)+S(y)$ (superposition principle)

  * •

Convolutive: $\forall v_{k}\neq v_{l},\quad(\hat{S})_{kl}=0$

  * •

Energy preserving (isometric): $E_{x}=E_{S(x)};E_{x}=x\times
x=\sum_{i\in\mathcal{V}}|x_{i}|^{2}$

Additionally, to ensure numeric stability, we assume that the spectral norm of
the shift operator is upper bounded, i.e., $\|S\|_{2}\leq\rho$.

In [heimowitzUnifiedViewDiffusion2017], the authors proposed using the Markov
matrix as a graph shift operator, which would redistribute the graph signal
according to the transition probabilities, so, for node $v_{i}$, the shift
operation would result in $\sum_{m=1}^{N}P(v_{i}v_{m})x_{m}$. Realizing that,
according to the diffusion map process, most of the graph’s geometry is
described by the largest eigenvalues and their associated eigenvectors of the
Markov matrix, we can take the next logical step and consider the diffusion
map embedding as a graph shift operator.

The diffusion map process can also act as a filter for graph signals, with the
benefit of reducing complexity [talmonDiffusionMapsSignal2013]. Take the data
set and compute the diffusion maps $\Psi_{t}(i)$ as before. Since the
eigenvectors are orthogonal, they form the complete basis for any actual
signal, allowing it to be viewed as $x=\sum_{j=0}^{M-1}\lambda_{j}\langle
x(j)\Psi_{j}\rangle\Psi_{j}$. Therefore, a subset of the eigenvectors $l<<N$
exists that conveys the signal information, with the remainder representing
noise. This can be expressed as a graph filter $h(x)_{l}$, a linear projection
onto the eigenvectors spanning the parameter subspace.

The primary benefit of this definition is that it allows us to define a graph
signal as a function of the diffusion embedding vectors instead of a function
of the nodes of a graph since the diffusion map embedding preserves the local
structure of the graph. Thus, we consider the diffusion map embedding of the
original data as a graph shift operator, $S_{DM}$, defined as
$S_{DM}=\Psi_{P}=\lambda_{P}\psi_{P}$. Intuitively, this embedding has the
same information as the transition matrix $P$ but in a lower dimension and can
detect temporal patterns in the graph. The diffusion map framework establishes
an embedding for the graph nodes contingent on the timescale parameter $t$.
Consequently, $S_{DM}$ is closely related to a time shift and can be employed
to identify temporal patterns in the graph, expressed as

$x=V\hat{x}=\Psi_{0}\hat{x},\quad
P^{\top}x=V\Lambda^{\top}\hat{x}=\Psi_{T}\hat{x}.$

This means that the graph shift operation can be interpreted as a shift in the
family of embedding, for example, from $\\{\Psi_{0}\\}$ to $\\{\Psi_{T}\\}$.
[heimowitzUnifiedViewDiffusion2017] One additional assumption is that the
graph signals are stationary with respect to a specific GSO $S$ if either the
signal can be expressed as the output of a graph filter with white inputs or
if its covariance matrix is simultaneously diagonalizable with S.
[giraultStationaryGraphSignals2015]

######  Theorem 0.1.

$S_{DM}$ is a GSO

######  Proof 0.2.

To determine if $S_{DM}$ is a GSO, apply the elements of the definition:

  * •

The symmetry of $S_{DM}$ depends on P. If P was generated using a symmetric
distance function, P is symmetric. Alternatively, use LABEL:DM-process to
construct $S_{DM}$ based on $L_{sym}$.

  * •

Symmetry implies real eigenvalues and orthogonal real eigenvectors, so
$S_{DM}$ is real.

  * •

$\|S_{DM}\|=\sqrt{\lambda_{1}}$ where $\lambda_{1}$ is the largest eigenvalue
of $S_{DM}$. $S_{DM}=\Psi_{P}=\lambda_{P}\psi_{P}$. Since P is a Markov
matrix, it is row stochastic, and $\rho(P)=1$, so $|\lambda|\leq 1$. Thus
$\lambda_{1}\leq 1$ so $\|S_{S}M\|$ is bounded.

  * •

For $x,y\in\mathcal{X}$, $S_{DM}(\alpha x+y)=\Psi_{i}(\alpha
x_{i}+y_{i})=\alpha\Psi_{i}x_{i}+\Psi_{i}y_{i}$ since
$\Psi_{i}=\lambda_{i}\psi_{i}=$ is constant.

  * •

The convolutive property refers to the ability of $\hat{S_{DM}}$ to have an
orthogonal diagonalization. Since $S_{DM}$ is real, symmetric, and isometric,
then $F=\Psi(S_{DM})$ is orthogonal, so $\hat{S_{DM}}=FS_{DM}F^{-1}$ and then
$\forall l:v_{k}\neq v_{i},\quad\hat{S_{DM}}_{ki}=0$

  * •

Isometric matrices have orthogonal columns, resulting from $S_{DM}$ being
real. Furthermore, since P is stochastic, we have $D(P)=1$. As
$S_{DM}=\Psi\Lambda\Psi^{-1}$, applying the graph shift operator $n$ times to
the graph signal $x$ results in
$x_{n}=S_{DM}^{n}x=\Psi\Lambda^{n}\Psi^{-1}x\approx\pi x$ because
$\pi=P\pi=\frac{\psi_{0}}{\sum_{i}\psi_{i}}$. Therefore, $S_{DM}$ preserves
the energy of the signal $x$.

The graph shift operator $S_{DM}$ can be approximated using diffusion maps as
$S_{DM}^{t}\approx I-P^{t}$, as the diffusion operator $P^{t}$ is related to
the graph Laplacian by formula $P^{t}=\exp(-tL)$, where $t$ is diffusion time.
Thus, $S_{DM}$ represents one discrete step in how information moves from one
node to its neighbors. In figure 1, different simple filters $H(S)$ are
applied to a simple signal $\delta(i)=1$ associated with an equally weighted
2D-lattice graph. Repeated iterations show how the signal spreads based on the
adjoining nodes. Since the kernel function ensures that closely connected
nodes in the graph have close diffusion embedding vectors, the graph signal on
the vertices is mapped onto the set of diffusion embedding vectors. This
suggests that band-limited graph signals can approximate the graph signal when
the graph shift operator is the Markov matrix. By enforcing the preservation
of the power spectrum (i.e., through isometry and convolution), we allowed the
graph translation to preserve the bandwidth of the signal through multiple
applications of the operator. [gaviliShiftOperatorGraph2017]

Figure 1: applying different $H(S)^{t}$ to $x$

Since the choice of the GSO determines the signal model, using $S_{DM}$, we
can now write the general graph filter as

$H(x)=h_{0}\Psi_{0}\hat{x}+h_{1}\Psi_{1}\hat{x}+\dots+h_{L}\Psi_{L}\hat{x}.$

Usually, $l<<N$, so $H(x)$ will be a low-degree matrix polynomial, and the
output will strictly depend on the connectivity structure and signal values in
the node’s local neighborhood. As $S_{DM}=\Psi\Lambda\Psi^{\top}$ The GFT is
now given by $\hat{\mathbf{x}}=\mathbf{\Psi}^{\top}\mathbf{x}$. The iGFT is
$\mathbf{x}=\mathbf{\Psi}\hat{\mathbf{x}}$. In general, filtering a graph
signal, $x\in\mathbb{R}^{m}$ by $H$ is defined as

$y=H(S_{DM})x=\sum_{i}u_{i}h(\lambda_{i})\hat{x}$

where $\\{u_{i},\lambda_{i}\\}$ are the pairs of eigenvalues-eigenvectors of
$S_{DM}$, and $\hat{x}\in\mathbb{R}^{m}$ is the graph Fourier representation
of $x$. The convolution operator is now
$(f*g)(i)=\sum_{l=0}^{N-1}\hat{f}(\lambda_{l})\hat{g}(\lambda_{l})u_{l}(i)$,
so convolution in the vertex domain is equal to multiplication in the graph
spectral domain.

Figure 2: Effect of a low-pass filter222Inspired from
https://pygsp.readthedocs.io/en/latest/examples/filtering.html

##  0.2 application of $S_{DM}$ in graph learning

The objective of graph learning is to derive information or relationships
about the underlying graph structure given a set of data. A reasonable
assumption would be that the data should be smooth so that small changes would
be local, based on some understanding of distance, while significant changes
would indicate a more considerable distance. Classically, the total variation
of a signal is defined as the sum of squared differences in consecutive signal
samples $\Sigma_{n}=(x_{n}-x_{n-1})^{2}$, and a signal is called smooth if the
quadratic form takes a small value relative to its norm. This concept can be
extended to graphs where the notion of neighborhood replaces that of
consecutive nodes to obtain

$tr(X^{\top}LX)=TV_{G}(x)=x\top
Lx=\sum_{n=1}^{N}\sum_{M\in\mathcal{N}(n)}w_{mn}(x_{n}-x_{m})^{2}.$

This definition of total variation allows us to interpret the ordering of the
Laplacian’s eigenvalues in terms of frequencies. The larger eigenvalues
correspond to higher frequencies (larger total variation). Since, for any
function $f\in\mathbb{R}^{N}$,
$f^{\top}Lf=\frac{1}{2}\sum_{i,j=1}^{N}w_{ij}(f_{i}-f_{j})^{2}$, if we let $f$
equal the eigenvector of the graph Laplacian $u_{k}$, the smoothness
$TV(u_{k})=u_{k}^{\top}Lu_{k}$ gives us the corresponding eigenvalue
$\lambda_{k}$.

Although there are various methods for learning a graph, we initially adopt
the approach proposed in [egilmezGraphLearningFiltered2018]. This method
assumes smooth signals and incorporates a pre-filtering step to enhance
convergence. The process involves eigendecomposing the sample covariance
matrix and applying a graph filter $h(x)$. The output is used to estimate a
$L$ that minimizes the total variance. This procedure is repeated iteratively
until the maximum likelihood error falls within a specified tolerance.

As an enhancement to the total variation approach,
[heimowitzMarkovVariationApproach2020] introduced a Markov measure of
variation defined as

$MV(x)=\|x-D^{-1}Wx\|=\|x-Px\|,$

where a smooth vector implies $\|x-Px\|_{1}<f(N)$, with $f$ denoting a
function related to the number of nodes. The $l_{1}$ norm is given by

$\|x-Px\|_{1}=\sum_{i=1}^{N}\frac{1}{d(v_{i})}\left|\sum_{m=1}^{N}W_{i,m}(x_{i}-x_{m})\right|.$

Since $P$ can be decomposed into $V\Lambda V^{-1}$ and
$\Psi_{t}=V\Lambda^{t}$, the Markov measure of variation is expressed as
$MV(x)=|\sum_{i=1}^{N}\ x_{i}-\Psi_{1}^{T}(i)\hat{x}|<f(N)$.

It is reasonable to assume that the measured signals include distortion or
other noise. If we assume that this noise is white, a low-pass filter would
reduce the noise while preserving the signal. After the noise is removed, we
can analyze the smoothness of the signal to derive information about the
underlying graph itself. The process first applies the graph Fourier transform
so that the signal is now in the spectral domain $x_{g}=U^{\top}x$. Next,
apply a convolution to the filter, remaining in the spectral domain
$(f*g)(i)=\sum_{l=0}^{N-1}\hat{f}(\lambda_{l})\hat{g}(\lambda_{l})u_{l}(i)=u_{i}h(\lambda_{i})\hat{x}$.
Last, use an inverse graph Fourier transform and bring the result back in the
vertex domain $x=u^{-1}\hat{x}$. For this paper, we choose a filter based on
Tikhanov Regularization, $h(x)=\frac{1}{1-\tau x}$, where $\tau\in[0,1]$.
Since a diffusion map does not uniquely define a graph, we will seek $L$ in
the graph learning process, which will require a different representation of
$S_{DM}$ as $L$ by $TV=x^{\top}Lx=x^{\top}D^{-1}(P-I)D^{-1}x\approx
x^{\top}D^{-1}(\Psi^{l}-I)D^{-1}x$.

Figure 3: Comparing graph filters and time filters

Since the choice of the GSO leads to a specific signal model, the filters and
transforms are all based on $S$. Using this framework, we can compare the
results of different GSOs. Now, since the filter $H(x)$ smooths the signals
and since $\Psi_{t}(i)$ separates the signal information from noise, using
$S_{DM}$ should minimize the oscillation, or total variation, of the signal
with respect to the learned graph. Therefore, we have the following:

######  Theorem 0.3.

For smooth signals on a connected graph satisfying the following assumptions,
the Markov Variation (MV) optimization process converges to a signal with
minimal variation and a distribution that approaches the Markov stable
distribution of the signals. This convergence is significantly faster than the
convergence of Total Variation (TV).

Assumptions:

  1. 1.

Connected Graph: The underlying graph $\mathcal{G}$ is connected. This ensures
that information can flow freely throughout the graph, which is crucial for
accurate diffusion map construction and convergence.

  2. 2.

Irreducible Markov Chain: The Markov chain associated with the diffusion map
process is irreducible. This means that any state can be reached from any
other state with a non-zero probability in a finite number of steps. This
property guarantees that the diffusion process can explore the entire graph
and capture the global structure of the data.

  3. 3.

Finite State Space: The number of states in the Markov chain is finite,
corresponding to the number of nodes in the graph. This allows for efficient
computation of the diffusion map and ensures that the process reaches a stable
state.

  4. 4.

Bounded Transition Probabilities: The transition probabilities in the
transition matrix $P$ are bounded between 0 and 1. This ensures that the
diffusion process reaches a stationary distribution and prevents the signal
from exploding or oscillating indefinitely.

  5. 5.

Signal Smoothness: The graph signal $x$ is bounded and sufficiently smooth,
meaning its second derivative is bounded by a constant $C$. This assumption
ensures that the signal variations are gradual and allows for the application
of the Rayleigh-Ritz theorem for bounding the difference between TV and MV.

  6. 6.

Hermitian Graph Laplacian: The graph Laplacian $L$ is Hermitian, meaning
$L=L^{\top}$. This property guarantees the symmetry of the diffusion process
and ensures the validity of certain analytical tools used in the proof.

######  Proof 0.4.

Spectral Properties: Let $\lambda_{1}$ and $\lambda_{2}$ be the first and
second smallest eigenvalues of $P$, respectively. Due to connectedness and
aperiodicity, $\lambda_{1}>0$ and there exists a positive spectral gap
$\lambda_{1}-\lambda_{2}>0$. This positive spectral gap is crucial, as it
ensures the ergodicity of the diffusion map process and guarantees its
convergence to the unique and stable Markov distribution.

Bound on the Difference between TV and MV: The Rayleigh-Ritz theorem states
that for a symmetric matrix $A$ and a real vector $x$, the following
inequality holds:

$\lambda_{min}\cdot\|x\|^{2}\leq x^{\top}Ax\leq\lambda_{max}\cdot\|x\|^{2}$

where $\lambda_{min}$ and $\lambda_{max}$ are the smallest and largest
eigenvalues of $A$, respectively. In the context of this proof, the matrix $A$
is the graph Laplacian $L$, and the vector $x$ represents the graph signal.
The leading eigenvector, $\psi_{1}$, corresponding to the smallest eigenvalue
$\lambda_{1}$, captures the most significant information about the graph
structure and encodes the diffusion process. We project the graph signal $x$
onto the leading eigenvector $\psi_{1}$, obtaining the projection
$\hat{x}=\psi_{1}^{\top}x$. The smoothness assumption of the signal, which
limits the magnitude of its second derivative, allows us to relate the
difference between the original signal and its projection to the distance from
the leading eigenvector. The spectral gap plays a crucial role in ensuring the
convergence and precision of the diffusion
map[vonluxburgTutorialSpectralClustering2007]. Using this gap and the
smoothness bound, we obtain:

$|TV(x)-MV(x)|\leq\sqrt{C}\cdot\sqrt{\lambda_{2}}\cdot\|x-\Psi(x)\|$

where $C$ is the smoothness of the signal and $\Psi=\lambda\psi$ is the
diffusion map generated from the signals $x$.

Convergence Rates: Both TV and MV aim to minimize $\|x-\Psi(x)\|$.

  * •

TV: Uses gradient descent, which has an average convergence rate of
$O(1/k^{2})$ for smooth signals, although the convergence can be slow in the
initial stages.

  * •

MV: Directly minimizes the distance to the diffusion map representation
through projection, resulting in a significantly faster $O(1/k)$ convergence
rate due to the spectral gap. This faster rate is achieved because the
projection directly aligns the signal with the leading eigenvector, which
captures the most significant information in the diffusion process.

In general, $O(1/k^{2})$ is faster than $O(1/k)$. However, in the context of
the proof, the statement is comparing the average convergence rates of two
different algorithms. Additionally, the dependence of the bound on the
spectral gap explains why MV converges faster than TV. The spectral gap
indicates the rate at which irrelevant information decays in the diffusion
process, and a larger gap implies faster convergence of MV towards the signal
with minimal variation.

Ergodicity and Stable Distribution: The positive spectral gap guarantees the
ergodicity of the diffusion map process. Ergodicity implies that the
distribution of the transformed signal converges to a unique and stable Markov
distribution, regardless of the initial state of the process. This stable
distribution captures the long-term behavior of the diffusion process and
reflects the underlying structure of the graph.

In conclusion, MV minimizes the distance to the diffusion map representation,
which aligns the signal with the underlying diffusion process. This process
inherently minimizes information loss and ensures that signal variations are
smooth and consistent with the global structure of the graph. Therefore, the
signal obtained through MV optimization can be considered to have minimal
variation within the context of the graph and its diffusion process and has a
distribution that approaches the Markov stable distribution of the signals.
This convergence is significantly faster than the convergence of Total
Variation (TV).

By minimizing the distance to the diffusion map representation, MV not only
reduces the overall variation in the signal, but also promotes its similarity
to a signal generated by the Markov process associated with the graph. This
means that the signal obtained through MV optimization effectively captures
the essential features and dynamics of the diffusion process, making it
particularly useful for tasks such as anomaly detection, community detection,
and graph signal denoising. We can now incorporate this approach into the
graph learning algorithm and measure the accuracy against a known graph and
actual data.


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Metastable supersolid in spin-orbit coupled Bose-Einstein condensates

Wei-Lei Xia  Guangdong Provincial Key Laboratory of Quantum Engineering and
Quantum Materials, School of Physics and Telecommunication Engineering, South
China Normal University, Guangzhou 510006, China  Guangdong-Hong Kong Joint
Laboratory of Quantum Matter, Frontier Research Institute for Physics, South
China Normal University, Guangzhou 510006, China  Lei Chen  Guangdong
Provincial Key Laboratory of Quantum Engineering and Quantum Materials, School
of Physics and Telecommunication Engineering, South China Normal University,
Guangzhou 510006, China  School of Physics and Electronic Science, Zunyi
Normal University, Zunyi 563006, China  Tian-Tian Li  Guangdong Provincial Key
Laboratory of Quantum Engineering and Quantum Materials, School of Physics and
Telecommunication Engineering, South China Normal University, Guangzhou
510006, China  Guangdong-Hong Kong Joint Laboratory of Quantum Matter,
Frontier Research Institute for Physics, South China Normal University,
Guangzhou 510006, China  Yongping Zhang  Department of Physics, Shanghai
University, Shanghai 200444, China  Qizhong Zhu<EMAIL_ADDRESS>Guangdong
Provincial Key Laboratory of Quantum Engineering and Quantum Materials, School
of Physics and Telecommunication Engineering, South China Normal University,
Guangzhou 510006, China  Guangdong-Hong Kong Joint Laboratory of Quantum
Matter, Frontier Research Institute for Physics, South China Normal
University, Guangzhou 510006, China

(September 3, 2024)

###### Abstract

Supersolid is a special state of matter with both superfluid properties and
spontaneous modulation of particle density. In this paper, we focus on the
supersolid stripe phase realized in a spin-orbit coupled Bose-Einstein
condensate and explore the properties of a class of metastable supersolids. In
particular, we study a one-dimensional supersolid whose characteristic wave
number $k$ (magnitude of wave vector) deviates from $k_{m}$, i.e., the one at
ground state. In other words, the period of density modulation is shorter or
longer than the one at ground state. We find that this class of supersolids
can still be stable if their wave numbers fall in the range $k_{c1}<k<k_{c2}$,
with two thresholds $k_{c1}$ and $k_{c2}$. Stripes with $k$ outside this range
suffer from dynamical instability with complex Bogoliubov excitation spectrum
at long wavelength. Experimentally, these stripes with $k$ away from $k_{m}$
are accessible by exciting the longitudinal spin dipole mode, resulting in
temporal oscillation of stripe period as well as $k$. Within the mean-field
Gross-Pitaevskii theory, we numerically confirm that for a large enough
amplitude of spin dipole oscillation, the stripe states become unstable
through breaking periodicity, in qualitative agreement with the existence of
thresholds of $k$ for stable stripes. Our work extends the concept of
supersolid and uncovers a new class of metastable supersolids to explore.

##  I introduction

Supersolid is an exotic state of matter, which features both superfluid
properties and spontaneous modulation of particle density. Supersolids are
characterized by two types of gapless Goldstone modes, associated with
spontaneous breaking of U(1) gauge symmetry and spatial translation symmetry,
respectively. Previously, it has been long speculated that solid helium-4 at
low temperature could be a supersolid [1, 2], which, however, is overturned by
subsequent experiments [3, 4]. In recent years, this concept has stimulated
renewed interest among researchers thanks to significant progress made in the
community of cold atomic gases. Till now, there are three successful schemes
to realize a supersolid, i.e., the stripe phase of spin-orbit coupled Bose-
Einstein condensate (BEC) [5, 6, 7, 8, 9], BEC coupled with modes of two
optical cavities [10], and recently in dipolar gas with roton excitation
spectrum [11, 12, 13, 14, 15, 16, 17].

Figure 1: (a) Non-interacting dispersion of spin-orbit coupled BEC. The green
(blue) dots denote the momentum of plane waves with momentum $k_{1}<k_{m}$
($k_{2}>k_{m}$), whose interference results in a stripe with period longer
(shorter) than the one at ground state, as illustrated by their density
modulations in (b). (c),(d) Nonlinear dispersions of chemical potential
$\mu(k)$ determined by Eq. 6. (e),(f) Dispersions of energy $E(k)$ defined in
Eq. 5 corresponding to (c) and (d), respectively. Line color in (c)-(f)
denotes the density contrast $\mathcal{C}$, as shown by the colorbar.
$\bar{n}g/k_{0}^{2}=0.484$, $\bar{n}g_{12}/k_{0}^{2}=0.44$,
$\Omega/k_{0}^{2}=0.09$ for (c),(e), and $\Omega/k_{0}^{2}=0.63$ for (a), (d)
and (f), respectively. $\bar{n}$ is the average condensate density.

A supersolid at ground state and its low energy excitations have been
extensively studied, where special attention is paid to the existence of two
types of gapless Goldstone modes [18, 19, 14, 20, 21]. These studies consider
the collective excitations in the linear response regime, which reflects the
properties of a supersolid close to ground state. In fact there are other
class of excited states of a supersolid. For example, in the plane wave phase
of spin-orbit coupled BEC, the wave vector does not necessarily coincide with
the one at ground state, but instead can take values away from the lowest
energy one. In particular, these states correspond to an excited state
carrying mass current, and under certain conditions they become metastable
superfluid [22, 23, 24]. Analogously, for the stripe phase, one can also vary
the wave vector, or equivalently, the period of stripe, to realize a
supersolid at excited state (see Figs. 1(a) and 1(b)). This scenario is also
experimentally relevant as the supersolid prepared may not definitely reach
the ground state during limited observation time. Will these supersolids be
metastable? If not, what is the underlying instability mechanism? Does there
exist a critical spin supercurrent for these metastable supersolids, in
analogy to the critical mass supercurrent in plane wave phase? These are the
interesting questions we aim to answer in this paper.

To be specific, we focus on the supersolid realized in the stripe phase of
spin-orbit coupled BEC, and consider a supersolid whose wave number $k$
(magnitude of wave vector) does not coincide with the one at ground state,
denoted as $k_{m}$, but instead a supersolid with general $k$. We find that
nonlinear dispersion of chemical potential with respect to $k$ can exhibit a
loop structure. For $k\neq k_{m}$, the supersolid is at excited state and may
not be stable. We study the stability of a general stripe state by calculating
its Bogoliubov excitation spectrum, and find that stripe states with $k$
within the range $k_{c1}<k<k_{c2}$ are metastable. When $k$ falls outside this
range, the supersolid suffers from dynamical instability. These stripes with
$k\neq k_{m}$ can be accessed by exciting the longitudinal spin dipole mode of
stripe [25], and the thresholds of $k$ are related with the upper limit of the
oscillation amplitude of spin dipole mode.

The rest of this paper is organized as follows. In Sec. II, we formulate the
spin-orbit coupled BEC within mean-field theory, numerically solve the wave
function of stripe state, and calculate the nonlinear dispersion $\mu(k)$ as a
function of wave number $k$. In Sec. III, we calculate the excitation spectrum
of stripe phase with Bogoliubov theory, from which we determine the stability
condition of stripes, and examine the instability mechanism. We also show that
there exists a critical spin supercurrent carried by those stable stripes. In
Sec. IV, it is demonstrated that stripes with $k\neq k_{m}$ can be
experimentally accessed by exciting the spin dipole mode, and the stability
condition of stripe obtained in Sec. III is related with the stability of
stripe in response to the spin dipole mode excitation within the framework of
Gross-Pitaevskii theory. Finally, we summarize our results and conclude in
Sec. V.

##  II Nonlinear dispersion

Spin-orbit coupled BEC with equal Rashba and Dresselhauss type can be
engineered by Raman coupling of two internal atomic states, as first
experimentally realized in NIST group [26]. Assume that the Raman lasers
introduce momentum transfer in $x$ direction, and the trap frequencies in $y$
and $z$ directions are high enough, so that the low energy degree of freedom
of this system can be described by the following one-dimensional effective
Hamiltonian with spin-orbit coupling,

$\hat{H}_{0}=\frac{\hbar^{2}}{2m}\left(-i\partial_{x}-k_{0}\sigma_{z}\right)^{2}+\frac{\Omega}{2}\sigma_{x},$
(1)

where $k_{0}$ is the strength of spin-orbit coupling, $\Omega$ is the Raman
coupling strength and Pauli matrices are defined in the pseudospin-$1/2$
Hilbert space. For simplicity, in the following calculations, we set
$\hbar=m=1$. $k_{0}$ and $k_{0}^{2}$ are chosen as momentum and energy units,
respectively. The dispersion of this Hamiltonian has two branches, as shown in
Fig. 1(a).

Now consider Bose-Einstein condensation occurs in this system, and the two-
component condensate can be described by the mean-field energy functional
[27],

$\displaystyle{E}\left(\Psi_{\uparrow},\Psi_{\downarrow}\right)=$
$\displaystyle\int
dx\left(\Psi_{\uparrow}^{*},\Psi_{\downarrow}^{*}\right)\hat{H}_{0}\left(\begin{array}[]{c}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{array}\right)+g_{\uparrow\downarrow}|\Psi_{\uparrow}|^{2}|\Psi_{\downarrow}|^{2}$
(4)
$\displaystyle+\frac{g}{2}\left(|\Psi_{\uparrow}|^{4}+|\Psi_{\downarrow}|^{4}\right).$
(5)

Here $\Psi_{\uparrow}$ ($\Psi_{\downarrow}$) is the condensate wave function
of pseudospin up (down), and $g$ ($g_{12}$) is the interaction strength
between same (different) spin species. The bare values of $g$ and $g_{12}$ are
related with the $s$-wave scattering lengths between atoms of same and
different hyperfine states. For 87Rb atom, $g$ and $g_{12}$ are very close
[26], making the density contrast of stripe state very low [27]. However,
their difference can be enhanced by various schemes, e.g., separating the two
components in a bilayer configuration [28].

It is predicted that the ground state of this model has three phases, with the
change of parameters such as $\Omega/k_{0}^{2}$, $g$ and $g_{12}$ [27]. For
atomic density below a critical value, with the increase of
$\Omega/k_{0}^{2}$, the ground state changes from the stripe phase to the
plane wave phase, and then to the zero momentum phase. The stripe phase
spontaneously breaks translation symmetry of the Hamiltonian and develops
periodic modulation of atomic density, which is the focus of this paper.

Minimizing the energy functional with respect to the spinor wave function, one
arrives at the stationary Gross-Pitaevskii (GP) equation,

$\mu\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix}=\begin{pmatrix}H_{\uparrow\uparrow}&\frac{\Omega}{2}\\\
\frac{\Omega}{2}&H_{\downarrow\downarrow}\end{pmatrix}\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix},$ (6)

where
$H_{\uparrow\uparrow}={\hbar^{2}(-i\partial_{x}-k_{0})^{2}}/{2m}+g|\Psi_{\uparrow}|^{2}+g_{12}|\Psi_{\downarrow}|^{2}$,
and
$H_{\downarrow\downarrow}={\hbar^{2}(-i\partial_{x}+k_{0})^{2}}/{2m}+g|\Psi_{\downarrow}|^{2}+g_{12}|\Psi_{\uparrow}|^{2}$.
The wave function of a general stripe characterized by $k$ has the following
form,

$\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix}=\sum_{n=-2L-1}^{2L+1}\begin{pmatrix}a_{n}\\\
b_{n}\end{pmatrix}e^{inkx},$ (7)

where $L$ is a positive integer and the spacing of integer $n$ is $2$. Note
that for the stripe phase, apart from an overall phase factor, the
coefficients satisfy the relation $b_{n}=a_{-n}^{*}$, i.e., symmetric under
the simultaneous operation of space inversion and spin flip. Since the density
now is a spatially periodic function, the exact stripe wave function should
have infinitely many plane waves similar to Bloch waves in periodic
potentials, although in realistic calculations a finite cutoff $L$ has to be
taken. Plugging the above expansion into the stationary GP equation, one
arrives at a set of nonlinear equations that the coefficients satisfy. By
solving the set of nonlinear equations, one can obtain the nonlinear
dispersion, i.e., $\mu$ as a function of $k$, and also the coefficients of
stripe wave function.

Figure 1 shows the nonlinear dispersion of both chemical potential $\mu(k)$
and energy ${E}(k)$. It is clear that, at an optimal wave number $k_{m}$, the
energy of stripe state has the lowest energy. $k_{m}$ is close to the one
determined in non-interacting model, i.e.,
$k_{m}=k_{0}\sqrt{1-\Omega^{2}/4k_{0}^{4}}$ [27]. $k_{m}$ also corresponds to
an optimal stripe period $d_{m}$, through the relation $d_{m}=\pi/k_{m}$. Away
from $k_{m}$, the energy of stripe state increases, and a loop or swallowtail
structure appears in $\mu(k)$ near $k=0$ for small $\Omega$, similar to the
plane wave case [29]. The density contrast of the stripe state, defined as
$\mathcal{C}\equiv(n_{\mathrm{max}}-n_{\mathrm{min}})/(n_{\mathrm{max}}+n_{\mathrm{min}})$
with $n_{\mathrm{max}}$ ($n_{\mathrm{min}}$) being local density maximum
(minimum), at different $k$ is also indicated by color. One can observe that
stripes with small $k$ or long period has larger density contrast. We also
find that the detailed structure of dispersion $\mu(k)$ depends on the cutoff
$L$, but the lower branch of dispersion near $k_{m}$ does not, which is the
range of $k$ we are interested in.

Previous studies have only considered the stripe state with this optimal wave
number $k_{m}$. Here, we have released this constrain and considered a general
stripe with arbitrary $k$, corresponding to ground state or excited states
depending on whether $k=k_{m}$. For $k>k_{m}$, the stripe is effectively
compressed with shorter period; otherwise, it is stretched. Since for $k\neq
k_{m}$, these stripes are at excited state, it is natural to ask whether these
stripes are stable or not. In other words, does there exist metastable
supersolid? This is the issue we will address in the following section.

##  III Stability of supersolid

The stability of stripes at excited state, or $k\neq k_{m}$, can be examined
by studying their Bogoliubov excitation spectra. If the excitation spectrum at
all excitation momentum is real and non-negative, these stripes are stable;
otherwise, they suffer from Landau instability or dynamical instability [30].

Figure 2: (a) The dispersion $\mu(k)$ with five characteristic momenta $k$,
denoted by red dots. From left to right, they correspond to $k$ in (b)-(f).
The dashed lines indicate the lower and upper thresholds of $k$ of stable
stripes. (b)-(f) Bogoliubov excitation spectra for stripes with five $k$
denoted in (a), corresponding to $k<k_{c1}$, $k_{c1}<k<k_{m}$, $k=k_{m}$,
$k_{m}<k<k_{c2}$, and $k>k_{c2}$, respectively. Excitation spectra in (c)-(e)
are all real and non-negative. Dashed lines in (b) and (f) denote the
imaginary parts of complex excitation spectra. $\bar{n}g/k_{0}^{2}=0.484$,
$\bar{n}g_{12}/k_{0}^{2}=0.44$, $\Omega/k_{0}^{2}=0.59$.

We start from the time-dependent GP equation,

$i\hbar\frac{\partial}{\partial t}\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix}=\begin{pmatrix}H_{\uparrow\uparrow}&\frac{\Omega}{2}\\\
\frac{\Omega}{2}&H_{\downarrow\downarrow}\end{pmatrix}\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix},$ (8)

and study the stability of these stripe states with respect to weak
perturbations. Due to periodic modulation of density, the eigenmodes of
perturbation take the form of Bloch waves,

$\displaystyle\begin{pmatrix}\delta\Psi_{\uparrow}\\\
\delta\Psi_{\downarrow}\end{pmatrix}=$ $\displaystyle
e^{-i\mu(k)t/\hbar}\sum_{n=-2L-1}^{2L+1}\begin{pmatrix}u_{n}^{\uparrow}\\\
u_{n}^{\downarrow}\end{pmatrix}e^{i(nk+q)x}e^{-i\varepsilon(q)t/\hbar}$
$\displaystyle+$ $\displaystyle
e^{-i\mu(k)t/\hbar}\sum_{n=-2L-1}^{2L+1}\begin{pmatrix}v_{n}^{\uparrow*}\\\
v_{n}^{\downarrow*}\end{pmatrix}e^{i(nk-q)x}e^{i\varepsilon(q)t/\hbar},$ (9)

where $u_{n}^{\uparrow(\downarrow)}$ and $v_{n}^{\uparrow(\downarrow)}$ are
perturbation amplitudes for spin up (down) component. $q$ is the excitation
momentum. Writing $\Psi_{\sigma}^{\prime}=\Psi_{\sigma}+\delta\Psi_{\sigma}$,
plugging $\Psi_{\sigma}^{\prime}$ into the time-dependent GP equation, and
only retaining linear term of perturbation, one arrives at the Bogoliubov
equation that $u_{n}^{\uparrow(\downarrow)}$ and
$v_{n}^{\uparrow(\downarrow)}$ satisfy, and obtains the excitation spectrum
$\varepsilon(q)$ as well.

The excitation spectra of stripes at different $k$ are depicted in Fig. 2. Due
to periodic density modulation, excitation spectrum $\varepsilon(q)$ has
similar feature with Bloch band of periodic potentials. For a stripe at ground
state, i.e., $k=k_{m}$, the excitation spectrum features two branches of real
and non-negative gapless excitatons (see Fig. 2(d)), corresponding to two
Goldstone modes associated with spontaneous breaking of U(1) gauge symmetry
and translation symmetry [18]. With the deviation of $k$ from $k_{m}$, the
lower branch of excitation spectrum becomes softer (Figs. 2(c) and 2(e)). When
$k<k_{c1}$ or $k>k_{c2}$, with $k_{c1}$ and $k_{c2}$ being the lower and upper
threshold of stable stripes, the excitation spectrum becomes complex at small
$q$, signalling the dynamical instability of those stripes (Figs. 2(b) and
2(f)). We find that the excitation spectrum always turns complex starting from
$q\rightarrow 0$, indicating that instability first appears in long wavelength
excitation.

Both $k_{c1}$ and $k_{c2}$ depend on the Raman coupling strength and the
interaction strengths $g$, $g_{12}$. Figure 3 shows the change of two
thresholds $k_{c1}$ and $k_{c2}$ as functions of $\Omega/k_{0}^{2}$, for two
sets of interaction strengths. One can find that, in both Figs. 3(a) and 3(b),
$k_{c1}$ and $k_{c2}$ are getting closer to each other with the increase of
$\Omega/k_{0}^{2}$, when approaching the stripe and plane wave phase boundary
[27]. Besides, $k_{c1}$ and $k_{c2}$ almost lie symmetrically on opposite
sides of $k_{m}$. The range of stable $k$ becomes narrower, consistent with
the expectation that the stripe states are more susceptible to dynamical
instability when approaching the phase boundary. Comparison between Figs. 3(a)
and 3(b) also shows that, the larger ratio $g/g_{12}$ leads to a wider range
of $k$ of stable stripes.

The stability condition also sets a limit on the magnitude of spin
supercurrent [31]. The stripe state carries a spin current with the $z$
component of spin current density tensor given by [32, 33]

$\boldsymbol{J}_{z}^{s}=\frac{1}{2}\left\\{\sum_{\sigma=\uparrow,\downarrow}\Psi_{\sigma}^{*}(s_{z})_{\sigma\sigma}\hat{v}_{\sigma\sigma}\Psi_{\sigma}+\mathrm{c.c.}\right\\},$
(10)

where the velocity operator is defined as
$\hat{v}=-i\hbar\partial_{x}/m-k_{0}\sigma_{z}$. For easier illustration, we
first take the stripe wave function with two modes as an example, in which
case the spin current reads

$\displaystyle\boldsymbol{J}_{z}^{s}=$
$\displaystyle\frac{k}{2}\left(|a_{1}|^{2}-|b_{1}|^{2}-|a_{-1}|^{2}+|b_{-1}|^{2}\right)-\frac{k_{0}}{2}$
$\displaystyle-2k_{0}\mathrm{Re}\left[\left(a_{1}^{*}a_{-1}+b_{1}^{*}b_{-1}\right)\exp(-2ikx)\right].$
(11)

After integrating the spin current over a period, the average spin current
becomes

$\boldsymbol{\bar{J}}_{z}^{s}=\frac{\pi}{2}\left(|a_{1}|^{2}-|b_{1}|^{2}-|a_{-1}|^{2}+|b_{-1}|^{2}\right)-\frac{\pi
k_{0}}{2k}.$ (12)

For a general stripe wave function with multiple momenta, the average spin
current is found to be

$\boldsymbol{\bar{J}}_{z}^{s}=\frac{\pi}{2}\sum_{m=1}^{2L+1}m\left(|a_{m}|^{2}-|b_{m}|^{2}-|a_{-m}|^{2}+|b_{-m}|^{2}\right)-\frac{\pi
k_{0}}{2k},$ (13)

where $m=1,3,\dots,2L+1$. So the thresholds of $k$ for stable stripes also set
a limit on the magnitude of spin supercurrent. The magnitude of spin
supercurrent lies within the range
$\boldsymbol{\bar{J}}_{c1}\leq\boldsymbol{\bar{J}}_{z}^{s}\leq\boldsymbol{\bar{J}}_{c2}$,
with
$\boldsymbol{\bar{J}}_{c1}={\pi}\sum_{m=1}^{2L+1}m\left(|a_{m}|^{2}-|b_{m}|^{2}-|a_{-m}|^{2}+|b_{-m}|^{2}\right)/2-{\pi
k_{0}}/{2k_{c1}}$, and
$\boldsymbol{\bar{J}}_{c2}={\pi}\sum_{m=1}^{2L+1}m\left(|a_{m}|^{2}-|b_{m}|^{2}-|a_{-m}|^{2}+|b_{-m}|^{2}\right)/2-{\pi
k_{0}}/{2k_{c2}}$. Note that the coefficients $a_{\pm m}$ and $b_{\pm m}$ are
also functions of $k$.

Figure 3: The change of upper and lower thresholds $k_{c1}$ and $k_{c2}$ with
respect to $\Omega$ are denoted by red and blue lines, respectively. The
change of $k_{m}$ is also shown in green line for comparison. Orange dashed
lines indicate the critical value of $\Omega$, beyond which the stripe phase
is not ground state. $\bar{n}g/k_{0}^{2}=0.484$,
$\bar{n}g_{12}/k_{0}^{2}=0.44$ for (a), and $\bar{n}g/k_{0}^{2}=0.57$,
$\bar{n}g_{12}/k_{0}^{2}=0.456$ for (b).

##  IV Experimental observation

A stripe with $k$ away from $k_{m}$ can be accessed by exciting the spin
dipole mode of spin-orbit coupled BEC, which is associated with the temporal
change of stripe period [25]. Analogously, to examine the stability of plane
wave phase in spin-orbit coupled BEC, it has been previously proposed to
excite the dipole oscillation of the condensate, whose center-of-mass momentum
sweeps over a finite range around the minimum of dispersion [34, 24].

Here the spin dipole mode can be excited by first preparing the stripe at
ground state with the perturbation $-m\omega^{2}x_{0}x\sigma_{z}$ added to the
Hamiltonian, and then suddenly releasing the perturbation and observing the
temporal evolution of stripe wave function as well as oscillation of density
modulation period or $k$ [25]. Numerically we first find the ground state of
the system with the perturbation $-m\omega^{2}x_{0}x\sigma_{z}$ and a harmonic
trap $V(x)=m\omega^{2}x^{2}/2$ by evolving the GP equation in imaginary time.
In this way, the spin up/down component of stripe is shifted in opposite
directions at ground state. Subsequently, at time $t=0$, the perturbation is
released and real time evolution of stripe is obtained with the time-dependent
GP equation,

$i\hbar\frac{\partial}{\partial t}\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix}=\begin{pmatrix}H_{\uparrow\uparrow}&\frac{\Omega}{2}\\\
\frac{\Omega}{2}&H_{\downarrow\downarrow}\end{pmatrix}\begin{pmatrix}\Psi_{\uparrow}\\\
\Psi_{\downarrow}\end{pmatrix},$ (14)

where
$H_{\uparrow\uparrow}={\hbar^{2}(-i\partial_{x}-k_{0})^{2}}/{2m}+g|\Psi_{\uparrow}|^{2}+g_{12}|\Psi_{\downarrow}|^{2}+V(x)$,
and
$H_{\downarrow\downarrow}={\hbar^{2}(-i\partial_{x}+k_{0})^{2}}/{2m}+g|\Psi_{\downarrow}|^{2}+g_{12}|\Psi_{\uparrow}|^{2}+V(x)$.

Figures 4(a) and 4(c) depict the overall density distribution $|\Psi|^{2}$ of
stripe at two different moments for two different initial displacements
$x_{0}$. Also shown are the momentum distribution obtained by Fourier
transform of overall density (Fig. 4(d)), and the oscillation of $k$ with time
(Fig. 4(b)). We find that, for a small initial displacement $x_{0}$, the
overall profile of density remains similar and the momentum distribution
always features a peak centred at $k_{x}=0$ and a pair of peaks centred at
$k_{x}=\pm k$ which is associated with stripe period, as shown in Figs. 4(a)
and 4(b). In this case, the stripe has well-defined periodicity, and the
oscillation amplitude of $k$ around $k_{m}$ is small, lying within the stable
regime in Fig. 3(a) (see the orange dashed line within stable range). The
oscillation behaviour of $k$ is consistent with recent studies [25]. On the
other hand, for large initial $x_{0}$, the momentum distribution develops
multiple peaks centred at nonzero momentum, which means the crucial feature of
stripe is destroyed with no definite period, as shown in Figs. 4(c) and 4(d).
This behaviour qualitatively confirms our expectation that, when the
oscillation magnitude of $k$ is large for large initial $x_{0}$, the stripe
can suffer from dynamical instability, with significant change in the feature
of density distribution. So the stability region of $k$ effectively
corresponds to a limit on the amplitude of spin dipole oscillation. Note also
that the agreement here is qualitative, as there is a harmonic trap included
in the time evolution of GP equation but not present in the calculation of
Bogoliubov excitation spectrum. The agreement should be better for a weaker
trap with harmonic oscillator length $a_{x}$ much larger than stripe period.

Figure 4: (a) Overall density at two moments specified in the inset, for a
small initial displacement $x_{0}=2.5$. The stripe has well-defined period and
$k$ in this case. The temporal oscillation of $k$ is shown in (b). (c) Overall
density for a large initial $x_{0}=25$, with its momentum distribution at
$\omega t/2\pi=1.273$ shown in (d). The period and $k$ of stripe are not well-
defined in this case. Momentum distribution is obtained by Fourier transform
of density, with momentum divided by $2$ for comparison with $k$ defined in
this work. The trap frequency is $\hbar\omega/k_{0}^{2}=0.002$. The other
parameters are $\bar{n}g/k_{0}^{2}=0.484$, $\bar{n}g_{12}/k_{0}^{2}=0.44$ and
$\Omega/k_{0}^{2}=0.647$, same as Fig. 3(a).

##  V Conclusion and Acknowledgement

In summary, we have studied the properties of a class of metastable supersolid
stripes in spin-orbit coupled BEC. We find that for a stripe with general wave
number $k$, the nonlinear dispersion $\mu(k)$ can exhibit a loop structure.
Besides, when the wave number $k$ of stripes deviates from the one at ground
state, supersolids can still be stable with respect to low energy excitations.
There exist two thresholds $k_{c1}$ and $k_{c2}$, and for stripes with
$k_{c1}<k<k_{c2}$, their Bogoliubov excitation spectra are always real and
non-negative. For stripes with $k$ outside this regime, the excitation
spectrum at long wavelength becomes complex, signalling the dynamical
instability of those stripes. We also point out that, the thresholds of $k$
for stable stripes are related with the stability criteria of stripes in
response to spin dipole oscillation, which is accompanied by temporal
oscillation of $k$. For spin dipole oscillation of large amplitude, the
dramatic change in the feature of stripe wave functions can be qualitatively
explained by the instability mechanism we propose here. With the recent
experimental progress on compressional oscillations in dipolar supersolid
[15], it will be an interesting subject to explore similar stability
properties there.

W.X., T.L., and Q.Z. are supported by the National Key Research and
Development Program of China (Grant No. 2022YFA1405304), the National Natural
Science Foundation of China (Grant No. 12004118), and the Guangdong Basic and
Applied Basic Research Foundation (Grants No. 2020A1515110228 and No.
2021A1515010212). L.C. is supported by the National Natural Science Foundation
of China (Grant No. 12264061) and the Science Foundation of Guizhou Science
and Technology Department (Grant No. QKHJZ[2021]033). Y.Z. is supported by the
National Natural Science Foundation of China with Grants No. 11974235 and
11774219.

## References

  * Kim and Chan [2004a] E. Kim and M. H. W. Chan, Nature 427, 225 (2004a).
  * Kim and Chan [2004b] E. Kim and M. H. W. Chan, Science 305, 1941 (2004b).
  * Day and Beamish [2007] J. Day and J. Beamish, Nature 450, 853 (2007).
  * Kim and Chan [2012] D. Y. Kim and M. H. W. Chan, Phys. Rev. Lett. 109, 155301 (2012).
  * Wang  _et al._ [2010] C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010).
  * Ho and Zhang [2011] T.-L. Ho and S. Zhang, Phys. Rev. Lett. 107, 150403 (2011).
  * Li  _et al._ [2017] J.-R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F. Ç. Top, A. O. Jamison, and W. Ketterle, Nature 543, 91 (2017).
  * Bersano  _et al._ [2019] T. M. Bersano, J. Hou, S. Mossman, V. Gokhroo, X.-W. Luo, K. Sun, C. Zhang, and P. Engels, Phys. Rev. A 99, 051602(R) (2019).
  * Putra  _et al._ [2020] A. Putra, F. Salces-Cárcoba, Y. Yue, S. Sugawa, and I. B. Spielman, Phys. Rev. Lett. 124, 053605 (2020).
  * Léonard  _et al._ [2017] J. Léonard, A. Morales, P. Zupancic, T. Esslinger, and T. Donner, Nature 543, 87 (2017).
  * Tanzi  _et al._ [2019a] L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019a).
  * Böttcher  _et al._ [2019] F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019).
  * Chomaz  _et al._ [2019] L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. X 9, 021012 (2019).
  * Guo  _et al._ [2019] M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, Nature 574, 386 (2019).
  * Tanzi  _et al._ [2019b] L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature 574, 382 (2019b).
  * Norcia  _et al._ [2021] M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and F. Ferlaino, Nature 596, 357 (2021).
  * Petter  _et al._ [2021] D. Petter, A. Patscheider, G. Natale, M. J. Mark, M. A. Baranov, R. van Bijnen, S. M. Roccuzzo, A. Recati, B. Blakie, D. Baillie, L. Chomaz, and F. Ferlaino, Phys. Rev. A 104, L011302 (2021).
  * Li  _et al._ [2013] Y. Li, G. I. Martone, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 110, 235302 (2013).
  * Natale  _et al._ [2019] G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 123, 050402 (2019).
  * Li  _et al._ [2021] G.-Q. Li, X.-W. Luo, J. Hou, and C. Zhang, Phys. Rev. A 104, 023311 (2021).
  * Geier  _et al._ [2021] K. T. Geier, G. I. Martone, P. Hauke, and S. Stringari, Phys. Rev. Lett. 127, 115301 (2021).
  * Zhu  _et al._ [2012] Q. Zhu, C. Zhang, and B. Wu, EPL 100, 50003 (2012).
  * Zheng  _et al._ [2013] W. Zheng, Z.-Q. Yu, X. Cui, and H. Zhai, J. Phys. B 46, 134007 (2013).
  * Ozawa  _et al._ [2013] T. Ozawa, L. P. Pitaevskii, and S. Stringari, Phys. Rev. A 87, 063610 (2013).
  * Geier  _et al._ [2022] K. T. Geier, G. I. Martone, P. Hauke, W. Ketterle, and S. Stringari, arXiv:2210.10064 .
  * Lin  _et al._ [2011] Y.-J. Lin, K. Jiménez-García, and I. B. Spielman, Nature 471, 83 (2011).
  * Li  _et al._ [2012] Y. Li, L. P. Pitaevskii, and S. Stringari, Phys. Rev. Lett. 108, 225301 (2012).
  * Martone  _et al._ [2014] G. I. Martone, Y. Li, and S. Stringari, Phys. Rev. A 90, 041604(R) (2014).
  * Zhang  _et al._ [2019] Y. Zhang, Z. Gui, and Y. Chen, Phys. Rev. A 99, 023616 (2019).
  * Wu and Niu [2003] B. Wu and Q. Niu, New J. Phys. 5, 104 (2003).
  * Li  _et al._ [2019] C.-H. Li, C. Qu, R. J. Niffenegger, S.-J. Wang, M. He, D. B. Blasing, A. J. Olson, C. H. Greene, Y. Lyanda-Geller, Q. Zhou, C. Zhang, and Y. P. Chen, Nat. Commun. 10, 375 (2019).
  * Sun and Xie [2005] Q.-f. Sun and X. C. Xie, Phys. Rev. B 72, 245305 (2005).
  * Zhu  _et al._ [2015] Q. Zhu, Q.-f. Sun, and B. Wu, Phys. Rev. A 91, 023633 (2015).
  * Zhang  _et al._ [2012] J.-Y. Zhang, S.-C. Ji, Z. Chen, L. Zhang, Z.-D. Du, B. Yan, G.-S. Pan, B. Zhao, Y.-J. Deng, H. Zhai, S. Chen, and J.-W. Pan, Phys. Rev. Lett. 109, 115301 (2012).


# From Inverse Optimization to Feasibility to ERM

Saurabh Mishra  Anant Raj  Sharan Vaswani

Inverse Optimization, Predict and Optimize, Convex Feasibility

## Impact Statement

This paper presents work whose goal is to advance the field of Machine
Learning. There are many potential societal consequences of our work, none of
which we feel must be specifically highlighted here.

## Acknowledgements

We would like to thank Michael Lu, Anh Dang, Reza Babanezhad, Vibhuti Dhingra,
Karan Desai, Tarannum Khan and Hardik Chauhan for helpful feedback on the
paper.

## References

  * Amos (2019) Amos, B.  Differentiable optimization-based modeling for machine learning.  _Ph. D. thesis_ , 2019.
  * Amos & Kolter (2017) Amos, B. and Kolter, J. Z.  Optnet: Differentiable optimization as a layer in neural networks.  In _International Conference on Machine Learning_ , pp. 136–145. PMLR, 2017.
  * Angalakudati et al. (2014) Angalakudati, M., Balwani, S., Calzada, J., Chatterjee, B., Perakis, G., Raad, N., and Uichanco, J.  Business analytics for flexible resource allocation under random emergencies.  _Management Science_ , 60:1552–1573, 2014.
  * Armijo (1966) Armijo, L.  Minimization of functions having Lipschitz continuous first partial derivatives.  _Pacific Journal of Mathematics_ , 16(1):1 – 3, 1966.
  * Bansal (2005) Bansal, R.  Optimization methods for electric power systems: An overview.  _International Journal of Emerging Electric Power Systems_ , 2, 2005.
  * Bassily et al. (2018) Bassily, R., Belkin, M., and Ma, S.  On exponential convergence of sgd in non-convex over-parametrized learning.  _arXiv preprint arXiv:1811.02564_ , 2018.
  * Bauschke & Borwein (1993) Bauschke, H. H. and Borwein, J. M.  On the convergence of von neumann’s alternating projection algorithm for two sets.  _Set-Valued Analysis_ , 1:185–212, 1993.
  * Bauschke & Borwein (1996) Bauschke, H. H. and Borwein, J. M.  On projection algorithms for solving convex feasibility problems.  _SIAM review_ , 38(3):367–426, 1996.
  * Berthet et al. (2020) Berthet, Q., Blondel, M., Teboul, O., Cuturi, M., Vert, J.-P., and Bach, F.  Learning with differentiable pertubed optimizers.  _Advances in neural information processing systems_ , 33:9508–9519, 2020.
  * Bertsimas & Kallus (2020) Bertsimas, D. and Kallus, N.  From predictive to prescriptive analytics.  _Management Science_ , 66:1025–1044, 2020.
  * Bertsimas et al. (2015) Bertsimas, D., Gupta, V., and Paschalidis, I. C.  Data-driven estimation in equilibrium using inverse optimization.  _Mathematical Programming_ , 153:595–633, 2015.
  * Besbes et al. (2023) Besbes, O., Fonseca, Y., and Lobel, I.  Contextual inverse optimization: Offline and online learning.  _Operations Research_ , 2023.
  * Birge et al. (2017) Birge, J. R., Hortaçsu, A., and Pavlin, J. M.  Inverse optimization for the recovery of market structure from market outcomes: An application to the miso electricity market.  _Operations Research_ , 65(4):837–855, 2017.
  * Bousquet & Elisseeff (2002) Bousquet, O. and Elisseeff, A.  Stability and generalization.  _The Journal of Machine Learning Research_ , 2:499–526, 2002.
  * Boyd & Vandenberghe (2004) Boyd, S. P. and Vandenberghe, L.  _Convex optimization_.  Cambridge university press, 2004.
  * Boyd et al. (2001) Boyd, S. P., Lee, T. H., et al.  Optimal design of a cmos op-amp via geometric programming.  _IEEE Transactions on Computer-aided design of integrated circuits and systems_ , 20:1–21, 2001.
  * Cameron et al. (2022) Cameron, C., Hartford, J., Lundy, T., and Leyton-Brown, K.  The perils of learning before optimizing.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pp. 3708–3715, 2022.
  * Chan et al. (2022) Chan, T. C., Eberg, M., Forster, K., Holloway, C., Ieraci, L., Shalaby, Y., and Yousefi, N.  An inverse optimization approach to measuring clinical pathway concordance.  _Management Science_ , 68(3):1882–1903, 2022.
  * Chan et al. (2023) Chan, T. C., Mahmood, R., and Zhu, I. Y.  Inverse optimization: Theory and applications.  _Operations Research_ , 2023.
  * Cornuejols & Tütüncü (2006) Cornuejols, G. and Tütüncü, R.  _Optimization methods in finance_ , volume 5.  Cambridge University Press, 2006.
  * Deutsch (1984) Deutsch, F.  Rate of convergence of the method of alternating projections.  _Parametric optimization and approximation (Oberwolfach, 1983)_ , 72:96–107, 1984.
  * Diamond & Boyd (2016) Diamond, S. and Boyd, S.  Cvxpy: A python-embedded modeling language for convex optimization.  _The Journal of Machine Learning Research_ , 17:2909–2913, 2016.
  * Domahidi et al. (2013) Domahidi, A., Chu, E., and Boyd, S.  Ecos: An socp solver for embedded systems.  In _2013 European control conference (ECC)_ , pp. 3071–3076. IEEE, 2013.
  * Duchi et al. (2011) Duchi, J., Hazan, E., and Singer, Y.  Adaptive subgradient methods for online learning and stochastic optimization.  _Journal of machine learning research_ , 12(7), 2011.
  * El Balghiti et al. (2019) El Balghiti, O., Elmachtoub, A. N., Grigas, P., and Tewari, A.  Generalization bounds in the predict-then-optimize framework.  _Advances in neural information processing systems_ , 32, 2019.
  * Elmachtoub & Grigas (2022) Elmachtoub, A. N. and Grigas, P.  Smart “predict, then optimize”.  _Management Science_ , 68(1):9–26, 2022.
  * Fabozzi et al. (2008) Fabozzi, F. J., Markowitz, H. M., and Gupta, F.  Portfolio selection.  _Handbook of finance_ , 2, 2008.
  * Garrigos (2023) Garrigos, G.  Square distance functions are polyak-$\\{$$\backslash$L$\\}$ ojasiewicz and vice-versa.  _arXiv preprint arXiv:2301.10332_ , 2023.
  * Geng et al. (2023) Geng, H., Ruan, H., Wang, R., Li, Y., Wang, Y., Chen, L., and Yan, J.  Rethinking and benchmarking predict-then-optimize paradigm for combinatorial optimization problems.  _arXiv preprint arXiv:2311.07633_ , 2023.
  * Gower et al. (2021) Gower, R., Sebbouh, O., and Loizou, N.  Sgd for structured nonconvex functions: Learning rates, minibatching and interpolation.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1315–1323. PMLR, 2021.
  * Guyomarch (2017) Guyomarch, J.  Warcraft ii open-source map editor.  _URL http://github. com/war2/war2edit_ , 2017.
  * Hardt et al. (2016) Hardt, M., Recht, B., and Singer, Y.  Train faster, generalize better: Stability of stochastic gradient descent.  In _International conference on machine learning_ , pp. 1225–1234. PMLR, 2016.
  * He et al. (2016) He, K., Zhang, X., Ren, S., and Sun, J.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 770–778, 2016.
  * Heuberger (2004) Heuberger, C.  Inverse combinatorial optimization: A survey on problems, methods, and results.  _Journal of combinatorial optimization_ , 8:329–361, 2004.
  * Iyengar & Kang (2005) Iyengar, G. and Kang, W.  Inverse conic programming with applications.  _Operations Research Letters_ , 33(3):319–330, 2005.
  * Karimi et al. (2016) Karimi, H., Nutini, J., and Schmidt, M.  Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition.  In _Machine Learning and Knowledge Discovery in Databases: European Conference, ECML PKDD 2016, Riva del Garda, Italy, September 19-23, 2016, Proceedings, Part I 16_ , pp. 795–811. Springer, 2016.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kuhn & Tucker (1951) Kuhn, H. W. and Tucker, A. W.  Nonlinear programming, paper presented at proceedings of the second berkeley symposium on mathematical statistics and probability, 1951.
  * Lavington et al. (2023) Lavington, J. W., Vaswani, S., Babanezhad, R., Schmidt, M., and Roux, N. L.  Target-based surrogates for stochastic optimization.  _arXiv preprint arXiv:2302.02607_ , 2023.
  * Lei & Ying (2020a) Lei, Y. and Ying, Y.  Fine-grained analysis of stability and generalization for stochastic gradient descent.  In _International Conference on Machine Learning_ , pp. 5809–5819. PMLR, 2020a.
  * Lei & Ying (2020b) Lei, Y. and Ying, Y.  Sharper generalization bounds for learning with gradient-dominated objective functions.  In _International Conference on Learning Representations_ , 2020b.
  * Lei & Ying (2021) Lei, Y. and Ying, Y.  Sharper generalization bounds for learning with gradient-dominated objective functions.  In _International Conference on Learning Representations_ , 2021.
  * Li et al. (2022) Li, B., Wu, G., He, Y., Fan, M., and Pedrycz, W.  An overview and experimental study of learning-based optimization algorithms for the vehicle routing problem.  _IEEE/CAA Journal of Automatica Sinica_ , 9:1115–1138, 2022.
  * Li et al. (2018) Li, J., Liu, F., Wang, Z., Low, S. H., and Mei, S.  Optimal power flow in stand-alone dc microgrids.  _IEEE Transactions on Power Systems_ , 33:5496–5506, 2018.
  * Liu et al. (2022) Liu, C., Zhu, L., and Belkin, M.  Loss landscapes and optimization in over-parameterized non-linear systems and neural networks.  _Applied and Computational Harmonic Analysis_ , 59:85–116, 2022.
  * Liu et al. (2023) Liu, C., Drusvyatskiy, D., Belkin, M., Davis, D., and Ma, Y.-A.  Aiming towards the minimizers: fast convergence of sgd for overparametrized problems.  _arXiv preprint arXiv:2306.02601_ , 2023.
  * Luenberger et al. (1984) Luenberger, D. G., Ye, Y., et al.  _Linear and nonlinear programming_ , volume 2.  Springer, 1984.
  * Ma et al. (2018) Ma, S., Bassily, R., and Belkin, M.  The power of interpolation: Understanding the effectiveness of sgd in modern over-parametrized learning.  In _International Conference on Machine Learning_ , pp. 3325–3334. PMLR, 2018.
  * Mohajerin Esfahani et al. (2018) Mohajerin Esfahani, P., Shafieezadeh-Abadeh, S., Hanasusanto, G. A., and Kuhn, D.  Data-driven inverse optimization with imperfect information.  _Mathematical Programming_ , 167:191–234, 2018.
  * Ng et al. (2000) Ng, A. Y., Russell, S., et al.  Algorithms for inverse reinforcement learning.  In _Icml_ , volume 1, pp. 2, 2000.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Polyak (1964) Polyak, B. T.  Gradient methods for solving equations and inequalities.  _USSR Computational Mathematics and Mathematical Physics_ , 4(6):17–32, 1964.
  * Raj & Bach (2021) Raj, A. and Bach, F.  Explicit regularization of stochastic gradient methods through duality.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 1882–1890. PMLR, 2021.
  * Raja & Pugazhenthi (2012) Raja, P. and Pugazhenthi, S.  Optimal path planning of mobile robots: A review.  _International journal of physical sciences_ , 7:1314–1320, 2012.
  * Robbins & Monro (1951) Robbins, H. and Monro, S.  A stochastic approximation method.  _The annals of mathematical statistics_ , pp. 400–407, 1951.
  * Sahoo et al. (2022) Sahoo, S. S., Paulus, A., Vlastelica, M., Musil, V., Kuleshov, V., and Martius, G.  Backpropagation through combinatorial algorithms: Identity with projection works.  _arXiv preprint arXiv:2205.15213_ , 2022.
  * Stellato et al. (2020) Stellato, B., Banjac, G., Goulart, P., Bemporad, A., and Boyd, S.  Osqp: An operator splitting solver for quadratic programs.  _Mathematical Programming Computation_ , 12(4):637–672, 2020.
  * Sun et al. (2023) Sun, C., Liu, S., and Li, X.  Maximum optimality margin: A unified approach for contextual linear programming and inverse linear programming.  _arXiv preprint arXiv:2301.11260_ , 2023.
  * Vaswani et al. (2019) Vaswani, S., Bach, F., and Schmidt, M.  Fast and faster convergence of sgd for over-parameterized models and an accelerated perceptron.  In _The 22nd international conference on artificial intelligence and statistics_ , pp. 1195–1204. PMLR, 2019.
  * Vlastelica et al. (2019) Vlastelica, M., Paulus, A., Musil, V., Martius, G., and Rolínek, M.  Differentiation of blackbox combinatorial solvers.  _arXiv preprint arXiv:1912.02175_ , 2019.
  * Von Neumann (1949) Von Neumann, J.  On rings of operators. reduction theory.  _Annals of Mathematics_ , pp. 401–485, 1949.
  * Wahdany et al. (2023) Wahdany, D., Schmitt, C., and Cremer, J. L.  More than accuracy: end-to-end wind power forecasting that optimises the energy system.  _Electric Power Systems Research_ , 221:109384, 2023.
  * Wilder et al. (2019) Wilder, B., Dilkina, B., and Tambe, M.  Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 33, pp. 1658–1665, 2019.
  * Zhang & Yin (2013) Zhang, H. and Yin, W.  Gradient methods for convex minimization: better rates under weaker conditions.  _arXiv preprint arXiv:1303.4645_ , 2013.


# Diversity Measures:
Domain-Independent Proxies for Failure in Language Model Queries

Noel Ngu, Nathaniel Lee, Paulo Shakarian

###### Abstract

Error prediction in large language models often relies on domain-specific
information. In this paper, we present measures for quantification of error in
the response of a large language model based on the diversity of responses to
a given prompt - hence independent of the underlying application. We describe
how three such measures - based on entropy, Gini impurity, and centroid
distance - can be employed. We perform a suite of experiments on multiple
datasets and temperature settings to demonstrate that these measures strongly
correlate with the probability of failure. Additionally, we present empirical
results demonstrating how these measures can be applied to few-shot prompting,
chain-of-thought reasoning, and error detection.

##  1 Introduction

Prompt engineering techniques for large language models (LLMs) have yielded
significant results in reducing errors. Advances such as self-consistency and
chain-of-thought prompting (Wei et al. 2022; Wang et al. 2022; Yao et al.
2023), in particular, have shown the ability to reduce errors for a given
problem. In this paper, we take a different, but complementary approach to
these methods. Specifically, we study diversity measures that can function as
a domain-independent proxy for failure probability for question-answer tasks
that rely on the prompting of the LLM. The intuition is that the diversity of
the response suggests uncertainty. In our approach, we avoid requirements of
direct access to the LLM (Kuhn, Gal, and Farquhar 2023) as well as semantic
understanding (which is often domain-specific) (Lin, Trivedi, and Sun 2023;
Kuhn, Gal, and Farquhar 2023; Shakarian et al. 2023). Our contributions are as
follows.

  * •

We formalize three diversity measures with respect to LLM prompting (Section 3
and show, across three datasets and five temperature settings that all three
measures correlate closely with the probability of failure with $R^{2}$ values
exceeding $0.8$ in most cases and often higher (Sections 5.2-5.3). We note
that these measures do not require access to the LLM or domain knowledge.

  * •

We introduce the concept of diversity-based prompt selection (Section 4.1).
Experimentally (Section 5.4), we show in two datasets, significant improvement
($14\%,45\%$) over baselines (and near the baseline in a third). In all cases,
this method significantly improved over worst-case prompt outcomes.

  * •

We also extend the results on prompt selection for use with chain-of-thought
(CoT) prompting (Wei et al. 2022). We show experimentally (Section 5.5) that
CoT and diversity-driven prompting converge to the same failure probability.
When combined, failure probability is reduced by $30\%$ as compared to when
either technique is used in isolation.

  * •

We describe how such measures can be used in machine learning models for
predicting failure (Section 4.2) and show experimentally (Section 5.6) that
they can be viable for error prediction. We show how diversity measures can be
used in a supervised model to predict errors significantly outperforming the
baseline, in the best case obtaining an AUPRC of $0.703$.

The rest of the paper is organized as follows. In Section 2 we review related
work. In Section 3 we discuss our approach, providing specifics on how we
compute entropy, Gini impurity, and centroid-based diversity measures with
respect to the output of self-consistency-style prompting (Wang et al. 2022).
This is followed by a discussion on how the approach is applied to diversity-
based prompt selection and error prediction in Section 4. We then present
experimental evaluations in Section 5 where we study how the diversity
measures correlated with failure, diversity-based prompt selection, synergies
with CoT prompting, and how such measures can be used with machine learning
approaches to predict errors. We provide an implementation available at
https://github.com/lab-v2/diversity˙measures

##  2 Related Work

There have been numerous recent studies evaluating the performance of LLMs in
various settings (Liu et al. 2023; Shakarian et al. 2023; Frieder et al. 2023;
Tu et al. 2023; Ray 2023; Kocoń et al. 2023; Yuan et al. 2023; Pardos and
Bhandari 2023; OpenAI 2023) these studies compare LLM performance across
tasks, models, and datasets. While these lead to important insights, they
generally do not look into the use of specific metrics of the output to
improve performance in the general case. One exception is (Shakarian et al.
2023) which studies indicators of ChatGPT failure with respect to math word
problems by counting the mathematical operators in the template equations
present in the metadata of the DRAW-1K dataset. Drawing on this idea, a simple
ML model was created based on counts of various mathematical operations. The
approach in this paper is not tied specifically to mathematical operations.

Examining the diversity of LLM outputs has also been a topic of previous
studies. The selection of tokens using various methods such as temperature,
top-$k$, and nucleus sampling (Ficler and Goldberg 2023; Fan, Lewis, and
Dauphin 2018; Holtzman et al. 2019) to create the output has been studied.
These results have generally been integrated into the system design of
contemporary LLMs. Essentially, these findings allow for the selection of
tokens that appear in the generated text of an LLM. As the methods are
probabilistic, multiple runs lead to slightly different results. The
aforementioned work provides hyper-parameters to adjust the variance of the
LLM response to a given prompt. Recent work on prompt engineering (Wei et al.
2022; Wang et al. 2022; Yao et al. 2023) has leveraged this stochastic
property to improve results by taking a sample of LLM responses to a given
prompt. However, that work focuses on leveraging a set of LLM outputs to
improve accuracy for a given problem, rather than on using such outputs to
assess the correctness of the result.

This work is complementary to other active areas of research such as meta-
learning (Hospedales et al. 2022; Zhou et al. 2022) and introspection (Daftry
et al. 2016; Ramanagopal et al. 2018). However, these methods often rely on
features to predict failure or make a determination of what model to use. The
diversity measures introduced in this paper may enable such techniques to be
leveraged in the context of LLMs.

Concurrent with this work there have been several other studies either
recently published or in preprint on language model uncertainty. For example,
(Zhou, Jurafsky, and Hashimoto 2023) found that the presence of words relating
to uncertainty in the result of an LLM lead is indicative of errors. Studies
such as (Tian et al. 2023) study how to leverage human calibration to create
confidence sources for language model results. In (Agrawal, Mackey, and Kalai
2023) the authors look to evaluate LLM results for hallucinated sources by
comparison with external sources. The work of (Kuhn, Gal, and Farquhar 2023)
also uses entropy measurements, but relies on access to the token
probabilities of the LLM (which we do not assume) and utilizes an additional
step of semantic comparison among pairs of language model results - which
requires a second supervised model. The recent preprint of (Lin, Trivedi, and
Sun 2023) builds on that work by creating a second trained model to
approximate token probabilities. In this work, we avoid the training of an
additional second model, pairwise comparisons, or access to the original LLM.

##  3 Approach

Technical Preliminaries. Our intuition is to build on the concept of self-
consistency (Wang et al. 2022). For a given prompt-question pair
$\mathcal{P},\mathcal{Q}$, the decoder returns a reasoning path, answer pair
$\mathbf{r},\mathbf{a}$ where $\mathbf{a}$ is an element of some (possibly
infinite sized) answer set $\mathbb{A}$. In this paradigm, the decoder is
queried $m$ times producing $m$ such $\mathbf{r},\mathbf{a}$ pairs. The final
answer is computed based on the answer $a\in\mathbb{A}$ that maximized the
following conditional probability.

$\displaystyle
P(a|\mathcal{P},\mathcal{Q})=\sum_{i=1}^{m}\mathbf{1}(\mathbf{a}_{i}=a)P(\mathbf{r}_{i},\mathbf{a}_{i}|\mathcal{P},\mathcal{Q})$
(1)

In practice, Equation 1 can be accurately approximated with a majority vote
over the $m$ answers (shown in Equation 2). This has the nice property of
avoiding the requirement of knowing
$P(\mathbf{r}_{i},\mathbf{a}_{i}|\mathcal{P},\mathcal{Q})$ \- which is not
available for most commercial models and is a core assumption in this paper.
Hence, we shall use Equation 2 as the definition of this probability.

$\displaystyle
P(a|\mathcal{P},\mathcal{Q})\approx\frac{1}{m}\sum_{i=1}^{m}\mathbf{1}(\mathbf{a}_{i}=a)$
(2)

Diversity Measures. Our intuition is that given answers
$\mathbf{a}_{1},\ldots,\mathbf{a}_{m}$ the diversity among the answers
provides insights into the correctness of the decoder’s response. We consider
two paradigms in which to view a given $\mathbf{a}_{i}$. First, for “set-
based” output we shall assume a universe of elements $\mathbb{U}$ and that
each $\mathbf{a}_{i}\subseteq\mathbb{U}$. The intuition here is that this view
would be suitable for applications such as math and logic problems where the
$\mathbf{a}_{i}$ is a set of natural numbers or symbols. In the second
paradigm, “vector-based” output, we use an embedding function $\mathit{embed}$
that maps each $\mathbf{a}_{i}$ into a $D$ dimensional space
($\mathit{embed}:\mathbb{A}\mapsto\mathbb{R}^{D}$). This is suitable where a
set representation of the output is too coarse-grain or where the reasoning
path ($\mathbf{r}_{i}$) is considered part of the answer itself.

Set-Based Output. With set-based output, each $\mathbf{a}_{i}$ is a set of
elements from some universe $\mathbb{U}$. For a given element
$e\in\mathbb{U}$, we will define the probability of $e$ appearing in the
result as:

$\displaystyle P(e)=\frac{1}{m}\sum_{i=1}^{m}\mathbf{1}(e\in\mathbf{a}_{i})$
(3)

With this quantity, we can directly apply Entropy and Gini impurity measures
to create diversity measurements as follows.

$\displaystyle H=-\sum_{e\in\mathbb{U}}P(e)\log(P(e))$ (4) $\displaystyle
G=1-\sum_{e\in\mathbb{U}}P(e)^{2}$ (5)

Vector-Based Output. Here we look to evaluate each $\mathbf{a}_{i}$ by
examining its vector embedding of dimension $D$. We will use a distance
function $d:\mathbb{R}^{D}\times\mathbb{R}^{D}\mapsto\mathbb{R}^{+}$ to
compare two vectors of size $D$. Given the vectors created by the
$\mathit{embed}$ function
$\mathit{embed}(\mathbf{a}_{1}),\ldots,\mathit{embed}(\mathbf{a}_{m})$ we can
calculate a $D$ dimension centroid $\mathit{cent}$. Note that we can vary the
calculation of $\mathit{embed},d,$ and $\mathit{cent}$ and explore such
variations empirically. Using this information, we propose the following
diversity measure.

$\displaystyle\frac{1}{m}\sum_{i}d(\mathit{embed}(a_{i}),\mathit{cent})$ (6)

##  4 Applications

In this section, we describe two applications of diversity measures introduced
in the previous section: diversity-based prompt selection (Section 4.1) and
error prediction (Section 4.2).

###  4.1 Diversity-Based Prompt Selection

A key problem in developing LLM-powered applications is the design of the
prompt. A diversity measure, such as entropy, can be used to evaluate prompts
during training but can also be used in a strategy that employs multiple
prompts in what we term diversity-based prompt selection. For example, suppose
we have $N$ different prompts, each queried $m$ times. Using majority voting,
we get $N$ responses. We can then measure each response by calculating the
diversity (e.g., entropy) over the $m$ queries associated with each response
and take the response with the lowest diversity.

###  4.2 Error Prediction

Following the intuition of work on introspection (Daftry et al. 2016;
Ramanagopal et al. 2018), we look to train a machine learning model in a
supervised fashion to predict failures for a given use case. The input to the
model consists of entropy, Gini impurity, and centroid-based features.

##  5 Evaluation

Using three datasets described in Section 5.1, we examine the entropy, Gini
impurity, and distance from centroid measures in each of the three datasets
and report on their correlation with failure probability in Sections 5.2 for
set-based measures and Section 5.3 for vector-based measures. We then describe
experiments on our application results focused on diversity-based prompt
selection (Section 5.4) and extend those results for CoT prompting (Section
5.5). We then look at the application to failure prediction by describing our
supervised model in Section 5.6.

###  5.1 Experimental Setup

In all of our experiments, we used OpenAI’s GPT 3.5 model for the LLM and
interfaced with the LLM through the OpenAI API using Python 3.x. For
experiments on vector-based output, the $\mathit{embed}$ function was
implemented using sentence-BERT (Reimers and Gurevych 2019).

We evaluated three datasets listed in Table 1. These datasets were chosen with
the intention of examining LLM performance across a diverse set of reasoning
tasks (commonsense, arithmetic, and symbolic reasoning.) CSQA consists of
$9,741$ samples, DRAW-1K consists of $1,000$ samples, and LL consists of
$3,000$ samples. Splits on training and testing as well as temperature
settings are described for individual experiments.

Dataset | Domain | Response
---|---|---
|  | Type
CSQA | Knowledge- | Multiple
(Saha et al. 2018) | based Q&A | choice
DRAW-1K | Math word | Numerical
(Upadhyay and Chang 2017) | problems | answers
LL | Last-letter | Textual
(Diao et al. 2023) | concatenation | concatenation
Table 1: Datasets used in evaluation.

###  5.2 Measure Evaluation: Set-Based Output

We evaluated entropy and Gini impurity measures on the entirety of each of the
three datasets with temperature settings of $\\{0.3,0.5,0.7,0.8,0.9\\}$. We
examine cumulative distributions (both for minimum and maximum values) for all
values of the measure with a minimum bucket size of 100 samples. Entropy
results are shown in Figure 1 while Gini impurity results are shown in Figure
2.

Figure 1: Entropy vs. failure probability (top row: CSQA, middle row: DRAW-1K,
bottom row: LL)

Entropy Results. Entropy measures appeared to be strongly related to the
probability of failure for cumulative plots in both directions. In most cases,
the relationship appeared linear with $R^{2}$ values (for cumulative minimum
values with temperature of $0.7$) at $0.808$ for CSQA, $0.973$ for DRAW-1K,
and $0.876$ for LL. Relationships had a tighter linear fit for cumulative
plots based on maximum values with $R^{2}$ values (temperature of $0.7$) of
$0.799$ for CSQA, $0.978$ for DRAW-1K, and $0.979$ for LL.

Figure 2: Gini impurity vs. failure probability (top row: CSQA, middle row:
DRAW-1K, bottom row: LL)

Gini Impurity Results. Gini impurity was also closely related to entropy in
this case and produced similar trends. We included a plot of the cumulative
relationships in Figure 2.

###  5.3 Measure Evaluation: Vector-Based Output

For vector-based output, we evaluated the average distance to the centroid
across the same datasets and temperature settings for set-based output (see
Section 5.2 for details). Cumulative plots are shown in Figure 3. In cases
where the failure probability is less than $0.5$, we observed a clear
monotonic relationship between the average distance to the centroid and
failure probability for cumulative plots in both directions. For cumulative
plots based on the maximum average distance from the centroid, the
relationship becomes less interesting in CSQA and LL for the temperature
settings of $0.3$ and $0.5$ \- exhibiting little slope. For cumulative plots
based on the minimum value, we see the trend reversing in the case of LL for
these temperature settings - and a similar reversal in the same plots of CSQA
across all temperature settings once the maximum failure probability is
crested. This contradicts our original intuition - that diversity will, in
general, lead to uncertainty in the result (and increase the chance of an
incorrect response being returned). However, in these cases of high failure
probability and extreme temperature, such diversity may actually lead to a
correct answer being produced by chance. It is notable that this phenomenon is
less pronounced in DRAW-1K, likely due to that problem not being multiple
choice (making guessing more prone to error).

Figure 3: Average distance to centroid vs. failure probability (top row: CSQA,
middle row: DRAW-1K, bottom row: LL)

###  5.4 Prompt Selection Experiments

As described earlier, we studied diversity measure’s applications to prompt
engineering, studying the effects of selecting the prompt providing the lowest
entropy response. Here we generated $20$ prompts from each dataset by randomly
selecting $30$ samples with ground-truth responses to serve as few-shot
samples for each prompt. We then randomly selected 100 samples from each
dataset, ensuring that these samples were drawn from a population not
previously used for any prompt, to create a test set. Results are shown in
Figure 4. The panels on the left show the failure probabilities for each of
the $20$ prompts used and compare them with the selection of a prompt based on
minimum entropy, Gini impurity, and average distance to the centroid (note for
each query, all prompts are used and the one providing the lowest diversity
measure is used to return the result).

Entropy and Gini Impurity Prompt Selection Results. Comparisons with
individual prompts are shown in the left-hand column of Figure 4. Entropy and
Gini impurity measures performed comparably. Selection by both measures
outperformed $19$ of the $20$ prompts and obtained the same probability of
failure on the one prompt it did not outperform. On average, it provided a
$14\%$ improvement (i.e., an average difference in failure probability of
$0.04$) over the CSQA prompts. Performance was more pronounced with DRAW-1K,
which also outperformed $19$ of the $20$ prompts, but on average provided
$45\%$ improvement (i.e., an average difference in probability of $0.17$) for
entropy, with a slightly better performance obtained with Gini impurity.
Performance with LL was less impressive. However, unlike the other datasets,
the failure probability for all prompts in this experiment exceeded $0.5$.
Here, minimum entropy selection outperformed $6$ of the $20$ prompts, but on
average did $1.6\%$ worse (i.e., average increase in probability by $0.01$).
However, in the cases where it did outperform the prompts, it provided a
$2.6\%$ improvement - so the approach at least shows it can mitigate against
poorly-performing prompts.

We note that the use of $20$ prompts per problem (and each used with $20$
queries) can become an expensive operation, especially considering API and
runtime costs. Therefore, we examined in the right-hand panels of Figure 4 the
failure probability as a function of the number of prompts used. Note that an
average in this case functions as an oracle, as operationally we would not
know which prompt to select. We note that our method (when using entropy or
Gini impurity) showed the ability to reduce failure probability below the
average without requiring all $20$ prompts, stably reaching a lower failure
probability with 5 prompts for CSQA, and 2 prompts for DRAW-1K. We did not see
consistent improvement over average with LL (the only case where failure
exceeded $0.5$) but did see the failure probability of minimum entropy
selection improve over the worst-case stay within $1.5\%$ of the average with
only $5$ prompts.

Centroid Distance Prompt Selection Results. Centroid-based prompt selection
consistently outperformed the worst-case scenario, but under-performs
(relative to entropy and Gini-based prompt selection) with respect to the
average case - in this regard, the result was poorest for DRAW-1K ($18\%$
worse than the average prompt failure probability). This is likely due to
vector representations of numeric answers near each other in the vector space
differing significantly in semantics (e.g., a missing negative sign may result
in such a phenomenon).

Figure 4: Prompt engineering experiment. Top row: CSQA, middle row: DRAW-1K,
bottom row: LL

###  5.5 CoT Prompting

Chain of thought (CoT) prompting (Wei et al. 2022) utilizes explanations in
the few-shot examples. One would suspect that the use of CoT prompting would
limit the variety of responses to a prompt, hence limiting the impact of a
diversity technique employed to select the correct prompt. As the DRAW-1K
dataset included the associated math equations with ground truth, we used this
information to replicate the experiment of Section 5.4 with CoT-style prompts.

Entropy and Gini Impurity Results (with CoT). Again, entropy and Gini impurity
performed comparably. First, we observe (in Figure 5, left) that both CoT
(without entropy/Gini) and minimum entropy/Gini selection (without CoT)
converge to a failure probability of $0.31$ which suggest that both are
reducing diversity.

However, we believe that few-shot prompting and CoT reduce diversity in
different ways, as there was a synergistic effect from combining them. Using
only $3$ prompts, minimum entropy/Gini selection with CoT consistently
provides over a $30\%$ reduction in failure probability over either technique
alone. We also observe the two techniques together provide improved failure
probability in $18$ of the $20$ CoT prompts (Figure 5, right).

Centroid Distance Results (with CoT). Again, centroid distance under-performs
when compared to entropy and Gini impurity, though combined with CoT provides
better results than without. We also note that DRAW-1K is the dataset for
which centroid distance performs the poorest with standard prompt selection
and underperformance is likely due to similar issues as described in Section
5.4.

Figure 5: Prompt engineering experiment using Chain of Thought (CoT) on
DRAW-1K.

###  5.6 Error Prediction Experiments

We created machine learning models based on the features described in Section
4.2 and trained models (specific to each dataset) with a training dataset
balanced between cases where the prompt failed and where it was successful. To
balance the training data, the majority class was randomly oversampled. We
then evaluated it on a separate test dataset that was not artificially
balanced. Precision-recall curves for a 10-layer multi-layer perceptron
(implemented in TensorFlow) are shown in Figure 6. We note that in all cases,
the models provide significant improvement in precision over a knowledgless
selection. Further, for low recall setting (under $0.2$) we obtain over $0.7$
precision for DRAW-1K and over $0.8$ precision in LL.

Figure 6: Precision-Recall Curves for the MLP-10 failure prediction model CSQA
(left), DRAW-1K (center), and LL (right) for temperature=0.7, dashed lines
show knowledgeless baseline

Ablation Study. We studied the effects of ablating different diversity
measures within the MLP-10 model. This model was trained to optimize accuracy
in predicting overall performance in the LLM for each of the three data sets.
Results are depicted in Figure 7. Ablating entropy or Gini impurity separately
generally had little effect, which was unsurprising as other experiments in
this paper demonstrated that they behave similarly. However, ablating both
results in a significant reduction in performance, most notably with DRAW-1K.

Figure 7: Ablation of various diversity measures for the failure prediction
model. Accuracy, precision, recall, F1, and AUPRC are all with respect to
predicting LLM failure (MLP-10 for classifier, temperature=0.7).

Alternative Classification Methods. We also examined precision and recall
results for a variety of classifiers including dense networks for various
depths, XGBoost, and AdaBoost. Results for CSQA, DRAW-1K, and LL are shown in
Tables 2, 3, and 4 respectively.

|  | Incorrect | Correct
---|---|---|---
Classifier | Acc. | Prec. | Rec. | Prec. | Rec.
MLP 10 | 0.667 | 0.592 | 0.461 | 0.699 | 0.798
MLP 15 | 0.667 | 0.593 | 0.461 | 0.699 | 0.799
XGB | 0.628 | 0.525 | 0.474 | 0.685 | 0.727
AdaBoost | 0.665 | 0.589 | 0.458 | 0.698 | 0.796
Table 2: CSQA failure prediction results |  | Incorrect | Correct
---|---|---|---
Classifier | Acc. | Prec. | Rec. | Prec. | Rec.
MLP 5 | 0.748 | 0.398 | 0.701 | 0.920 | 0.759
MLP 10 | 0.762 | 0.414 | 0.690 | 0.919 | 0.778
MLP 15 | 0.757 | 0.412 | 0.696 | 0.919 | 0.771
XGB | 0.769 | 0.379 | 0.381 | 0.860 | 0.857
AdaBoost | 0.741 | 0.380 | 0.636 | 0.904 | 0.765
Table 3: DRAW-1K failure prediction results |  | Incorrect | Correct
---|---|---|---
Classifier | Acc. | Prec. | Rec. | Prec. | Rec.
MLP 5 | 0.657 | 0.596 | 0.706 | 0.733 | 0.620
MLP 10 | 0.657 | 0.613 | 0.738 | 0.766 | 0.613
MLP 15 | 0.680 | 0.637 | 0.656 | 0.730 | 0.698
XGB | 0.669 | 0.619 | 0.644 | 0.714 | 0.688
AdaBoost | 0.678 | 0.614 | 0.734 | 0.756 | 0.635
Table 4: LL failure prediction results

##  6 Conclusion

In this paper, we presented methods that can be used as a proxy for failure
probability in LLM responses to QA tasks without domain-specific or access to
the LLM. The three measures introduced were closely correlated with failure
probability across three datasets and we demonstrated several applications
including error prediction, prompt engineering, and error correction. These
applications demonstrated how, often without the requirement for supervised
data, we could obtain improved results over traditional prompting methods.
Further integration with other techniques such as the recently introduced
“tree of thought” prompting (Yao et al. 2023) is a direction for future work.
We also believe the application to automatic code generation (Vaithilingam,
Zhang, and Glassman ), and the creation of logical sentences (Liu et al. 2022)
for reasoning tasks (e.g., the rules used in (Aditya et al. 2023)) is an
important application due to the precision required in the output of an LLM
when translating natural language to computer code and/or formal logic.
Additionally, diversity measures may also have applications to creative use
cases where there is no “correct” answer. In some of these use cases (e.g.,
for brainstorming in an artistic project) an increase in diversity of response
may be desirable.

##  7 Acknowledgments

Some of the authors are supported by internal funding from ASU Fulton Schools
of Engineering.

## References

  * Aditya et al. (2023) Aditya, D.; Mukherji, K.; Balasubramanian, S.; Chaudhary, A.; and Shakarian, P. 2023\.  PyReason: Software for Open World Temporal Logic.  In _AAAI Spring Symposium_.
  * Agrawal, Mackey, and Kalai (2023) Agrawal, A.; Mackey, L.; and Kalai, A. T. 2023.  Do Language Models Know When They’re Hallucinating References?  arXiv:2305.18248.
  * Daftry et al. (2016) Daftry, S.; Zeng, S.; Bagnell, J. A.; and Hebert, M. 2016.  Introspective Perception: Learning to Predict Failures in Vision Systems.
  * Diao et al. (2023) Diao, S.; Wang, P.; Lin, Y.; and Zhang, T. 2023.  Active Prompting with Chain-of-Thought for Large Language Models.  arXiv:2302.12246.
  * Fan, Lewis, and Dauphin (2018) Fan, A.; Lewis, M.; and Dauphin, Y. 2018.  Hierarchical Neural Story Generation.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , 889–898. Association for Computational Linguistics.
  * Ficler and Goldberg (2023) Ficler, J.; and Goldberg, Y. 2023.  Controlling Linguistic Style Aspects in Neural Language Generation.  In _Proceedings of the Workshop on pages = 94–104, Stylistic Variation_. Association for Computational Linguistics.
  * Frieder et al. (2023) Frieder, S.; Pinchetti, L.; Griffiths, R.-R.; Salvatori, T.; Lukasiewicz, T.; Petersen, P. C.; Chevalier, A.; and Berner, J. 2023.  Mathematical Capabilities of ChatGPT.
  * Holtzman et al. (2019) Holtzman, A.; Buys, J.; Du, L.; Forbes, M.; and Choi, Y. 2019.  The Curious Case of Neural Text Degeneration.
  * Hospedales et al. (2022) Hospedales, T.; Antoniou, A.; Micaelli, P.; and Storkey, A. 2022.  Meta-Learning in Neural Networks: A Survey.  44(9): 5149–5169.  Publisher: IEEE Computer Society.
  * Kocoń et al. (2023) Kocoń, J.; Cichecki, I.; Kaszyca, O.; Kochanek, M.; Szydło, D.; Baran, J.; Bielaniewicz, J.; Gruza, M.; Janz, A.; Kanclerz, K.; Kocoń, A.; Koptyra, B.; Mieleszczenko-Kowszewicz, W.; Miłkowski, P.; Oleksy, M.; Piasecki, M.; Radliński, .; Wojtasik, K.; Woźniak, S.; and Kazienko, P. 2023.  ChatGPT: Jack of all trades, master of none.
  * Kuhn, Gal, and Farquhar (2023) Kuhn, L.; Gal, Y.; and Farquhar, S. 2023.  Semantic Uncertainty: Linguistic Invariances for Uncertainty Estimation in Natural Language Generation.  In _The Eleventh International Conference on Learning Representations_.
  * Lin, Trivedi, and Sun (2023) Lin, Z.; Trivedi, S.; and Sun, J. 2023.  Generating with Confidence: Uncertainty Quantification for Black-box Large Language Models.  arXiv:2305.19187.
  * Liu et al. (2023) Liu, H.; Ning, R.; Teng, Z.; Liu, J.; Zhou, Q.; and Zhang, Y. 2023.  Evaluating the Logical Reasoning Ability of ChatGPT and GPT-4.
  * Liu et al. (2022) Liu, J. X.; Yang, Z.; Schornstein, B.; Liang, S.; Idrees, I.; Tellex, S.; and Shah, A. 2022.  Lang2LTL: Translating Natural Language Commands to Temporal Specification with Large Language Models.  In _Workshop on Language and Robotics at CoRL 2022_.
  * OpenAI (2023) OpenAI. 2023.  GPT-4 Technical Report.
  * Pardos and Bhandari (2023) Pardos, Z. A.; and Bhandari, S. 2023.  Learning gain differences between ChatGPT and human tutor generated algebra hints.
  * Ramanagopal et al. (2018) Ramanagopal, M. S.; Anderson, C.; Vasudevan, R.; and Johnson-Roberson, M. 2018.  Failing to Learn: Autonomously Identifying Perception Failures for Self-driving Cars.  3(4): 3860–3867.
  * Ray (2023) Ray, P. P. 2023.  ChatGPT: A comprehensive review on background, applications, key challenges, bias, ethics, limitations and future scope.  3: 121–154.
  * Reimers and Gurevych (2019) Reimers, N.; and Gurevych, I. 2019.  Sentence-BERT: Sentence Embeddings using Siamese BERT-Networks.  In _Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP)_ , 3982–3992. Hong Kong, China: Association for Computational Linguistics.
  * Saha et al. (2018) Saha, A.; Pahuja, V.; Khapra, M. M.; Sankaranarayanan, K.; and Chandar, S. 2018\.  Complex Sequential Question Answering: Towards Learning to Converse over Linked Question Answer Pairs with a Knowledge Graph.  In _Proceedings of the Thirty-Second AAAI Conference on Artificial Intelligence and Thirtieth Innovative Applications of Artificial Intelligence Conference and Eighth AAAI Symposium on Educational Advances in Artificial Intelligence_ , AAAI’18/IAAI’18/EAAI’18. AAAI Press.  ISBN 978-1-57735-800-8.
  * Shakarian et al. (2023) Shakarian, P.; Koyyalamudi, A.; Ngu, N.; and Mareedu, L. 2023.  An Independent Evaluation of ChatGPT on Mathematical Word Problems (MWP).
  * Tian et al. (2023) Tian, K.; Mitchell, E.; Zhou, A.; Sharma, A.; Rafailov, R.; Yao, H.; Finn, C.; and Manning, C. D. 2023.  Just Ask for Calibration: Strategies for Eliciting Calibrated Confidence Scores from Language Models Fine-Tuned with Human Feedback.  arXiv:2305.14975.
  * Tu et al. (2023) Tu, S.; Li, C.; Yu, J.; Wang, X.; Hou, L.; and Li, J. 2023.  ChatLog: Recording and Analyzing ChatGPT Across Time.
  * Upadhyay and Chang (2017) Upadhyay, S.; and Chang, M.-W. 2017.  Annotating Derivations: A New Evaluation Strategy and Dataset for Algebra Word Problems.
  * (25) Vaithilingam, P.; Zhang, T.; and Glassman, E. L. ????  Expectation vs. Experience: Evaluating the Usability of Code Generation Tools Powered by Large Language Models.  _CHI Conference on Human Factors in Computing Systems Extended Abstracts (CHI ’22 Extended Abstracts)_.
  * Wang et al. (2022) Wang, X.; Wei, J.; Schuurmans, D.; Le, Q.; Chi, E.; Narang, S.; Chowdhery, A.; and Zhou, D. 2022.  Self-Consistency Improves Chain of Thought Reasoning in Language Models.
  * Wei et al. (2022) Wei, J.; Wang, X.; Schuurmans, D.; Bosma, M.; Ichter, B.; Xia, F.; Chi, E.; Le, Q.; and Zhou, D. 2022.  Chain-of-Thought Prompting Elicits Reasoning in Large Language Models.
  * Yao et al. (2023) Yao, S.; Yu, D.; Zhao, J.; Shafran, I.; Griffiths, T. L.; Cao, Y.; and Narasimhan, K. 2023.  Tree of Thoughts: Deliberate Problem Solving with Large Language Models.
  * Yuan et al. (2023) Yuan, Z.; Yuan, H.; Tan, C.; Wang, W.; and Huang, S. 2023.  How well do Large Language Models perform in Arithmetic tasks?
  * Zhou, Jurafsky, and Hashimoto (2023) Zhou, K.; Jurafsky, D.; and Hashimoto, T. 2023.  Navigating the Grey Area: Expressions of Overconfidence and Uncertainty in Language Models.  arXiv:2302.13439.
  * Zhou et al. (2022) Zhou, K.; Liu, Z.; Qiao, Y.; Xiang, T.; and Loy, C. C. 2022.  Domain Generalization: A Survey.  1–20.  Conference Name: IEEE Transactions on Pattern Analysis and Machine Intelligence.