# Triglobal resolvent analysis of swept-wing wakes

Jean Hélder Marques Ribeiro1<EMAIL_ADDRESS>Chi-An Yeh1,2  Kunihiko
Taira1  1 Department of Mechanical and Aerospace Engineering, University of
California, Los Angeles, CA 90095, USA 2 Department of Mechanical and
Aerospace Engineering, North Carolina State University, Raleigh, NC 27695, USA

###### Abstract

Through triglobal resolvent analysis, we reveal the effects of wing tip and
sweep angle on laminar separated wakes over swept wings. For the present
study, we consider wings with semi-aspect ratios from $1$ to $4$, sweep angles
from $0^{\circ}$ to $45^{\circ}$, and angles of attack of $20^{\circ}$ and
$30^{\circ}$ at a chord-based Reynolds number of $400$ and a Mach number of
$0.1$. Using direct numerical simulations, we observe that unswept wings
develop vortex shedding near the wing root with a quasi-steady tip vortex. For
swept wings, vortex shedding is seen near the wing tip for low sweep angles,
while the wakes are steady for wings with high sweep angles. To gain further
insights into the mechanisms of flow unsteadiness, triglobal resolvent
analysis is used to identify the optimal spatial input-output mode pairs and
the associated gains over a range of frequencies. The three-dimensional
forcing and response modes reveal that harmonic fluctuations are directed
towards the root for unswept wings and towards the wing tip for swept wings.
The overlapping region of the forcing-response mode pairs uncovers triglobal
resolvent wavemakers associated with self-sustained unsteady wakes of swept
wings. Furthermore, we show that for low aspect ratio wings optimal
perturbations develop globally over the entire wingspan. The present study
uncovers physical insights on the effects of tip and sweep on the growth of
optimal harmonic perturbations and the wake dynamics of separated flows over
swept wings.

##  1 Introduction

Understanding flow separation over finite swept wings is essential to the
study of aircraft and biological flight (Anderson, 2010; Videler et al., 2004;
Lentink et al., 2007). The aspect ratio, angle of attack, and sweep play
important roles in influencing stall and wake characteristics (Zhang et al.,
2020a, b). Although a number of studies have deepened our knowledge of laminar
separated wakes around swept wings, coherent flow structures associated with
the three-dimensional ($3$-D) flow separation have not been characterized in a
comprehensive manner. Such findings would be crucial to explain the role
played by the perturbations in characterizing the wakes and support efforts to
control flow separation around finite wings.

Previous studies have shown the effect of sweep on post-stall wake
characteristics with focus on the role of spanwise flow over wings (Harper &
Maki, 1964). The spanwise flow induced by sweep delays the emergence of stall
(Yen & Hsu, 2007; Yen & Huang, 2009) and reduces wake oscillations, as shown
for high-Reynolds number flows over transonic buffets in biglobal (Crouch et
al., 2019; Paladini et al., 2019; Plante et al., 2021) and triglobal linear
stability analysis (Timme, 2020; He & Timme, 2021). Similar observations have
been made for flows around aircraft models in experiments (Masini et al.,
2020) and computations (Houtman et al., 2022).

At a low Reynolds number, direct numerical simulations (DNS) from Zhang et al.
(2020a) showed that sweep angle can significantly alter the wake patterns. For
wings with low sweep angles, vortex shedding develops near the wing tip, while
unsteadiness is suppressed for flows over highly swept wings. Similar
attenuation of flow unsteadiness was further studied for forward-swept wings
(Zhang & Taira, 2022) revealing that wing sweep has a strong effect on
attenuating wake oscillations in laminar flows. Furthermore, linear
instabilities around swept wings were examined for a variety of swept and
unswept wings, showing that the sweep angle suppresses the emergence of wake
modes (Burtsev et al., 2022; Ribeiro et al., 2022a).

The aspect ratio of the wing also affects the wake dynamics on separated flows
due to the wing tip vortex in steady (Devenport et al., 1996; Torres &
Mueller, 2004; Taira & Colonius, 2009) and unsteady wing motion (Buchholz &
Smits, 2006; Yilmaz & Rockwell, 2012). For low-aspect-ratio wings, the tip
vortex may suppress the leading-edge vortex formation, reducing the wake
unsteadiness (Taira & Colonius, 2009). Tip vortices can also produce adverse
effects on the wing, with induced drag and a reduced lift.

To alter the wake dynamics with a proper actuation input, we need to identify
the optimal forcing structures that can be amplified in the flow field
(Edstrand et al., 2018a, b). For this task, we may use modal analysis
techniques (Taira et al., 2017, 2020) to study the dynamics of flow
oscillations. Resolvent analysis is an attractive tool for the present study
because it identifies the optimal input perturbations in the flow field, their
energy amplification, and the characteristics of their unsteady response
(Trefethen et al., 1993; Jovanović & Bamieh, 2005). Furthermore, with the
diverse steady and unsteady wakes observed around swept wings, resolvent
analysis can provide a comprehensive study of the input-output dynamics around
wings with different aspect ratios, angles of attack, and sweep.

Resolvent analysis has been used to study a broad range of fluid flows
(Moarref et al., 2013; Thomareis & Papadakis, 2018; Schmidt et al., 2018;
Skene & Schmid, 2019; Yeh et al., 2020; Ricciardi et al., 2022). This approach
was initially formulated for steady base flows, to identify modal structures
that can be amplified in stable flow regimes (Trefethen et al., 1993). This
perspective on fluid dynamics was later extended to unstable systems by
Jovanović & Bamieh (2005) and to unsteady and turbulent flows by McKeon &
Sharma (2010). In these formulations, a time-averaged flow is used as a base
state and nonlinear terms act as sustained forcing in the flow field. In both
steady and unsteady flows, resolvent analysis identifies harmonic forcings
that produce an amplified response in the flow.

In this study, we identify the optimal spatial input-output modes around the
wing through a $3$-D global (triglobal) resolvent analysis, that assumes no
spatial homogeneity. Moreover, we gain insights into the self-sustained
fluctuations that support unsteadiness on laminar separated flows using
resolvent wavemakers, which are similar in spirit to the eigenvector-based
wavemakers (Giannetti & Luchini, 2007; Giannetti et al., 2010). The resolvent
wavemakers, also named as structural sensitivity, are obtained from the
overlap of forcing and response modes (Qadri & Schmid, 2017; Skene et al.,
2022b). These findings provide a comprehensive analysis of the energy
amplification mechanisms in flows around swept wings through an input-output
process, identifying the optimal locations where perturbations can be
introduced to alter the wake behavior. Therefore, these findings are crucial
for the development of efficient flow control strategies (Yeh & Taira, 2019;
Liu et al., 2021) that aim to improve the aerodynamic performance of swept
wings experiencing massive flow separation.

The present paper on triglobal resolvent analysis is organized as follows. In
section 2, we describe the problem setup for the current work. In section 3,
we discuss our main findings from triglobal resolvent analysis. We identify
the emergence of wake unsteadiness caused by the overlap of optimal forcing
and response modes in the near wake. Perturbations are directed towards the
region where vortex shedding takes place. The locations of the optimal forcing
and response modes over the wingspan also suggest that wakes of highly swept
wings are more resilient to external perturbations. Furthermore, we find that
low-aspect-ratio wings limit the growth of perturbations to global modes
extending over the entire wingspan. Finally, our conclusions are presented in
section 4.

##  2 Problem setup

$\Lambda$

$x$

$y$

$z$

$y$

$x$

$U_{\infty}$

$c$

$\alpha$

$x$

$z$

$b$

DNS domain($x$,$y$) plane($x$,$z$) planeresolvent analysis

$-0.2$

$1.0$

$u_{x}$

Figure 1: Setup for finite swept wing simulation. In the gray boxes, the
instantaneous flow field for $\alpha=20^{\circ}$, $\Lambda=15^{\circ}$, and
$sAR=b/c=4$, with $Q=2$ isosurfaces colored by instantaneous $u_{x}$. Mesh
colored by time-averaged $\bar{u}_{x}$. In the blue box, isosurfaces of the
primary response mode with mesh in light gray.

We consider laminar flows over untapered swept wings with NACA 0015 cross-
sectional profile, as shown in figure 1. The spatial coordinates are defined
with $(x,y,z)$ being the streamwise, transverse, and spanwise directions,
respectively, with the origin placed at the leading edge of the wing root. The
NACA 0015 airfoil geometry is defined on the $(x,y)$ plane. The wingspan is
formed by extruding the airfoil profile in the spanwise direction. The semi-
aspect ratio is defined through the half-span length $b$ and the chord-length
$c$ as $sAR=b/c$, with values set between $1\leq sAR\leq 4$. For swept wings,
the $3$-D computational setup is sheared in the $x$-direction and the sweep
angle is defined between the $z$-direction and the leading edge. In the
present work, we consider sweep angles $0\leq\Lambda\leq 45^{\circ}$. The
angle of attack, $\alpha=20^{\circ}$ and $30^{\circ}$, is defined between the
streamwise direction and the airfoil chord line. To focus on the effects of
wing tip and sweep in the wake dynamics, we analyze a half-span model with
symmetry boundary conditions imposed at the root plane. The wings have a
straight-cut tip and sharp trailing edge. For all flows analyzed herein, we
define the chord-based Reynolds number $Re_{c}=U_{\infty}c/\nu=400$, where
$U_{\infty}$ is the freestream velocity and $\nu$ is the kinematic viscosity.
The freestream Mach number is set to $M_{\infty}=U_{\infty}/a_{\infty}=0.1$,
where $a_{\infty}$ is the freestream speed of sound.

###  2.1 Direct numerical simulation

We perform DNS with the compressible flow solver CharLES (Khalighi et al.,
2011; Brès et al., 2017), which uses a second-order accurate finite-volume
method in space with a third-order accurate scheme in time. With the origin at
the leading edge of the airfoil $(x/c,\ y/c,\ z/c)=(0,0,0)$, the computational
domain extends over $(x/c,y/c,z/c)\in[-20,25]\times[-20,20]\times[0,20]$. We
build a C-type grid for each angle of attack with $\min$ $(\Delta x,\Delta
y,\Delta z)/c=(0.005,0.005,0.04)$ applying mesh refinement near the airfoil
and in the wake, as shown in figure 1.

Inlet and farfield boundaries are prescribed with Dirichlet boundary
conditions
$(\rho,u_{x},u_{y},u_{z},p)=(\rho_{\infty},U_{\infty},0,0,p_{\infty})$, where
$\rho$ is density, $p$ is pressure, and $u_{x}$, $u_{y}$, and $u_{z}$ are
velocity components in the $x$, $y$, and $z$ directions, respectively.
Variables with subscript $\infty$ denote freestream values. For all setups
considered herein, the velocity boundary conditions applied at the inlet and
farfield are aligned with the $x$-direction, which enforces the same
streamwise flow over wings with different sweep angles. The airfoil surface is
provided with adiabatic no-slip boundary condition. To simulate a half-wing
model, we prescribe the symmetry boundary condition along the root plane. A
sponge layer is applied at the outlet over $x/c\in[15,25]$ with the target
state being the running-averaged flow over $5$ convective time units $t\equiv
c/U_{\infty}$ (Freund, 1997). Simulations start with uniform flow and time
integration is performed with a constant acoustic Courant–Friedrichs–Lewy
(CFL) number of $1$. After transients are flushed out of the computational
domain, the time-averaged base flow $\bar{\mathbi{q}}$ is determined over $50$
convective time units. The present results were carefully verified and
validated. Close agreement for instantaneous and time-averaged velocity
components was achieved with those from Zhang et al. (2020a). We have further
validated our computations for time-averaged drag and lift coefficients,

