# Towards Stable Interstellar Flight:
Levitation of a Laser-Propelled Sailcraft*

Afroza Shirin1, Edl Schamiloglu2, and Rafael Fierro2 *This work is supported
by Breakthrough Starshot Foundation LLC.1Aerospace and Mechanical Engineering,
University of Texas at El Paso, El Paso, Texas, USA. E-mail:
<EMAIL_ADDRESS>& Computer Engineering, University of New Mexico,
Albuquerque, NM, USA. E-mail: {edls<EMAIL_ADDRESS>
###### Abstract

Exploring and traveling to distant stars has long fascinated humanity but has
been limited due to the vast distances. The Breakthrough Starshot Program aims
at eliminating this limitation by traveling to Alpha Centauri, which is 4.37
light-years away. This is only possible if a vehicle travels at a substantial
fraction of the speed of light. The Breakthrough Starshot Program initiative
is to develop a proof-of-concept that is accelerating a sail to relativistic
speeds using a laser beam aimed at the sail. At this high speed, while stable
beam riding is one of the crucial concerns of this concept, the dynamic
stability analysis of a sail is hardly present in the previous literature.
Furthermore, it is important to investigate the dynamic stability in the
experiment before driving the sail to relativistic speeds. As a proof-of-
concept, we study the dynamic stability of the sail levitated at a certain
height by a laser beam. The sail’s dynamics are modeled as a rigid body whose
shape is parameterized by a sweep function. We estimate the region of
attraction (ROA) for dynamic stability analysis using Lyapunov theory and sum-
of-square (SOS) programming. The ROA confirms how many transverse and angular
perturbations a levitated sail can tolerate. We also conclude on some of the
important parameters of the sail that affects the dynamic stability.
Simulation results validate our theoretical analysis.

###### Index Terms:

Levitation, dynamic stability, Lyapunov stability, sailcraft, laser
propulsion, region of attraction, sum of squares, positive semidefinite
programming.

##  I Introduction

The Breakthrough Starshot project aims to demonstrate a proof-of-concept for
ultra-fast laser beam-driven sail and lays the foundations for the first
launch to Alpha Centauri within the next generation [1, 2, 3]. The idea is to
propel a sail, which will carry ultra-low-mass “starchips” to 20% of the speed
of light. The spacecraft would spend 20 years getting to Alpha Centauri and
another 4.25 years sending back its data so that they can be analyzed within
next generation. The stability of the sail is a key consideration in the
success of this concept and is affected by sail shape, laser beam profile,
beam jitter, mass distribution, payload, spinning dynamics, flexibility,
material properties, etc. There are several investigations on the stability of
a sail based on the sail configuration [4, 5, 6, 7, 8]. There are some other
works where the material and optical properties have been included to
investigate the stability [9, 10, 11, 12, 13].

Although these initiatives towards stability analysis of the sail are
promising, they are based on static stability analysis and carry some
limitations. First, static stability, which is the linearization around an
equilibrium, can be used to indicate stability/instability for the particular
parameters they used, but the conclusions may not apply to the parameters or
timescales other than those considered in the study. Second, the precise nano-
patterning and handling of a lightweight macroscopic substrate is a delicate
task. Also, the loading of the object in the laser beam will be difficult as
it will require the system to be released in a very controlled way.

As the sail itself does not have any propeller or control and the laser will
be ground based or in a low earth orbit, the perturbations in translational
and angular positions with respect to the laser beam are very critical and
measuring the maneuverability is very important. This analysis can be done by
using dynamic stability or nonlinear Lyapunov stability [14], which has barely
made an appearance in the current sail literature. One of the aspects of
dynamic stability analysis is to estimate the region of attraction (ROA). Any
perturbation of the states within the ROA will converge to the equilibrium
point. The ROA can be used to quantify the maneuverability with respect to the
laser beam. However, estimation of the ROA of an equilibrium point for
uncertain systems is itself a challenging problem that requires the
development of computationally efficient methods. To address the dynamic
stability analysis of a levitated sail, we estimate the ROA. In this article,
we consider the sail dynamics as a 3D rigid body. We parameterize the sail
shape. Then we use Lyapunov theory and sum-of-squares (SOS) programming to
estimate the ROA [15]. Finally, we analyze how the ROA is affected as a sail’s
parameter changes.

##  II Parameterized Sail

The whole vehicle consists of a sail, a mast, and a payload. The payload is
assumed to be attached by the mast to the sail. To generate the sail surface a
parameterized sweep curve in the x-z plane is rotated about the z-axis [16].
In this article, two types of sail shapes are considered: 1) conical sail and
2) spherical-cap sail. The conical sail is studied in [4, 7] and a spherical-
cap sail in [6]. For the conical sail, the following curve is used

$g_{s}(x)=c_{0}+c_{1}\left(\frac{x}{R}\right)+c_{2}\left(\frac{x}{R}\right)^{2}+c_{3}\left(\frac{x}{R}\right)^{3}+c_{4}\left(\frac{x}{R}\right)^{4}.$
(1)

Here $c_{0},c_{1},c_{2},c_{3}$, and $c_{4}$ are coefficients of the
polynomials and define different shapes of the sail with base radius $R$. A
conical sail with base radius $R$ and cone angle $\alpha$ can be obtained by
setting $c_{0}=R\tan(\alpha),c_{1}=-R\tan(\alpha),c_{2}=c_{3}=c_{4}=0$ in Eq.
(1). Here the angle $\alpha$ is between the conical sail surface and the
horizontal plane. For the conical sail, the cone angle $\alpha$ is considered
an important parameter and it is analyzed how the dynamical stability is
affected as the parameter $\alpha$ changes. To generate a spherical-cap sail,
the following function is used

$g_{s}(x)=\sqrt{R^{2}-x^{2}}-\sqrt{R^{2}-a^{2}}.$ (2)

Notice that the mast length is an important parameter to analyze the dynamic
stability. A conical sail and a spherical-cap sail are presented in Fig 1.

