# On the Role of Tikhonov Regularizations in Standard Optimization Problems

J. Adriazola  Department of Mathematics, University of California, Santa
Barbara, CA, USA, 93106

###### Abstract.

Tikhonov regularization is a common technique used when solving poorly behaved
optimization problems. Often, and with good reason, this technique is applied
by practitioners in an ad hoc fashion. In this note, we systematically
illustrate the role of Tikhonov regularizations in two simple, yet instructive
examples. In one example, we use regular perturbation theory to predict the
impact Tikhonov regularizations have on condition numbers of symmetric,
positive semi-definite matrices. We then use a numerical example to confirm
our result. In another example, we construct an exactly solvable optimal
control problem that exhibits a boundary layer phenomena. Since optimal
control problems are rarely exactly solvable, this brings clarity to how vital
Tikhonov regularizations are for the class of problems this example
represents. We solve the problem numerically using MATLAB’s built-in
optimization software and compare results with a Galerkin-type/genetic
algorithm. We find the second method generally outperforms MATLAB’s
optimization routines and we provide the MATLAB code used to generate the
numerics.

##  1\. Introduction

Ill-posed and ill-conditioned problems frequently arise in applied and
computational mathematics. To overcome this common challenge, one typically
introduces a type of regularization that either converts the problem to a
well-posed one or better conditions the resulting computational problem. The
Tikhonov regularization technique is often chosen as the first strategy toward
this goal and has been successfully applied in several contexts ranging from
ill-posed optimal control problems [9, 11, 13, 2, 1] to the numerical solution
of poorly-conditioned integral equations [21, 18, 19, 20, 15]. In statistics,
the method is better known as ridge regression [10].

To illustrate how Tikhonov regularizations apply in a general optimization
setting, consider the problem of finding $x_{*}$ such that

(1) $J(x_{*})=\min_{x\in X}J(x),$

subject to

(2) $f(x)=0,$

where $J:X\to\mathbb{R}$ is called an objective functional, $f:X\to Y$ is a
constraint function, and $X,\ Y$ are general Banach spaces. Although the
problem may be solvable, it may either be too computationally intensive to
solve or fail to admit solutions which remain feasible to implement in a
scientific or technological application. Thus, one often seeks to define a
modified problem

(3) $J_{\varepsilon}(x)=J(x)+\varepsilon K(x)$

while still remaining faithful to the originally intended goal.

Of the many possible pitfalls that may be encountered when attempting to solve
an ill-posed problem, one important example is where the optimal solution
$x_{*}\in X$ lacks a desired degree of regularity or smoothness. A Tikhonov
regularization enforces smoothness of solutions to this problem by modifying
the objective to include a penalization. A typical choice is setting

(4) $K(x)=\norm*{x}^{2}_{\dot{H}^{1}(\Omega)}$

in Equation (3), where $\Omega\subset\mathbb{C}^{n},\ $ $n\in\mathbb{Z}^{+},$
and $\varepsilon\in\mathbb{R}^{+}.$ The notation $\dot{H}^{1}(\Omega)$
indicates a homogeneous Sobolev space with norm given by
$\left\|Dx\right\|_{L^{2}(\Omega)},$ where $L^{2}(\Omega)$ is the standard
Lebesgue (Hilbert) space [12] and $D$ indicates differentiation with respect
to elements from $\Omega.$ Simply put, the introduction of the regularization
penalizes solutions that rapidly vary over $\Omega.$

The problem of finding an optimal value of $\varepsilon$ is known in
statistics as the so-called bias-variance dilemma [14]. Although an optimal
choice of the Tikhonov parameter $\varepsilon$ is known in certain
circumstances [23], the parameter is typically chosen in an ad hoc manner.
Through the simple and instructive examples studied here, we reveal suitable
regularization parameters in simpler settings which may guide choices in more
complex ones. Of course, choosing suitable regularization parameters in more
complex, real-world applications remains to be done on a case-by-case basis.

This work is organized as follows. Section 2 shows how the Tikhonov technique
better conditions a symmetric, positive semi-definite matrix whose inverse is
desired. There, we quantitatively show the type of impact a small Tikhonov
regularization can have on an extremely ill-conditioned matrix. In Section 3,
we show an example of an optimal control problem admitting a closed-form
solution. Since the problem is exactly solvable, this provides a testbed for
standard optimization methods of which we compare two; an interior-point
method [4] implemented by MATLAB’s fmincon and a method, due to Calarco, et
al. [6, 7], called the Chopped Random Basis (CRAB) method which we also
implement in MATLAB. In Appendix A, we provide MATLAB code reflecting our
implementation of the methods used to solve the optimal control problem
discussed in Section 3.

##  2\. Regularizing an Ill-Conditioned Matrix

The first example we consider is matrix inversion of an ill-conditioned
matrix. Let $n\in\mathbb{Z}^{+}$ and consider the following optimization
problem:

(5)
$\min_{x\in\mathbb{R}^{n}}J=\min_{x\in\mathbb{R}^{n}}\left\\{\frac{1}{2}x^{\mathsf{T}}Ax-x^{\mathsf{T}}b+\frac{\varepsilon}{2}x^{\mathsf{T}}D^{\mathsf{T}}Dx\right\\}$

where $A$ is some real, symmetric, and positive semi-definite $n\times n$
matrix, $b\in\mathbb{R}^{n},$ $D$ is some, as of yet, unspecified $n\times m$
matrix, $m\in\mathbb{Z}^{+},$ and $\varepsilon$ is a small, real number. The
restrictions placed on $A$ are done so its spectrum can be analyzed in detail.

By taking the gradient of the objective $J$ in Equation (5) with respect to
the variable $x$, we see that the optimal vector $x^{*}$ satisfies

(6) $\left(A+\varepsilon D^{\mathsf{T}}D\right)x^{*}=b.$

Note that in the case $\varepsilon=0$, the solution to Equation (6) is given
by $x^{*}=A^{-1}b.$ If $A$ is singular, then the task of finding $x^{*}$ is
ill-posed, and if $A$ is nearly-singular, computing $x^{*}$ is challenging and
often requires substantial computational resources for large $n$.

It is well-known that the condition number of a symmetric, positive semi-
definite matrix $A$, with condition number denoted by $c(A),$ is given by the
ratio of its largest eigenvalue to its smallest eigenvalue [17]. For this
reason, we aim to better understand the effect of the Tikhonov regularizer
$\varepsilon D^{\mathsf{T}}D$ on the spectrum of $A$ for small enough
$\varepsilon$ since this informs us of the condition number of the regularized
matrix $A+\varepsilon D^{\mathsf{T}}D$.

To this end, further assume the eigenvalues of $A$ are simple, and consider
the asymptotic expansions

(7a)
$\displaystyle\lambda(\varepsilon)\sim\sum_{j=0}^{\infty}\lambda_{j}\varepsilon^{j}\qquad\mathrm{as}\
\varepsilon\to 0,$ (7b)
$\displaystyle\nu(\varepsilon)\sim\sum_{j=0}^{\infty}\nu_{j}\varepsilon^{j}\qquad\mathrm{as}\
\varepsilon\to 0,$

where $\nu_{0}$ and $\lambda_{0}$ are defined as the respective
eigenvector/eigenvalue pair of the matrix $A$. To leading order, we have

(8) $A\nu_{0}=\lambda_{0}\nu_{0},$

which is automatically satisfied by the definition of $\nu_{0}$ and
$\lambda_{0}$.

To next order, i.e., $\mathcal{O}(\varepsilon),$ we have

(9) $(A-\lambda_{0}I)\nu_{1}=(\lambda_{1}-D^{\mathsf{T}}D)\nu_{0},$

where $I$ is the $n\times n$ identity matrix. The Fredholm alternative [17]
requires the right-hand side to be orthogonal to the left-eigenvector of $A$
so that Equation (9) is solvable. Since $A$ is symmetric, the left eigenvector
of $A$ is simply $\nu_{0}^{\mathsf{T}}.$ The solvability condition of Equation
(9) is then given by

(10)
$\lambda_{1}=\frac{\nu_{0}^{\mathsf{T}}D^{\mathsf{T}}D\nu_{0}}{\nu_{0}^{\mathsf{T}}\nu_{0}}.$

Thus, the resulting condition number is

(11) $c(A+\varepsilon
D^{\mathsf{T}}D)=\frac{\lambda^{\mathrm{max}}}{\lambda^{\mathrm{min}}}\sim\frac{\lambda_{0}^{\mathrm{max}}+\varepsilon\lambda_{1}^{\mathrm{max}}}{\lambda_{0}^{\mathrm{min}}+\varepsilon\lambda_{1}^{\mathrm{min}}}+\mathcal{O}(\varepsilon^{2}),\qquad\mathrm{as}\
\varepsilon\to 0,$

where $\lambda_{0}^{\rm min/max}$ are, respectively, the smallest and largest
eigenvalues satisfying Equation (8) while $\lambda_{1}^{\rm min/max}$ are the
smallest and largest numbers satisfying Equation (10). Presumably,
$\lambda_{0}^{\mathrm{min}}$ is nearly zero, so long as
$\lambda_{1}^{\mathrm{min}}$ is not, the regularization has a high likelihood
of being effective.

Indeed, the asymptotic formula (11) gives a criteria for finding effective
$A-$specific Tikhonov matrices, i.e., the Tikhonov matrix $D$ solves the
following max-min problem

(12)
$\max_{D\in\mathcal{D}}\min_{\nu_{0}\in\mathcal{V}}\left\\{\frac{\nu_{0}^{\mathsf{T}}D^{\mathsf{T}}D\nu_{0}}{\nu_{0}^{\mathsf{T}}\nu_{0}}\right\\},$

where $\mathcal{D}$ is the space of all $n\times m$ matrices and $\mathcal{V}$
is the set of eigenvectors of the ill-conditioned matrix $A.$ We do not take
up a study of Equation (12) here, and instead leave further investigation for
future study.

To better illustrate the effect a Tikhonov regularization has on a nearly
singular operator, take as an example the symmetric, positive-definite matrix

(13) $A=\begin{pmatrix}1&1\\\ 1&1+\mu\end{pmatrix},\qquad\mu>0$

which has eigenvalues

(14) $\lambda_{0}^{\rm
min}=\frac{1}{2}\left(\mu+2-\sqrt{\mu^{2}+4}\right),\quad\lambda_{0}^{\rm
max}=\frac{1}{2}\left(\mu+2+\sqrt{\mu^{2}+4}\right).$

and corresponding eigenvectors

(15) $\nu_{0}^{\rm
min}=\begin{pmatrix}\frac{1}{2}\left(-\mu-\sqrt{\mu^{2}+4}\right)\\\
1\end{pmatrix},\quad\nu_{0}^{\rm
max}=\begin{pmatrix}\frac{1}{2}\left(-\mu+\sqrt{\mu^{2}+4}\right)\\\
1\end{pmatrix},$

Clearly, $A$ is ill-conditioned since $c(A)\to\infty$ as $\mu\to 0$ from
above. Indeed, for $\mu=10^{-6},$ the condition number $c(A)=4\times 10^{6}$
to seven digits of precision. However, the asymptotic result (11) predicts
that simply choosing $D$ to be the forward difference operator

(16) $D=\begin{pmatrix}1&0\\\ -1&1\\\ 0&-1\end{pmatrix},$

a choice consistent with Equation (4) in this finite-dimensional setting, and
a small regularization parameter of $\varepsilon=0.01,$ the regularized
condition number $c(A+\varepsilon D^{\mathsf{T}}D)=203$ to three digits of
precision; a noteworthy improvement. In fact, the asymptotic result is well
within the same order of magnitude of the true value of the condition number

(17) $c(A+\varepsilon
D^{\mathsf{T}}D)=\frac{2+4\varepsilon+\mu+\sqrt{4-8\varepsilon+4\varepsilon^{2}+\mu^{2}}}{2+4\varepsilon+\mu-\sqrt{4-8\varepsilon+4\varepsilon^{2}+\mu^{2}}}$