$C_{D}=\frac{F_{x}}{\frac{1}{2}\rho U_{\infty}^{2}bc}\quad\text{and}\quad
C_{L}=\frac{F_{y}}{\frac{1}{2}\rho U_{\infty}^{2}bc}\mbox{ ,}$ (1)

respectively, where $F_{x}$ is the drag and $F_{y}$ is the lift over the wing,
as reported in table 1.

| $\Lambda=0^{\circ}$ |  | $\Lambda=15^{\circ}$ |  | $\Lambda=30^{\circ}$ |  | $\Lambda=45^{\circ}$
---|---|---|---|---|---|---|---
| $\overline{C_{L}}$ | $\overline{C_{D}}$ |  | $\overline{C_{L}}$ | $\overline{C_{D}}$ |  | $\overline{C_{L}}$ | $\overline{C_{D}}$ |  | $\overline{C_{L}}$ | $\overline{C_{D}}$
Present study |  $0.53$ |  $0.35$ |  |  $0.50$ |  $0.34$ |  |  $0.45$ |  $0.31$ |  |  $0.40$ |  $0.29$
Zhang et al. (2020a) |  $0.53$ |  $0.35$ |  |  $0.51$ |  $0.33$ |  |  $0.44$ |  $0.30$ |  |  $0.40$ |  $0.29$
Table 1: Time-averaged lift and drag coefficients ($\overline{C_{L}}$ and
$\overline{C_{D}}$) compared to Zhang et al. (2020a) for laminar separated
flow over NACA 0015 wings with $sAR=4$, $\alpha=20^{\circ}$, and
$\Lambda=0^{\circ}$, $15^{\circ}$, $30^{\circ}$, and $45^{\circ}$.

$\Lambda$

$\Lambda$

$\alpha=20^{\circ}$

$\alpha=30^{\circ}$

$sAR$

$-0.2$

$1.0$

$u_{x}$

Figure 2: Instanteneous isosurfaces of $Q=2$ colored by $u_{x}$. Unsteady
shedding near wing root ( $\bm{\diamond}$), unsteady shedding near wing tip (
$\bm{\square}$), steady flow with root structures (
$\bm{\blacktriangleright}$), steady flow with streamwise vortices (
$\bm{\blacktriangledown}$).

A variety of wake patterns can be observed for different $\alpha$, $\Lambda$,
and $sAR$, as summarized in figure 2. In the bottom plot, the flow over the
wing with $(sAR,\alpha,\Lambda)=(2,30^{\circ},0^{\circ})$ exhibits a quasi-
steady streamwise oriented tip vortex. This structure is characteristic of
flows over unswept wings and also appears around wings with different $\alpha$
and $sAR$. For such wings, unsteady spanwise vortices develop at the root
plane. Between the root and the wing tip, there is an intermediate zone with
braid-like vortices.

Wing sweep affects the wake dynamics and structures. At low sweep angles, a
spanwise flow develops over the wing and advects unsteady vortices towards the
wing tip. For instance, for $(sAR,\alpha,\Lambda)=(4,20^{\circ},15^{\circ})$,
spanwise vortices still appear. These structures are similar to the ones
observed over unswept wings, although they form closer to the wing tip, and
break into helical structures in the wake. For the flow around this wing,
streamwise-oriented tip vortices are absent.

For higher $\Lambda$, wing sweep can attenuate wake oscillations. For
$(sAR,\alpha,\Lambda)=(4,30^{\circ},30^{\circ})$ near-wake unsteadiness is
reduced and unsteady vortices appear further downstream in the wake. We notice
that these structures are absent when the angle of attack is lowered to
$20^{\circ}$. Similarly, increasing the sweep angle to $\Lambda=45^{\circ}$
suppresses unsteady vortices on both angles of attack and the wake becomes
steady. On such highly swept wings, ram-horn shaped streamwise-oriented
vortices develop from the root plane and extend into the wake.

For each $(\alpha,\Lambda)$ pair, the wake exhibits similar characteristics
for wings with $sAR=4$ and $2$. Reducing the semi-aspect ratio to $sAR=1$ has
a strong influence on the wake dynamics, as shown in figure 2. For such wings,
tip effects can suppress the formation of leading edge vortices at lower
angles of attack. For instance, at $\alpha=20^{\circ}$, wake unsteadiness is
reduced and swept wings exhibit steady flows with root structures.

Unsteady vortices are observed in flows over $sAR=1$ wings at
$\alpha=30^{\circ}$ for all considered sweep angles. The unsteadiness appears
near the root for lower $\Lambda$. Further downstream, unsteady vortices
appear over the entire wingspan. For higher $\Lambda$, vortices are generated
near the wing tip and helical structures are observed in the wake. These
observations agree with the characterizations by Zhang et al. (2020a). To
deepen our insights into swept-wing wake dynamics, we now call for the
triglobal resolvent analysis.

###  2.2 Resolvent analysis

Let us consider the Reynolds decomposition of state variable
$\mathbi{q}=\bar{\mathbi{q}}+\mathbi{q}^{\prime}$, where $\bar{\mathbi{q}}$ is
the time-averaged flow and $\mathbi{q}^{\prime}$ is the statistically
stationary fluctuation component (McKeon & Sharma, 2010). This decomposition
along with spatial discretization is used to linearize the compressible
Navier–Stokes equations about $\bar{\bm{q}}$ to yield

$\frac{\partial\mathbi{q}^{\prime}}{\partial
t}=\mathsfbi{L}_{\bar{\mathbi{q}}}\mathbi{q}^{\prime}+\mathbi{f}^{\prime}\mbox{
,}$ (2)

where $\mathsfbi{L}_{\bar{\mathbi{q}}}$ is the discrete linearized
Navier–Stokes operator (Sun et al., 2017) and $\mathbi{f}^{\prime}$ accounts
for the external forcing and nonlinear terms. With the Fourier representation

$\left[\mathbi{q}^{\prime}(\mathbi{x},t),\mathbi{f}^{\prime}(\mathbi{x},t)\right]=\int_{-\infty}^{\infty}\left[\hat{\mathbi{q}}_{\omega}(\mathbi{x}),\hat{\mathbi{f}}_{\omega}(\mathbi{x})\right]e^{-i\omega
t}{\rm d}\omega\mbox{ ,}$ (3)

we obtain

$-i\omega\hat{\mathbi{q}}_{\omega}=\mathsfbi{L}_{\bar{\mathbi{q}}}\hat{\mathbi{q}}_{\omega}+\hat{\mathbi{f}}_{\omega}\mbox{
,}$ (4)

where $\mathbi{x}=(x,y,z)$ and the triglobal response and forcing modes are
$\hat{\mathbi{q}}_{\omega}$ and $\hat{\mathbi{f}}_{\omega}$, respectively, for
a temporal frequency $\omega$. This expression leads to

$\hat{\mathbi{q}}_{\omega}=\mathsfbi{H}_{\bar{\mathbi{q}},\omega}\hat{\mathbi{f}}_{\omega}\mbox{
,}$ (5)

in which the resolvent operator
$\mathsfbi{H}_{\bar{\mathbi{q}},\omega}\in\mathbb{C}^{m\times m}$, with $m$
defined by the product of the number of state variables and the number of
spatial grid points. For the present triglobal base flows, the linear
operators have size $m$ between $3$ and $5\times 10^{6}$. We analyze the
resolvent operator with the singular value decomposition (SVD)

$\mathsfbi{H}_{\bar{\mathbi{q}}}=\left[-i\omega\mathsfbi{I}-\mathsfbi{L}_{\bar{\mathbi{q}}}\right]^{-1}=\mathsfbi{Q}\mathbi{\Sigma}\mathsfbi{F}^{*},$
(6)

where
$\mathsfbi{F}=[\hat{\mathbi{f}}_{1},\hat{\mathbi{f}}_{2},\dots,\hat{\mathbi{f}}_{m}]$
is an orthonormal matrix holding the forcing modes, $\bm{\Sigma}={\rm
diag}[\sigma_{1},\sigma_{2},\dots,\sigma_{m}]$ is the diagonal matrix with
singular values (gain) in descending order, and
$\mathsfbi{Q}=[\hat{\mathbi{q}}_{1},\hat{\mathbi{q}}_{2},\dots,\hat{\mathbi{q}}_{m}]$
is the orthonormal matrix comprised of the response modes (Trefethen et al.,
1993; Jovanović & Bamieh, 2005). For visualization purposes, we show only the
real part of the complex-valued resolvent modes. Here, we employ the Chu norm
(Chu, 1965) incorporating it within the resolvent operator through a
similarity transform
$\mathsfbi{H}_{\bar{\mathbi{q}}}\rightarrow\mathsfbi{W}^{\frac{1}{2}}\mathsfbi{H}_{\bar{\mathbi{q}}}\mathsfbi{W}^{-\frac{1}{2}}$,
where $\mathsfbi{W}$ is the weight matrix that accounts for numerical
quadrature and energy weights.

Resolvent analysis requires careful consideration of the eigenvalues of
$\mathsfbi{L}_{\bar{\mathbi{q}}}$. In the presence of unstable modes in the
linear operator eigenspectrum, the asymptotic input-output relationship is
buried under the unstable dynamics behavior. The present resolvent analysis
utilizes a time-averaged flow as the base state. Since such flow is not the
equilibrium state, stability characterization cannot be performed in a strict
sense. However, it is important to check the location of the eigenvalues in
the complex plane to capture the growth rate of the most unstable modes of
$\mathsfbi{L}_{\bar{\mathbi{q}}}$.

To use resolvent analysis to study the wake dynamics of unstable base flows,
we examine the dynamics through the lens of temporal discounting (Jovanović,
2004). Discounting applies a temporal damping on forcing and response modes as
$[\hat{\mathbi{q}}_{\omega},\hat{\mathbi{f}}_{\omega}]e^{-\beta t}$, where
$\beta>0$ is a time-discounting parameter defined within the discounted
resolvent operator (Jovanović, 2004). With the discounted resolvent analysis,
we can examine amplification dynamics that takes place on a time scale shorter
than that of the most unstable mode. Detailed discussions on our choice of
$\beta$ are provided in appendix A. Through the discounted resolvent analysis,
valuable insights have been provided in past studies for the dynamics and
control of flows over airfoils (Yeh & Taira, 2019; Yeh et al., 2020; Ricciardi
et al., 2022; Ribeiro et al., 2022b).