Figure 1: a) A conical sail with base radius $R=1$, $h=1$ and
$\alpha=\frac{\pi}{4}$. b) A portion of a spherical-cap sail with base radius
$a=0.5$, $h=0.13$ and with a radius of curvature $R=1$.

##  III Modeling the Dynamics of the Sail

For mathematical modeling, we assume that the whole vehicle, which includes a
sail, a mast, and a payload, is a rigid body. The sail shape is chosen from
Section II, the length of the mast is $l$, and the total mass of the vehicle
is $m$. The sail is fully reflective and not a multi-reflector 111multi-
reflector: A surface that reflects the light more than once.. The vehicle
motion is represented as a translation of the vehicle’s center of mass (CM)
and a rotation of the body frame ($x_{b}-y_{b}-z_{b}$) with respect to the
inertial frame ($x-y-z$). The vehicle’s CM position is described by its
inertial Cartesian coordinates, $\textbf{r}=(x,y,z)$, and its orientation with
respect to the inertial reference frame by three Euler angles
$\boldsymbol{\alpha}_{E}=(\psi,\theta,\phi)$. The convention of rotations is
considered as $x-y^{\prime}-z^{\prime\prime}$ by following the sequence: (1)
rotation about the x-axis by angle $\phi$; (2) rotation about the new position
of the y-axis by angle $\theta$; (3) rotation about the new position of the
z-axis by angle $\psi$. Now the vehicle dynamics is described by 12 ODEs [17]:

$\displaystyle\dot{\textbf{r}}$ $\displaystyle=\textbf{v},$ (3a)
$\displaystyle\dot{\boldsymbol{\alpha}}_{E}$
$\displaystyle=L^{I}_{B}\boldsymbol{\omega}_{b},$ (3b)
$\displaystyle\dot{\textbf{v}}$
$\displaystyle=\frac{\textbf{F}}{M}+\textbf{g},$ (3c)
$\displaystyle\dot{\boldsymbol{\omega}}_{b}$
$\displaystyle=J^{-1}[-\boldsymbol{\omega}_{b}^{X}\times
J\boldsymbol{\omega}_{b}+\boldsymbol{\tau}_{b}],$ (3d)

where v is the inertial velocity of the vehicle, $\boldsymbol{\omega}_{b}^{X}$
is the skew-symmetric matrix associated with the vehicle’s angular velocity
vector, $\boldsymbol{\omega}_{b}$ is the angular velocity of the vehicle in
the inertial frame and expressed in the body-frame, $J$ is the inertia matrix
with respect to the vehicle’s CM, $\textbf{g}=[0,0,-9.8]$ m/s is the gravity
vector, $M$ is the total mass of the vehicle. In this formulation, we assume
that the origin of the body axis moves along the inertial $z$-axis. The force
F and the torque $\boldsymbol{\tau}_{b}$ are the net external force and torque
acting on the body, respectively. When the laser beam is radiated onto the
sail, the force and the torque are given by

$\displaystyle\textbf{F}(\textbf{x},\Theta)=\left(R_{I}^{B}\right)^{-1}$
$\displaystyle\left[\iint_{S(\Theta_{1})}2\frac{P(\textbf{x},\Theta_{2})\hat{\textbf{b}}\cdot\hat{\textbf{n}}(\textbf{x})}{c}\hat{\textbf{n}_{b}}dS\right],$
(4)
$\displaystyle\quad\hat{\textbf{n}}(\textbf{x})=\left(R_{I}^{B}\right)^{-1}\hat{\textbf{n}_{b}}.$

$\boldsymbol{\tau}_{b}(\textbf{x},\Theta)=\iint_{S(\Theta_{1})}2\frac{P(\textbf{x},\Theta_{2})\hat{\textbf{b}}\cdot\hat{\textbf{n}}(\textbf{x})}{c}(\textbf{r}_{b}(\textbf{x})\times\hat{\textbf{n}}_{b})dS,$
(5)

where $c$ is the speed of light, $S$ is the domain of integration of the
surface of the sail, $\hat{\textbf{n}}_{b}$ is the unit vector normal to the
sail surface at the generic point on the sail, $\hat{\textbf{b}}$ is a unit
vector parallel to the beam axis, $\textbf{r}_{b}$ is the position vector from
the vehicle CM to the point x in the body frame. Here
$\Theta=\\{\Theta_{1},\Theta_{2}\\}$, $\Theta_{1}$ is the set of parameters
used to describe the sail geometry, and $\Theta_{2}$ is the set of parameters
to express the beam profile. The laser beam is assumed to have a radially
symmetric Gaussian power distribution with full-width-at-half-maximum (FWHM)
$w=2\sigma\sqrt{2\ln(2)}$, where $\sigma$ is the standard deviation of the
Gaussian beam. The beam power flux $P$ at the point x is given as

$P(\textbf{x},\Theta_{2})=\frac{P_{0}}{2\pi\sigma^{2}}e^{-\frac{(x^{2}+y^{2})}{2\sigma^{2}}},$
(6)

where $P_{0}$ (in GW) is the total beam power. In Eq. (3b), the matrix
$L_{B}^{I}$ relates the time derivative of the Euler angles to the body frame
components of the angular velocity and derived as [16]. In Eq. (4),
$R_{B}^{I}$ is the rotation matrix associated with the considered rotation
convention [16].

We now assume the material of the sail has a damping effect. Common qualified
dampers may include viscous and magnetic induction devices and physical
implementation can be obtained [18]. These effects will improve the stability
performance of the beam-driven sail [19]. We model the damping effect in the
equations of motion, and the dynamics in (3c) and (3d) can be modified as

$\displaystyle\dot{\textbf{v}}$
$\displaystyle=\frac{\textbf{F}}{M}+\textbf{g}-\frac{D}{M}\textbf{v}.$

The matrix $D$ captures the damping effect and satisfies the condition:
$Dij=-Dji$. In Fig. 2, a schematic illustration of a sail riding on a Gaussian
beam is presented.

Figure 2: Schematic illustration of a sail riding on a Gaussian laser beam.

##  IV Levitation of the Laser Beam-driven Sail

###  IV-A Hypothetical Setup of the Experiment