which is 67 to two digits of precision.

Rounding out the example, we choose an example vector
$b=\begin{pmatrix}\frac{1}{2},\ \frac{1}{2}\end{pmatrix}^{\mathsf{T}}$ to
demonstrate that we can still solve the matrix inversion in a meaningful way.
When $\varepsilon=0$, the exact solution is given by $x_{\rm
unreg}=\begin{pmatrix}\frac{1}{2},\ 0\end{pmatrix}^{\mathsf{T}},$ while, to
three digits of precision and when $\varepsilon=0.01,$ the solution $x_{\rm
reg}=\begin{pmatrix}0.248,\ 0.248\end{pmatrix}^{\mathsf{T}}.$ The regularized
solution has a residual, i.e., $J(x_{\rm reg})$ when $\varepsilon=0,$ of
$1.24\times 10^{-5}.$ The vast improvement on the condition number of $A,$
while still remaining faithful to the unregularized numerical computation,
clearly demonstrates that even a small Tikhonov regularization can greatly
mitigate an extremely ill-conditioned numerical computation.

##  3\. Regularizing an Exactly Solvable Optimal Control Problem

The second example we consider is an optimal control problem representative of
optimization problems that frequently arise in controlling dispersive waves,
e.g. [2, 1, 11]. This optimal control problem, by virtue of its exact
solution, helps reveal some of the typical numerical difficulties exhibited by
this class of problems. Consider the following Lagrange control problem

(18)
$\min_{u\in\mathcal{U}}J=\min_{u\in\mathcal{U}}\int_{0}^{\pi}\left(x+\frac{1}{2}u^{2}+\frac{\varepsilon}{2}\left(\frac{\mathrm{d}u}{\mathrm{d}t}\right)^{2}\right)\mathrm{d}t,$

where the admissible class $\mathcal{U}=\left\\{u\in
H^{1}([0,\pi]):u(0)=u(\pi)=0\right\\},$ subject to the forced harmonic
oscillator

(19a) $\displaystyle\frac{\mathrm{d}x}{\mathrm{d}t}$ $\displaystyle=p,\qquad$
$\displaystyle x(0)=0,$ (19b) $\displaystyle\frac{\mathrm{d}p}{\mathrm{d}t}$
$\displaystyle=u-x,\qquad$ $\displaystyle p(0)=0.$

The prescription of boundary conditions in defining the class of admissible
controls $\mathcal{U}$ will soon play a critical role when considering the
unregularized case, i.e., $\varepsilon=0$. For ease of notation in what
follows, we use overhead dots to denote differentiation with respect to time,
e.g., $\dot{x}=\frac{\mathrm{d}x}{\mathrm{d}t}.$

We recast the problem into a form more amenable to standard tools from the
calculus of variations. Using Lagrange multipliers, or costates in the
language of optimal control theory and denoted here by $\lambda(t)$ and
$\mu(t)$, the equivalent unconstrained saddle-point problem is given by

(20)
$\displaystyle\begin{split}&\min_{u\in\mathcal{U}}\int_{0}^{\pi}\mathcal{L}\left(x,\dot{x},p,\dot{p},\lambda,\mu,u,\dot{u}\right)\mathrm{d}t\\\
=&\min_{u\in\mathcal{U}}\int_{0}^{\pi}\left(x+\frac{1}{2}u^{2}+\frac{\varepsilon}{2}\dot{u}^{2}+\lambda(\dot{x}-p)+\mu(\dot{p}+x-u)\right)\mathrm{d}t.\end{split}$

The Euler-Lagrange equations (see, for example, [8, 24, 5]) are given by

(21a) $\displaystyle\frac{\delta\mathcal{L}}{\delta x}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial{x}}-\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial{\dot{x}}}=1-\dot{\lambda}+\mu=0,\qquad$
$\displaystyle\frac{\partial\mathcal{L}}{\partial{\dot{x}}}\bigg{|}_{t=\pi}=\lambda(\pi)=0,$
(21b) $\displaystyle\frac{\delta\mathcal{L}}{\delta p}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial{p}}-\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial{\dot{p}}}=\lambda-\dot{\mu}=0,\qquad$
$\displaystyle\frac{\partial\mathcal{L}}{\partial{\dot{p}}}\bigg{|}_{t=\pi}=\mu(\pi)=0,$
(21c) $\displaystyle\frac{\delta\mathcal{L}}{\delta\lambda}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial{\lambda}}=\dot{x}-p=0,\qquad$
$\displaystyle x(0)=0,$ (21d)
$\displaystyle\frac{\delta\mathcal{L}}{\delta\mu}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial{\mu}}=\dot{p}+x-u=0,\qquad$
$\displaystyle p(0)=0,$ (21e) $\displaystyle\frac{\delta\mathcal{L}}{\delta
u}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial{u}}-\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial{\dot{u}}}=u-\mu-\varepsilon\ddot{u}=0,\qquad$
$\displaystyle u(0)=u(\pi)=0.$

Note that the costates’ terminal conditions, present in Equations (21a) and
(21b), are a result of an integration by parts again test functions
$\varphi\in C^{\infty}([0,T]),$ which do not necessarily vanish at $t=\pi,$ in
the derivation of the Euler-Lagrange equations.

By direct integration, it’s easy to see the following functions

(22a) $\displaystyle\lambda_{*}$ $\displaystyle=-\sin(t),$ (22b)
$\displaystyle\mu_{*}$ $\displaystyle=-\cos(t)-1$

are the optimal costates. This implies, after using the method of variation of
parameters, see e.g. [3], to solve Equation (21e), the optimal control is

(23)
$u_{*}=\frac{\text{csch}\left(\frac{\pi}{\sqrt{\varepsilon}}\right)\left((\varepsilon+2)\sinh\left(\frac{\pi-t}{\sqrt{\varepsilon}}\right)+\varepsilon\sinh\left(\frac{t}{\sqrt{\varepsilon}}\right)\right)-\cos(t)-\varepsilon-1}{\varepsilon+1}.$

Additionally, we solve for the optimal state $x_{*}$ in a similar fashion
which results in

(24)
$x_{*}=\frac{\sqrt{\varepsilon}\text{csch}\left(\frac{\pi}{\varepsilon}\right)}{(\varepsilon+1)^{2}}\left(\varepsilon^{3/2}\sinh\left(\frac{t}{\sqrt{\varepsilon}}\right)-\varepsilon\sin(t)\right)-\frac{t\sin(t)}{2(\varepsilon+1)}+\cos(t)-1\\\
+\frac{\sqrt{\varepsilon}(\varepsilon+2)\text{csch}\left(\frac{\pi}{\varepsilon}\right)}{(\varepsilon+1)^{2}}\left(\sqrt{\varepsilon}\left(\sinh\left(\frac{\pi-t}{\sqrt{\varepsilon}}\right)-\sinh\left(\frac{\pi}{\sqrt{\varepsilon}}\right)\cos(t)\right)+\cosh\left(\frac{\pi}{\sqrt{\varepsilon}}\right)\sin(t)\right),$

while the optimal objective $J^{*}$ is found to be

(25)
$J_{*}=\sqrt{\varepsilon}\tanh\left(\frac{\pi}{2\sqrt{\varepsilon}}\right)\left(\frac{\coth^{2}\left(\frac{\pi}{2\sqrt{\varepsilon}}\right)}{(\varepsilon+1)^{2}}+1\right)-\frac{\pi(2\varepsilon+3)}{4(\varepsilon+1)}.$

With these calculations in hand and further considering any time $0<t\leq T$,
we may take these optimal functions’ limits as $\varepsilon$ approaches zero
from above. We find that

(26a) $\displaystyle\lim_{\varepsilon\to 0^{+}}u_{*}$ $\displaystyle:=u_{\rm
limit}=-1-\cos(t)$ (26b) $\displaystyle\lim_{\varepsilon\to 0^{+}}x_{*}$
$\displaystyle:=x_{\rm limit}=-\frac{1}{2}t\sin(t)+\cos(t)-1\,$ (26c)
$\displaystyle\lim_{\varepsilon\to 0^{+}}J_{*}$ $\displaystyle:=J_{\rm
limit}=-\frac{3\pi}{4}.$

Although the unregularized problem exhibits no features of an ill-posed
problem, note that the limiting optimal control fails to remain in the
admissible class $\mathcal{U}$ since $u_{\rm limit}$ does not vanish at $t=0$.
However, from the explicit calculation, we have that for each $\varepsilon>0,$
there exists a $u\in\mathcal{U}$ which solves the optimal control problem
(18). Thus, taking small values of the regularization parameter $\varepsilon$
manifests in a type of boundary layer in the vicinity of $t=0$ for optimal
control $u_{*}$. That is, the optimal control develops a steep gradient at
$t=0$ as $\varepsilon\to 0^{+},$ and, consequently, $\lim_{\varepsilon\to
0^{+}}\norm*{u_{*}-u_{\rm limit}}_{H^{1}([0,T])}\neq 0.$ The effect this has
on numerical computations will be discussed shortly.

We numerically solve optimal control problem (18) in two ways. The first,
performed pointwise in time, is to use MATLAB’s fmincon which implements an
interior-point method [4]. The second is to use a method due to Calarco, et
al. [6, 7], called the $\mathbf{C}$hopped $\mathbf{Ra}$ndom $\mathbf{B}$asis
(CRAB) method. It relies on choosing controls from the span of an
appropriately chosen finite set of basis functions so that the optimization is
performed over a set of unknown coefficients. The choice of basis is problem-
specific but is always chosen so that controls remain in the appropriate
admissible space. The choice of representation we make here is the following
sine series:

(27)
$u_{r}(t)=\sum_{j=1}^{N}\frac{r_{j}}{j^{2}}\sin\left(\frac{2j\pi}{T}\right),\quad
t\in[0,T],\\\ $

where the coefficients $r_{j}\in[-1,1]$ are parameters to be optimized over.
Note that we choose the amplitudes in Equation (27) to decay quadratically
because the Fourier series of absolutely continuous functions exhibits this
type of decay [22]. In this way, the search space of amplitudes is not
severely restricted, yet early iterations of the controls generated by the
CRAB method are not highly oscillatory.

Figure 1. A numerical solution of optimal control problem (18), with
$\varepsilon$=1, via the CRAB method and MATLAB’s fmincon. Left Panel: The
controls found via both numerical methods and the exact optimal control given
in closed-form by Equation (23). Right Panel: The states found via both
numerical methods and the exact state given in closed-form by Equation (24).

The CRAB method is of Galerkin-type. Recall that the intent behind any
Galerkin method is to choose the number of basis functions $N$ simultaneously
large enough to define an accurate approximation, yet small enough so that the
overall procedure remains computationally inexpensive. For this problem, we
find that just 12 modes are sufficient. The optimization problem is now a
small-scale nonlinear programming (NLP) problem and can be solved using any
number of industry-standard techniques. The technique we use to solve the
resulting NLP problem is Differential Evolution (DE) [16].

In our numerical study, we focus on investigating the role of the Tikhonov
parameter $\varepsilon.$ To this end, we set $\varepsilon=1$ and show the
numerical solution of Problem (18) in Figure 1. We solve the state equations
(19) with MATLAB’s ODE45 solver. The numerics shown there demonstrate the
adequacy of both methods used to solve Problem (18). The numerical
implementation of these methods is provided in the form of MATLAB code in
Appendix A.

By setting $\varepsilon=1,$ we find that we have strongly regularized the
problem. This is seen in a pointwise comparison between $x_{*},$ shown in
Figure 1, and $x_{\rm limit}$, shown in Figure 2, especially for later times
$t$. Yet, by taking $\varepsilon=0.001,$ we find that both optimization
methods struggle to resolve the sharp transition of the control in the
boundary layer near $t=0.$ We find $\varepsilon=0.04$ to be a good
intermediate value under the tension of under- and over-regularizing the
control problem (the bias-variance dilemma). In this case, we see that the
CRAB method successfully finds the optimal control $u_{*},$ while fmincon
develops undesirable oscillations in the control near the boundary layer.