The $\mathsfbi{H}_{\bar{\mathbi{q}},\omega}$ operators were discretized over
$3$-D structured grids with the leading edge at the root positioned at $(x/c,\
y/c,\ z/c)=(0,0,0)$, extending over
$(x/c,y/c,z/c)\in[-10,15]\times[-10,10]\times[0,10]$ with near wake grids
shown at the bottom left of figure 1. The computational grids used for
resolvent analysis have a smaller domain size than those used for DNS. For the
base flow on the mesh for resolvent analysis, we perform a linear
interpolation of the flow field from DNS mesh to the resolvent mesh. We
prescribe homogeneous Neumann boundary conditions for $T^{\prime}$ and
homogeneous Dirichlet boundary conditions for the fluctuating variables
$\rho^{\prime}$ and $u^{\prime}$ along the farfield, airfoil surface, and
outlet. Sponges are applied far from the airfoil and in conjunction with the
boundary conditions (Freund, 1997).

For the large linear operators in the present work, efficient numerical tools
are needed for SVD (Halko et al., 2011). We use the randomized resolvent
analysis algorithm from Ribeiro et al. (2020), sketching
$\mathsfbi{H}_{\bar{\mathbi{q}},\omega}$ with $10$ random test vectors. Each
entry of the test vectors is associated with a particular grid point and the
five state variables, scaled by $[\|\nabla\rho\|,\|\nabla u_{x}\|,\|\nabla
u_{y}\|,\|\nabla u_{z}\|,\|\nabla T\|]$ at each spatial location for each
state variable (Ribeiro et al., 2020; House et al., 2022). A convergence
analysis of the randomized resolvent algorithm is provided in appendix B.

The computation of resolvent modes for large linear operators can be
challenging for the resolvent analysis of high-Reynolds number flows that
require a large grid. The bottleneck is related to the time and memory
requirements of the linear systems solvers within the SVD. Building an optimal
basis to avoid linear system solvers is possible (Barthel et al., 2022),
although a generalization for complex geometries is still challenging. It is
possible, however, to obtain accurate resolvent modes with time-stepping
instead of direct solvers. Those methods tend to penalize the computational
time costs, although a considerable reduction in memory requirements can be
achieved (Barkley et al., 2008; Monokrousos et al., 2010; Gómez et al., 2016).
The computational time required by time-steppers can also be reduced by
incorporating streaming discrete Fourier transforms (Martini et al., 2021;
Farghadan et al., 2021). The use of iterative solvers has shown promising
results to compute resolvent modes around a commercial aircraft model (Houtman
et al., 2022).

In the present work, the direct and adjoint linear systems were directly
solved using the MUMPS (multifrontal massively parallel sparse direct solver)
package (Amestoy et al., 2001). Moreover, we incorporate the adjoint-based
sensitivity analysis to interpolate the resolvent norm over frequencies
$\omega$ (Schmid & Brandt, 2014; Fosas de Pando & Schmid, 2017). This approach
is used to calculate the gradient of $\sigma$ with respect to $\omega$,
allowing an accurate interpolation among frequencies (Skene & Schmid, 2019).
The codes used to compute the resolvent modes are part of the ‘linear analysis
package’ made available by Skene et al. (2022a).

##  3 Triglobal resolvent analysis

###  3.1 Forcing and response modes structures

Let us first examine the dominant gains, forcing, and response modes for
$(sAR,\alpha,\Lambda)=(4,20^{\circ},0^{\circ})$, as shown in figure 3. The
dominant resolvent modes are observed at $St=0.14$, where

$St=\frac{\omega}{2\pi}\frac{c\sin\alpha}{U_{\infty}\cos\Lambda}\mbox{ }$ (7)

is the Fage–Johansen Strouhal number (Fage & Johansen, 1927) with a
$1/\cos\Lambda$ scaling that incorporates the influence of the sweep angle.
This frequency scaling is inspired by the independence principle (Wygnanski et
al., 2011) and collapses the spectral behavior of the resolvent modes over
different sweep angles (Ribeiro et al., 2022b). This frequency matches the
peak frequency for the lift coefficient shown in figure 3 (bottom left). The
dominant frequency for $\sigma_{1}$ and $\hat{C}_{L}$ agrees for all unsteady
flows presented herein.

The spatial structures of forcing-response mode pairs are shown in figure 3
(right) for representative frequencies. For $St=0.14$, primary modes exhibit
modal structures near the root plane. The forcing mode appears near and
upstream of the wing, while the response mode develops downstream in the wake.
The modal structures for the primary forcing and response modes are aligned
with the wingspan, with the response mode similar to the unsteady vortices
revealed from DNS. At this frequency, the secondary modes are comprised of
spanwise-aligned vortices near the root plane, similar to the primary modes.

$\sigma_{j}$

$\hat{C}_{L}$

$St=(\omega/2\pi)(c\sin\alpha/U_{\infty}\cos\Lambda)$

$St=0.14$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{2}$

$\hat{q}_{2}$

$\hat{f}_{2}$

$\hat{q}_{2}$

$\hat{f}_{2}$

$\hat{q}_{2}$

$\hat{f}_{2}$

$\hat{q}_{2}$

$St=0.14$

$St=0.16$

$St=0.18$

$St=0.20$

Figure 3: Resolvent gains and forcing-response mode pairs for
$(sAR,\alpha,\Lambda)=(4,20^{\circ},0^{\circ})$. For each mode, forcing
($\hat{\mathbi{f}}$) is the top-half while response ($\hat{\mathbi{q}}$) it
the bottom-half with isosurfaces of velocity $\hat{u}_{y}\in[-0.2,0.2]$, with
freestream directed to the right. On bottom left, power spectrum density of
lift coefficient $\hat{C}_{L}$.

As we increase the frequency, the resolvent gains decay in magnitude and
$\sigma_{1}$ decays faster than $\sigma_{2}$. Their magnitudes become
approximately the same at $St=0.16$. At this frequency, the spatial
characteristics of the primary and secondary forcing-response mode pairs
exhibit distinct behavior. The primary forcing and response modes are aligned
with the wingspan and near the root plane, similar to those at lower $St$. The
secondary modes, however, exhibit modal structures near the wing tip, in
contrast to the secondary modes at lower frequencies which reside near the
wing root.

For $St=0.18$, the primary forcing-response mode pair appears near the wing
tip, while the secondary mode pair develops at the root plane. Such behavior
persists as we increase the frequency to $St=0.20$. For $St\geq 0.18$, primary
modes are tip-dominated while secondary modes are root-dominated around this
wing. This means that root and wing tip modes switch their order of
amplification at $St\approx 0.18$, i.e., mode switching.

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$St$

$St$

$St$

$\sigma$

$\Lambda=15^{\circ}$

$\Lambda=30^{\circ}$

$\Lambda=45^{\circ}$

Figure 4: Resolvent gain distribution for the top three mode pairs and
forcing-response mode pairs for selected frequencies for
$(sAR,\alpha)=(4,20^{\circ})$ and $\Lambda=15^{\circ}$, $30^{\circ}$, and
$45^{\circ}$. Primary forcing ($\hat{f}$) and response ($\hat{q}$) modes shown
with isosurfaces of velocity $\hat{u}_{y}\in[-0.2,0.2]$.

The mode switching phenomenon is also observed for swept wings with
$\Lambda=15^{\circ}$. For such wings, root-supported structures appear as the
primary forcing-response pairs at $St=0.14$, as shown in figure 4 (left). A
distinct mode switching is observed over this wing, as the forcing-response
pairs gradually transition towards the root at $z/c\approx 0$ with the
increase in $St$. This type of concentrated resolvent mode at the wing root
also appears for the unswept wings at $St=0.20$ as a secondary mode, as shown
in figure 3 (bottom, right).

For $\Lambda=30^{\circ}$, mode switching also occurs towards the root with the
increase in $St$, in an opposite trend to the unswept wings. The dominant
response modes at lower frequencies appear at the wing tip, as shown in figure
4 (middle). There is a gradual transition to root-supported modes as the
frequency increases. At a higher sweep angle, $\Lambda=45^{\circ}$, no mode
switching occurs. The region of dominance of the forcing and response modes is
slightly invariant for the frequencies shown herein.

In contrast with the lower sweep angle wings, for $\Lambda=45^{\circ}$,
forcing and response modes are dominant at distinct wingspan locations, as
shown in figure 4 (right). Response modes are tip-dominated while forcing
structures appear upstream near the root plane, extending over the wingspan
aligned with the sweep angle. For all $(sAR,\alpha)=(4,20^{\circ})$ wings, the
highest amplification is found for $\Lambda=30^{\circ}$, at $St\approx 0.12$.
At $\Lambda=45^{\circ}$, the dominant gain is an order of magnitude lower.
This finding suggests that it is challenging to perturb flows over
$\Lambda=45^{\circ}$ wings. These wings are steady because self-sustained flow
disturbances cannot introduce sufficient energy into the wake to generate
vortex shedding.

###  3.2 Resolvent wavemakers

To characterize the self-sustained unsteadiness in the flows over swept wings,
we study the spatial overlap between the forcing and response modes that
supports the continuous formation of vortical structures. Since the forcing
modes show regions receptive to external perturbations and the response modes
reveal the structures being excited due to the forcing, the region over which
forcing and response modes overlap can be interpreted as a mechanism for self-
sustained oscillations in the flow. This idea is similar to the wavemaker
concept deduced from direct and adjoint eigenmodes presented in Giannetti &
Luchini (2007).

Through the wavemaker analysis, previous studies identified critical points
responsible for sustaining wake shedding on laminar wakes around cylinders
(Strykowski & Sreenivasan, 1990; Hill, 1992) and regions associated with their
primary and secondary instability modes (Giannetti & Luchini, 2007; Giannetti
et al., 2010). Moreover, wavemakers revealed the physical mechanisms
responsible for tonal noise generation in high-Reynolds number flows over
airfoils (Fosas de Pando et al., 2017) and self-sustained flow instabilities
in transonic buffet regimes (Paladini et al., 2019).

In the aforementioned studies, wavemakers were derived from direct and adjoint
global stability eigenmodes. Our formulation derives wavemakers from global
resolvent modes and it is closely related to the structural sensitivity
devised by Qadri & Schmid (2017) and to the resolvent wavemaker studied by
Skene et al. (2022b). The present resolvent wavemaker is not identical to the
eigenvector-based wavemaker. Using the time-averaged base flow, the present
forcing terms encapsulate nonlinear effects as a internal feedback mechanism
within the flow. Hence, the spatial overlap between forcing and response
identifies regions responsible for self-sustained wake oscillations.

Herein, the resolvent wavemaker modes are directly obtained from the resolvent
modes, as the Hadamard product of forcing and response modes

$\hat{\mathbi{w}}=\hat{\mathbi{f}}\circ\hat{\mathbi{q}}\mbox{ ,}$ (8)

where $\hat{\mathbi{w}}$ is the resolvent wavemaker mode. The resolvent modes
presented herein are defined with the five state variables,
$\hat{\mathbi{f}}=[\hat{\mathbi{f}}_{\rho},\hat{\mathbi{f}}_{u_{x}},\hat{\mathbi{f}}_{u_{y}},\hat{\mathbi{f}}_{u_{z}},\hat{\mathbi{f}}_{T}]$
and
$\hat{\mathbi{q}}=[\hat{\mathbi{q}}_{\rho},\hat{\mathbi{q}}_{u_{x}},\hat{\mathbi{q}}_{u_{y}},\hat{\mathbi{q}}_{u_{z}},\hat{\mathbi{q}}_{T}]$.
We define our resolvent wavemaker gain $\xi$ as