For an experimental setup, we assume that a sail is placed in a vacuum chamber
and is driven by a Gaussian laser beam as discussed previously. We also assume
that the velocometer acts like a sensor and can send the feedback signal to
the actuator (laser) to produce the force from the beam to levitate the sail.
A hypothetical representation of the experiment is shown in Fig. 3.

Figure 3: Schematic of a sub-scale experimental laser-driven sail. Credit:
Breakthrough Foundation.

###  IV-B Lyapunov Based Design

It turns out that we need to design a nonlinear closed-loop control law that
not only can levitate the sail at an assigned height but also guarantees the
closed-loop stability of the sail. In this case, the force directly drives the
states $z$ and $v_{z}$ and the other states indirectly. Thus, the states $z$
and $v_{z}$ are considered as actuated states. The dynamics of the actuated
system can be written as

$\displaystyle\dot{z}$ $\displaystyle=v_{z},$ (8a) $\displaystyle\dot{v}_{z}$
$\displaystyle=g+\frac{F_{z}}{M},$ (8b)

where $F_{z}$ is the $z$ component of the force F given in Eq. (4). We can
rewrite Eq. (8) as,

$\begin{bmatrix}\dot{z}\\\
\dot{v}_{z}\end{bmatrix}=\left[\begin{array}[]{c}v_{z}\\\
g\end{array}\right]+\begin{bmatrix}0\\\ \frac{G_{z}}{M}\end{bmatrix}u,$ (9)

where the vector quantity,

G
$\displaystyle=\frac{1}{c\pi\sigma^{2}}(R_{I}^{B})^{-1}\iint_{S(\Theta_{1})}e^{-\frac{(x^{2}+y^{2})}{2\sigma^{2}}}\hat{\textbf{b}}.\hat{\textbf{n}}(\textbf{x})\textbf{dS},$
$\displaystyle u$ $\displaystyle=P_{0},$

and $G_{z}$ is the $z$ component of G. Note $F_{x}=G_{x}$, $F_{y}=G_{y}$, and
$F_{z}=G_{z}u$. The actuator only can produce force along the $z$-direction
but not any applied torque. But the torque due to the force F is still
considered. If there is any angular perturbation, i.e., $(R_{I}^{B})^{-1}\neq
I_{3}$, non-zero $x$ and $y$ components of F will be produced. Moreover, we
assume that spinning of the sail around the $z$-axis is allowed and will not
affect the dynamic stability, but the spinning around $x$ and $y$ will as
these spinnings will tumble the sail. We separate the actuated and
underactuated states as follows

$\displaystyle\textbf{x}_{a}$
$\displaystyle=\begin{bmatrix}z&v_{z}\end{bmatrix}^{T},$
$\displaystyle\textbf{x}_{u}$
$\displaystyle=\begin{bmatrix}x&y&\psi&\theta&\phi&v_{x}&v_{y}&\omega_{x}&\omega_{y}&\omega_{z}\end{bmatrix}^{T}.$

Consider the Lyapunov function for the actuated states to design a closed-loop
controller

$V_{a}(\textbf{x}_{a})=\frac{1}{2}z^{2}+\frac{1}{2}v_{z}^{2}$ (11)

of which the time derivative is

$\displaystyle\dot{V}_{a}$
$\displaystyle=zv_{z}+v_{z}\left(g+\frac{G}{M}u\right).$ (12)

The control input

$u=P_{0}=\frac{-M}{G}\left[g+z+v_{z}\right]$ (13)

leads to $\dot{V}_{a}=-v_{z}^{2}<0$. If we want to levitate the sail at
$z=z_{d}$, the transformation $e=z-z_{d}$ accommodates the case and we obtain

$u=P_{0}=\frac{-M}{G}[g+e+v_{z}].$ (14)

This feedback control law $u$ levitates the sail at $z_{d}$. Now we need to
analyze the local stability of the internal dynamics. The internal dynamics
can be written as

$\begin{bmatrix}\dot{x}\\\ \dot{y}\\\ \dot{\boldsymbol{\alpha}}_{E}\\\
\dot{v}_{x}\\\ \dot{v}_{y}\\\
\dot{\boldsymbol{\omega}}_{b}\end{bmatrix}=\underbrace{\begin{bmatrix}v_{x}\\\
v_{y}\\\ L^{I}_{B}\boldsymbol{\omega}_{b}\\\ 0\\\ 0\\\
J^{-1}[-\boldsymbol{\omega}_{b}^{X}\times
J\boldsymbol{\omega}_{b}]\end{bmatrix}}_{\bar{\textbf{f}}}+\underbrace{\begin{bmatrix}0\\\
0\\\ 0\\\ \frac{F_{x}}{M}\\\ \frac{F_{y}}{M}\\\
J^{-1}\boldsymbol{\tau}_{b}\end{bmatrix}}_{\bar{\textbf{g}}}.$ (15)

The states’ associated with the internal dynamics can be written as,

$\textbf{x}_{u}=\begin{bmatrix}x&y&\psi&\theta&\phi&v_{x}&v_{y}&\omega_{x}&\omega_{y}&\omega_{z}\end{bmatrix}^{T}.$

The internal dynamics can be considered as a perturbation of a nominal system
$\dot{\textbf{x}}_{u}(t)=\bar{\textbf{f}}(\textbf{x}_{u}(t))$ and the
perturbed system is given in the form,

$\dot{\textbf{x}}_{u}(t)=\bar{\textbf{f}}(\textbf{x}_{u}(t))+\bar{\textbf{g}}(\textbf{x}_{u}).$
(16)

With the control input $u$ in Eq. (13), the full dynamics is (locally)
exponentially stable. We assume that $\bar{\textbf{g}}(\textbf{x}_{u})$ is
bounded on some set $\Gamma$. By considering the stability of perturbed
systems, more specifically vanishing perturbation, [14, Section 9.1], the
internal dynamics are stable and there exists a Lyapunov function
$V(\textbf{x}_{u}(t))$ in $[0,\inf)\times\Gamma$, where
$\Gamma=\\{\textbf{x}_{u}\in\mathbb{R}^{10}:\|\textbf{x}_{u}\|<\bar{c}\\}$
such that $\dot{V}(\textbf{x}_{u})\leq 0$. Here $\bar{c}$ is a positive
constant. Such a set $\Gamma$ suggests that there exists an invariant subset
of $\Gamma$ which can be estimated as a ROA described in the following
subsection. Let