Figure 2. Examples of when the Tikhonov parameter $\varepsilon$ is small. The
conventions here are consistent with Figure 1. The top panels correspond to
$\varepsilon=0.001$ while the bottom panels have $\varepsilon=0.04$. The
bottom panels also include the limiting control and state functions given by
Equation (26).

##  4\. Concluding Remarks

In this work, we provide a pair of standard optimization problems where an
application of the Tikhonov regularization technique is crucial. Moreover,
they provide instructive examples which can be used in basic courses on
numerical linear algebra and optimization. We also demonstrate the utility of
the CRAB method, over more standard numerical optimization methods, when
trying to solve control problems that exhibit sharp transition regions in the
optimal control.

##  5\. Acknowledgments

The author gratefully acknowledges helpful discussions with T. Witelski, R.H.
Goodman, and M. Porter.

## References

  * [1] J. Adriazola and R. H. Goodman.  Optimal control approach to gradient-index design for beam reshaping.  Journal of the Optical Society of America A, 39(5):907, apr 2022\.
  * [2] J. Adriazola and R. H. Goodman.  Reduction-based strategy for optimal control of bose-einstein condensates.  Physical Review E, 105(2), feb 2022.
  * [3] C. S. Bender and S. Orszag.  Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory.  Springer-Verlag New York Inc., New York, NY, 2010.
  * [4] S. Boyd and L. Vandenberghe.  Convex optimization, volume 1.  Cambridge University Press, Cambridge, 2004.
  * [5] A. E. Bryson and Y.-C. Ho.  Applied optimal control: optimization, estimation, and control.  Hemisphere Publishing Corporation, 1975.
  * [6] T. Caneva, T. Calarco, and S. Montangero.  Chopped random-basis quantum optimization.  Phys. Rev. A, 84:022326, Aug 2011.
  * [7] P. Doria, T. Calarco, and S. Montangero.  Optimal control technique for many-body quantum dynamics.  Phys. Rev. Lett., 106:190501, May 2011.
  * [8] I. Gelfand and S. Fomin.  Calculus of Variations.  Prentice-Hall, Englewood Cliffs, New Jersey, 1963.
  * [9] M. Hintermuller, D. Marahrens, P. A. Markowich, and C. Sparber.  Optimal bilinear control of Gross-Pitaevskii equations.  SIAM (Society for Industrial and Applied Mathematics) Journal on Control and Optimization, 51(3):2509–2543, 2013.
  * [10] A. E. Hoerl and R. W. Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  Soviet Mathematics, 12(1):55–67, 1970.
  * [11] U. Hohenester, P. K. Rekdal, A. Borzi, and J. Schmiedmayer.  Optimal quantum control of Bose-Einstein condensates in magnetic microtraps.  Physical Review A, 75(2):023602, 2007.
  * [12] E. Lieb and M. Loss.  Analysis.  American Mathematical Society (AMS), Philadelphia, PA, second edition, 2010.
  * [13] J. Mennemann, D. Matthes, R. Weishaupl, and T. Langen.  Optimal control of Bose-Einstein condensates in three dimensions.  New Journal of Physics, 17(11):113027, 2015.
  * [14] B. Neal.  On the bias-variance tradeoff: Textbooks need an update.  CoRR, abs/1912.08286, 2019.
  * [15] D. L. Phillips.  A technique for the numerical solution of certain integral equations of the first kind.  Journal of the ACM, 9:84–97, 1963.
  * [16] R. Storn and K. Price.  Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces.  Journal of Global Optimization, 11(4):341–359, 1997.
  * [17] G. Strang.  Introduction to linear algebra, volume 3.  Wellesley-Cambridge Press, Wellesley, MA, 1993.
  * [18] A. N. Tikhonov.  On the stability of inverse problems.  Doklady Akademii Nauk SSSR, 39(5):195–198, 1943.
  * [19] A. N. Tikhonov.  Solution of incorrectly formulated problems and the regularization method.  Soviet Mathematics, 39(5):1035–1038, 1963.
  * [20] A. N. Tikhonov and V. Y. Arsenin.  Solution of Ill-posed Problems.  Winston and Sons, 1977.
  * [21] A. N. Tikhonov, A. V. Goncharsky, V. V. Stepanov, and A. G. Yagola.  Numerical Methods for the Solution of Ill-Posed Problems.  Springer, Dordrecht, 1995.
  * [22] L. N. Trefethen.  Spectral methods in MATLAB.  Society for Industrial and Applied Mathematics (SIAM), New York, NY, 2000\.
  * [23] G. Wahba.  Spline Models for Observational Data.  Society for Industrial and Applied Mathematics, 1990.
  * [24] T. Witelski and M. Bowen.  Methods of Mathematical Modeling.  Spring Undergraduate Mathematics Series, Cham, 2016.

##  Appendix A Main Script

    %% Solves an Optimal Control Problem via CRAB and fmincon
    clear, close all

    % Tikhonov weight g and Optimal control in closed-form
    g=.04;
    uExact=@(x)(-1).*(1+g).^(-1).*(1+g+cos(x)+(-1).*csch(g.^(-1/2).*pi).*((2+g).* ...
      sinh(g.^(-1/2).*(pi+(-1).*x))+g.*sinh(g.^(-1/2).*x)));

    % Initial conditions for the state and control boundary conditions
    x0=[0 0];
    u0=0;uT=0;t0=0;T=pi;

    %Define the objectives for each case and number of optimization points
    JCRAB=@(u)getObjective(u,g,false,x0,t0,T);
    Jfmin=@(u)getObjective(u,g,true,x0,t0,T);

    % Differential Evolution Parameters
    DEParams.F=0.8;DEParams.CR=0.7;
    DEParams.NP=120;DEParams.Nmax=300;
    DEParams.ND=12;

    % Amplitude hyperrectangle for the CRAB method
    l=zeros(1,DEParams.ND);
    for j=1:DEParams.ND
        l(j)=(2*rand-1)/j^2;
    end

    % Execute differential evolution
    optimalamps=diffevoND(JCRAB,l,DEParams);

    % Reconstruct the control from the CRAB amplitudes and compute resulting state
    uCRAB=recon(linspace(t0,T,1000),optimalamps(:,end));

    % Prepare for a blackbox optimization
    options = optimoptions(@fmincon,’Display’,’iter’,’Algorithm’,’interior-point’);

    % Keep number of function evaluations consistent with CRAB
    options.MaxFunctionEvaluations=DEParams.NP*DEParams.Nmax;
    fminPts=15;
    t=linspace(t0,T,fminPts);

    % initial guess on the control
    u00=t.*(t-T);

    % set equality constraints on the control
    Aeq = zeros(fminPts);
    Aeq(1,1)=1;
    Aeq(end,end)=1;
    beq=zeros(fminPts,1);

    % Execute an interior-point method using MATLAB’s fmincon
    ufmin=fmincon(Jfmin,u00,[],[],Aeq,beq,[],[],[],options);
    ufmin=spline(t,ufmin);
    ufmin=@(t)ppval(t,ufmin);

    % Solve for the states
    xCRAB =SolveStates(t0,T,uCRAB,x0);
    xfmin =SolveStates(t0,T,ufmin,x0);
    xExact=SolveStates(t0,T,uExact,x0);

##  Appendix B Functions Used

    function x=SolveStates(t0,T,u,x0)
    % Solve for states using ODE45
    [~,states]=ode45(@(t,x)StateRHS(t,x,u),[t0 T],x0);
    x=states(:,1);x=spline(t,x);x=@(t)ppval(t,x);
    end

    function dxdt=StateRHS(t,x,u)
    dxdt=[x(2)
        u(t)-x(1)];
    end

    function J=getObjective(u,g,blackbox,x0,t0,T)
    if blackbox
        % this is coming from fmincon
        t=linspace(t0,T,length(u));
        dt=t(2)-t(1);
        u=spline(t,u);
        u=@(t)ppval(t,u);
    else
        % this is coming from CRAB
        t=linspace(t0,T,1000);
        dt=t(2)-t(1);
        u=recon(t,u);
    end
    % solve for the states
    x=SolveStates(t0,T,u,x0);
    % compute the derivative of the control
    udot=gradient(u(t),dt);
    % the objective functional
    J=trapz(t,x(t)+.5*u(t).^2+.5*g*udot.^2);
    end

    function u=recon(t,params)
    % reconstruct the control from the CRAB amplitudes
    t0=t(1);T=t(end)-t0;
    sinepart=0;
    for i=1:length(params)
        sinepart=sinepart+params(i)*sin(i*pi*(t-t0)’/T);
    end
    u=sinepart;
    u=spline(t,u);
    u=@(t)ppval(t,u);
    end

    function optimalamps=diffevoND(f,l,DEparams)
    %Minimizes f(x) using the differential evolution algorithm
    %of Storn, Price
    %Inputs:
    %f is the function handle we want to minimize
    %l is the length of the interval for the defined function
    %DEparams.CR,DEparams.NP,DEparams.F are the diffevo
    %parameters discussed in the literature
    %Output: the possible global minimizer
    %Initialize a population of candidate minimizers from the intervals [-l,l]
    pop=zeros(DEparams.ND,DEparams.NP);
    for j=1:DEparams.ND
        for i=1:DEparams.NP
            pop(j,i)=(l(j)*(2*rand-1));
        end
    end
    fpop=zeros(DEparams.NP,1);
    z=zeros(DEparams.ND,1);
    for i=1:DEparams.NP
        fpop(i)=f(pop(:,i));
    end

    %Set DEparams.tolerance and maximum number of iterations
    counter=0;tic
    while counter<DEparams.Nmax
        newpop=pop;
        newfpop=fpop;
        for i=1:DEparams.NP
            for j=1:DEparams.ND
                currentmember=pop(j,i);
                %Choose 3 distinct elements different from current member
                randind=randsample(DEparams.NP,3);
                while length(unique(randind))<3||~isempty(find(randind==i,1))
                    randind=randsample(DEparams.NP,3);
                end
                a=pop(j,randind(1));
                b=pop(j,randind(2));
                c=pop(j,randind(3));

                %Compute currentmember’s potentially new location y
                r=rand;
                %Potential Cross Over
                if r<DEparams.CR
                    y=a+DEparams.F*(b-c);
                else
                    y=currentmember;
                end
                z(j)=y;
            end
            fCandidate=f(z);
            if fCandidate<fpop(i)
                newpop(:,i)=z;
                newfpop(i)=fCandidate;
            end
        end
        %Update
        counter=counter+1;
        pop=newpop;
        fpop=newfpop;
        fprintf(’%f cost after %i iterations of DE complete in %f seconds.\n’,...
            min(fpop),counter,toc)
        [~,optindex]=min(fpop);
        optimalamps(:,counter)=pop(:,optindex);
    end
    end


# Mixing time of fractional random walk on finite fields

Jimmy He  Department of Mathematics, Stanford University, Stanford, CA, USA
<EMAIL_ADDRESS>,  Huy Tuan Pham  Department of Mathematics, Stanford
University, Stanford, CA, USA<EMAIL_ADDRESS>and  Max Wenqiang Xu
Department of Mathematics, Stanford University, Stanford, CA, USA
<EMAIL_ADDRESS>
###### Abstract.

We study a random walk on $\mathbb{F}_{p}$ defined by
$X_{n+1}=1/X_{n}+\varepsilon_{n+1}$ if $X_{n}\neq 0$, and
$X_{n+1}=\varepsilon_{n+1}$ if $X_{n}=0$, where $\varepsilon_{n+1}$ are
independent and identically distributed. This can be seen as a non-linear
analogue of the Chung–Diaconis–Graham process. We show that the mixing time is
of order $\log p$, answering a question of Chatterjee and Diaconis [12].

##  1\. Introduction

This paper studies a non-linear analogue of the Chung–Diaconis–Graham process,
defined on $\mathbb{F}_{p}$ by

$X_{n+1}=aX_{n}+\varepsilon_{n+1},$