$\xi=\sigma^{2}\int_{S}|\hat{\mathbi{w}}(\mathbi{x})|\ {\rm d}S\mbox{ ,}$ (9)

which follows
$\xi=\sigma^{2}|\langle\hat{\mathbi{f}},\hat{\mathbi{q}}\rangle|$, derived by
Skene et al. (2022b). This is similar to the one presented in Ribeiro et al.
(2022b). Qualitatively, both definitions of $\xi$ result in similar
discussions and interpretations. The present expression for the resolvent
wavemaker gain provides a proper quantitative definition (Skene et al.,
2022b). The resolvent wavemaker gain $\xi$ can also be computed for each
spanwise slice and each frequency. To this end, we consider $S=S(x,y)$, as
$z$-normal planes at different spanwise locations, to build the $\xi$-contours
shown in figure 5. Through this analysis, we highlight the spatial support of
the resolvent wavemaker over the wingspan.

$\xi$

$St$

$St$

$z/c$

$z/c$

$\Lambda=0^{\circ}$

$\Lambda=15^{\circ}$

$\Lambda=30^{\circ}$

$\Lambda=45^{\circ}$

Figure 5: Wingspan location of primary resolvent wavemakers with isocontours
of $\xi$ for $0.05\leq St\leq 0.25$ for $(sAR,\alpha)=(4,20^{\circ})$ wings
with $0^{\circ}\leq\Lambda\leq 45^{\circ}$. Resolvent wavemaker modes at
$St=0.14$ are shown with isosurfaces of
$\hat{u}_{y}/\|\hat{u}_{y}\|_{\infty}=\pm 0.1$ and instantaneous $Q=1$ are
gray-colored isosurfaces.

Let us focus our resolvent wavemaker analysis on the flow over the unswept
wing with $(sAR,\alpha,\Lambda)=(4,20^{\circ},0^{\circ})$, as shown in figure
5 (top, left). At $St=0.14$, triglobal resolvent wavemakers with high $\xi$
appear between $2\lesssim z/c\lesssim 3$ in the near wake. The resolvent
wavemakers at this region support the formation of unsteady root vortices that
propagate downstream in the wake. This resolvent wavemaker region is also
characterized by the formation of braid-like structures that connect to the
root shedding as vortex loops (Zhang et al., 2020a). Resolvent wavemakers for
$(sAR,\alpha,\Lambda)=(4,20^{\circ},15^{\circ})$ also show similar shedding
behavior, as seen in figure 5 (top, right).

Resolvent wavemakers are also revealed for steady flows. The overlap of
forcing and response modes for flows over wings with
$(sAR,\alpha,\Lambda)=(4,20^{\circ},30^{\circ})$, as shown in figure 5
(bottom, left), develops over the wing and extends into the wake aligned with
the wing tip. These resolvent wavemakers extend over the entire wingspan,
being stronger and larger than the ones exhibited around wings with lower
sweep angles. The $\xi$ peak appears at the wing tip at $St=0.12$, indicating
that the tip region is more susceptible to develop unsteadiness around this
wing.

For $\Lambda=45^{\circ}$, shown in figure 5 (bottom, right), we reveal that
resolvent wavemakers emerge from the leading edge near the root plane towards
the wing tip and downstream at the wake, overlapping the region where steady
ram-horn-shaped vortices appear in the DNS. These resolvent wavemakers exhibit
a region of the flow field with high receptiveness to amplify forcing
structures and disturb the steady ram-horn vortex. Because the dominant
resolvent wavemaker around the $\Lambda=45^{\circ}$ wing have a low $\xi$, in
spite of occupying a large region of the wake, the energy they introduce to
the flow field is insufficient to disturb the wake.

The resolvent wavemakers further exhibit the root- and tip-dominated modal
characteristics and the mode switching phenomenon in figure 5, in agreement
with the forcing-response modal behavior shown in figures 3 and 4. For
instance, the resolvent wavemaker modes at the peak $\xi$ values for the
unswept wing appear near $z/c\approx 2$, with a gradual transition from root-
supported to tip-dominated modes as $St$ increases. Moreover, for the
$\Lambda=15^{\circ}$ wing, there is a transition in the dominant region of
resolvent wavemaker support from $z/c\approx 1.5$ at lower frequencies to
$z/c\approx 0$ at higher frequencies, as shown in figure 5 (top, right).
Lastly, for the $\Lambda=30^{\circ}$ wing, there is a tip-to-root transition
with the increase in $St$ while the peak resolvent wavemakers for
$\Lambda=45^{\circ}$ are invariant over the frequencies, appearing near
$z/c\approx 2$.

###  3.3 Forcing-to-response dynamics

Let us further explain how perturbations emerge around swept wings, by
analyzing the overlap between the forcing and response modes in the spanwise
direction. To this end, we integrate the norm of $\hat{\mathbi{f}}$ and
$\hat{\mathbi{q}}$ over $z$-normal planes, analogous to the resolvent
wavemaker modes analysis, as

$\mathbi{\Omega}_{\hat{\mathbi{f}}}(z)=\int_{S(x,y)}\|\hat{\mathbi{f}}\|_{2}\
{\rm
d}S\quad\text{and}\quad\mathbi{\Omega}_{\hat{\mathbi{q}}}(z)=\int_{S(x,y)}\|\hat{\mathbi{q}}\|_{2}\
{\rm d}S\mbox{ ,}$ (10)

where $\|\hat{\mathbi{f}}\|_{2}$ and $\|\hat{\mathbi{q}}\|_{2}$ are the
$2$-norm of
$[\hat{\mathbi{f}}_{\rho},\hat{\mathbi{f}}_{u_{x}},\hat{\mathbi{f}}_{u_{y}},\hat{\mathbi{f}}_{u_{z}},\hat{\mathbi{f}}_{T}]$
and
$[\hat{\mathbi{q}}_{\rho},\hat{\mathbi{q}}_{u_{x}},\hat{\mathbi{q}}_{u_{y}},\hat{\mathbi{q}}_{u_{z}},\hat{\mathbi{q}}_{T}]$,
respectively, at each grid point of the computational domain. By performing
the integral over $S(x,y)$, we obtain $\mathbi{\Omega}_{\hat{\mathbi{f}}}$ and
$\mathbi{\Omega}_{\hat{\mathbi{q}}}$ computed for each spanwise slice and for
each frequency. Here, we plot their contours normalized by the maximum
$\mathbi{\Omega}_{\hat{\mathbi{f}}}$ and $\mathbi{\Omega}_{\hat{\mathbi{q}}}$
at each $St$, to emphasize the spatial support of forcing and response over
the wingspan, as shown in figure 6 for wings at $\alpha=30^{\circ}$ with
$0^{\circ}\leq\Lambda\leq 45^{\circ}$ . The locations of the maximum strength
of forcing and response modes are shown by the dot-dashed lines. Black arrows
indicate the direction from the maximum forcing to the maximum response at
$St=0.15$. This analysis depicts the preferential direction in which optimal
forcing is transferred to optimal response over the wingspan at each
frequency.

For unswept wings, shown in figure 6 (left), the optimal forcing structures
appear closer to the wing tip than the response modes, which are slightly
shifted towards the root, suggesting that fluctuations are directed towards
the root. Indeed, as seen in the DNS, unsteadiness is concentrated towards the
root, as evident from figure 2, also in agreement with the results reported by
Zhang et al. (2020b). In addition, the flow around the wing tip for unswept
wings is characterized by an almost steady tip vortex, suggesting that it is
likely hard to amplify flow oscillations near the tip.

$St$

$z/c$

$z/c$

$z/c$

$z/c$

$\Lambda=0^{\circ}$

$\Lambda=15^{\circ}$

$\Lambda=30^{\circ}$

$\Lambda=45^{\circ}$

Figure 6: Wingspan locations of dominant forcing (red) and response (blue)
with contours of $\mathbi{\Omega}_{\hat{\mathbi{f}}}$ and
$\mathbi{\Omega}_{\hat{\mathbi{q}}}\in[0.4,1.0]$ for $0.05\leq St\leq 0.35$
for $sAR=4$ wings with $\alpha=30^{\circ}$ with $0^{\circ}\leq\Lambda\leq
45^{\circ}$. Dot-dashed lines are polynomial fit of maximum $z/c$ of forcing
and response at each $St$. Arrows show direction of optimal forcing-to-
response at $St=0.15$.

For swept wings, fluctuations are directed towards the wing tip. For
$\Lambda=15^{\circ}$, both forcing and response modes appear near the wing
root. At the vortex shedding frequency for this wing, $St\approx 0.15$, we
observe forcing and response modes to be dominant at $z/c\approx 1$, with the
forcing mode supported closer to the wing root than the response mode. This
concurs with the flow field we observe in the DNS, as vortices are formed near
the wing root and evolve towards the wing tip where spanwise vortices appear
and propagate in the wake. For the $\Lambda=30^{\circ}$ wing, the dominant
forcing-response mode pair emerges near the wing tip at low $St$, as seen in
figure 6. For this wing, low-frequency vortical structures emerge downstream
in the wake aligned at the tip, as shown in figure 2.

For the $\Lambda=45^{\circ}$ wing, the distance between the maximum forcing
and response mode locations significantly increases. For this sweep angle, the
region of forcing is centered at $z/c\approx 1$, while the response is
supported mostly at $z/c\approx 3$. As the peak $\sigma_{1}$ is lower for this
wing compared to lower sweep angle planforms, we argue that a significant
amount of energy is required for an external forcing to perturb the wakes of
highly swept wings. For all wings with $sAR=4$, this distance between the
dominant forcing-response mode pairs is strongly associated with the sweep
angle, while having a minor dependency on the angle of attack and presenting a
gradual decrease with the frequency.

The direction from forcing-to-response revealed by the optimal triglobal
resolvent modes suggests a spanwise advection of flow structures associated
with the sweep angle. As shown previously, we can relate the forcing-to-
response characteristics to the vortical fluctuations observed in the DNS. We
can further relate these findings to the modal convective speed from biglobal
stability analysis over swept wings (Crouch et al., 2019; Paladini et al.,
2019; Plante et al., 2021). Triglobal resolvent modes also reveal the
advection of perturbations over the wingspan related to the sweep angle, the
attenuation of flow unsteadiness, and the resilience to amplify perturbations
at high sweep angles. Even for unswept wings, the triglobal analysis uncovers
a preferential root-direction for advection of oscillations.

###  3.4 Influence of the aspect ratio

High sweep angle and low aspect ratio restrict the emergence of fluctuations
in flows over finite wings. As shown in figure 6, tip- and root-dominated
modes may extend over $1$ or $2$ chord-lengths over the wingspan. For this
reason, for flows over wings with $sAR<2$, the dominance of the global modes
may not be associated with root or tip regions, as they extend over the entire
wingspan.