$V(\textbf{x}_{u})=\textbf{x}_{u}^{T}P\textbf{x}_{u}(t)$ (17)

be a Local Lyapunov Function (LLF) for the system in Eq. (15), where $P$ is a
$10\times 10$ positive definite matrix. This LLF is used to compute the ROA to
study the dynamic stability of the sail.

##  V ROA Estimation

To address the dynamic stability of the sail, it is sufficient to analyze the
dynamic stability of the internal dynamics. Specifically, we estimate the ROA
of the internal dynamics by using Lyapunov theory and SOS programming. In SOS
programming, the SOS decomposition technique is used to convert the cost and
constraints such that the ultimate optimization problem becomes a Semidefinite
Programming Problem (SDP) [15]. In [20], it has been proved that locally
stable polynomial systems admit polynomial LLFs on compact sets. For
polynomial systems and polynomial LLFs, one can maximize the size of the ROA
using SOS programming (e.g., [21, 22]). Here the procedure to estimate the ROA
is presented for a general nonlinear system of the form
$\dot{\textbf{x}}(t)=\textbf{f}(\textbf{x})$.

###  V-A Procedure to Estimate the ROA

Consider that $V(\textbf{x})$ is an LLF for the equilibrium
$\textbf{x}=\textbf{0}$ of a nonlinear system
$\dot{\textbf{x}}(t)=\textbf{f}(\textbf{x})$. If $\textbf{x}=\textbf{0}$ is
asymptotically stable such that the following condition holds

$V(\textbf{x})\leq\rho\implies\dot{V}(\textbf{x})<0,$ (18)

then the $\rho$ sublevel set of $V(\textbf{x})$ is an inner subset of the ROA.
This condition can be solved by using the S-procedure [23, Chapter 2.6.3]. For
a general polynomial, S-procedure is stated by the following lemma.

######  Lemma 1.

Let $q_{0},q_{1},\cdots,q_{m}$ be quadratic functions on
$\textbf{x}\in\mathbb{R}^{n}$. If it is required that

$q_{0}(\textbf{x})\geq 0\text{ such that }q_{i}(\textbf{x})\geq
0,i=1,2,\cdots,m$ (19)

then there exist polynomials $s_{i}(\textbf{x})\geq 0,i=1,\cdots,m$ such that

$q_{0}(\textbf{x})-\sum_{i=1}^{i=m}s_{i}(\textbf{x})q_{i}(\textbf{x})\geq 0.$
(20)

Here $s(\textbf{x})$ is a multiplier polynomial.

By using the above lemma, a sufficient condition is,

$p(\textbf{x})(V(\textbf{x})-\rho)+s(\textbf{x})\dot{V}(\textbf{x})\geq 0$
(21)

to hold the condition in Eq. (18) to be true. By replacing the non-negativity
condition in Eq. (21) with an SOS constraint, the following SOS optimization
problem can be formulated to estimate the ROA

$\displaystyle\max_{\rho\in\mathbb{R},s\in\mathcal{S}}\;$ $\displaystyle\rho$
(22a) subject to $\displaystyle
p(\textbf{x})(V(\textbf{x})-\rho)+s(\textbf{x})\dot{V}(\textbf{x})\in SOS.$
(22b)

Here $\mathcal{S}$ is a given subspace of polynomials in Eq. (22a), for
example, one consisting of all quadratic polynomials or all quartic
polynomials, $p(\textbf{x})$ is a deterministic positive definite polynomial,
for example, $p(\textbf{x})=\textbf{x}^{T}\textbf{x}$. The optimization
problem in Eq. (22) is an SDP and can be solved numerically. For polynomial
systems and polynomial LLF one can maximize the invariant subset of the ROA by
using SOS programming. Therefore, the step-by-step procedure is as follows,

  * •

Express the nonlinear system $\dot{\textbf{x}}(t)=\textbf{f}(\textbf{x})$ in
the form of polynomials by series expansion. Although the dynamical system
might not be directly expressed in the form of polynomials, it can be
approximated by series expansion [22].

  * •

Find an LLF $V(x)$ for the equilibrium point $\textbf{x}=\textbf{0}$ of the
nonlinear system. An LLF can be found by using Eq. (17) and then compute
$\dot{V}=\frac{\partial V}{\partial\textbf{x}}\dot{\textbf{x}}$.

  * •

Solve the SOS optimization problem in Eq. (22).

###  V-B SOS Programming Solver

A software package using the Systems Polynomial Optimization Toolbox (SPOT)
[24] features is used to set up the polynomial algebra and the large-scale SOS
programming arising from our problem. We use MOSEK as the optimization solver.
A complete implementation is available online
111https://github.com/anirudhamajumdar/spotless/tree/spotless_isos.

##  VI Numerical Results

###  VI-A Parameter Setting for the Dynamics

For simulations, we consider the proposed Starshot design parameters: the beam
power $P_{0}$ is measured in GW, the FWHM of the beam $w$ is set to the base
radius of the sail, a sail with mass $m_{s}=10$ g, the mass of the payload
$m_{p}=10$ g (if any, otherwise, $m_{p}=0$). The integrals in Eqs. (4) and (5)
were approximated by discretizing the beam into a grid of 100 × 100 rays. The
path of each ray is then traced as it intersects the sail and is reflected off
of its surface. The net change in momentum of each ray is calculated, and the
resulting forces and torques are applied to the sail. For matrix $D$ we set
$D_{11}=D_{12}=D_{22}=0.01$, $D_{13}=D_{23}=D_{33}=0.0$.

###  VI-B Numerical Computation of the ROA