where $a\in\mathbb{F}_{p}^{\times}$ is fixed and the $\varepsilon_{i}$ are
independent and identically distributed. The mixing time of this Markov chain
has been extensively studied, and it is now known that for certain $a$ (such
as $a=2$) and almost all $p$, cutoff occurs at time $c\log p$ for an explicit
constant $c$ [18].

Since simple random walk on $\mathbb{F}_{p}$ requires order $p^{2}$ steps to
mix (see [28] for example), the Chung–Diaconis–Graham process provides an
explicit example of a random walk where applying a deterministic bijection (in
this case, $x\mapsto ax$) between steps of the walk exponentially speeds up
mixing. This was studied in [12], where it was asked whether other explicit
examples could be provided.

In this paper, we consider the random walk on $\mathbb{F}_{p}$ defined by

$X_{n+1}=\iota(X_{n})+\varepsilon_{n+1},$

where $\iota(x)=1/x$ if $x\neq 0$ and $\iota(0)=0$, and the $\varepsilon_{i}$
are independent and identically distributed. Chatterjee and Diaconis asked
about the order of the mixing time for this walk [12], which was originally
suggested by Soundararajan. We solve this problem by showing that this random
walk mixes in order $\log p$ steps. Since simple random walk on
$\mathbb{F}_{p}$ mixes in order $p^{2}$ steps, this gives an exponential
speedup for the mixing time, and provides an explicit example of when adding a
deterministic function between steps of a Markov chain gives an exponential
speedup.

This model is similar in spirit to random walk on $\mathbb{F}_{p}$ given by
moving from $x$ to one of $x+1$, $x-1$ or $x^{-1}$, which is known to be an
expander (see Theorem 8.8 of [36] for example). Our results are not directly
implied by this, but similar ideas underlie the results for both models. We
discuss the relationship with this model in more detail in Remark 1.3.

The proof uses comparison theory developed by Smith [34] to relate the random
walk on $\mathbb{F}_{p}$ to a random walk on the projective line
$\mathbb{P}^{1}(\mathbb{F}_{p})$. The walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$
is the quotient of a random walk on $\operatorname{SL}_{2}(\mathbb{F}_{p})$
which is known to have a constant order spectral gap by results of Bourgain
and Gamburd [6]. While our methods identify the order of the mixing time, they
are not strong enough to identify the constant or obtain cutoff, and we leave
this as an open problem.

Our result also implies bounds on the number of solutions to the congruence
$xy\equiv 1\pmod{p}$ for $x\in I$ and $y\in J$ intervals of the same length,
although these bounds are weaker than the ones established in [14].

###  1.1. Main results

We now formally state our main results. For two probability measures $\mu$ and
$\nu$, let

$\|\mu-\nu\|_{TV}=\sup_{A}|\mu(A)-\nu(A)|$

denote the _total variation distance_. Let $X_{n}$ be an ergodic Markov chain
with transition matrix $P$ and stationary distribution $\pi$, . Let
$P^{n}(x,y)$ denote the probability that $X_{n}=y$ given that $X_{0}=x$. Then
the _mixing time_ is defined to be

$t_{mix}(\varepsilon)=\inf\left\\{n\colon\sup_{x}\|P^{n}(x,\cdot)-\pi\|_{TV}\leq\varepsilon\right\\}.$

Let $\iota:\mathbb{F}_{p}\to\mathbb{F}_{p}$ be the function defined by
$\iota(x)=1/x$ if $x\neq 0$ and $\iota(0)=0$. Let $\varepsilon_{i}$ be
independent and identically distributed on $\mathbb{Z}$. Let $K$ denote the
transition matrix for the Markov chain on $\mathbb{F}_{p}$ defined by

(1.1) $X_{n+1}=\iota(X_{n})+\varepsilon_{n+1}.$

We show that for this Markov chain, $t_{mix}(\varepsilon)$ is $\Theta(\log
p)$, where the implied constants are allowed to depend on both $\mu$ and
$\varepsilon$. Informally, this means that order $\log p$ steps are necessary
and sufficient for the Markov chain to converge to its stationary
distribution. The following theorems give the desired lower and upper bounds
on $t_{mix}$ respectively .

######  Theorem 1.1 (Lower bound).

Let $\mu$ denote a probability measure on $\mathbb{Z}$, and assume that

$H(\mu)=\sum-\mu(x)\log\mu(x)<\infty.$

Let $K$ be the transition matrix for the Markov chain on $\mathbb{F}_{p}$
defined by

$X_{n+1}=\iota(X_{n})+\varepsilon_{n+1},$

where the $\varepsilon_{n}$ are independent and identically distributed
according to $\mu$. Let $\pi$ denote the uniform distribution on
$\mathbb{F}_{p}$. Then for all $x$ in $\mathbb{F}_{p}$

$\|K^{n}(x,\cdot)-\pi\|_{TV}\geq 1-\frac{nH(\mu)+\log 2}{\log p}.$

######  Theorem 1.2 (Upper bound).

Let $\mu$ denote a probability measure on $\mathbb{Z}$ whose support contains
at least two values. Let $K$ be the transition matrix for the Markov chain on
$\mathbb{F}_{p}$ defined by

$X_{n+1}=\iota(X_{n})+\varepsilon_{n+1},$

where the $\varepsilon_{n}$ are independent and identically distributed
according to $\mu$. Let $\pi$ denote the uniform distribution on
$\mathbb{F}_{p}$. Then there exists a constant $C>0$ depending only on $\mu$
such that for all $x$ in $\mathbb{F}_{p}$ and sufficiently large $p$,

$\|K^{n}(x,\cdot)-\pi\|_{TV}\leq\frac{\sqrt{p}}{2}e^{-Cn}.$

The lower bound follows easily from entropy considerations. We make an
assumption that the entropy of $\mu$ is finite, which cannot be removed
entirely (see Remark 2.2). However, this assumption is relatively weak, and
indeed all finitely supported distributions have finite entropy.

We now describe our plan to prove the upper bound. We establish the upper
bound by connecting the Markov chain $X_{n}$ with a Markov chain on
$\mathbb{P}^{1}(\mathbb{F}_{p})$, which in turn is a projection of a Markov
chain on a Cayley graph of $\operatorname{SL}_{2}(\mathbb{F}_{p})$. We next
sketch this connection.

Recall that given a group $G$ and a generating set $S$, we define the _Cayley
graph_ $\mathcal{G}(G,S)$ as the graph with vertex set $G$ and two elements
$x$ and $y$ are connected if and only if $x=\sigma y$ for some $\sigma\in S$.
We use $A(\mathcal{G})$ to denote the normalised adjacency matrix of
$\mathcal{G}$, whose largest eigenvalue is 1. Let $S$ be a subset of
$\operatorname{SL}_{2}(\mathbb{Z})$ so that $S_{p}$ generates
$\operatorname{SL}_{2}(\mathbb{F}_{p})$, where $S_{p}$ is the set $S$ mod $p$.
Results of Bourgain and Gamburd [6] give conditions for when
$\mathcal{G}(\operatorname{SL}_{2}(\mathbb{F}_{p}),S_{p})$ has a constant
spectral gap, i.e.,
$\lambda_{2}(A(\mathcal{G}(\operatorname{SL}_{2}(\mathbb{F}_{p}),S_{p})))\leq
1-c$ for some constant $c>0$ independent of $p$, and thus also a $O(\log p)$
mixing time.

The random walk on $\mathcal{G}(\operatorname{SL}_{2}(\mathbb{F}_{p}),S_{p})$
has a natural projection on $\mathbb{P}^{1}(\mathbb{F}_{p})$, defined by the
action

$\begin{pmatrix}a&b\\\ c&d\\\ \end{pmatrix}\cdot x=\frac{ax+b}{cx+d}.$

This can also be seen as a random walk on the Schreier graph corresponding to
this action. Since the spectral gaps of the walks on
$\operatorname{SL}_{2}(\mathbb{F}_{p})$ and $\mathbb{P}^{1}(\mathbb{F}_{p})$
are related, this gives an $O(\log p)$ mixing time for the walk on
$\mathbb{P}^{1}(\mathbb{F}_{p}).$

Finally, we connect the Markov chain on $\mathbb{P}^{1}(\mathbb{F}_{p})$ with
the fractional Markov chain $X_{n}$ on $\mathbb{F}_{p}$. With an appropriate
choice of the set $S$, the Markov chain on $\mathbb{P}^{1}(\mathbb{F}_{p})$
induced by the random walk on
$\mathcal{G}(\operatorname{SL}_{2}(\mathbb{F}_{p}),S_{p})$ and $X_{n}$ have
very similar transition probabilities. The mixing time bound of $X_{n}$ can
thus be obtained from that of the chain on $\mathbb{P}^{1}(\mathbb{F}_{p})$
via comparison theory of Markov chains on different state spaces.

A key input to the above plan is the spectral gap of random walks on Cayley
graphs of $\operatorname{SL}_{2}(\mathbb{F}_{p})$, which we discuss in further
detail in Section 1.4.

######  Remark 1.3.

The random walk on $\mathbb{F}_{p}$ obtained by moving from $x$ to $x+1$,
$x-1$ or $x^{-1}$ with equal probability has been previously studied. In
Section 3.3 of [31], the walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$ is shown to
have a constant order spectral gap, and Theorem 8.8 of [36] states that the
same is true on $\mathbb{F}_{p}$. This model can be seen as an explicit
version of the speedup proposed in [24], who study random walks on graphs with
an additional random perfect matching.

The proofs in both models also follow similar ideas, although the random walk
we study is more challenging due to the lack of reversibility. In particular,
our results imply that the walks above have a constant spectral gap, while the
reverse is not true. Theorem 8.8 of [36] is stated without proof, and we could
not find an explanation for the reduction between the study of the random walk
on $\mathbb{F}_{p}$ to $\mathbb{P}^{1}(\mathbb{F}_{p})$ in the literature, and
so we note that the obvious analogue of Proposition 4.4 in this setting would
give a proof.

###  1.2. Bounds on solutions to a congruence equation

We also give an application of these ideas to bounding the number of solutions
to the congruence

$xy\equiv 1\pmod{p}$

with $x\in I$ and $y\in J$ for intervals $I$ and $J$ of the same length $m\leq
p/2$. In particular, for large enough $p$ and $m$, we establish that

$|\\{(x,y)\in I\times J\mid xy\equiv 1\pmod{p}\\}|\leq(1-\delta)m$

for some absolute constant $\delta>0$ (see Theorem 5.1). The proof is simple
given the spectral gap estimates in the proof of Theorem 1.2.

Our bound on the number of solutions to the congruence $xy\equiv 1\pmod{p}$ is
weaker than the bounds obtained in [11] and [14] for intervals that are not
too large and weaker than the standard ones coming from estimates on
incomplete Kloosterman sums for intervals that are not too small. The number
of solutions has also been estimated in [19, 15, 33]. In particular, the
method used in [33] also relies on $\operatorname{SL}_{2}$ action. While this
result is not new, the proof is straightforward and we hope that our ideas can
lead to further applications.

###  1.3. Related work

The random walk we study may be viewed as a non-linear version of the
Chung–Diaconis–Graham process, or the $ax+b$ process, which is the random walk
on $\mathbb{F}_{p}$ (or more generally $\mathbb{Z}/n\mathbb{Z}$ for composite
$n$) defined by

$X_{n+1}=aX_{n}+\varepsilon_{n+1},$

where $a\in\mathbb{F}_{p}^{\times}$ is fixed and the $\varepsilon_{i}$ are
independent and identically distributed, drawn from some distribution. It was
introduced in [13], where the case of $a=2$ and $\varepsilon_{i}$ distributed
uniformly on $\\{-1,0,1\\}$ was studied in detail. They showed that the mixing
time was at most order $\log p\log\log p$, and that for almost all odd $p$,
order $\log p$ steps were sufficient, although for infinitely many odd $p$,
order $\log p\log\log p$ steps were necessary. This process was subsequently
studied in [25, 26, 29, 27, 10, 18]. In particular, Eberhard and Varjú
recently showed in [18] that for almost all $p$, the walk exhibits cutoff at
$c\log p$ for some explicit constant $c\approx 1.01136$.