$St$

$z/c$

$z/c$

$z/c$

$z/c$

$sAR=1$, $\Lambda=15^{\circ}$

$sAR=1$, $\Lambda=30^{\circ}$

$sAR=2$, $\Lambda=15^{\circ}$

$sAR=2$, $\Lambda=30^{\circ}$

Figure 7: Wingspan locations of dominant forcing (red) and response (blue)
with contours of $\mathbi{\Omega}_{\hat{\mathbi{f}}}$ and
$\mathbi{\Omega}_{\hat{\mathbi{q}}}\in[0.4,1.0]$ for $0.05\leq St\leq 0.35$
for wings at $\alpha=30^{\circ}$, $sAR=2$ and $1$, and $\Lambda=15^{\circ}$
and $30^{\circ}$.

For flows over $sAR=2$ wings, we observe a gradual transition between root-
dominated and tip-dominated forcing and response modes, as shown in figure 7.
For $\Lambda=15^{\circ}$, the optimal forcing-response mode pair appears near
the root for lower frequencies and at the wing tip for higher frequencies,
characterizing a root-to-tip mode switching. For $\Lambda=30^{\circ}$, the
trend is opposite, with wing tip modes at lower frequencies and root modes at
higher frequencies, characterizing a tip-to-root mode switching. These
features are similar to the mode switching observed for these sweep angles
with $sAR=4$, as shown in figure 6.

For wings with a low aspect ratio, the growth of root-dominated and tip-
dominated perturbations is constrained and mode switching does not occur for
$sAR=1$, as shown in figure 8. Distinguishing between root-dominated and tip-
dominated modes may be challenging for flows over $sAR=1$ wings as forcing and
response mode pairs appear globally, extending over the entire wingspan,
independent of the sweep angle. Therefore, flows around wings with $sAR=1$
tend to exhibit similar wake characteristics over different sweep angles.
Indeed, the wake patterns for flows over $sAR=1$ wings at a particular angle
of attack and sweep exhibit different characteristics from the flows over
higher aspect ratio wings, e.g., $sAR=2$ and $4$.

For high-aspect-ratio wings, for instance, the flow around
$(sAR,\alpha,\Lambda)=(4,30^{\circ},15^{\circ})$ wings, we observe in the DNS
that the wake shedding structures appear over the entire wingspan. The
resolvent modes depict these structures in three different flow mechanisms. As
shown in figure 8 (right) for $sAR=4$, there are two types of root-dominated
modes, which were also previously identified for this wing at
$\alpha=20^{\circ}$, shown in figure 4 (left). The first one is characterized
by root-dominated structures and appears at $St=0.15$, while the second type,
with a high $\sigma_{1}$, develops at $St=0.25$ with compact root-concentrated
modes. The third type is comprised of tip-dominated modes that become primary
as the frequency increases to $St=0.28$. These modes were primary at
$\alpha=20^{\circ}$ and $St\approx 0.40$, as shown in figure 4, although for
$\alpha=30^{\circ}$ they present a higher amplification gain.

For $sAR=2$, root-dominated modes are primary for $St<0.20$. Root-concentrated
modes are absent and tip-dominated modes are the primary perturbations for
$St\geq 0.20$, as shown in figure 8 (middle), characterizing a root-to-tip
mode switching. The overall mode switching for $sAR=2$ is the same, with root-
to-tip transition. For $sAR=1$, shown in figure 8 (left), mode switching is
absent. Both primary and secondary modes develop over the entire wingspan for
all frequencies, as shown for the primary modes at $St=0.14$. Although mode
switching is absent we can still reveal two distinct root- and tip-dominated
mechanisms on a single mode over low-aspect-ratio wings. For instance, at
$St=0.24$, modes emerge from the leading edge at the root and from the
trailing edge near the tip. Combined, these two types of flow unsteadiness
yield a global mode that appears over the entire wingspan.

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$\hat{f}_{1}$

$\hat{q}_{1}$

$sAR=1$

$sAR=2$

$sAR=4$

$\sigma$

$St$

$St$

$St$

Figure 8: Resolvent gain distribution and forcing-response mode pairs over
frequency for $(\alpha,\Lambda)=(30^{\circ},15^{\circ})$ and $1\leq sAR\leq
4$. Forcing ($\hat{f}$) and response ($\hat{q}$) modes shown with isosurfaces
of velocity $\hat{u}_{y}\in[-0.2,0.2]$. Mode switching is absent for $sAR=1$
due to merging of root and wingtip perturbations on the wake.

##  4 Conclusions

We presented the triglobal resolvent analysis of laminar separated flows over
swept wings and characterized the effects of wing tip and sweep angle on the
wake dynamics. We revealed the forcing and response structures that can be
amplified from harmonic oscillations or external actuation over finite wings.
In the present triglobal analysis, we have identified the wingspan locations
where forcing structures can be amplified near the wing and the regions where
the unsteady response develops. We have further characterized the region of
dominance of modal structures over the wingspan, with forcing-response mode
pairs appearing near the wing root or tip as a function of their
characteristic frequency.

Through resolvent wavemakers, we studied the steady to unsteady flow
characteristics over swept wings. We showed the regions where self-sustained
unsteadiness appear over swept wings and related those to the vortex shedding
structures observed in the DNS. We also revealed the most sensitive regions
for perturbation amplification in steady wakes over highly swept wings. The
forcing-response mode pairs further revealed the mechanisms of spanwise
advection of flow structures, which is further related to the spanwise
convective speed found for two- and three-dimensional resolvent analysis over
swept wings and also associated with the nonlinear flow characteristics
observed in DNS.

At last, we showed for low-aspect-ratio wings that localized perturbations
with root- or tip-dominant characteristics are limited as modes evolve
globally over the entire wingspan. In fact, we have shown that root-and tip-
dominant structures can appear over the wing in a single mode for low-aspect-
ratio wings. This behavior explains the characteristics of the laminar flows
around these wings, as observed in DNS, to be different from the flows over
$sAR=2$ and $4$ wings. These findings provide fundamental insights into future
studies on flow separation over swept wings at higher Reynolds numbers, in
which a wider spectrum of fluctuations is present.

##  Appendix A Choice of resolvent discounting parameter

$St_{r}$

$St_{i}$

$St_{i}$

$\alpha=20^{\circ}$

$\alpha=30^{\circ}$

Figure 9: Eigenvalues of $\mathsfbi{L}_{\bar{\mathbi{q}}}$ for $sAR=4$ wings
with $\alpha=20^{\circ}$ and $30^{\circ}$, and $0\leq\Lambda\leq 45^{\circ}$.
Here $St_{r}$ and $St_{i}$ represent the growth rate and temporal frequency,
respectively. Green-dashed line shows $St_{r}=0$. Cyan-dashed line shows where
discounting parameter is set.

Prior to performing the resolvent analysis, we examine the stability
characteristics of $\mathsfbi{L}_{\bar{\mathbi{q}}}$ and the need of a
discounting parameter. For this task we analyze the eigenspectrum of the
linear operator with respect to the time-average base flow. Since this base
flow is not the equilibrium point, the eigenvalues of
$\mathsfbi{L}_{\bar{\mathbi{q}}}$ do not necessarily present stability
characteristics. However, eigenvalue properties are needed to enable the
examination of the amplification mechanisms of perturbation over the
appropriate time scale. The eigenvalues of $\mathsfbi{L}_{\bar{\mathbi{q}}}$
are computed using the Krylov–Schur method (Stewart, 2002) with 128 vectors
for the Krylov subspace and a tolerance residual of $10^{-10}$. This analysis
reveals eigenvalues $-i\omega=-i\omega_{r}+\omega_{i}$, where $\omega_{i}$ is
the growth rate and $\omega_{r}$ is the temporal frequency. The Strouhal
number scaling (equation 7) is used throughout this study to report
$\omega_{i}$ and $\omega_{r}$ as $St_{i}$ and $St_{r}$, respectively, as shown
in figure 9. The discounting parameter is defined in a similar manner as
$St_{\beta}=(\beta/2\pi)(c\sin\alpha/U_{\infty}\cos\Lambda)$. In this manner,
$St_{\beta}$ is directly associated with a physical time window
$t_{\beta}=(2\pi/\beta)$, which is chosen to be shorter than the the time
scale associated with the largest $St_{r}$. As shown in figure 9 for flows
over $sAR=4$ wings with $\alpha=20^{\circ}$ and $30^{\circ}$, and
$0\leq\Lambda\leq 45^{\circ}$, modes which appear with positive $St_{r}$
(above green-dashed line) are so-called unstable.

Through the global stability analysis of all $\mathsfbi{L}_{\bar{\mathbi{q}}}$
operators examined in the present study, we observe that eigenvalues of the
unstable modes with higher growth rate lie near the stability margin. The
discounting parameter $St_{\beta}$ must be chosen such that
$St_{\beta}>\max(St_{r})$ for all linear operators considered. For
$St_{\beta}=0.01$, the resolvent discounting is able to encompass all the
unstable modes in the present study, as shown in figure 9 with a cyan-dashed
line. This discounting corresponds to a physical time window of
$t_{\beta}(U_{\infty}\cos\Lambda/c\sin\alpha)=100$.

##  Appendix B Convergence test for resolvent analysis

We document the randomized resolvent analysis computations in table 2, for a
selected case of $(sAR,\alpha,\Lambda)=(2,30^{\circ},15^{\circ})$. The
convergence of randomized SVD algorithm depends on the number of test vectors
$k$ used for sketching (Halko et al., 2011; Ribeiro et al., 2020). As our
discussions focus on the modal characteristics of the dominant and sub-
dominant resolvent modes, we have analyzed their convergence using $k=5$,
$10$, and $20$. The use of $k=10$ test vectors was shown to be sufficient to
guarantee converged leading modes with error smaller than $1\%$ for the
leading $5$ resolvent gains.

|  | $k=5$ |  | $k=10$ |  | $k=20$
---|---|---|---|---|---|---
$\sigma_{1}$ |  |  $3058.119140625$ |  |  $3058.1840820312$ |  |  $3058.0112304688$
$\sigma_{2}$ |  |  $395.9677734375$ |  |  $397.50375366211$ |  |  $397.54058837891$
$\sigma_{3}$ |  |  $283.04278564453$ |  |  $285.14834594727$ |  |  $285.30682373047$
$\sigma_{4}$ |  |  $172.99363708496$ |  |  $178.36595153809$ |  |  $179.13203430176$
$\sigma_{5}$ |  |  $131.17066955566$ |  |  $155.90252685547$ |  |  $156.57176208496$
Table 2: Convergence of the randomized resolvent analysis using $k=5$, $10$,
and $20$ test vectors shown for the five leading singular values with
$(sAR,\alpha,\Lambda)=(2,30^{\circ},15^{\circ})$ at $St=0.14$.

## Acknowledgments