In this section, we numerically compute the ROA of the sail and conduct its
dynamic stability analysis. For the ROA computation, first, express the
internal dynamics in Eq. (15) in the form of polynomials by Taylor expansion
of degree 3. We find a LLF $V(\textbf{x}_{u})$ as we discussed previously.
Then, compute $\dot{V}=\frac{\partial V}{\partial\textbf{x}}\dot{\textbf{x}}$.
The ROA can be approximated by solving the optimization problem in Eq. (22).
For simulations, we choose $\mathcal{S}$ as a set of quadratic polynomials.

####  VI-B1 The ROA of a Conical Sail

We consider a conical sail with $c_{0}=R\tan(\alpha)$, $c_{1}=-R\tan(\alpha)$,
$c_{2}=c_{3}=c_{4}=0$ in Eq. (1). The base radius is set to $R=1$, the mast
length $l=2$ m, and the levitation height $z_{d}=10$ m. We conduct the
stability analysis with respect to the cone angle $\alpha$. It is found that
the $\alpha$ is one of the crucial parameters for the conical sail [4]. We
compute the ROA for $\alpha=40^{\circ}$ and $\alpha=45^{\circ}$ and plot the
projected regions in the $x-y$, $x-\phi$, $y-\theta$, and $\theta-\phi$
planes. The other projections of the ROA are omitted due to symmetry or
redundancy. In Fig. 4(a)-(d), we plot the projected region of the ROA on the
$x-y$, $x-\phi$, $\theta-\phi$, and $y-\theta$ planes, respectively. The
projected regions are shaded with colors.

Figure 4: The ROA of the conical sail projected in different planes (a)-(d).
The pink region is associated with the cone angle $\alpha=40^{\circ}$ and the
light blue region is associated with $\alpha=45^{\circ}$. The projected region
in the (a) $x-y$ plane, (b) $x-\phi$ plane, (c) $\theta-\phi$ plane, and (d)
$y-\theta$ plane. Figure 5: State histories for the conical sail under
actuation. The cone angle of the sail is $\alpha=40^{\circ}$ and the base
radius $R=1$ m. The initial condition is chosen such that the roll and pitch
angle perturbations are about $10^{\circ}$ (0.1745 rad).

The sail with $\alpha=40^{\circ}$ can tolerate transverse perturbations of at
least 8% of the base radius $R$ (Fig. 4(a)), and roll and pitch perturbations
of about $10^{\circ}$ (0.17 rad) (Fig. 4(c)). The same conical sail with
$\alpha=45^{\circ}$ can tolerate transverse perturbations of at least 9% of
the base radius, and roll and pitch perturbations of about $8^{\circ}$ (0.13
rad). Upon comparison with the regions, we observe that the ROA of the conical
sail with $\alpha=40^{\circ}$ can tolerate more pitch and roll perturbations
than the other one. To check the dynamic stability of the conical sail, the
state trajectories of the full system in Eq. (3) are computed with the initial
roll and pitch angles perturbed about $10^{\circ}$ (0.17 rad). In Fig. 5, the
state time histories of the dynamics with closed-loop control $u$ are plotted
for the conical sail with $\alpha=40^{\circ}$. The sail is levitated at 10 m.
Also, in each case, the states can tolerate the initial perturbations and
becomes steady after a certain amount of time.

Figure 6: The ROA of the spherical-cap sail projected in different planes
(a)-(d). The pink region is associated with mast length $l=2.3$ m and the
light blue region is associated with $l=2.5$. The projected region in the (a)
$x-y$ plane, (b) $x-\phi$ plane, (c) $\theta-\phi$ plane, and (d) $y-\theta$
plane. Figure 7: State histories for the spherical-cap sail under actuation.
The length of the payload is $l=2.3$ m and the base radius $R=1$ m. The
initial condition is chosen such that the roll and pitch angle perturbations
are about $10^{\circ}$ (0.1745 rad).

####  VI-B2 The ROA of a Spherical-cap Sail

We next consider a part of a spherical sail with base radius $a=0.5$ m and
with a radius of curvature $R=1$ m in Eq. (2). The levitated height is set at
$10$ m. Now we analyze the stability of the spherical-cap sail. The
interesting parameter, in this case, is the mast length $l$. In [8], it has
been shown that for a spherical-cap sail with base radius $a$ and a radius of
curvature $R$, the condition to be stable is $l>2R$. Now we compute the ROA
for $l=2.3$ m and $l=2.5$ m and plot the projected region on different planes
as we did for the conical sail. It is important to mention that the spherical-
cap sail with $l=2.1$ m is hardly stable. In Fig. 6(a)-(d), we plot the
projected region in the $x-y$, $x-\phi$, $\theta-\phi$, and $y-\theta$ planes,
respectively. The sail with $l=2.3$ m can tolerate transverse perturbations of
at least $20\%$ of the sail radius and roll and pitch perturbation of at least
$6^{\circ}$ (0.1 rad). The same sail with $l=2.5$ m can tolerate transverse
perturbations of 10% and roll and pitch perturbations of $5^{\circ}$. Upon
comparison with two shaded regions, it is observed that the spherical-cap sail
with $l=2.3$ m can tolerate more transverse and angular perturbations than the
sail with $l=2.5$ m. To check the dynamic stability of the spherical-cap sail,
the state trajectories of the full system in Eq. (3) are computed with the
initial roll and pitch angles perturbed of about $10^{\circ}$ (0.17 rad). For
this computation, the differential equations of motion described in Section
III are integrated using the standard fourth-order Runge–Kutta method. In Fig.
7, the state time histories of the dynamics with closed-loop control are
plotted for the spherical-cap sail with $l=2.3$ m. The sail is levitated at 10
m. Also, in each case, the states can tolerate the initial perturbations and
become steady after a certain amount of time.

##  VII Discussion