The Chung–Diaconis–Graham process gives an example of a speedup phenomenon
that occurs when applying bijections between steps of a Markov chain. This was
studied in [12], where it was shown that for a Markov chain on a state space
with $n$ elements, almost all bijections, when applied between steps of a
Markov chain, cause the chain to mix in $O(\log n)$ steps. Our work gives
another example of this phenomenon, which fits into a broader theme of
studying how convergence can be sped up, via either deterministic or random
means [24, 5, 2].

In [22], the first author studied a more general non-linear version of the
Chung–Diaconis–Graham process, defined by

$X_{n+1}=f(X_{n})+\varepsilon_{n+1},$

where $f$ is a bijection on $\mathbb{F}_{p}$. It was shown that for functions
$f(x)$ which were extensions of rational functions, the mixing time was of
order at most $p^{1+\varepsilon}$ for any $\varepsilon>0$, where the implicit
constant depends only on $\varepsilon$ and the degree of the polynomials
appearing in $f(x)$. In particular, this applies for the function $\iota$ that
we consider. However, the only known lower bound on the mixing time is of
order $\log p$. Our work closes the large gap between the upper and lower
bounds in this special case.

###  1.4. Random walk on Cayley graphs of
$\operatorname{SL}_{2}(\mathbb{F}_{p})$

Selberg’s Theorem [32] shows that if $S$ is a subset of
$\operatorname{SL}_{2}(\mathbb{Z})$ such that $S$ generates a subgroup of
$\operatorname{SL}_{2}(\mathbb{Z})$ of finite index, then

$\lim\sup_{p\to\infty}\lambda_{2}(A(\mathcal{G}(\operatorname{SL}_{2}(\mathbb{F}_{p}),S_{p})))<1.$

Bourgain and Gamburd [6] strengthen this result and show that the above
property holds if and only if $S$ generates a non-elementary subgroup of
$\operatorname{SL}_{2}(\mathbb{Z})$. Weigel [35] gives a convenient
characterization that $\langle S\rangle$ is non-elementary if and only if
$\langle S_{p}\rangle=\operatorname{SL}_{2}(\mathbb{F}_{p})$ for some prime
$p\geq 5$ (see also [30]).

We remark that a simple reduction allows us to deduce Theorem 1.2 from the
special case where the distribution of $\varepsilon_{n}$ is supported on two
points. In this case, we can construct the set $S$ corresponding to the
fractional Markov chain and easily verify that $\langle S\rangle$ is non-
elementary.

Golumbev and Kamper [21] show that under more stringent assumptions on the set
of generators $S$, the Cayley graphs
$\mathcal{G}(\operatorname{SL}_{2}(\mathbb{F}_{p}),S_{p})$ have a stronger
property that implies cut-off for the associated projected random walk on
$\mathbb{P}^{1}(\mathbb{F}_{p})$. This property also holds with high
probability for random generating sets $S_{p}$. We also remark that the result
of Bourgain-Gamburd has been further developed in various aspects, e.g., [8,
9, 20, 21]. These results also have many applications in number theory,
including sum product problems on finite fields [7, 23].

###  1.5. Outline

In Section 2, we prove Theorem 1.1. In Section 3, we explain the comparison
theory for Markov chains used to relate the spectral gaps for the walks on
$\mathbb{F}_{p}$ and $\mathbb{P}^{1}(\mathbb{F}_{p})$. In Section 4, we give a
proof of Theorem 1.2. Finally, in Section 5, we give a bound for the number of
solutions to $xy\equiv 1\pmod{p}$ in a square $I_{1}\times
I_{2}\subseteq\mathbb{F}_{p}^{2}$.

###  1.6. Notations

Throughout, we let $\iota:\mathbb{F}_{p}\to\mathbb{F}_{p}$ be the function
defined by $\iota(x)=1/x$ if $x\neq 0$ and $\iota(x)=0$ if $x=0$. We view
$\mathbb{P}^{1}(\mathbb{F}_{p})$ as a superset of $\mathbb{F}_{p}$ with one
extra element $\infty$. We define the function
$\overline{\iota}:\mathbb{P}^{1}(\mathbb{F}_{p})\to\mathbb{P}^{1}(\mathbb{F}_{p})$
by $\overline{\iota}(x)=1/x$ if $x\neq 0,\infty$, and
$\overline{\iota}(0)=\infty$ and $\overline{\iota}(\infty)=0$.

If $P$ is a symmetric matrix, we let $\lambda_{i}(P)$ denote the $i$th largest
eigenvalue of $P$.

##  2\. Lower bound

In this section, we prove Theorem 1.1, which follows from entropy
considerations. For a discrete probability measure $\mu$, we let

$H(\mu)=\sum-\mu(x)\log\mu(x)$

denote the _entropy_ of $\mu$, with the convention that $0\log 0=0$. If $X$ is
a random variable, we let $H(X)$ denote the entropy of the law of $X$. We need
two basic properties of entropy. The first is that it is subadditive, in the
sense that if $X$ and $Y$ are independent, then $H(X+Y)\leq H(X)+H(Y)$. The
second is that for any function $f$, $H(f(X))\leq H(X)$, with equality if $f$
is bijective.

The lower bound follows by noting that the entropy of the random walk at time
$n$ is at most $nH(\mu)$, which is too low to be close to uniform if $n$ is
too small. This idea is certainly not new, and similar arguments have appeared
before, see [18, 1, 4] for example. Our argument follows [1], although we fill
in some details and fix a minor error.

The following lemma shows that measures close to uniform in total variation
must have large entropy. It was stated without the $\log 2$ term in [1], but
this cannot be correct, as can be seen if $\mu$ is uniform and $\nu$ is
concentrated at a single point. It can be seen as a special case of Theorem 1
of [3] (in fact, it is established as part of the proof).

######  Lemma 2.1.

Let $X$ be a finite set of size $n$, and let $\pi$ denote the uniform measure
on $X$. Let $\nu$ be any probability measure on $X$, and let
$\delta=\|\nu-\pi\|_{TV}$. Then

$\begin{split}|H(\mu)-H(\nu)|&\leq\delta\log(n-1)-\delta\log\delta-(1-\delta)\log(1-\delta)\\\
&\leq\delta\log n+\log(2).\end{split}$

###### Proof.

The first inequality is Equation 11 of [3] and the second inequality is clear.
∎

###### Proof of Theorem 1.1.

Because $\iota$ is a bijection, $H(X)=H(\iota(X))$. Furthermore, if $\mu_{p}$
denotes the distribution given by reducing $\mu$ mod $p$, then $H(\mu_{p})\leq
H(\mu)$. Then by subadditivity of entropy, $H(K^{n}(x,\cdot))\leq nH(\mu)$.
Thus, Lemma 2.1 gives

$\|K^{n}(x,\cdot)-\pi\|_{TV}\geq 1-\frac{nH(\mu)+\log 2}{\log p}.$

∎

######  Remark 2.2.

We note that the finite entropy assumption cannot be removed entirely. Fix any
constant $\delta\in(0,1)$. Consider the measure on $\mathbb{Z}$ given by

$\mu(i)\propto\frac{1}{i(\log i)^{2-\delta}}.$

There is some constant $c=c(\delta)>0$ such that for all $p$ and $i\leq p$,

$\sum_{k\in\mathbb{Z}}\mu(i+kp)\geq\sum_{k\geq
p}\frac{1}{kp(\log(kp))^{2-\delta}}\geq\frac{c}{p(\log p)^{1-\delta}}.$

Then on $\mathbb{F}_{p}$, the walk defined by (1.1) satisfies

$\mathbb{P}(X_{n+1}=x|X_{n}=y)\geq\frac{c}{p(\log p)^{1-\delta}}$

for all $x$ and $y$, and so we can couple two copies of the walk, $X_{n}$ and
$X_{n}^{\prime}$, so that independent of the initial states,

$\mathbb{P}(X_{n}\neq X^{\prime}_{n})\leq\left(1-\frac{c}{(\log
p)^{1-\delta}}\right)^{n}.$

Since the distance to stationarity can be controlled by this probability (see
[28, Corollary 5.5]), this implies that

$\|K^{n}(x,\cdot)-\pi\|_{TV}\leq\left(1-\frac{c}{(\log
p)^{1-\delta}}\right)^{n},$

and so $t_{mix}(\varepsilon)=O_{\varepsilon}((\log p)^{1-\delta})$, which is
much smaller than $\log p$.

For positive integers $k$, define $\log^{(k)}$ iteratively by
$\log^{(1)}(x)=\log x$ and $\log^{(k)}(x)=\log(\log^{(k-1)}(x))$ for $k\geq
2$. The measure

$\mu(i)\propto\frac{1}{i(\log
i)(\log^{(2)}i)\cdots(\log^{(k-1)}i)(\log^{(k)}i)^{2}}$

gives a probability distribution on $\mathbb{Z}$ with infinite entropy for
which the mixing time $t_{mix}(\varepsilon)=O_{\varepsilon}(\log^{(k)}p)$.

##  3\. Comparison theory

Comparison theory for Markov chains was introduced by Diaconis and Saloff-
Coste in [17, 16]. This theory allows the spectral gaps (and also log-Sobolev
constants) of different Markov chains on the same state space to be compared.
Unfortunately, since we wish to compare a Markov chain on $\mathbb{F}_{p}$
with one on $\mathbb{P}^{1}(\mathbb{F}_{p})$, this theory does not immediately
apply. Smith extended these ideas in [34] to the case of random walks on state
spaces $X_{0}\subseteq X$.

We first recall the relationship between the Dirichlet form and the spectral
gap, and then explain the comparison theory developed by Smith for Markov
chains on different state spaces. This is used to compare the spectral gaps of
the random walks on $\mathbb{P}^{1}(\mathbb{F}_{p})$ and $\mathbb{F}_{p}$. For
a more thorough treatment and proofs, we refer the reader to [34].

###  3.1. The Dirichlet form

Recall the following standard notions. Let $P$ be a reversible Markov chain on
a finite space $X$ with stationary distribution $\pi$. Define the function

$V_{\pi}(f)=\frac{1}{2}\sum_{x,y\in X}|f(x)-f(y)|^{2}\pi(x)\pi(y)$

and the _Dirichlet form_

$\mathcal{E}_{P}(f,f)=\frac{1}{2}\sum_{x,y\in X}|f(x)-f(y)|^{2}P(x,y)\pi(x)$

for $f:X\to\mathbb{R}$. Then the spectral gap can be computed by

(3.1) $1-\lambda_{2}(P)=\inf_{\begin{subarray}{c}f:X\to\mathbb{R}\\\ f\text{
not constant}\end{subarray}}\frac{\mathcal{E}_{P}(f,f)}{V_{\pi}(f)}.$

###  3.2. Comparison theory for Markov chains on different state spaces

Let $X_{0}\subseteq X$, and let $P_{0}$ be a Markov chain on $X_{0}$ and $P$
be a Markov chain on $X$, with stationary distributions $\pi_{0}$ and $\pi$
respectively. Assume that both are ergodic and reversible. The following
results of Smith compare $V_{\pi_{0}}$ and $\mathcal{E}_{P_{0}}$ with
$V_{\pi}$ and $\mathcal{E}_{P}$.

######  Lemma 3.1 ([34, Lemma 2]).

Let $f_{0}:X_{0}\to\mathbb{R}$ and let $f:X\to\mathbb{R}$ be any extension of
$f_{0}$. Then

$V_{\pi_{0}}(f_{0})\leq CV_{\pi}(f),$

where

$C=\sup_{x\in X_{0}}\frac{\pi_{0}(x)}{\pi(x)}.$

The analogous comparison result for the Dirichlet form requires further setup.
For each $x\in X$, fix a probability measure $Q_{x}$ on $X_{0}$ such that
$Q_{x}=\delta_{x}$ if $x\in X_{0}$. For any $f_{0}:X_{0}\to\mathbb{R}$, this
defines an extension $f:X\to\mathbb{R}$ by