We acknowledge support from the US Air Force Office of Scientific Research
(program manager: Dr. G. Abate, grant: FA9550-21-1-0174) and the US Army
Research Office (program managers: Drs. M. J. Munson and R. Anthenien, grant:
W911NF-21-1-0060). We thank C. S. Skene, T. R. Ricciardi, M. Amitay, V.
Theofilis, and P. J. Schmid for enlightening discussions. Computational
resources were provided by the High Performance Computing Modernization
Program at the US Department of Defense and the Texas Advanced Computing
Center.

## Declaration of interest

The authors report no conflict of interest.

## References

  * Amestoy et al. (2001) Amestoy, P.R., Duff, I. S., L’Excellent, J.-Y. & Koster, J. 2001 A fully asynchronous multifrontal solver using distributed dynamic scheduling. SIAM J. Mat. Anal. App. 23 (1), 15–41.
  * Anderson (2010) Anderson, J. D. 2010 Fundamentals of aerodynamics. McGraw-Hill.
  * Barkley et al. (2008) Barkley, D., Blackburn, H. M. & Sherwin, S. J. 2008 Direct optimal growth analysis for timesteppers. Int. J. Numer. Meth. Fluids 57 (9), 1435–1458.
  * Barthel et al. (2022) Barthel, B., Gomez, S. & McKeon, B. J. 2022 Variational formulation of resolvent analysis. Phys. Rev. Fluids 7 (1), 013905.
  * Brès et al. (2017) Brès, G. A., Ham, F. E., Nichols, J. W. & Lele, S. K. 2017 Unstructured large-eddy simulations of supersonic jets. AIAA J. 55 (4), 1164–1184.
  * Buchholz & Smits (2006) Buchholz, J. H. J. & Smits, A. J. 2006 On the evolution of the wake structure produced by a low-aspect-ratio pitching panel. J. Fluid Mech. 546, 433–443.
  * Burtsev et al. (2022) Burtsev, A., He, W., Hayostek, S., Zhang, K., Theofilis, V., Taira, K. & Amitay, M. 2022 Linear modal instabilities around post-stall swept finite wings at low Reynolds numbers. J. Fluid Mech. 944, A6.
  * Chu (1965) Chu, B.-T. 1965 On the energy transfer to small disturbances in fluid flow (part I). Acta Mechanica 1 (3), 215–234.
  * Crouch et al. (2019) Crouch, J. D., Garbaruk, A. & Strelets, M. 2019 Global instability in the onset of transonic-wing buffet. J. Fluid Mech. 881, 3–22.
  * Devenport et al. (1996) Devenport, W. J., Rife, M. C., Liapis, S. I. & Follin, G. J. 1996 The structure and development of a wing-tip vortex. J. Fluid Mech. 312, 67–106.
  * Edstrand et al. (2018a) Edstrand, A. M., Schmid, P. J., Taira, K. & Cattafesta III, L. N. 2018a A parallel stability analysis of a trailing vortex wake. J. Fluid Mech. 837, 858.
  * Edstrand et al. (2018b) Edstrand, A. M., Sun, Y., Schmid, P. J., Taira, K. & Cattafesta, L. N. 2018b Active attenuation of a trailing vortex inspired by a parabolized stability analysis. J. Fluid Mech. 855, R2.
  * Fage & Johansen (1927) Fage, A. & Johansen, F. C. 1927 On the flow of air behind an inclined flat plate of infinite span. Proc. R. Soc. Lond. A 116 (773), 170–197.
  * Farghadan et al. (2021) Farghadan, A., Towne, A., Martini, E. & Cavalieri, A. 2021 A randomized time-domain algorithm for efficiently computing resolvent modes. AIAA Paper 2021-2896.
  * Fosas de Pando & Schmid (2017) Fosas de Pando, M. & Schmid, P. J. 2017 Optimal frequency-response sensitivity of compressible flow over roughness elements. J. Turb. 18 (4), 338–351.
  * Fosas de Pando et al. (2017) Fosas de Pando, M., Schmid, P. J. & Sipp, D. 2017 On the receptivity of aerofoil tonal noise: an adjoint analysis. J. Fluid Mech. 812, 771–791.
  * Freund (1997) Freund, J. B. 1997 Proposed inflow/outflow boundary condition for direct computation of aerodynamic sound. AIAA J. 35 (4), 740–742.
  * Giannetti et al. (2010) Giannetti, F., Camarri, S. & Luchini, P. 2010 Structural sensitivity of the secondary instability in the wake of a circular cylinder. J. Fluid Mech. 651, 319–337.
  * Giannetti & Luchini (2007) Giannetti, F. & Luchini, P. 2007 Structural sensitivity of the first instability of the cylinder wake. J. Fluid Mech. 581 (1), 167–197.
  * Gómez et al. (2016) Gómez, F., Blackburn, H. M., Rudman, M., Sharma, A. S. & McKeon, B. J. 2016 A reduced-order model of three-dimensional unsteady flow in a cavity based on the resolvent operator. J. Fluid Mech. 798, R2.
  * Halko et al. (2011) Halko, N., Martinsson, P.-G. & Tropp, J. A. 2011 Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. SIAM review 53 (2), 217–288.
  * Harper & Maki (1964) Harper, C. W. & Maki, R. L. 1964 A review of the stall characteristics of swept wings. Tech. Rep. NASA/TN D-2373. NASA, Washington, DC.
  * He & Timme (2021) He, W. & Timme, S. 2021 Triglobal infinite-wing shock-buffet study. J. Fluid Mech. 925, A27.
  * Hill (1992) Hill, D. 1992 A theoretical approach for analyzing the restabilization of wakes. AIAA Paper 1992-67.
  * House et al. (2022) House, D. C., Skene, C. S., Ribeiro, J. H. M., Yeh, C.-A. & Taira, K. 2022 Sketch-based resolvent analysis. AIAA Paper 2022-3335.
  * Houtman et al. (2022) Houtman, J., Timme, S. & Sharma, A. 2022 Resolvent analysis of large aircraft wings in edge-of-the-envelope transonic flow. AIAA Paper 2022-1329.
  * Jovanović (2004) Jovanović, M. R. 2004 Modeling, analysis, and control of spatially distributed systems. University of California at Santa Barbara, Dept. of Mechanical Engineering.
  * Jovanović & Bamieh (2005) Jovanović, M. R. & Bamieh, B. 2005 Componentwise energy amplification in channel flows. J. Fluid Mech. 534, 145–183.
  * Khalighi et al. (2011) Khalighi, Y., Ham, F., Nichols, J., Lele, S. K. & Moin, P. 2011 Unstructured large eddy simulation for prediction of noise issued from turbulent jets in various configurations. AIAA Paper 2011–2886.
  * Lentink et al. (2007) Lentink, D., Müller, U. K., Stamhuis, E. J., De Kat, R., Van Gestel, W., Veldhuis, L. L. M., Henningsson, P., Hedenström, A., Videler, J. J. & Van Leeuwen, J. L. 2007 How swifts control their glide performance with morphing wings. Nature 446 (7139), 1082–1085.
  * Liu et al. (2021) Liu, Q., Sun, Y., Yeh, C.-A., Ukeiley, L. S., Cattafesta, L. N. & Taira, K. 2021 Unsteady control of supersonic turbulent cavity flow based on resolvent analysis. J. Fluid Mech. 925, A5.
  * Martini et al. (2021) Martini, E., Rodríguez, D., Towne, A. & Cavalieri, A. V. G. 2021 Efficient computation of global resolvent modes. J. Fluid Mech. 919, A3.
  * Masini et al. (2020) Masini, L., Timme, S. & Peace, A. J. 2020 Analysis of a civil aircraft wing transonic shock buffet experiment. J. Fluid Mech. 884, A1.
  * McKeon & Sharma (2010) McKeon, B. J. & Sharma, A. S. 2010 A critical-layer framework for turbulent pipe flow. J. Fluid Mech. 658, 336–382.
  * Moarref et al. (2013) Moarref, R., Sharma, A. S., Tropp, J. A. & McKeon, B. J. 2013 Model-based scaling of the streamwise energy density in high-reynolds-number turbulent channels. J. Fluid Mech. 734, 275–316.
  * Monokrousos et al. (2010) Monokrousos, A., Åkervik, E., Brandt, L. & Henningson, D. S. 2010 Global three-dimensional optimal disturbances in the blasius boundary-layer flow using time-steppers. J. Fluid Mech. 650, 181–214.
  * Paladini et al. (2019) Paladini, E., Beneddine, S., Dandois, J., Sipp, D. & Robinet, J.-C. 2019 Transonic buffet instability: From two-dimensional airfoils to three-dimensional swept wings. Phys. Rev. Fluids 4 (10), 103906.
  * Plante et al. (2021) Plante, F., Dandois, J., Beneddine, S., Laurendeau, É. & Sipp, D. 2021 Link between subsonic stall and transonic buffet on swept and unswept wings: from global stability analysis to nonlinear dynamics. J. Fluid Mech. 908, A16.
  * Qadri & Schmid (2017) Qadri, U. A. & Schmid, P. J. 2017 Frequency selection mechanisms in the flow of a laminar boundary layer over a shallow cavity. Phys. Rev. Fluids 2, 043902.
  * Ribeiro et al. (2020) Ribeiro, J. H. M., Yeh, C.-A. & Taira, K. 2020 Randomized resolvent analysis. Phys. Rev. Fluids 5 (3), 033902\.
  * Ribeiro et al. (2022a) Ribeiro, J. H. M., Yeh, C.-A., Zhang, K. & Taira, K. 2022a From biglobal to triglobal resolvent analysis: laminar separated flows over swept wings. AIAA Paper 2022-2428.
  * Ribeiro et al. (2022b) Ribeiro, J. H. M., Yeh, C.-A., Zhang, K. & Taira, K. 2022b Wing sweep effects on laminar separated flows. J. Fluid Mech. 950, A23.
  * Ricciardi et al. (2022) Ricciardi, T. R., Wolf, W. R. & Taira, K. 2022 Transition, intermittency and phase interference effects in airfoil secondary tones and acoustic feedback loop. J. Fluid Mech. 937.
  * Schmid & Brandt (2014) Schmid, P. J. & Brandt, L. 2014 Analysis of fluid systems: Stability, receptivity, sensitivity. Applied Mechanics Reviews 66 (2), 024803.
  * Schmidt et al. (2018) Schmidt, O. T., Towne, A., Rigas, G., Colonius, T. & Brès, G. A. 2018 Spectral analysis of jet turbulence. J. Fluid Mech. 855, 953–982.
  * Skene et al. (2022a) Skene, C. S., Ribeiro, J. H. M. & Taira, K. 2022a csskene/linear-analysis-tools: Initial release. https://doi.org/10.5281/zenodo.6550726.
  * Skene & Schmid (2019) Skene, C. S. & Schmid, P. J. 2019 Adjoint-based parametric sensitivity analysis for swirling M-flames. J. Fluid Mech. 859, 516–542.
  * Skene et al. (2022b) Skene, C. S., Yeh, C.-A., Schmid, P. J. & Taira, K. 2022b Sparsifying the resolvent forcing mode via gradient-based optimisation. J. Fluid Mech. 944, A52.
  * Stewart (2002) Stewart, G. W. 2002 A Krylov–Schur algorithm for large eigenproblems. SIAM J. Mat. Anal. App. 23 (3), 601–614.
  * Strykowski & Sreenivasan (1990) Strykowski, P. J. & Sreenivasan, K. R. 1990 On the formation and suppression of vortex ‘shedding’at low reynolds numbers. J. Fluid Mech. 218, 71–107.
  * Sun et al. (2017) Sun, Y., Taira, K., III, L. N. Cattafesta & Ukeiley, L. S. 2017 Biglobal instabilities of compressible open-cavity flows. J. Fluid Mech. 826, 270–301.
  * Taira et al. (2017) Taira, K., Brunton, S. L., Dawson, S. T. M., Rowley, C. W., Colonius, T., McKeon, B. J., Schmidt, O. T., Gordeyev, S., Theofilis, V. & Ukeiley, L. S. 2017 Modal analysis of fluid flows: An overview. AIAA J. 55 (12), 4013–4041.
  * Taira & Colonius (2009) Taira, K. & Colonius, T. 2009 Three-dimensional flows around low-aspect-ratio flat-plate wings at low Reynolds numbers. J. Fluid Mech. 623, 187–207.
  * Taira et al. (2020) Taira, K., Hemati, M. S., Brunton, S. L., Sun, Y., Duraisamy, K., Bagheri, S., Dawson, S. T. M. & Yeh, C.-A. 2020 Modal analysis of fluid flows: Applications and outlook. AIAA J. 58 (3), 998–1022.
  * Thomareis & Papadakis (2018) Thomareis, N. & Papadakis, G. 2018 Resolvent analysis of separated and attached flows around an airfoil at transitional reynolds number. Phys. Rev. Fluids 3 (7), 073901.
  * Timme (2020) Timme, S. 2020 Global instability of wing shock-buffet onset. J. Fluid Mech. 885, A37.
  * Torres & Mueller (2004) Torres, G. E. & Mueller, T. J. 2004 Low-aspect-ratio aerodynamics at low Reynolds numbers. AIAA J. 42 (5), 865–873.
  * Trefethen et al. (1993) Trefethen, L. N., Trefethen, A. E., Reddy, S. C. & Driscoll, T. A. 1993 Hydrodynamic stability without eigenvalues. Science 261 (5121), 578–584.
  * Videler et al. (2004) Videler, J. J., Stamhuis, E. J. & Povel, G. D. E. 2004 Leading-edge vortex lifts swifts. Science 306 (5703), 1960–1962.
  * Wygnanski et al. (2011) Wygnanski, I., Tewes, P., Kurz, H., Taubert, L. & Chen, C. 2011 The application of boundary layer independence principle to three-dimensional turbulent mixing layers. J. Fluid Mech. 675, 336–346.
  * Yeh et al. (2020) Yeh, C.-A., Benton, S. I., Taira, K. & Garmann, D. J. 2020 Resolvent analysis of an airfoil laminar separation bubble at Re = 500 000. Phys. Rev. Fluids 5 (8), 083906\.
  * Yeh & Taira (2019) Yeh, C.-A. & Taira, K. 2019 Resolvent-analysis-based design of airfoil separation control. J. Fluid Mech. 867, 572–610.
  * Yen & Hsu (2007) Yen, S.-C. & Hsu, C. M. 2007 Flow patterns and wake structure of a swept-back wing. AIAA J. 45 (1), 228–236.
  * Yen & Huang (2009) Yen, S.-C. & Huang, L.-C. 2009 Flow patterns and aerodynamic performance of unswept and swept-back wings. J. Fluids Eng. 131 (11).
  * Yilmaz & Rockwell (2012) Yilmaz, T. O. & Rockwell, D. 2012 Flow structure on finite-span wings due to pitch-up motion. J. Fluid Mech. 691, 518–545.
  * Zhang et al. (2020a) Zhang, K., Hayostek, S., Amitay, M., Burstev, A., Theofilis, V. & Taira, K. 2020a Laminar separated flows over finite-aspect-ratio swept wings. J. Fluid Mech. 905, R1.
  * Zhang et al. (2020b) Zhang, K., Hayostek, S., Amitay, M., He, W., Theofilis, V. & Taira, K. 2020b On the formation of three-dimensional separated flows over wings under tip effects. J. Fluid Mech. 895, A9.
  * Zhang & Taira (2022) Zhang, K. & Taira, K. 2022 Laminar vortex dynamics around forward-swept wings. Phys. Rev. Fluids 7 (2), 024704\.