The dynamics of the beam-driven sail are modeled as a rigid body. The sail
shape parameterized as a sweep function has been presented to define the sail
shape. A Lyapunov-based closed-loop controller which levitates a sail at an
assigned height, is designed. The dynamic stability analysis of a levitated
sail is conducted through the estimation of the ROA. The ROA is computed using
Lyapunov theory and SOS. Moreover, the ROA indicates how much perturbation a
sail can tolerate. A full simulation of the sail system with the control law
confirms these regions. The dynamic stability analysis conveys that the system
with closed-loop control is not only levitated, but also confirms its dynamic
stability. It is also noticeable that with the control law, the sail is not
only levitated but also stable even for large perturbations of the states.
Upon comparison with the ROA for the conical sail with $\alpha=40^{\circ}$ and
$45^{\circ}$, this analysis suggests that the conical sail with angle
$\alpha=40^{\circ}$ performs better than the sail with $\alpha=45^{\circ}$. A
similar analysis is also performed for the spherical-cap sail for which the
parameter of interest is the mast length $l$. The ROA for the spherical-cap
sail with mast lengths $l=2.3$ m and $2.5$ m. This dynamic stability analysis
indicates that the spherical-cap sail with $l=2.3$ m can tolerate more
perturbations than that of $l=2.5$ m. By comparing the ROA of a conical sail
with a spherical-cap sail, we conclude that in general, a conical sail can
tolerate more angular perturbations while the spherical-cap sail can tolerate
more transverse perturbations.

## References

  * [1] A. Finkbeiner, “Near-light-speed mission to alpha centauri,” _Scientific American_ , vol. 316, no. 3, pp. 30–37, 2017.
  * [2] P. Daukantas, “Breakthrough starshot,” _Optics and Photonics News_ , vol. 28, no. 5, pp. 26–33, 2017.
  * [3] K. L. Parkin, “The Breakthrough Starshot system model,” _Acta Astronautica_ , vol. 152, pp. 370–384, Nov 2018.
  * [4] E. Chahine, C. Abdallah, E. Schamiloglu, and D. Georgiev, “Dynamics and Control of Microwave-propelled Sails,” in _AIAA Space 2003 Conference & Exposition_. Reston, Virginia: American Institute of Aeronautics and Astronautics, Sep 2003.
  * [5] J. Benford, G. Benford, O. Gornostaeva, E. Garate, M. Anderson, A. Prichard, and H. Harris, “Experimental tests of beam-riding sail dynamics,” in _AIP Conference Proceedings_ , vol. 608, 2002, pp. 457–461.
  * [6] E. Popova, M. Efendiev, and I. Gabitov, “On the stability of a space vehicle riding on an intense laser beam,” _Mathematical Methods in the Applied Sciences_ , vol. 40, no. 4, pp. 1346–1354, March 2017.
  * [7] Z. Manchester and A. Loeb, “Stability of a Light Sail Riding on a Laser Beam,” _The Astrophysical Journal_ , vol. 837, no. 2, p. L20, March 2017\.
  * [8] G. Genta and D. Riccobono, “Preliminary analysis of the deformed shape of an inflated light sail,” in _5th Tennessee Valley Interstellar Workshop (TVIW)_ , Huntsville, AL, Oct 4-6 2017.
  * [9] J. Siegel, A. Y. Wang, S. G. Menabde, M. A. Kats, M. S. Jang, and V. W. Brar, “Self-Stabilizing Laser Sails Based on Optical Metasurfaces,” _ACS Photonics_ , vol. 6, no. 8, pp. 2032–2040, Aug 2019.
  * [10] O. Ilic and H. A. Atwater, “Self-stabilizing photonic levitation and propulsion of nanostructured macroscopic objects,” _Nature Photonics_ , vol. 13, no. 4, pp. 289–295, April 2019.
  * [11] Y.-J. L. Chu, M. Meem, P. R. Srivastava, R. Menon, and G. A. Swartzlander, “Parametric control of a diffractive axicon beam rider,” _Optics Letters_ , vol. 46, no. 20, pp. 5141–5144, 2021.
  * [12] M. M. Salary and H. Mosallaei, “Photonic metasurfaces as relativistic light sails for doppler-broadened stable beam-riding and radiative cooling,” _Laser & Photonics Reviews_, vol. 14, no. 8, p. 1900311, 2020.
  * [13] R. Gao, M. D. Kelzenberg, Y. Kim, A. Merkt, O. Ilic, and H. A. Atwater, “Restoring optical forces and torques for passive stabilization of microscopic lightsails,” in _Optical Trapping and Optical Micromanipulation XVIII_ , vol. 11798. International Society for Optics and Photonics, 2021, p. 117980Y.
  * [14] H. K. Khalil, _Nonlinear Systems_ , 3rd ed. Prentice hall, 2002.
  * [15] P. A. Parrilo, “Structured semidefinite programs and semialgebraic geometry methods in robustness and optimization,” Ph.D. dissertation, California Institute of Technology, 2000.
  * [16] A. Shirin, E. Schamiloglu, C. Sultan, Y. Yang, J. Benford, and R. Fierro, “Modeling and stability of a laser beam-driven sail,” in _2021 American Control Conference (ACC)_. IEEE, 2021, pp. 4269–4275.
  * [17] P. C. Hughes, _Spacecraft attitude dynamics_. Courier Corporation, 2012.
  * [18] M. Rafat, H. R. Dullin, B. T. Kuhlmey, A. Tuniz, H. Luo, D. Roy, S. Skinner, T. J. Alexander, M. S. Wheatland, and C. M. de Sterke, “Self-stabilization of light sails by damped internal degrees of freedom,” _Physical Review Applied_ , vol. 17, no. 2, p. 024016, 2022.
  * [19] P. R. Srivastava, Y.-J. L. Chu, and G. A. Swartzlander, “Stable diffractive beam rider,” _Optics Letters_ , vol. 44, no. 12, p. 3082, June 2019.
  * [20] M. M. Peet, “Exponentially stable nonlinear systems have polynomial Lyapunov functions on bounded regions,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 5, pp. 979–987, 2009.
  * [21] S. Prajna, A. Papachristodoulou, P. Seiler, and P. A. Parrilo, “SOSTOOLS and its control applications,” in _Positive polynomials in control_. Springer, 2005, pp. 273–292.
  * [22] A. Majumdar, A. A. Ahmadi, and R. Tedrake, “Control and verification of high-dimensional systems with DSOS and SDSOS programming,” in _53rd IEEE Conference on Decision and Control_ , 2014, pp. 394–401.
  * [23] S. Boyd, L. El Ghaoui, E. Feron, and V. Balakrishnan, _Linear matrix inequalities in system and control theory_. SIAM, 1994.
  * [24] A. Megretski, “Systems polynomial optimization tools (SPOT),” _Massachusetts Inst. Technol., Cambridge, MA, USA_ , 2010.