$f(x)=\sum_{y\in X_{0}}Q_{x}(y)f(y).$

Next, choose a coupling of $Q_{x}$ and $Q_{y}$ for all $x,y\in X$ for which
$P(x,y)>0$ (note this is a symmetric relation). Denote these couplings by
$Q_{x,y}$.

A _path_ in $X_{0}$ from $x$ to $y$ is a sequence of $x_{i}\in X_{0}$ for
$i=0,\dotsc,k$ for which $x_{0}=x$, $x_{k}=y$ and $P_{0}(x_{i},x_{i+1})>0$ for
all $i$. For a path $\gamma$ from $x$ to $y$, call $x$ the _initial vertex_
and $y$ the _final vertex_ , and denote them by $i(\gamma)$ and $o(\gamma)$
respectively. Let $|\gamma|$ denote the length of the path.

A _flow_ on $X_{0}$ is a function $F$ on the set of paths in $X_{0}$ whose
restriction to paths from $x$ to $y$ gives a probability measure for all $x$,
$y$. We require a choice of flow on $X_{0}$. Actually it suffices to define
the flow only for paths from $a$ to $b$ such that there exists $x,y\in X$ with
$Q_{x,y}(a,b)>0$. For our purposes, it will suffice to choose the flow such
that when restricted to paths from $a$ to $b$, it concentrates on a single
path.

The following theorem compares $\mathcal{E}_{P_{0}}$ with $\mathcal{E}_{P}$ in
terms of the chosen data.

######  Theorem 3.2 ([34, Theorem 4]).

Consider the setup defined above. Then for $f_{0}:X_{0}\to\mathbb{R}$ and $f$
the extension to $X$ defined by the measures $Q_{x}$, we have

$\mathcal{E}_{P}(f,f)\leq\mathcal{A}\mathcal{E}_{P_{0}}(f_{0},f_{0}),$

where

$\begin{split}\mathcal{A}=&\sup_{P_{0}(x,y)>0}\frac{1}{P_{0}(x,y)\pi_{0}(x)}\bigg{(}\sum_{\gamma\ni(x,y)}|\gamma|F(\gamma)P(i(\gamma),o(\gamma))\pi(i(\gamma))\\\
&\qquad+2\sum_{\gamma\ni(x,y)}|\gamma|F(\gamma)\sum_{b\not\in
X_{0}}Q_{b}(o(\gamma))P(i(\gamma),b)\pi(i(\gamma))\\\
&\qquad+\sum_{\gamma\ni(x,y)}|\gamma|F(\gamma)\sum_{\begin{subarray}{c}a,b\not\in
X_{0}\\\
P(x,y)>0\end{subarray}}Q_{a,b}(i(\gamma),o(\gamma))P(a,b)\pi(a)\bigg{)}.\end{split}$

Note that by (3.1), an immediate consequence of these comparison results is
that

$1-\lambda_{2}(P_{0})\geq\frac{1}{C\mathcal{A}}(1-\lambda_{2}(P)).$

##  4\. Upper bound

In this section, we prove Theorem 1.2. The proof of the upper bound consists
of three main steps. First, we control the mixing time by the spectral gap of
a symmetrized random walk (Corollary 4.2). Next, we relate the spectral gap of
the symmetrized walk on $\mathbb{F}_{p}$ with a walk on
$\mathbb{P}^{1}(\mathbb{F}_{p})$ (Proposition 4.4). Finally, we show that the
spectral gap of the walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$ is of constant
order (Proposition 4.8).

###### Proof of Theorem 1.2.

Let $P$ denote the Markov matrix encoding the step $X\mapsto X+\varepsilon$
and let $\Pi$ denote the Markov matrix encoding $X\mapsto\iota(X)$. Note that
$P^{T}$ encodes the step $X\mapsto X-\varepsilon$. By Corollary 4.2, it
suffices to show that

$\lambda_{2}(P^{T}\Pi P^{T}P\Pi P)\leq 1-c$

for some constant $c>0$, independent of $p$.

Let $a_{1},a_{2}\in\mathbb{Z}$ be two distinct elements in the support of
$\mu$ and write $b=a_{1}-a_{2}$. Then in the Markov chain given by the
transition matrix $P^{T}\Pi P^{T}P\Pi P$, there is $u>0$ depending only on
$\mu$ such that we transition from $x$ to one of $x+b$, $x-b$,
$\iota(\iota(x+a_{1})+b)-a_{1}$ and $\iota(\iota(x+a_{1})-b)-a_{1}$ with
probability at least $u$. For example, note that
$x+b=\iota(\iota(x+a_{1})+a_{1}-a_{1})-a_{2}$ and so there is a positive
probability of moving from $x$ to $x+b$.

Then we can write $P^{T}\Pi P^{T}P\Pi P=uL_{0}+(1-u)L^{\prime}$ where $L_{0}$
is the transition matrix for the random walk going from $x$ to one of $x+b$,
$x-b$, $\iota(\iota(x+a_{1})+b)-a_{1}$ and $\iota(\iota(x+a_{1})-b)-a_{1}$
with equal probability, and $L^{\prime}$ a symmetric stochastic matrix. We
have

$\lambda_{2}(P^{T}\Pi P^{T}P\Pi P)\leq u\lambda_{2}(L_{0})+1-u,$

and so it suffices to show that $\lambda_{2}(L_{0})\leq 1-c$ for some constant
$c>0$.

Let
$\overline{\iota}:\mathbb{P}^{1}(\mathbb{F}_{p})\to\mathbb{P}^{1}(\mathbb{F}_{p})$
be the function defined by $\overline{\iota}(x)=1/x$ if $x\neq 0,\infty$, and
$\overline{\iota}(0)=\infty$, and $\overline{\iota}(\infty)=0$. Let $L$ denote
the transition matrix for the random walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$
going from $x$ to one of $x+b$, $x-b$,
$\overline{\iota}(\overline{\iota}(x+a_{1})+b)-a_{1}$ and
$\overline{\iota}(\overline{\iota}(x+a_{1})-b)-a_{1}$ with equal probability.
By Proposition 4.4, we can instead show

$\lambda_{2}(L)\leq 1-c$

for the walk defined by $L$ instead of the walk defined by $L_{0}$, and this
is exactly given by Proposition 4.8. ∎

###  4.1. Reduction to the symmetric walk

Recall that we are interested in studying the walk
$X_{n+1}=\iota(X_{n})+\varepsilon_{n+1}$ on $\mathbb{F}_{p}$, where
$\varepsilon_{n}$ are independent and identically distributed according to
$\mu$. This random walk is non-reversible, so the first step is to relate it
to a suitable symmetrization.

Let $\Pi$ denote the transition matrix for the (deterministic) walk
$X\mapsto\iota(X)$ on $\mathbb{F}_{p}$ and let $P$ denote the walk $X\mapsto
X+\varepsilon$. Note that $P^{T}$ is also a Markov matrix and encodes the walk
$X\mapsto X-\varepsilon$. Let $\lambda_{2}$ be the second largest eigenvalue
of $P^{T}\Pi P^{T}P\Pi P$. Note that all eigenvalues of $P^{T}\Pi P^{T}P\Pi P$
are real and non-negative.

The following result from [12] relates the original walk $K=P\Pi$ to the
symmetrized walk $P^{T}\Pi P^{T}P\Pi P$.

######  Lemma 4.1 ([12, Corollary 5.4]).

Let $k\geq 2$ and let $x\in\mathbb{R}^{p}$ such that $\sum x_{i}=0$. Then

$\|(K^{T})^{k}x\|_{2}\leq\lambda_{2}^{(k-2)/4}\|x\|_{2}.$

This result means that if $P^{T}\Pi P^{T}P\Pi P$ can be shown to have a
spectral gap of constant order, then the walk $K$ mixes in order $\log p$
steps. The following is an easy consequence of Lemma 4.1 (see Equation 5.5 in
[12]).

######  Corollary 4.2.

Let $\pi$ denote the uniform measure on $\mathbb{F}_{p}$. Then for all
$x\in\mathbb{F}_{p}$,

$\|K^{k}(x,\cdot)-\pi\|_{TV}\leq\frac{\sqrt{p}}{2}\lambda_{2}^{(k-2)/4}.$

###  4.2. Comparison of the random walks on $\mathbb{P}^{1}(\mathbb{F}_{p})$
and $\mathbb{F}_{p}$

We now explain how to apply the comparison results to bound the spectral gap
of the walk on $\mathbb{F}_{p}$ in terms of the spectral gap of the walk on
$\mathbb{P}^{1}(\mathbb{F}_{p})$.

Recall from the proof of Theorem 1.2 that we denote by $L_{0}$ the transition
matrix of the random walk on $\mathbb{F}_{p}$ which moves from
$x\in\mathbb{F}_{p}$ to one of $x+b$, $x-b$, $\iota(\iota(x+a_{1})+b)-a_{1}$
and $\iota(\iota(x+a_{1})-b)-a_{1}$ with equal probability. We denote by $L$
the transition matrix of the walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$ which
moves from $x\in\mathbb{F}_{p}$ to one of $x+b$, $x-b$,
$\overline{\iota}(\overline{\iota}(x+a_{1})+b)-a_{1}$ and
$\overline{\iota}(\overline{\iota}(x+a_{1})-b)-a_{1}$. Here, recall that
$\overline{\iota}(x)=1/x$ if $x\neq 0,\infty$, and
$\overline{\iota}(0)=\infty$, and $\overline{\iota}(\infty)=0$. The following
lemma is useful in constructing the data required in Theorem 3.2.

######  Lemma 4.3.

For all $x,y\in\mathbb{F}_{p}$ except $x=-a_{1}$ and $y=-a_{1}$, if
$L(x,y)>0$, then $L_{0}(x,y)>0$.

###### Proof.

First, note that if $x\neq-a_{1},b^{-1}-a_{1},-b^{-1}-a_{1}$, then the
statement is clear since in this case the functions $\iota$ and
$\overline{\iota}$ are identical for the transitions involved. It can be
checked that if $x=b^{-1}-a_{1}$ or $x=-b^{-1}-a_{1}$, and $y\neq\infty$ and
$L(x,y)>0$, then $L_{0}(x,y)>0$. Similarly, we can check that if $x=-a_{1}$
and $L(x,y)>0$ then $L_{0}(x,y)>0$ unless $y=-a_{1}$. ∎

######  Proposition 4.4.

Let $L_{0}$ and $L$ denote the transition matrices for the random walk on
$\mathbb{F}_{p}$ and $\mathbb{P}^{1}(\mathbb{F}_{p})$ respectively. Then there
exists an absolute constant $c>0$ such that

$1-\lambda_{2}(L_{0})\geq c(1-\lambda_{2}(L)).$

###### Proof.

We begin by defining the data needed to apply Theorem 3.2. Take
$X_{0}=\mathbb{F}_{p}$ and $X=\mathbb{P}^{1}(\mathbb{F}_{p})$. Define for
$x\in\mathbb{F}_{p}$, $Q_{x}=\delta_{x}$ and $Q_{\infty}$ is uniform on the
set of $y\in\mathbb{F}_{p}$ such that $L(\infty,y)>0$. Take all couplings to
be independent, so $Q_{x,y}(x_{0},y_{0})=Q_{x}(x_{0})Q_{y}(y_{0})$. Then note
that the condition that $Q_{x,y}(x_{0},y_{0})>0$ for some $x,y\in X$ is
exactly that either $L(x_{0},y_{0})>0$, or $L(x_{0},\infty)>0$ and
$L(y_{0},\infty)>0$.

By Lemma 4.3, if $L(x_{0},y_{0})>0$, then $L_{0}(x_{0},y_{0})>0$ except when
$x_{0}=-a_{1}$ and $y_{0}=-a_{1}$. Thus, if $L(x_{0},y_{0})>0$, we pick our
flow $F$ to be concentrated on the path consisting of the single edge
$(x_{0},y_{0})$, except when $x_{0}=y_{0}=-a_{1}$ where we take the path of
length 2 given by $-a_{1}\mapsto b^{-1}-a_{1}\mapsto-a_{1}$.