# Partial Tensorized Transformers for Natural Language Processing

Subhadra Vadlamannati1, Ryan Solgi1
1Mercer Island High School, 9100 SE 42nd St, Mercer Island, USA
2Department of Electrical and Computer Engineering, University of California
Santa Barbara, Santa Barbara, USA
{subhadra<EMAIL_ADDRESS>

###### Abstract

The transformer architecture has revolutionized Natural Language Processing
(NLP) and other machine-learning tasks, due to its unprecedented accuracy.
However, their extensive memory and parameter requirements often hinder their
practical applications. In this work, we study the effect of tensor-train
decomposition to improve the accuracy and compress transformer vision-language
neural networks, namely BERT and ViT. We focus both on embedding-layer
compression and partial tensorization of neural networks (PTNN) through an
algorithmic approach. Our novel PTNN approach significantly improves the
accuracy of existing models by up to 5%, all without the need for post-
training adjustments, breaking new ground in the field of tensor
decomposition.

##  1 INTRODUCTION

The transformer architecture has been extremely influential in the field of
Natural Language Processing (NLP) by achieving state-of-the-art performance in
multiple downstream tasks [Brown et al., 2020]. Specifically, transformer-
based models have achieved state-of-the-art performance in machine
translation, sentence representation, and language generation [Brown et al.,
2020]. However, due to the multi-head self-attention layers present in such
architectures and the sheer number of parameters present in each layer
(upwards of a million parameters during each forward pass of the model),
transformers are computationally demanding in terms of both memory and
computation time [Gu et al., 2022].

Effective network compression for convolutional neural networks with fully-
connected networks is still an open research question. Therefore, compressing
Transformer networks is essential for their use on edge devices and in low-
resource environments. Tensor Decomposition works by representing higher-order
tensors as a series of low-rank factors[Oseledets, 2011].

Tensor Decomposition (TD) works by representing higher-order tensors as a
series of simpler, lower-order tensors via elementary operations [Oseledets,
2011]. There are a variety of methods proposed to express tensors in a low-
rank format factored form, including CANDECOMP PARAFAC (CP), tucker, and
tensor-train decomposition [Bacciu and Mandic, 2020].

Several studies have explored the application of tensor decomposition methods
to reduce the size of neural networks and improve their efficiency. For
example, CP decomposition has been utilized to factorize weight matrices as
the sum of rank-1 tensors, reducing the model’s memory requirements and
speeding up inference time 8x [Maldonado-Chan et al., 2021].

Techniques like tensorized training can be applied to neural networks to
achieve both faster model inference and also a large reduction in the number
of parameters. More recently, tensorized training approaches have explored the
possibility of updating the factored tensors of the weights of deep neural
networks during the forward pass instead of the model weights directly [Bacciu
and Mandic, 2020].

In the field of natural language processing, while tensor decomposition
methods have been extensively studied for text-based neural networks, their
application to multimodal neural networks has received limited attention.
Multimodal neural networks such as ALIGN [Jia et al., 2021b] and VisualBERT
[Li et al., 2019], have gained significant attention in recent years due to
their remarkable ability to process both text and image data. These dual-
encoder architectures require a larger number of parameters compared to
traditional text-based models due to their need to capture information from
many modalities.

While there has been existing research on applying traditional tensor
decomposition methods to reduce the size of neural networks, their application
on multimodal neural networks has been severely limited [Maldonado-Chan et
al., 2021]. Existing methods have demonstrated that tensor decomposition
methods can effectively reduce memory time and increase efficiency on neural
networks while maintaining minimal accuracy loss (Zhong 2019; Bacciu and
Mandic 2020).

In this paper, we focus on vision-language transformer-based models as prime
candidates for models that require compression, as these models employ a dual-
encoder architecture that requires double the amount of parameters traditional
text-based models use. These models must encode both text and image data and
then generate a latent space that can be employed on a variety of downstream
tasks including visual question-answering, visual common-sense reasoning,
natural language for visual reasoning, and region-to-phrase grounding. Vision-
language tasks have a wide range of practical applications, including closed-
captioning for the hard-of-hearing and screen reading capabilities [Morris et
al., 2018].

This paper aims to advance the field of multimodal language processing by
applying state-of-the-art tensor decomposition techniques to compress popular
networks such as ViT and BERT while maintaining model accuracy. We
additionally study the effect of tensor size and hyperparameter values during
tensorized training.

Figure 1: Bimodal Space-saving distribution for the two models.

##  2 METHODS

###  2.1 Data

We use the ImageNet dataset to train ViT. ImageNet is a collection of 80
million tiny images for image classification. Instead of using the entire
dataset, we use the standard 1000 image test set of ImageNet known as CIFAR10
[Krizhevsky, 2009] to evaluate the compressed model.

###  2.2 Tensor-Train Decomposition

We use tensor-train (TT) decomposition throughout this paper to compress
weight tensors. In the tensor train (TT) format [Oseledets, 2011], a $d$-way
tensor $\textbf{$\mathcal{W}$}\in\mathbb{R}^{n_{1}\times....\times n_{d}}$ is
approximated with a set of $d$ cores
$\bar{\textbf{$\mathcal{G}$}}=\\{\textbf{$\mathcal{G}$}_{1},\textbf{$\mathcal{G}$}_{2},...,\textbf{$\mathcal{G}$}_{d}\\}$
where $\textbf{$\mathcal{G}$}_{j}\in\mathbb{R}^{r_{j-1}\times n_{j}\times
r_{j}}$ , $r_{j}$’s for $j=1,...,d-1$ are the ranks, $r_{0}=r_{d}=1$, and each
element of $\mathcal{Y}$ is approximated by the following equation:

$\displaystyle\hat{\textbf{{$\mathcal{W}$}}}[i_{1},...,i_{d}]=\sum_{l_{0},...,l_{d}}{\textbf{$\mathcal{G}$}_{1}[l_{0},i_{1},l_{1}]...\textbf{$\mathcal{G}$}_{d}[l_{d-1},i_{d},l_{d}]}.$

With a prescribed error tolerance ($\epsilon$), the fundamental components,
denoted as $\mathcal{G}_{j}$, are obtained through a sequence of ($d-1$)
consecutive Singular Value Decompositions (SVDs) performed on auxiliary
matrices derived from the unfolding of the tensor $\mathcal{Y}$ along various
axes. This decomposition method, referred to as TT-SVD, is detailed in
Algorithm 1.

Input: $d$-way tensor $\mathcal{Y}$, error bound $\epsilon$.

Output:
$\bar{\textbf{$\mathcal{G}$}}=\\{\textbf{$\mathcal{G}$}_{1},\textbf{$\mathcal{G}$}_{2},...,\textbf{$\mathcal{G}$}_{d}\\}$

$\sigma=\frac{\epsilon}{d-1}\|\textbf{$\mathcal{Y}$}\|_{F}$

$r_{0}=1$, $r_{d}=1$, $\textbf{W}=$
reshape$(\mathcal{Y},(n_{1},\frac{|\mathcal{Y}|}{n_{1}}))$

for  _$j=1$ to $j=d-1$_ do

$\textbf{W}=$
reshape$(\textbf{W},(r_{j-1}n_{j},\frac{|\textbf{W}|}{r_{j-1}n_{j}}))$

Compute $\sigma$-truncated SVD: $\textbf{W}=\textbf{USV}^{T}+\textbf{E}$,
where $\|\textbf{E}\|_{F}\leq\sigma$

$r_{j}$ = the rank of matrix W based on $\sigma$-truncated SVD

$\textbf{$\mathcal{G}$}_{j}=$ reshape$(\textbf{U},(r_{j-1},n_{j},r_{j}))$

$\textbf{W}=\textbf{SV}^{T}$

end for

$\textbf{$\mathcal{G}$}_{d}=$ reshape$(\textbf{W},(r_{d-1},n_{d},r_{d}))$

Return
$\bar{\textbf{$\mathcal{G}$}}=\\{\textbf{$\mathcal{G}$}_{1},\textbf{$\mathcal{G}$}_{2},...,\textbf{$\mathcal{G}$}_{d}\\}$

Algorithm 1 TT-SVD

###  2.3 Models and Pretraining

####  2.3.1 Models

We perform an analysis on the effects of various models to encode both visual
and textual data. For the text encoder, we focus on BERT [Liu et al., 2019]
which applies the bi-directional training of the Transformer to create in-
context word embeddings. To test BERT on image classification, we pair it with
EfficientNet for the visual encoding, which is consistent with the
implementation in the dual-encoder network ALIGN [Jia et al., 2021a]. BERT has
12 layers, a hidden layer size of 768, and around 120 million parameters.

For the visual encoder we focus on ViT [Dosovitskiy et al., 2021]. ViT is an
image processor based on the Transformer architecture.

####  2.3.2 Pre-Compression Training

We fine-tune the image encoder ViT as it’s trained on a much larger, more
general dataset ImageNet. Additionally, this fine-tuning improves baseline
model accuracy which ensures that our post-training adjustments truly improve
model accuracy. Specifically, we fine-tune ViT on the CIFAR10 train set with
50,000 images and hyperparameters as outlined in Section 4.1. Post fine-
tuning, we achieve an accuracy of 97.89%, a validation loss of 0.2564 and a
training loss of 0.4291.

###  2.4 Embedding Layer Tensorization

Due to the relatively large size of the embedding layer, comprising around 28%
of BERT’s total parameters, and around 40% of ViT’s total parameters, we first
explore the effects of embedding layer tensorization. Following the
aforementioned approach in 2.2, we compress the embedding layers of the model
and study the effects on model accuracy.

###  2.5 Partially Tensorized Neural Network (PTNN)

As the compression of the embedding layer caused the accuracy of the network
to increase, as mentioned in 3.0.1, we continue exploring this accuracy
improvement for the rest of the model using our method in 2.4. In this case,
we do not perform post-training adjustments (like retraining). Initially, we
extract model weights from all relevant encoders. We compress every layer by
initially transforming the weight’s 2D matrix into a tensor of higher
dimension with the same volume as the original matrix. Tensor-train
decomposition is then used to decompose model weights with a given error
bound. Lastly, we reshape the tensor back into the original dimensions of the
2D weight matrix. To perform this iteratively as described in Algorithm 2, we
re-initialize the model with the decomposed weights and perform the same
process to compress the remaining model weights.

while  _not at the end of the model_ do

compress layers $0$ to $n$ using TT decomp;

if  _model accuracy $\geq$ 5% of the original accuracy_ then

proceed to layer $n+1$;

else

leave layer $n$ uncompressed;

move to layer $n+1$;

end if

end while

Algorithm 2 Iterative model tensorization. Figure 2: BERT and ViT model
iterative compression with exclusion of low-accuracy layers. Lighter color
corresponding dashed lines indicate baseline accuracy.

###  2.6 Decomposition Evaluation

We evaluate the effectiveness of the decomposition using commonly known
metrics like compression ratio and space saving, as described in Equation (1).

We additionally evaluate the memory saved and time complexity of our
decomposed models. For memory reduction, we simply multiply the percentage
space saving with the total number of parameters in the non-decomposed matrix.
Then, we divide this over the total number of parameters in the model.

For time complexity, we use a standard approach to find the O(n) based on the
rank of the Tensor-Train decomposed matrices as outlined in [Oseledets, 2011].

##  3 RESULTS

We evaluate both models on both the train (50,000) and test (10,000) image
split of the CIFAR10 dataset. The ViT model has 98 weight tensors while BERT
contains 99 weight tensors. All layers are compressed with an epsilon of 0.5
and 3D projections are constructed with lowest norms.

####  3.0.1 Embedding Layer Decomposition

Model | Train | Test | Space-saving
---|---|---|---
Decomposed ViT | 0.986 | 0.966 | 0.267
Base ViT | 0.960 | 0.960 | NA
Tensorized BERT (ALIGN) | 0.767 | 0.833 | 0.269
Base BERT (ALIGN) | 0.750 | 0.746 | NA

Table 1: Embedding Layer Decomposition Results

The results of embedding layer decomposition are presented in Table 1.
Embedding layer decomposition results in an accuracy increase of around 2% and
still provides us with a noticeable parameter reduction. Even though we are
not able to achieve high space-saving on these layers (i.e. only around 26%),
the success of this approach in maintaining or even improving accuracy leads
us to apply this technique to compress more layers of the model.

As seen in 3.0.1, decomposing the embedding layer of BERT provides an accuracy
gain over baseline comparisons. Following these promising results, we follow
the methodology in Algorithm 2 to determine the amount of BERT that can be
compressed whilst still maintaining overall accuracy.

####  3.0.2 Iterative and Individual Decomposition

We first individually compressed ViT’s weight tensors, but did not saw
decreases in model accuracy after the embedding layer (results found in
APPENDIX). Therefore, individual decomposition was not tested extensively on
BERT due to computational limitations. Instead, we rather focused on the
synergistic effects between layers of the model as mentioned in 2.5.

After compressing both models utilizing the process described in Algorithm 2,
the results are seen in 2. In general, we see accuracy improvements even while
compressing around half of the model’s parameters (53% in BERT and around 49%
in ViT). This is achieved all without post-training adjustments. Most notably,
more layers in BERT tend to have an improved accuracy, however in ViT,
compression of many layers does not change the baseline accuracy, while still
resulting in a parameter reduction.

TT-decomposition can be an effective way to improve the accuracy of
transformer-based models. While some layers vastly decrease the accuracy of
the model, we can still achieve a high compression even when not including
such layers in the overall decomposition.

##  4 CONCLUSIONS

This paper introduces a novel approach to compressing transformer-based
vision-language neural networks using TensorTrain decomposition. By applying
this method to popular models like ViT and BERT, we achieved significant
improvements in accuracy, up to 5%, without the need for post-training
adjustments. The iterative compression of model layers, coupled with
retraining, enables us to preserve model accuracy while reducing up to 53% of
the model’s parameters. In the future, we would like to generalize our
approach to other dual encoder models and test our approach on other
multimodal tasks like visual question answering and caption generation. This
work represents a valuable advancement in the field of multimodal language
processing and contributes to our broader goal of making transformer-based
models more efficient and practical for real-world applications.

## ACKNOWLEDGEMENTS

The first author would like to extend her gratitude to Dr. Lina Kim and the
Research Mentorship Program Cohort for their support throughout the research
process.

## REFERENCES

  * Bacciu and Mandic, 2020 Bacciu, D. and Mandic, D. P. (2020).  Tensor decompositions in deep learning.  Computational Intelligence.
  * Brown et al., 2020 Brown, T. B., Mann, B., Ryder, N., Subbiah, M., Kaplan, J., Dhariwal, P., Neelakantan, A., Shyam, P., Sastry, G., Askell, A., Agarwal, S., Herbert-Voss, A., Krueger, G., Henighan, T., Child, R., Ramesh, A., Ziegler, D. M., Wu, J., Winter, C., and Amodei, D. (2020).  Language models are few-shot learners.  arXiv.
  * Dosovitskiy et al., 2021 Dosovitskiy, A., Beyer, L., Kolesnikov, A., Weissenborn, D., Zhai, X., Unterthiner, T., Dehghani, M., Minderer, M., Heigold, G., Gelly, S., Uszkoreit, J., and Houlsby, N. (2021).  An image is worth 16x16 words: Transformers for image recognition at scale.  arXiv.
  * Gu et al., 2022 Gu, J., Keller, B., Kossaifi, J., Anandkumar, A., Khailany, B., and Pan, D. Z. (2022).  Heat: Hardware-efficient automatic tensor decomposition for transformer compression.  arXiv.
  * Jia et al., 2021a Jia, C., Yang, Y., Xia, Y., Chen, Y., Parekh, Z., Pham, H., Le, Q. V., Sung, Y., Li, Z., and Duerig, T. (2021a).  Scaling up visual and vision-language representation learning with noisy text supervision.  CoRR, abs/2102.05918.
  * Jia et al., 2021b Jia, C., Yang, Y., Xia, Y., Chen, Y.-T., Parekh, Z., Pham, H., Le, Q. V., Sung, Y., Li, Z., and Duerig, T. (2021b).  Scaling up visual and vision-language representation learning with noisy text supervision.  arXiv.
  * Krizhevsky, 2009 Krizhevsky, A. (2009).  Learning multiple layers of features from tiny images.
  * Li et al., 2019 Li, L. H., Yatskar, M., Yin, D., Hsieh, C.-J., and Chang, K.-W. (2019).  Visualbert: A simple and performant baseline for vision and language.  arXiv.
  * Liu et al., 2019 Liu, Y., Ott, M., Goyal, N., Du, J., Joshi, M., Chen, D., Levy, O., Lewis, M., Zettlemoyer, L., and Stoyanov, V. (2019).  Roberta: A robustly optimized bert pretraining approach.  arXiv.
  * Maldonado-Chan et al., 2021 Maldonado-Chan, M., Mendez-Vazquez, A., and Guardado-Medina, R. O. (2021).  Multimodal tucker decomposition for gated rbm inference.  Applied Sciences, 11(16):16.
  * Morris et al., 2018 Morris, M. R., Johnson, J., Bennett, C. L., and Cutrell, E. (2018).  Rich representations of visual content for screen reader users.  In Proceedings of the 2018 CHI Conference on Human Factors in Computing Systems, pages 1–11.
  * Oseledets, 2011 Oseledets, I. V. (2011).  Tensor-train decomposition.  SIAM Journal on Scientific Computing, 33(5):2295–2317.

## APPENDIX

###  4.1 ViT Fine-tuning Hyperparameters

    Learning Rate: 5
    Train/Evaluation Batch Size: 32
    Optimizer: Adam, linear learning rate
    Epoch(s): 1

###  4.2 Individual Compression

Figure 3: ViT Individual Compression