# Manipulating Sparse Double Descent

Ya Shi Zhang
Department of Pure Mathematics and Mathematical Statistics
University of Cambridge
<EMAIL_ADDRESS>
https://yashizhang.github.io/

###### Abstract

This paper investigates the double descent phenomenon in two-layer neural
networks, focusing on the role of L1 regularization and representation
dimensions. It explores an alternative double descent phenomenon, named
’sparse double descent’. The study emphasizes the complex relationship between
model complexity, sparsity, and generalization, and suggests further research
into more diverse models and datasets. The findings contribute to a deeper
understanding of neural network training and optimization. The code is
available at https://github.com/yashizhang/sparsedoubledescent.

##  1 Introduction

In the modern era of deep learning, it is commonplace for practitioners to
specify models that have _many more_ parameters than the number of training
data. This is due to a statistically non-intuitive phenomenon known as _double
descent_ (Belkin et al., 2019). In a first course on regression, we are often
taught the dangers of _overfitting_ , when the model fits to the noise in the
data and severely deteriorates performance across the entire population. For
example, consider polynomial regression, where the degree of the polynomial is
significantly higher than the number of training data points. This leads to
the notion of a trade-off between the bias and variance of the model, inducing
a U-shaped curve (depicted in fig. 1) where the vertical axis depicts the
population risk and the horizontal axis depicts the ratio between the number
of model parameters and number of data points.

Recently, Curth et al. (2023) has attempted to bridge the gap of our
understanding on the prevalence of the double descent phenomenon in non-deep
learning models, such as trees, boosting, and linear regression — methods
known as _smoothers_. Notably, analysis of these model classes affirmed the
wisdom that model complexity and capacity is not necessarily an affine
function of the number of model parameters. Rather, there are higher
dimensional substructures in the set of model parameters such that increasing
model parameters in certain directions do not cause overfitting.

In this paper, we examine the possibility of controlling the double descent
phenomenon on simple two layer neural networks with varying representation
dimension. We view the first layer as a learned approximation of the ‘ground
truth kernel,’ and the second layer as a linear classifier. Particularly, we
extend Curth et al. (2023) by first considering LASSO regressors, then
implementing LASSO via $L_{1}$ regularization in neural network training as a
sparsifier. This is motivated by the perspective of the $L_{1}$ regularizer as
a convex surrogate for the “$L_{0}$” objective, aiming to prune down the
number of parameters active in the network. Benefits of pruning and/or
sparsification include (1) computational and storage savings, (2) model
capacity regularization motivated by Occam’s razor (LeCun et al., 1989), and
(3) potential lottery ticket subspaces (Frankle and Carbin, 2019).

Furthermore, the prevalence of a ‘sparse double descent’ first noted by He et
al. (2022) — as model sparsity increases, test performance first decreases
then increases (before finally decreasing as sparsity approaches 100%) — opens
the possibility of model sparsity also exhibiting many intrinsic complexity
axes.

Figure 1: Simple depiction of the double descent phenomenon. $d$ refers to the
number of the parameters in the model and $n$ refers to the number of training
data points.

##  2 Background

Consider the supervised learning problem, we are given training data
$(x_{i},y_{i})_{i=1}^{n}$. In a family of functions
$\mathcal{F}:=\\{f_{\theta}(\cdot)\mid\theta\in\mathbb{R}^{d}\\}$
parameterized by $\theta\in\mathbb{R}^{d}$, we would like to choose a
$\theta^{*}$ such that the function induced $f_{\theta^{*}}(\cdot)$ minimizes
some loss function over the training data, usually of the form
$\mathcal{L}(\theta):=\sum_{i=1}^{n}\ell(f_{\theta}(x_{i}),y_{i})$. Depending
on whether the task is classification or regression, common choices for
$\ell(\cdot,\cdot)$ would be the cross entropy loss or quadratic loss,
respectively.

We now modify our loss function by adding a _regularizer_ , typically of the
form $\alpha||\theta||_{p}$ for some $p\in(0,\infty]$ and
$\alpha\in\mathbb{R}^{+}$. In this paper we focus exclusively on the case of
$p=1$. This $L_{1}$ penalty — inspired by the LASSO regressor in linear
regression — inductively biases the network to learn parametrizations that are
sparse and well-generalizing. The former is due to the $L_{1}$ penalty being
the projection of the $L_{0}$ penalty (counting the number of non-zero
elements in a vector) to the space of convex functions (Ramirez et al., 2013).

If our family of functions are neural networks with $L>1$ layers, then we can
view the first $L-1$ layers as a learnable approximation to some unknown,
‘ground truth’ kernel. Through this perspective, we can view the training of
neural networks as first learning an appropriate kernel function, then
performing linear regression on the kernelized inputs. This is theoretically
supported in the lazy training regime (Atanasov et al., 2022; Geiger et al.,
2020; Jacot et al., 2021), and empirically supported via many observed
phenomena such as the prevalence of Neural Collapse observed in Papyan et al.
(2020).

##  3 Experiments

For our experiments, we will be training a two-layer multilayer perceptron
with ReLU non-linearity on the MNIST dataset (LeCun et al., 2010) using
stochastic gradient descent. The empirical experiments comprise of modifying
two components of the training procedure to observe a few interesting
phenomena. The first is the regularization coefficient $\alpha$. As we
increase $\alpha$, we place heavier emphasis on sparsity when training the
neural network. The second is the intermediate layers’ dimension, which we
refer to as the _kernel dimension_. By varying this, we are essentially
changing the learned kernel’s intrinsic dimension.

###  3.1 Results