If $L(x_{0},\infty)>0$ and $L(y_{0},\infty)>0$, then there must be a path of
length $2$ connecting $x_{0}$ to $y_{0}$. This is because
$x_{0},y_{0}\in\\{b^{-1}-a_{1},-b^{-1}-a_{1}\\}$ and so we can take our flow
to be concentrated on the path $x_{0}\mapsto-a_{1}\mapsto y_{0}$. Thus, our
flow is concentrated entirely on paths of lengths $1$ and $2$.

We can now compute the constants $C$ and $\mathcal{A}$. It’s easy to see that
$C\leq 2$. To bound $\mathcal{A}$, note that for each $(x,y)$ for which
$L_{0}(x,y)>0$, there is at most one path $\gamma$ of length $1$ on which
$F(\gamma)\neq 0$ containing $(x,y)$, namely $\gamma=(x,y)$, and there is a
bounded number of paths $\gamma$ of length $2$ containing $(x,y)$ for which
$F(\gamma)>0$. Thus, for any $(x,y)$ for which $L_{0}(x,y)>0$, there are a
bounded number of terms in the summation in the definition of $\mathcal{A}$,
each of bounded size. Hence, $\mathcal{A}$ is bounded above by a constant.
This implies that

$1-\lambda_{2}(L_{0})\geq c(1-\lambda_{2}(L))$

for some absolute constant $c>0$. ∎

######  Remark 4.5.

For a reversible Markov chain $P$ with stationary distribution $\pi$, let

$\phi(P)=\min_{\pi(S)\leq 1/2}\frac{\sum_{x\in S,y\in
S^{c}}\pi(x)P(x,y)}{\pi(S)}$

denote the _bottleneck ratio_ (sometimes also called the _Cheeger constant_ or
_conductance_). We note that a constant order spectral gap for $L_{0}$ could
also be derived using Cheeger’s inequality (see [28, Theorem 13.10] for
example), which states that

$\frac{\phi(P)^{2}}{2}\leq 1-\lambda_{2}(P)\leq 2\phi(P).$

The bottleneck ratios of $L_{0}$ and $L$ can be easily compared, which would
give a uniform lower bound for $1-\lambda_{2}(L_{0})$ using the one for
$1-\lambda_{2}(L)$. We prefer to use comparison theory because the argument
can be more easily generalized. In particular, it is more robust and would
give a sharper result if $L_{0}$ did not have a constant order spectral gap.

###  4.3. Spectral gap for random walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$

In this section, we prove that the random walk on
$\mathbb{P}^{1}(\mathbb{F}_{p})$ has a constant order spectral gap.

Recall that $\operatorname{SL}_{2}(\mathbb{F}_{p})$ has a transitive action on
$\mathbb{P}^{1}(\mathbb{F}_{p})$, viewed as lines in $\mathbb{F}_{p}^{2}$.
This is through linear fractional transformations, and may be formally defined
by

(4.1) $\begin{pmatrix}a&b\\\ c&d\\\ \end{pmatrix}\cdot
x=\begin{cases}\frac{ax+b}{cx+d}&\text{if $cx+d\neq 0$ and
$x\in\mathbb{F}_{p}$}\\\ \infty&\text{if $cx+d=0$ and $x\in\mathbb{F}_{p}$}\\\
\frac{a}{c}&\text{if $x=\infty$ and $c\neq 0$}\\\ \infty&\text{if $x=\infty$
and $c=0$}.\end{cases}$

The random walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$ given by moving from $x$
to one of $x+b$, $x-b$, $\overline{\iota}(\overline{\iota}(x+a_{1})+b)-a_{1}$
and $\overline{\iota}(\overline{\iota}(x+a_{1})-b)-a_{1}$ uniformly at random
can be viewed as the quotient of a random walk on
$\operatorname{SL}_{2}(\mathbb{F}_{p})$. The following result of Bourgain and
Gamburd gives a constant order spectral gap for random walks on
$\operatorname{SL}_{2}(\mathbb{F}_{p})$.

######  Theorem 4.6 ([6, Theorem 1]).

Let $S\subseteq\operatorname{SL}_{2}(\mathbb{Z})$ be a symmetric set which
generates a non-elementary subgroup of $\operatorname{SL}_{2}(\mathbb{Z})$.
Let $S_{p}$ denote the set of generators mod $p$ for a prime $p$. Let $P$
denote the transition matrix for the random walk on
$\operatorname{SL}_{2}(\mathbb{F}_{p})$ defined by

$X_{n+1}=X_{n}\varepsilon_{n+1},$

where the $\varepsilon_{i}$ are independent and uniformly distributed on
$S_{p}$. Then for all primes $p$ large enough,

$\lambda_{2}(P)\leq 1-c$

for some constant $c>0$ independent of $p$.

To show that this applies to the random walk we wish to study, we need the
following lemma.

######  Lemma 4.7.

Let

$S=\left\\{\begin{pmatrix}1&b\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}1-a_{1}b&-a_{1}^{2}b\\\ b&a_{1}b+1\\\
\end{pmatrix},\begin{pmatrix}1&-b\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}1+a_{1}b&a_{1}^{2}b\\\ -b&1-a_{1}b\\\
\end{pmatrix}\right\\}.$

Then $S$ generates a non-elementary subgroup of
$\operatorname{SL}_{2}(\mathbb{Z})$.

###### Proof.

By Theorem 2.5 of [30], non-elementary subgroups are the same as Zariski-dense
subgroups in $\operatorname{SL}_{2}(\mathbb{Z})$. By a result of Weigel [35]
(see also Theorem 2.2 of [30]), it suffices to show that after reducing mod
$p$ for some $p\geq 5$, $S_{p}$ generates
$\operatorname{SL}_{2}(\mathbb{F}_{p})$.

If $b\neq 0\pmod{p}$, then the matrices $\begin{pmatrix}1&b\\\ 0&1\\\
\end{pmatrix}$ and $\begin{pmatrix}1&-b\\\ 0&1\\\ \end{pmatrix}$ generate all
matrices of the form $\begin{pmatrix}1&t\\\ 0&1\\\ \end{pmatrix}$ for
$t\in\mathbb{F}_{p}$.

Using $\begin{pmatrix}1&t\\\ 0&1\\\ \end{pmatrix}$ and
$\begin{pmatrix}1-a_{1}b&-a_{1}^{2}b\\\ b&a_{1}b+1\\\ \end{pmatrix}$, we can
generate all matrices of the form

$\begin{pmatrix}1&t^{\prime}\\\ 0&1\\\
\end{pmatrix}\cdot\begin{pmatrix}1-a_{1}b&-a_{1}^{2}b\\\ b&a_{1}b+1\\\
\end{pmatrix}\cdot\begin{pmatrix}1&t\\\ 0&1\\\ \end{pmatrix}.$

In particular, we can generate all matrices in the subset $X_{b}$ of
$\operatorname{SL}_{2}(\mathbb{F}_{p})$ of matrices whose lower left corner
equal to $b$.

Similarly, $S_{p}$ generates matrices in the subset $X_{-b}$ of
$\operatorname{SL}_{2}(\mathbb{F}_{p})$ of matrices with lower left corner
equal to $-b$, using $\begin{pmatrix}1+a_{1}b&a_{1}^{2}b\\\ -b&1-a_{1}b\\\
\end{pmatrix}$ instead of $\begin{pmatrix}1-a_{1}b&-a_{1}^{2}b\\\
b&a_{1}b+1\\\ \end{pmatrix}$.

Finally, we can easily check that $X_{b}$ and $X_{-b}$ generate
$\operatorname{SL}_{2}(\mathbb{F}_{p})$. Thus, $S_{p}$ generates
$\operatorname{SL}_{2}(\mathbb{F}_{p})$ for all $p>|b|$, and hence, $S$
generates a non-elementary subgroup of $\operatorname{SL}_{2}(\mathbb{Z})$. ∎

Since the random walk on $\mathbb{P}^{1}(\mathbb{F}_{p})$ is a quotient of a
random walk on $\operatorname{SL}_{2}(\mathbb{F}_{p})$, we can obtain the
desired bound on the spectral gap.

######  Proposition 4.8.

Let $L$ denote the transition matrix for the random walk on
$\mathbb{P}^{1}(\mathbb{F}_{p})$ which moves from $x\in\mathbb{F}_{p}$ to one
of $x+b$, $x-b$, $\overline{\iota}(\overline{\iota}(x+a_{1})+b)-a_{1}$ and
$\overline{\iota}(\overline{\iota}(x+a_{1})-b)-a_{1}$ uniformly at random.
Then

$\lambda_{2}(L)\leq 1-c$

for some constant $c>0$ independent of $p$.

###### Proof.

Note that the formulas given in 4.1 actually define a
$\operatorname{GL}_{2}(\mathbb{F}_{p})$ action. It’s easy to see that $\iota$
and addition by $a\in\mathbb{F}_{p}$ are both linear frational
transformations, represented by the matrices

$\begin{pmatrix}1&a\\\ 0&1\\\ \end{pmatrix},\begin{pmatrix}0&1\\\ 1&0\\\
\end{pmatrix}$

in $\operatorname{GL}_{2}(\mathbb{F}_{p})$. Then we may write the function
$x\mapsto x+b$ as a product of these matrices, and similarly for $x-b$,
$\overline{\iota}(\overline{\iota}(x+a_{1})+b)-a_{1}$ and
$\overline{\iota}(\overline{\iota}(x+a_{1})-b)-a_{1}$. This allows us to view
the walk defined by $L$ as the quotient of the walk on
$\operatorname{SL}_{2}(\mathbb{F}_{p})$ generated by the set

$S=\left\\{\begin{pmatrix}1&b\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}1-a_{1}b&-a_{1}^{2}b\\\ b&a_{1}b+1\\\
\end{pmatrix},\begin{pmatrix}1&-b\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}1+a_{1}b&a_{1}^{2}b\\\ -b&1-a_{1}b\\\
\end{pmatrix}\right\\},$

with transition matrix $\widetilde{L}$. This means that the spectrum of $L$ is
contained in the spectrum of $\widetilde{L}$, and in particular
$\lambda_{2}(L)\leq\lambda_{2}(\widetilde{L})$. But since $S$ generates a non-
elementary subgroup of $\operatorname{SL}_{2}(\mathbb{Z})$ by Lemma 4.7,
Theorem 4.6 implies that

$\lambda_{2}(\widetilde{L})\leq 1-c$

for some constant $c>0$ independent of $p$, giving us the desired bound. ∎

##  5\. Points on modular hyperbolas

In this section, we prove the following theorem stating that for all subsets
$A\subseteq\mathbb{F}_{p}$ of size at most $p/2$, the sets $A$ and $\iota(A)$
cannot both be close to intervals. As a corollary, we obtain bounds on the
number of solutions to the congruence $xy=1\pmod{p}$ lying in $I\times
J\subseteq\mathbb{F}_{p}^{2}$ for $I$ and $J$ intervals of the same length.

######  Theorem 5.1.

Let $I$ and $J$ be two intervals in $\mathbb{F}_{p}$, each of length $m\leq
p/2$, and let $A\subseteq\mathbb{F}_{p}$ with $A\subseteq I$ and
$\iota(A)\subseteq J$. There exists an absolute constant $\delta>0$ such that
for all $p$ and $m$ sufficiently large, $|A|\leq(1-\delta)m$.

###### Proof.

Throughout the proof, we use the notation $[-k,k]$ to denote the set of
integers $\\{-k,-k+1,\dots,k-1,k\\}$. Let $Q=P^{T}\Pi P^{T}P\Pi P$, where $P$
denotes the transition matrix for the random walk on $\mathbb{F}_{p}$
generated by the uniform measure on $[-1,1]$ and $\Pi$ encodes the bijection
$\iota$. Then $Q$ is the transition matrix of a reversible Markov chain on
$\mathbb{F}_{p}$, and from the proof of Theorem 1.2 it has a constant order
spectral gap, i.e. for large enough $p$, $1-\lambda_{2}(Q)\geq\gamma$ for some
absolute constant $\gamma>0$.

Assume that there exists intervals $I$, $J$ such that $A\subseteq I$,
$\iota(A)\subseteq J$, and $|A|\geq(1-\delta)m$. We will show using Cheeger’s
inequality that this implies the spectral gap cannot be of constant order,
which gives a contradiction.

Observe that $Q(x,y)\geq c>0$ for some constant $c$ independent of $p$. Since
the stationary distribution is uniform, the bottleneck ratio of the chain
$X_{n}$ is given by

$\Phi=\min_{S\subseteq\mathbb{F}_{p},|S|\leq p/2}\frac{\sum_{x\in S,y\notin
S}Q(x,y)}{|S|}.$

By Cheeger’s inequality, we have $\Phi$ is bounded below by a positive
constant $\gamma^{\prime}$ depending only on $\gamma$. Thus,

$\frac{\sum_{x\in A,y\notin A}Q(x,y)}{|A|}\geq\gamma^{\prime}.$

For $x\in A$, we consider $y\in\mathbb{F}_{p}$ for which $Q(x,y)>0$. Then

$y\in\iota(\iota(A+[-1,1])+[-2,2])+[-1,1].$

Note that $\iota(A+[-1,1])\subseteq\iota(I+[-1,1])$, and

$\displaystyle\iota(I+[-1,1])$
$\displaystyle=\iota(A\cup((I+\\{-1,+1\\})\setminus A))$
$\displaystyle=\iota(A)\cup\iota((I+\\{-1,+1\\})\setminus A).$

Let $S_{1}=\iota((I+\\{-1,+1\\})\setminus A)$. Then
$\iota(A+[-1,1])\subseteq\iota(A)\cup S_{1}$, and $|S_{1}|\leq 2+\delta m$.

Similarly, since $\iota(A)\cup S_{1}\subseteq J\cup S_{1}$, we have

$\displaystyle\iota((\iota(A)\cup S_{1})+[-2,2])$
$\displaystyle\subseteq\iota(J+[-2,2])\cup\iota(S_{1}+[-2,2])$
$\displaystyle\subseteq
A\cup\iota((J+[-2,2])\setminus\iota(A))\cup\iota(S_{1}+[-2,2]).$

Here, $|\iota((J+[-2,2])\setminus\iota(A))|\leq 4+\delta m$, and
$|\iota(S_{1}+[-2,2])|\leq 5|S_{1}|$. Letting

$S_{2}=\iota((J+[-2,2])\setminus\iota(A))\cup\iota(S_{1}+[-2,2]),$

we have

$\iota(\iota(A+[-1,1])+[-2,2])\subseteq A\cup S_{2},$

where $|S_{2}|\leq 14+6\delta m$.

Finally, we have $y\in(A\cup S_{2})+[-1,1]$, and

$\displaystyle(A\cup S_{2})+[-1,1]$ $\displaystyle\subseteq(I\cup
S_{2})+[-1,1]$ $\displaystyle\subseteq(I+[-1,1])\cup(S_{2}+[-1,1])$
$\displaystyle\subseteq A\cup((I+[-1,1])\setminus A)\cup(S_{2}+[-1,1]).$

Let $S_{3}=((I+[-1,1])\setminus A)\cup(S_{2}+[-1,1])$, and note that
$|S_{3}|\leq 44+19\delta m$.

Thus, if $x\in A$, $y\notin A$, and $Q(x,y)>0$, then we must have $y\in
S_{3}$. Hence,

$\frac{\sum_{x\in A,y\notin A}Q(x,y)}{|A|}\leq\frac{\sum_{x\in A,y\in
S_{3}}Q(x,y)}{|A|}\leq\frac{44+19\delta
m}{(1-\delta)m}\leq\frac{20\delta}{1-\delta},$

where we assumed that $m$ is sufficiently large (in $\delta$), and used that

$\sum_{x\in A,y\in S_{3}}Q(x,y)\leq\sum_{x\in\mathbb{F}_{p},y\in
S_{3}}Q(y,x)=|S_{3}|.$

For $\delta$ sufficiently small, $20\delta/(1-\delta)<\gamma^{\prime}$, and we
have a contradiction. ∎

The following corollary follows easily from Theorem 5.1 by taking $A=\\{x\in
I\mid\iota(x)\in J\\}$, which away from $0$ counts the solutions of interest.

######  Corollary 5.2.

Let $I$ and $J$ be two intervals in $\mathbb{F}_{p}$ of the same length $m\leq
p/2$. There exists an absolute constant $\delta>0$ such that for sufficiently
large $p$ and $m$,

$\left|\\{(x,y)\in I\times J\mid xy\equiv 1\pmod{p}\\}\right|\leq(1-\delta)m.$

## Acknowledgements

The authors thank Sourav Chatterjee, Persi Diaconis, Jacob Fox, Sean Eberhard,
Ilya Shkredov, Kannan Soundararajan, Péter Varjú, Thuy Duong Vuong, and Yuval
Wigderson for their help and comments on earlier drafts. We are greatly
indebted to Péter Varjú for pointing out the many useful references to the
study of spectral gap of Cayley graphs of
$\operatorname{SL}_{2}(\mathbb{F}_{p})$, in particular, [6].

## References

  * [1] David Aldous.  Random walks on finite groups and rapidly mixing Markov chains.  In Seminar on probability, XVII, volume 986 of Lecture Notes in Math., pages 243–297. Springer, Berlin, 1983.
  * [2] Noga Alon, Itai Benjamini, Eyal Lubetzky, and Sasha Sodin.  Non-backtracking random walks mix faster.  Commun. Contemp. Math., 9(4):585–603, 2007.
  * [3] Koenraad M. R. Audenaert.  A sharp continuity estimate for the von Neumann entropy.  J. Phys. A, 40(28):8127–8136, 2007.
  * [4] Charles Bordenave, Pietro Caputo, and Justin Salez.  Cutoff at the “entropic time” for sparse Markov chains.  Probab. Theory Related Fields, 173(1-2):261–292, 2019.
  * [5] Charles Bordenave, Yanqi Qiu, and Yiwei Zhang.  Spectral gap of sparse bistochastic matrices with exchangeable rows.  Ann. Inst. Henri Poincaré Probab. Stat., 56(4):2971–2995, 2020\.
  * [6] Jean Bourgain and Alex Gamburd.  Uniform expansion bounds for Cayley graphs of ${\rm SL}_{2}(\mathbb{F}_{p})$.  Ann. of Math. (2), 167(2):625–642, 2008.
  * [7] Jean Bourgain, Alex Gamburd, and Peter Sarnak.  Affine linear sieve, expanders, and sum-product.  Invent. Math., 179(3):559–644, 2010.
  * [8] Jean Bourgain and Péter P. Varjú.  Expansion in $SL_{d}({\bf Z}/q{\bf Z}),\,q$ arbitrary.  Invent. Math., 188(1):151–173, 2012.
  * [9] Emmanuel Breuillard and Alex Gamburd.  Strong uniform expansion in ${\rm SL}(2,p)$.  Geom. Funct. Anal., 20(5):1201–1209, 2010.
  * [10] Emmanuel Breuillard and Péter P. Varjú.  Cut-off phenomenon for the ax+b Markov chain over a finite field, 2019\.
  * [11] Tsz Ho Chan and Igor E. Shparlinski.  On the concentration of points on modular hyperbolas and exponential curves.  Acta Arith., 142(1):59–66, 2010.
  * [12] Sourav Chatterjee and Persi Diaconis.  Speeding up Markov chains with deterministic jumps.  Probab. Theory Related Fields, 178(3-4):1193–1214, 2020.
  * [13] F. R. K. Chung, Persi Diaconis, and R. L. Graham.  Random walks arising in random number generation.  Ann. Probab., 15(3):1148–1165, 1987.
  * [14] Javier Cilleruelo and Moubariz Z. Garaev.  Concentration of points on two and three dimensional modular hyperbolas and applications.  Geom. Funct. Anal., 21(4):892–904, 2011.
  * [15] Javier Cilleruelo and Ana Zumalacárregui.  Saving the logarithmic factor in the error term estimates of some congruence problems.  Math. Z., 286(1-2):545–558, 2017.
  * [16] Persi Diaconis and Laurent Saloff-Coste.  Comparison techniques for random walk on finite groups.  Ann. Probab., 21(4):2131–2156, 1993.
  * [17] Persi Diaconis and Laurent Saloff-Coste.  Comparison theorems for reversible Markov chains.  Ann. Appl. Probab., 3(3):696–730, 1993.
  * [18] Sean Eberhard and Péter P. Varjú.  Mixing time of the Chung–Diaconis–Graham random process.  Probab. Theory Related Fields, 2020.
  * [19] M. Z. Garaev.  On the logarithmic factor in error term estimates in certain additive congruence problems.  Acta Arith., 124(1):27–39, 2006.
  * [20] A. Salehi Golsefidy and Péter P. Varjú.  Expansion in perfect groups.  Geom. Funct. Anal., 22(6):1832–1891, 2012.
  * [21] Konstantin Golubev and Amitay Kamber.  Cutoff on graphs and the sarnak-xue density of eigenvalues, 2019.
  * [22] Jimmy He.  Markov chains on finite fields with deterministic jumps, 2020.
  * [23] H. A. Helfgott.  Growth and generation in ${\rm SL}_{2}(\mathbb{Z}/p\mathbb{Z})$.  Ann. of Math. (2), 167(2):601–623, 2008.
  * [24] Jonathan Hermon, Allan Sly, and Perla Sousi.  Universality of cutoff for graphs with an added random matching, 2020\.
  * [25] Martin Hildebrand.  On the Chung-Diaconis-Graham random process.  Electron. Comm. Probab., 11:347–356, 2006.
  * [26] Martin Hildebrand.  A lower bound for the Chung-Diaconis-Graham random process.  Proc. Amer. Math. Soc., 137(4):1479–1487, 2009.
  * [27] Martin Hildebrand.  On a lower bound for the Chung-Diaconis-Graham random process.  Statist. Probab. Lett., 152:121–125, 2019.
  * [28] David A. Levin and Yuval Peres.  Markov chains and mixing times.  American Mathematical Society, Providence, RI, 2017.  Second edition, With contributions by Elizabeth L. Wilmer, With a chapter on “Coupling from the past” by James G. Propp and David B. Wilson.
  * [29] Richard Neville, III.  On lower bounds of the Chung-Diaconis-Graham random process.  ProQuest LLC, Ann Arbor, MI, 2011.  Thesis (Ph.D.)–State University of New York at Albany.
  * [30] Igor Rivin.  Zariski density and genericity.  Int. Math. Res. Not. IMRN, (19):3649–3657, 2010.
  * [31] Peter Sarnak.  Some applications of modular forms, volume 99 of Cambridge Tracts in Mathematics.  Cambridge University Press, Cambridge, 1990.
  * [32] Atle Selberg.  On the estimation of Fourier coefficients of modular forms.  In Proc. Sympos. Pure Math., Vol. VIII, pages 1–15. Amer. Math. Soc., Providence, R.I., 1965.
  * [33] I. D. Shkredov.  Modular hyperbolas and bilinear forms of Kloosterman sums.  J. Number Theory, 220:182–211, 2021.
  * [34] Aaron Smith.  Comparison theory for Markov chains on different state spaces and application to random walk on derangements.  J. Theoret. Probab., 28(4):1406–1430, 2015.
  * [35] Thomas Weigel.  On the profinite completion of arithmetic groups of split type.  In Lois d’algèbres et variétés algébriques (Colmar, 1991), volume 50 of Travaux en Cours, pages 79–101. Hermann, Paris, 1996.
  * [36] Avi Wigderson.  Mathematics and computation.  Princeton University Press, Princeton, NJ, 2019.  A theory revolutionizing technology and science.