As seen in table 1, for each kernel dimension, we train separate neural
networks with varying $\alpha$. In the top right corner, we have set the
middle layer of our two-layer multi-layer perceptron to have $5$ neurons. In
the top middle, we have 10 neurons. Following this, we also show results for
25, 50, 75, and 100 intermediate neurons.

We vary the coefficient of the $L_{1}$ regularizer to first observe, for all
choices of number of intermediate neurons, a version of the sparse double
descent phenomenon. We also observe the existence of many ascents and descents
as we decrease the number of intermediate neurons. Finally, we see that the
location of the minima is invariant with respect to $\alpha$. This seems to
suggest the independence of the sparse double descent phenomena with respect
to each layer’s width in a fixed neural network architecture.

Table 1: The experiment results as detailed in section 3. |
---|---
|
|

##  4 Discussion

Our exploration of the double descent phenomenon in simple two-layer neural
networks has yielded some novel insights into model capacity, sparsity, and
double descent phenomenona. The observed behavior in networks under varying
representation dimensions and L1 regularization strenghts underscores a
nuanced understanding of the balance between model complexity and overfitting.

However, this study’s focus on specific network architectures and datasets
suggests that broader research is required to generalize these findings.
Future work should aim to validate these phenomena in more complex models and
diverse datasets. Ultimately, this study contributes to the evolving discourse
on neural network training, positioning sparse double descent as a key
consideration in the pursuit of optimal model performance. This understanding
aligns our work with the broader trajectory of machine learning research,
highlighting the dynamic and often counter-intuitive nature of model training
and generalization.

###  4.1 Future Directions

A future direction for this research could be to expand the experiments to
include pretrained vision models fine-tuned on the MNIST dataset. This
approach would allow for an exploration of the neural network up to the
penultimate layer as a learned kernel, with the retraining of the final layer
classifier under different objectives or settings providing insights into
various forms of kernel regressions, such as ridge and LASSO.

Additionally, exploring other sparsification and pruning methods could be
highly valuable. Theoretical underpinnings for these methods could be provided
through Neural Tangent Kernel theory or other kernel theories. Integrating
intuitive explanations using concepts like Kolmogorov complexity and Minimum
Description Length (MDL) principles would further enhance the understanding.
Examining the role of data equivariances and symmetries in these contexts
would also contribute to a more comprehensive understanding of the dynamics
influencing neural network performance and complexity.

## References

  * Atanasov et al. (2022) Alexander Atanasov, Blake Bordelon, and Cengiz Pehlevan.  Neural networks as kernel learners: The silent alignment effect.  In _The Tenth International Conference on Learning Representations, ICLR 2022, Virtual Event, April 25-29, 2022_. OpenReview.net, 2022.  URL https://openreview.net/forum?id=1NvflqAdoom.
  * Belkin et al. (2019) Mikhail Belkin, Daniel Hsu, Siyuan Ma, and Soumik Mandal.  Reconciling modern machine-learning practice and the classical bias–variance trade-off.  _Proceedings of the National Academy of Sciences_ , 116(32):15849–15854, August 2019.  doi: 10.1073/pnas.1903070116.  URL https://www.pnas.org/doi/10.1073/pnas.1903070116.  Publisher: Proceedings of the National Academy of Sciences.
  * Curth et al. (2023) Alicia Curth, Alan Jeffares, and Mihaela van der Schaar.  A u-turn on double descent: Rethinking parameter counting in statistical learning.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.  URL https://openreview.net/forum?id=O0Lz8XZT2b.
  * Frankle and Carbin (2019) Jonathan Frankle and Michael Carbin.  The lottery ticket hypothesis: Finding sparse, trainable neural networks.  In _International Conference on Learning Representations_ , 2019.  URL https://openreview.net/forum?id=rJl-b3RcF7.
  * Geiger et al. (2020) Mario Geiger, Stefano Spigler, Arthur Jacot, and Matthieu Wyart.  Disentangling feature and lazy training in deep neural networks.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2020(11):113301, nov 2020.  doi: 10.1088/1742-5468/abc4de.  URL https://dx.doi.org/10.1088/1742-5468/abc4de.
  * He et al. (2022) Zheng He, Zeke Xie, Quanzhi Zhu, and Zengchang Qin.  Sparse double descent: Where network pruning aggravates overfitting.  In Kamalika Chaudhuri, Stefanie Jegelka, Le Song, Csaba Szepesvari, Gang Niu, and Sivan Sabato, editors, _Proceedings of the 39th International Conference on Machine Learning_ , volume 162 of _Proceedings of Machine Learning Research_ , pages 8635–8659. PMLR, 17–23 Jul 2022.  URL https://proceedings.mlr.press/v162/he22d.html.
  * Jacot et al. (2021) Arthur Jacot, Franck Gabriel, and Clément Hongler.  Neural tangent kernel: Convergence and generalization in neural networks (invited paper).  In _Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing_ , STOC 2021, page 6, New York, NY, USA, 2021. Association for Computing Machinery.  ISBN 9781450380539.  doi: 10.1145/3406325.3465355.  URL https://doi.org/10.1145/3406325.3465355.
  * LeCun et al. (1989) Yann LeCun, John Denker, and Sara Solla.  Optimal brain damage.  In D. Touretzky, editor, _Advances in Neural Information Processing Systems_ , volume 2. Morgan-Kaufmann, 1989.  URL https://proceedings.neurips.cc/paper_files/paper/1989/file/6c9882bbac1c7093bd25041881277658-Paper.pdf.
  * LeCun et al. (2010) Yann LeCun, Corinna Cortes, and CJ Burges.  Mnist handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , 2, 2010.
  * Papyan et al. (2020) Vardan Papyan, X. Y. Han, and David L. Donoho.  Prevalence of neural collapse during the terminal phase of deep learning training.  _Proceedings of the National Academy of Sciences_ , 117(40):24652–24663, 2020.  doi: 10.1073/pnas.2015509117.
  * Ramirez et al. (2013) Carlos Ramirez, Vladik Kreinovich, and Miguel Argaez.  Why l1 is a good approximation to l0: A geometric explanation.  _Journal of Uncertain Systems_ , 2013.