# Real-time stokes polarimetry based on a polarisation camera

Mitchell A. Cox  School of Electrical and Information Engineering, University
of the Witwatersrand, Johannesburg, South Africa  Carmelo Rosales-Guzmán
Centro de Investigaciones en Óptica, A. C., Loma del Bosque 115, Col. Lomas
del Campestre, 37150, León, Gto., México<EMAIL_ADDRESS>https://www.wits.ac.za/oclab

###### Abstract

This lab note introduces the “Stokes Camera,” a simple and novel experimental
arrangement for real-time measurement of spatial amplitude and polarisation
and thus spatially resolved Stokes parameters. It uses a polarisation
sensitive camera and a fixed quarter-wave plate, providing a one-shot, digital
solution for polarisation measurement that is only limited by the frame rate
of the camera and the computation speed of the provided code. The note also
provides background information on relevant polarisation theory and vector
vortex beams, which are used as a demonstration of the device.

##  I Introduction

Structured light research involves the use of optical modes with diverse
spatial amplitude, phase, and polarization profiles. The creation and
detection of such modes are fundamental to the field, and various experimental
techniques (often with the goal of being simpler, more accurate or more
sensitive) have been developed to accomplish these tasks Pinnell  _et al._
(2020); Tidwell  _et al._ (1990); Niziev  _et al._ (2006); Passilly  _et al._
(2005); Mendoza-Hernández  _et al._ (2019); Marrucci  _et al._ (2006); Naidoo
_et al._ (2016); Devlin  _et al._ (2017); Cox  _et al._ (2022).

In the last decade or two, the use of liquid crystal spatial light modulators
has transformed the field of structured light because of the ease and
efficiency of which the spatial amplitude and phase can be modulated Rosales-
Guzmán and Forbes (2017); Davis  _et al._ (2000); Maurer  _et al._ (2007);
Moreno  _et al._ (2012); Mitchell  _et al._ (2017); Rosales-Guzmán  _et al._
(2017); Rong  _et al._ (2014); Liu  _et al._ (2018a). The polarisation state
of a laser beam is an essential property in various scientific and engineering
applications, and control of this degree of freedom has opened many new
avenues in structured light research in the form of so-called vector modes
Zhan (2009); Rosales-Guzmán  _et al._ (2018); Yao-Li  _et al._ (2020); Bo  _et
al._ (2022); Rosales-Guzmán  _et al._ (2022); Hu  _et al._ (2021a), which are
useful in optical trappingYang  _et al._ (2021); Michihata  _et al._ (2009);
Kozawa and Sato (2010); Donato  _et al._ (2012); Zhao  _et al._ (2005),
microscopyDonato  _et al._ (2021); Hao  _et al._ (2010); Török and Munro
(2004); Segawa  _et al._ (2014a, b), metrologyNdagano  _et al._ (2017); Hu
_et al._ (2019); Töppel  _et al._ (2014); Berg-Johansen  _et al._ (2015) and
even optical communications Trichili  _et al._ (2020); Cox  _et al._ (2021);
Barreiro  _et al._ (2008); D’Ambrosio  _et al._ (2016); Otte  _et al._ (2020);
Ndagano  _et al._ (2018); Milione  _et al._ (2015a, b); Wang and Song (2017);
Liu  _et al._ (2018b); Willner (2018)

Control of the spatial polarisation of a beam is significantly more
challenging but recent progress is in part due to renewed interest in digital
micro-mirror devices (DMDs) Lerner  _et al._ (2012); Ren  _et al._ (2015);
Mitchell  _et al._ (2016); Gong  _et al._ (2014); Manthalkar  _et al._ (2020);
Zhao  _et al._ (2020); Cox and Drozdov (2021). This renewed interest in DMDs
is likely due to the ease of which the polarisation degree of freedom can be
used and manipulated Rosales-Guzmán  _et al._ (2020); Hu  _et al._ (2021b);
Scholes  _et al._ (2019); Zhao  _et al._ (2019); Hu and Rosales-Guzmán (2022).

Stokes polarimetry is a well-established technique in optics that has been
used for over a century to study the polarisation properties of light. The
technique was first introduced by Sir George Stokes in 1851 in his study of
the composition of light. Since then, it has become an essential tool in many
fields of optics, including laser physics, spectroscopy, and microscopy, among
others.

Stokes parameters can be measured in various ways but the simplest and most
intuitive is called the rotating wave-plate method Toninelli  _et al._ (2019).
This method is based on the principle of passing the laser beam through a
wave-plate at a specific angle and rotating the wave-plate. As the wave-plate
rotates, the polarisation state of the laser beam changes, and the intensity
of the laser beam is measured at several angles using a photo-detector. If a
spatially resolved measurement is desired (for example in structured light
applications) then a camera can be used. By analysing the measured
intensities, the Stokes parameters can be calculated. The rotating wave-plate
method is relatively simple, low cost, and provides high accuracy. However,
this method requires mechanical rotation, which can be challenging to
implement in some applications. An alternative method uses DMDs to “digitally”
measure the Stokes parameters of a beam, i.e. without mechanically moving
parts (such as wave plates) Zhao  _et al._ (2019); Manthalkar  _et al._
(2020).

In this lab note, we present a simple yet novel static experimental
arrangement and associated MATLAB code which we call a “Stokes Camera”. The
Stokes Camera enables full real-time measurement of spatial polarisation and
thus spatially resolved Stokes parameters. This arrangement employs a
polarisation sensitive camera coupled with a fixed quarter-wave plate and
software processing, offering a simple yet effective solution for one-shot,
digital polarization measurement.

To make this lab note as useful as possible to the reader, we provide an
expanded background in Sec. II on relevant polarisation theory as well as
vector vortex beams (a form of structured light making use of all three
degrees of freedom) which we use to demonstrate the operation of the Stokes
Camera. Boradly speaking, we present this background in a “top-down” manner.
We describe the experimental arrangement and required processing code in Sec.
III. In Sec. IV we provide several examples of output of the provided codes,
which are open source and available on Github Cox and Rosales-Guzmán (2023).

##  II Background

###  II.1 The Polarisation Ellipse and the Poincaré Sphere

In this section a summary of the various techniques and associated theory will
be provided in a tutorial-like manner. For the readers convenience as a lab
note we provide more detail than is strictly required to describe our “Stokes
Camera”. A comprehensive treatment of all the topics and theory herein would
be formidable, and so, the interested reader is invited to read more in the
references provided. We begin by writing the polarisation components of an
electric field in the form of column vector as,

$\mathbf{E}=\begin{bmatrix}E_{x}\\\ E_{y}\end{bmatrix},$ (1)

where $E_{x}$ is conventionally the complex amplitude (i.e. it has a magnitude
and phase) of the field on the horizontal axis and $E_{y}$ on the vertical. In
this form, this vector is known as the Jones vector and it can be conveniently
transformed using various Jones matrices which describe different optical
elements such as wave-plates and polarisers, thus allowing computation of the
resulting polarisation of a beam, discussed further in Sec. II.4. As a side
note, two polarisation vectors are orthogonal if the inner product between
them is zero.

Visualisation of the polarisation of a field is highly convenient, as dealing
purely with complex numbers is not always intuitive to the lay scientist or
engineer. Since polarisation can be linear, circular or somewhere in between
(elliptical), a useful representation of each polarisation state is through
the polarisation ellipse, shown in Fig. 1 (a), which captures the essence of
each polarisation state, the mathematical formalism of this will be derived
below. A more useful representation relies on representing each of these
ellipses in three dimensions using the surface of a unitary sphere, known as
the Poincaré sphere, which shown in Fig. 1 (b), in honour of Henri Poincaré,
the inventor of this representation Poincaré’ (1892). In this representation,
the right- and left-handed polarisation states are mapped to the north and
south pole, respectively, the linear polarisation states along the equator,
whereas all the elliptical states on the rest of the sphere. Representing each
polarisation state through a specific polarisation can be specially useful to
represent structured light beams with none homogeneous polarisation. More
importantly such polarisation structure can be obtained from intensirty
images, as will be explained below. These ellipses are often overlaid on the
intensity profile of a beam that have spatial interest - we call this a Stokes
plot for brevity.

Figure 1: (a) The Poincaré sphere showing all the possible states of
polarisation of light indicated in Dirac notation with relevant Stokes
parameters ($S_{1,2,3}$) also visualised. An alternative visualisation is the
polarisation ellipse in (b).

To understand how polarisation can be represented either as an ellipse or on
the surface of a sphere, it is convenient to represent each component of the
electric field as a periodic function of the form

$\begin{split}&E_{x}(t)=E_{0x}\cos(\tau+\delta_{x})\\\
&E_{y}(t)=E_{0y}\cos(\tau+\delta_{y})\end{split}$ (2)

where $\tau=\omega t-kz$. In what follows we will perform a series of
mathematical steps to derive the polarisation ellipse equation. To begin with,
we can apply the trigonometric identity
$\cos{(\theta_{1}+\theta_{2})}=\cos\theta_{1}\cos\theta_{2}-\sin\theta_{1}\sin\theta_{2}$
to both expressions in Eq. LABEL:eq:EF and divide by $E_{0x}$ and $E_{0y}$,
respectively, to obtain,

$\begin{split}&\frac{E_{x}(t)}{E_{0x}}=\cos\tau\cos\delta_{x}-\sin\tau\sin\delta_{x}\\\
&\frac{E_{y}(t)}{E_{0y}}=\cos\tau\cos\delta_{y}-\sin\tau\sin\delta_{y}.\end{split}$
(3)

After multiplying the first equation by $\sin\delta_{y}$ and the second by
$\sin\delta_{x}$ and subtracting the second from the first we obtain

$\frac{E_{x}(t)}{E_{0x}}\sin\delta_{y}-\frac{E_{y}(t)}{E_{0y}}\sin\delta_{x}=\cos\tau\sin(\delta_{x}-\delta_{y}).$
(4)

Squaring both sides of the equation then yields

$\begin{split}&\frac{E^{2}_{x}(t)}{E^{2}_{0x}}\sin^{2}\delta_{y}-2\frac{E_{x}(t)}{E_{0x}}\frac{E_{y}(t)}{E_{0y}}\sin\delta_{x}\sin\delta_{y}+\\\
&\frac{E^{2}_{y}(t)}{E^{2}_{0y}}\sin^{2}\delta_{x}=\cos^{2}\tau\sin^{2}(\delta_{x}-\delta_{y}).\end{split}$
(5)

In a similar way, we can now multiply Eq. 3 by $\cos\delta_{x}$ and
$\cos\delta_{y}$, subtract both equations and square the result to obtain

$\begin{split}&\frac{E^{2}_{x}(t)}{E^{2}_{0x}}\cos^{2}\delta_{y}-2\frac{E_{x}(t)}{E_{0x}}\frac{E_{y}(t)}{E_{0y}}\cos\delta_{x}\cos\delta_{y}+\\\
&\frac{E^{2}_{y}(t)}{E^{2}_{0y}}\cos^{2}\delta_{x}=\sin^{2}\tau\sin^{2}(\delta_{x}-\delta_{y}).\end{split}$
(6)

Finally, adding both equations yields the final result:

$\frac{E^{2}_{x}(t)}{E^{2}_{0x}}+\frac{E^{2}_{y}(t)}{E^{2}_{0y}}-2\frac{E_{x}(t)}{E_{0x}}\frac{E_{y}(t)}{E_{0y}}\cos\delta=\sin^{2}\delta,$
(7)

where, $\delta=\delta_{x}-\delta_{y}$. This equation is known as the
polarisation ellipse because it has the general form of a rotated ellipse, as
illustrate in Fig. 1, where $A$ and $B$ represent the major and minor axis and
$\alpha$ the angle of the major axis with respect to the component $E_{x}$ of
the electric field. From this equation any polarisation state can be obtained,
for example, for $\delta=0,\pi$ the above equation transforms into,

$E_{x}(t)=\pm\frac{E_{0x}(t)}{E_{0y}(t)}E_{y}(t),$ (8)

which is the equation of a line with slope $\pm E_{0x}(t)/E_{0y}(t)$, and
therefore, it is associated to linear polarisation states. As another example,
for $\delta=\pi/2,3\pi/2$, Eq. 7 transforms into,

$\frac{E^{2}_{x}(t)}{E^{2}_{0x}}+\frac{E^{2}_{y}(t)}{E^{2}_{0y}}=1,$ (9)

which is the equation of a circle and therefore it describes circular
polarisation states. The equation of the polarisation ellipse (Eq. 7) can be
transformed into its standard form (not rotated), from which we can obtain the
relations,

$\begin{split}&\tan(2\alpha)=\tan(2\phi)\cos\delta\qquad 0\leq\alpha\leq\pi\\\
&\sin(2\chi)=\sin(2\phi)\sin\delta,\qquad-\pi/4\leq\chi\leq\pi/4\end{split}$
(10)

from which, any polarisation state can be obtained in terms of the the
parameters $\alpha$ and $\chi$, which are associated to the orientation angle
(including direction) and ellipticity angle of polarisation, respectively
Goldstein (2011). Here, the parameter $\phi$ is related to the polarisation
components of the electric field through the relation,

$\tan(\phi)=\frac{E_{0y}}{E_{0x}}\qquad 0\leq\phi\leq\pi/2.$ (11)

Notice that, for $-\pi/4\leq\chi\leq 0$ (see Eq. 10) we have left handed
circular polarisation, whereas for $0\leq\chi\leq\pi/4$ right handed. As a
final comment, the major and minor axis of the polarisation ellipse, A and B,
respectively, are related to the ellipticity angle $\chi$ through the
relation,

$\tan(\chi)=\frac{\pm B}{A}.$ (12)

###  II.2 The Stokes parameters

In the polarisation ellipse equation, $E_{x}$ and $E_{y}$ are both implicitly
dependent on time, hence, it predicts a temporal evolution of polarisation,
but the fast oscillations of the electromagnetic field does not make possible
an experimental observation. Hence, in order to observe the effect of
polarisation, we must take an average over the time of observation of the
previous equation,

$\frac{\langle E^{2}_{x}(t)\rangle}{E^{2}_{0x}}+\frac{\langle
E^{2}_{y}(t)\rangle}{E^{2}_{0y}}-2\frac{\langle
E_{x}(t)E_{y}(t)\rangle}{E_{0xE_{0y}}}\cos\delta=\sin^{2}\delta,$ (13)

Using the time average for each term, whic are given by,

$\begin{split}&\langle E^{2}_{x}(t)\rangle=\frac{E^{2}_{0x}}{2}\\\ &\langle
E^{2}_{y}(t)\rangle=\frac{E^{2}_{0y}}{2}\\\ &\langle
E_{x}(t)E_{y}(t)\rangle=\frac{1}{2}E_{0x}E_{0y},\end{split}$ (14)

the previous equation can be written as

$\begin{split}&(E^{2}_{0x}+E^{2}_{0y})^{2}-(E^{2}_{0x}-E^{2}_{0y})^{2}-(2E_{0x}E_{0y}\cos\delta)^{2}=\\\
&(2E_{0x}E_{0y}\sin\delta)^{2},\end{split}$ (15)

At this point, it is useful to introduce the following definitions,

$\begin{split}&S_{0}=E^{2}_{0x}+E^{2}_{0y},\hskip
28.45274ptS_{2}=2E_{0x}E_{0y}\cos\delta\\\ &S_{1}=E^{2}_{0x}-E^{2}_{0y},\hskip
22.76219ptS_{3}=2E_{0x}E_{0y}\sin\delta,\end{split}$ (16)

which are known as the Stokes parameters. In terms of these parameters Eq. 15
can be written as,

$S_{0}^{2}=S_{1}^{2}+S_{2}^{2}+S_{3}^{2}.$ (17)

The Stokes parameters can also be written using a complex representation for
the electric field as,

$\begin{split}&S_{0}=E_{x}E^{*}_{x}+E_{y}E^{*}_{y},\hskip
28.45274ptS_{2}=E_{x}E^{*}_{y}+E_{y}E^{*}_{x}\\\
&S_{1}=E_{x}E^{*}_{x}-E_{y}E^{*}_{y},\hskip
22.76219ptS_{3}=i(E_{x}E^{*}_{y}-E_{y}E^{*}_{x}),\end{split}$ (18)

where, $E^{*}_{x}$ and $E^{*}_{y}$ represent the complex conjugated of $E_{x}$
and $E_{y}$, respectively. A more useful representation based only on
intensity measurements will be derived below, since this allows a
reconstruction of any polarisation state.

In summary, $S_{0}$ represents the total intensity of the light field, while
$S_{1}$ describes the amount of linear polarisation in the horizontal and
vertical directions. $S_{2}$ describes the amount of diagonal polarisation in
the $45\degree$ and $135\degree$ directions, while $S_{3}$ describes the
amount of circular polarisation in the left and right directions.

###  II.3 Stokes Polarimetry

The primary goal of Stokes polarimetry is to measure the polarisation
structure of an optical field (see Sec. II.2), and thus describe the intensity
and polarisation properties of the light field.

For the sake of clarity, we will explain how the Stokes parameters can be
obtained from intensity measurements performed with a polariser and a phase
retarder. To this end, we write the components of the electric field using
complex notation, namely,

$\begin{split}&E_{x}(t)=E_{0x}\exp{(i\delta_{x})}\exp{(i\omega t)}\\\
&E_{y}(t)=E_{0y}\exp{(i\delta_{y})}\exp{(i\omega t)}\end{split}$ (19)

As a first step, the beam is sent through a Quarter-Wave Plate (QWP), which
produces a phase delay $\theta_{1}$ between both components. As result, the
components of the new electric field take the form

$\begin{split}&E^{\prime}_{x}(t)=E_{x}(t)\exp{(i\theta_{1}/2)}\\\
&E^{\prime}_{y}(t)=E_{y}(t)\exp{(-i\theta_{1}/2)}\end{split}$ (20)

The beam that exits the QWP is then sent through a Linear Polariser (LP),
which transmits light only along the transmission axis of the same.
Mathematically, the total light transmitted through the LP, considering its
transmission axis at an angle $\theta_{2}$ from the horizontal, is described
by,

$\begin{split}&E_{f}(\theta_{1},\theta_{2})=\\\
&E_{x}(t)\exp{(i\theta_{1}/2)}\cos\theta_{2}+E_{y}(t)\exp{(-i\theta_{1}/2)}\sin\theta_{2}.\end{split}$
(21)

Since we can only measure intensities experimentally, we have to compute the
intensity through the relation $E_{f}(\theta_{1},\theta_{2})\cdot
E^{*}_{f}(\theta_{1},\theta_{2})$,

$I(\theta_{1},\theta_{2})=E_{x}E^{*}_{x}\cos^{2}\theta_{2}+E_{y}E^{*}_{y}\sin^{2}\theta_{2}+E_{y}E^{*}_{x}\exp{(-i\theta_{1})}\sin\theta_{2}\cos\theta_{2}+E^{*}_{y}E_{x}\exp{(i\theta_{1})}\sin\theta_{2}\cos\theta_{2},$
(22)

which can be transformed, using trigonometric identities into,

$\begin{split}I(\theta_{1},\theta_{2})&=\frac{1}{2}[E_{x}E^{*}_{x}+E_{y}E^{*}_{y}+(E_{x}E^{*}_{x}-E_{y}E^{*}_{y})\cos
2\theta_{2}\\\ &+(E_{y}E^{*}_{x}+E^{*}_{y}E_{x})\cos\theta_{1}\sin
2\theta_{2}+\\\ &i(E_{y}E^{*}_{x}-E^{*}_{y}E_{x})\sin\theta_{1}\sin
2\theta_{2}],\end{split}$ (23)

We can now use Eq. 18 to rewrite Eq. 23 in terms of the Stokes parameters as,

$\begin{split}&I(\theta_{1},\theta_{2})=\\\ &\frac{1}{2}[S_{0}+S_{1}\cos
2\theta_{2}+S_{2}\cos\theta_{1}\sin 2\theta_{2}+S_{3}\sin\theta_{1}\sin
2\theta_{2}].\end{split}$ (24)

Specific settings for the angles of the QWP and LP, $\theta_{1}$ and
$\theta_{2}$, respectively, give the relations,

$\begin{split}&I(0,0)=\frac{1}{2}(S_{0}+S_{1})\qquad
I(0,90\degree)=\frac{1}{2}(S_{0}+S_{1})\\\
&I(0,45\degree)=\frac{1}{2}(S_{0}+S_{1})\qquad
I(90\degree,45\degree)=\frac{1}{2}(S_{0}+S_{1}),\end{split}$ (25)

from which we can express the Stokes parameters in terms of intensity values
as,

$\begin{split}&S_{0}=I(0,0)+I(0,90\degree),\\\
&S_{1}=I(0,0)-I(0,90\degree),\\\
&S_{2}=2I(0,45\degree)-[I(0,0)+I(0,90\degree)]\\\ &=2I(0,45\degree)-S_{0},\\\
&S_{3}=2I(90\degree,45\degree)-[I(0,0)+I(0,90\degree)\\\
&=2I(90\degree,45\degree)-S_{0}.\end{split}$ (26)

The above equations can be interpreted as follows, the parameter $S_{0}$, the
total intensity of the beam, can be measured as the addition of the intensity
when both, the QWP and the LP are set to $0\degree$, which is the intensity
corresponding to the horizontal polarisation $I_{H}$, plus the intensity when
the QWP is set to $0\degree,$ and the LP to $90\degree$, which corresponds to
the vertical polarisation component $I_{v}$. The Stokes parameter $S_{1}$ is
simple the subtraction of the horizontal and vertical polarisation components.
The intensity $I(0,45\degree)$ corresponds to diagonal polarisation $I_{D}$,
whereas the intensity $I(90\degree,45\degree)$ to right circular polarisation
$I_{R}$. Hence, the Stokes parameters can be written simply as,

$\begin{split}&S_{0}=I_{H}+I_{V},\\\ &S_{1}=I_{H}-I_{V},\\\
&S_{2}=2I_{D}-S_{0},\\\ &S_{3}=2I_{R}-S_{0},\\\ \end{split}$ (27)

Using that the Stokes parameter $S_{0}$ can also be written as the addition of
diagonal $I_{D}$ plus antidiagonal $I_{A}$ polarisation, i.e.
$S_{0}=I_{D}+I_{A}$ and similarly as the addition of right- plus left-handed
circular polarisation, i.e. $S_{0}=I_{R}+I_{L}$, the above equations can also
be written as,

$\begin{split}&S_{0}=I_{H}+I_{V}\\\ &S_{1}=I_{H}-I_{V},\\\
&S_{2}=I_{D}-I_{A},\\\ &S_{3}=I_{R}-I_{L}.\\\ \end{split}$ (28)

Table 1 summarises the different configurations of the QWP and LP to measure
the required intensities to compute the four Stokes parameters.

Table 1: Configuration of the QWP and LP to measure the required intensities to compute the four Stokes parameters  Intensity | QWP | LP
---|---|---
$I_{H}$ | $0\degree$ | $0\degree$
$I_{V}$ | $0\degree$ | $90\degree$
$I_{D}$ | $0\degree$ | $45\degree$
$I_{R}$ | $90\degree$ | $45\degree$

The spatial distribution of the polarisation properties can be obtained by
measuring these parameters at multiple points across the beam using a camera
or other imaging device, which we call a Stokes plot. This distribution can be
plotted as a polarisation ellipse, which provides an intuitive way to
visualise the polarisation properties of the light field in two dimensions.
The major and minor axes of the ellipse represent the polarisation directions,
while the eccentricity of the ellipse represents the ellipticity or degree of
circular polarisation.

###  II.4 Jones Calculus

While the theory of Jones calculus and Vector Vortex Modes (Sections II.4 and
II.5) is not strictly necessary to understand the work presented in this lab
note, we present it nonetheless as we believe it is very useful for the
interested reader to form a well-rounded understanding of polarisation.

It is convenient to express the polarisation states as a column vector and
optical elements as a 4-entry matrix. This “Jones calculus” makes it easy to
obtain the polarisation state of a polarised beam after traversing a single or
a set of optical elements. In terms of a normalised Jones vector (where the
intensity is one) we rewrite Eq. 1 as

$\mathbf{E}_{\mathrm{linear}}=\begin{bmatrix}\cos\theta\\\
\sin\theta\end{bmatrix},$ (29a)

$\mathbf{E}_{\mathrm{circular}}=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ \pm
i\end{bmatrix}.$ (29b)

Note that we use $\theta$ to indicate the rotation of the linear polarisation
from the $x$-axis (horizontal), and $i$ is the complex number. For left-handed
circular polarisation the helicity is positive. For quick reference, and in
Dirac (“ket”) notation for brevity, these six standard polarisation’s are
given by

$\ket{H}=\begin{bmatrix}1\\\ 0\end{bmatrix},$ (30a)

$\ket{V}=\begin{bmatrix}0\\\ 1\end{bmatrix},$ (30b)

$\ket{D}=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ 1\end{bmatrix},$ (30c)

$\ket{A}=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ -1\end{bmatrix},$ (30d)

$\ket{R}=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ -i\end{bmatrix},$ (30e)

$\ket{L}=\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\ i\end{bmatrix}.$ (30f)

Of relevance to this lab note are quarter-wave plates (QWP), half-wave plates
(HWP), linear polarisers and Q-plates Marrucci  _et al._ (2006), with their
respective Jones matrices given by

$\mathrm{QWP}:\frac{1}{\sqrt{2}}\begin{bmatrix}1+i\cos 2\theta&i\sin
2\theta\\\ i\sin 2\theta&1-i\cos 2\theta\end{bmatrix},$ (31a)

$\mathrm{HWP}:\begin{bmatrix}\cos 2\theta&\sin 2\theta\\\ \sin 2\theta&-\cos
2\theta\end{bmatrix},$ (31b)

$\mathrm{Linear~{}Polariser}:\begin{bmatrix}\cos^{2}\theta&\sin\theta\cos\theta\\\
\sin\theta\cos\theta&\sin^{2}\theta\end{bmatrix}.$ (31c)

$\mathrm{Q-plate}:\begin{bmatrix}\cos 2q\psi&\sin 2q\psi\\\ \sin 2q\psi&-\cos
2q\psi\end{bmatrix},$ (31d)

In the special case of the Q-plate, $\psi$ is a spatially resolved azimuthal
angle which describes the azimuthally varying polarisation (i.e. vector beam)
that is produced by the Q-plate.

In order to calculate the final polarisation of a laser beam that has
propagated through several elements, $\mathbf{E}^{{}^{\prime}}$, we simply
write the matrices for the various elements in reverse propagation order (i.e.
the last element is written first), with the incoming beam Jones vector last,
and then perform the multiplication. For example, assuming a diagonally
polarised source is passed through a horizontal polariser and then a quarter-
wave plate at an angle of $45\degree$, the calculation is as follows

$\begin{bmatrix}E_{x}^{{}^{\prime}}\\\
E_{y}^{{}^{\prime}}\end{bmatrix}=\underbracket{\frac{1}{\sqrt{2}}\begin{bmatrix}1+i\cos
2\theta&i\sin 2\theta\\\ i\sin 2\theta&1-i\cos
2\theta\end{bmatrix}}_{\mathrm{QWP,}~{}\theta=45\degree}\underbracket{\begin{bmatrix}\cos^{2}\theta&\sin\theta\cos\theta\\\
\sin\theta\cos\theta&\sin^{2}\theta\end{bmatrix}}_{\mathrm{Polariser,}~{}\theta=0\degree}\underbracket{\frac{1}{\sqrt{2}}\begin{bmatrix}1\\\
1\end{bmatrix}}_{\mathrm{Beam,}~{}\theta=45\degree}=\frac{1}{2}\begin{bmatrix}1\\\
i\end{bmatrix},$ (32)

which is, as we would expect, a left-polarised beam with reduced intensity.

###  II.5 Vector Vortex Modes

While the focus of this lab note is not on vector modes or $q$-plates, we will
briefly discuss few examples of vector modes, as they are the most general
state on polarisation, non-homogeneous, and have attracted a great amount of
attention in fields such as optical communications, optical tweezers optical
metrology, amongst others Rosales-Guzmán  _et al._ (2018); Yang  _et al._
(2021); Ndagano  _et al._ (2018). Given the interest in vector beams, we will
use cylindrical vector vortex (CVV) modes to illustrate the operation of the
device. CVV modes have a donut-like intensity profile, which are generated
from the weighted superposition of light beams with azimuthally varying
spatial phase in opposite direction, which carry orbital angular momentum
(OAM) Cox  _et al._ (2016). Crucially, they also have an azimuthally varying
polarisation, the direction of which depends on the type of mode as described
below.

The generation of both scalar and CVV modes can be accomplished using a
$q$-plate in combination with quarter- and half- wave plates according to the
unitary operation $\hat{Q}$ of a $q$-plateMarrucci  _et al._ (2011):

$\hat{Q}\begin{bmatrix}\ket{R}\\\
\ket{L}\end{bmatrix}=\begin{bmatrix}\ket{+2q}\ket{L}\\\
\ket{-2q}\ket{R}\end{bmatrix},$ (33a) $\hat{Q}\begin{bmatrix}\ket{H}\\\
\ket{V}\end{bmatrix}=\begin{bmatrix}\ket{+2q}\ket{R}+\ket{-2q}\ket{L}\\\
\ket{+2q}\ket{R}-\ket{-2q}\ket{L}\end{bmatrix},$ (33b)

where $q$ is the charge of the topological singularity of the $q$-plate (in
this work $q=1/2$), related to the OAM topological charge, $\ell=2q$. When the
input beam is circularly polarized, the output is a scalar beam of the
opposite polarisation embedded with an OAM of $\pm\ell$. Alternatively, when
the input is linearly polarised the output is a CVV in which polarisation and
OAM are coupled. As such, there are four main types of CVV mode with specific
names as a legacy from radio frequency wave-guide theory, as follows:

$\ket{TM}=\frac{1}{\sqrt{2}}(\ket{\ell}\ket{R}+\ket{-\ell}\ket{L}),$ (34a)

$\ket{TE}=\frac{1}{\sqrt{2}}(\ket{\ell}\ket{R}-\ket{-\ell}\ket{L}),$ (34b)

$\ket{HE^{e}}=\frac{1}{\sqrt{2}}(\ket{\ell}\ket{L}+\ket{-\ell}\ket{R}),$ (34c)

$\ket{HE^{o}}=\frac{1}{\sqrt{2}}(\ket{\ell}\ket{L}-\ket{-\ell}\ket{R}).$ (34d)

###  II.6 Polarisation Cameras

The polarisation camera used in this work makes use of the Sony IMX250 sensor
Yamazaki  _et al._ (2016); Sony Semiconductor Solutions Group . The sensor is
relatively large, with an 11.1 mm diagonal (2/3 inch), a resolution of
2464$\times$2056 3.45 $\mu$m pixels and a maximum frame rate of about 163
frames per second at full resolution (although this would be dependant on the
specific camera, drivers, region of interest, etc.). This sensor is readily
available in several popular brands of camera.

Rather than red-green-blue colour filters in front of the pixels in the
camera, a special polarising wire grid is used which is situated between an
array of micro-lenses and the photodiode layer. The polarisation filter
repeats in a $2\times 2$ grid with linear polarisation’s [H D; A V], shown in
Fig 2 (c). According to the sensor datasheet, the extinction ratios at 532 nm
and 633 nm are approximately 300:1 and 130:1, respectively. For comparison, a
typical Thorlabs polariser has an extinction ratio of about 1000:1, so very
high precision measurements are unfortunately not possible. Perhaps newer
versions of the sensor will improve on this, however, for an 8 bit image the
error due to the extinction ratio is on the same order of magnitude as the
quantisation error.

##  III Experimental Setup and Analysis Code

###  III.1 Experimental Setup

Figure 2: The structure of the polarisation camera pixels with an inset
photograph of the wire-grid (credit Yamazaki  _et al._ (2016)).

Since the camera is only able to measure linear polarisation states, we must
use a quarter wave plate to transform the circular components of the incoming
beam into linear polarisations. We then measure the original beam as well as
the transformed beam simultaneously, and use some software processing to
extract the relevant information.

A photo and drawing of the setup is shown in Fig 2 (a) and (b), respectively.
The setup is quite tight, but easy to align. A 50:50 beam splitter (not
polarising) creates the aforementioned split, and several mirrors are used to
direct the two resulting beams onto halves of the camera sensor. It is
important that these two beams travel in parallel. Any tilt on one of the
beams with respect to the other will result in some error. As such, do not be
tempted to direct the beam from the beam splitter directly to the camera. One
should also be aware that a mirror flips the circular polarisation of a beam,
and so it is important to have an even number of mirrors after the quarter
wave plate.

###  III.2 Analysis Code

The code to calculate the Stokes images, representing $S_{0}$ to $S_{3}$, is
straightforward but not trivial and is provided on GitHub Cox and Rosales-
Guzmán (2023). Users who make enhancements and fixes to the code are
encouraged to submit pull requests for the good of the community.

In essence, on the first run, the code performs a cross-correlation between
the two halves of the camera image. This is cached for subsequent runs as it
takes a few minutes, depending on the resolution of the camera image or region
of interest. After this is done, the measurements are robust to misalignments
of the incoming beam, as long as the components of the Stokes camera (the beam
splitter, wave plate, and mirrors) are not moved. The offset between the two
halves is found using this cross-correlation and is used to digitally align
them for the next step. For each frame, the different polarisation components
are separated (see Sec. II.6), offset accordingly, and the arithmetic (see Eq.
28) is performed to calculate the four Stokes images.

The code is designed to run continuously, as a camera preview window. We
recommend running this preview in a dedicated MATLAB instance so that the
window can remain open while other commands are executed. In addition, if the
user presses ‘p’ while the window is in focus, all of the latest frames are
stored to a PNG image and a MATLAB file for later processing. The GitHub
repository contains additional information and tips, which are expected to
accumulate over time in the spirit of open source.

In addition, we provide a code that can be used separately to plot so-called
Stokes plots. These plots show $S_{0}$ (the intensity of the beam) with a
number of polarisation ellipses overlaid, representing the state of
polarisation at those points of the beam.

##  IV Results and Discussion

Figure 3: Stokes plots of example beams (from left to right: approximately
vertically polarised, circularly polarised and a vector vortex mode) measured
with the stokes camera and plotted showing the polarisation ellipses over the
beam which are generated using the four stokes images. The polarisation of
these beams was controlled using two fibre paddles and in (c) a Q-plate.

We tested the Stokes Camera using a 520 nm, green single-mode fibre-coupled
laser diode. Two fibre paddles were used to control the state of polarisation
of the output Gaussian-like beam. Using the Stokes camera preview window, we
were able to see the polarisation of the beam in real-time, much like
commercially available polarimeters. Two example Gaussian beams are shown in
Fig. 3 (a) and (b).

For our second test, we inserted a Q-plate in the path of the beam Rosales-
Guzmán  _et al._ (2018); Yang  _et al._ (2021); Ndagano  _et al._ (2018); Cox
_et al._ (2016). Figure 3 (c) shows one example of the resulting vector vortex
mode.

The Stokes camera works well, is extremely convenient compared to manual
methods, and is capable of real-time observation and measurement. With a
limited region of interest and a GigE version of the polarisation camera, we
achieved 11 frames per second. This could be improved by using a USB3 camera,
which has a faster inherent frame rate (5 Gbps versus 1 Gbps). The
computational requirements of the code have been optimised significantly (for
example, caching the cross-correlation), but are still formidable. We believe
there is further room for optimisation, but a Python or C++ implementation
would likely enable faster frame rates.

If fast measurements are required, then a video of the camera output (at
hundreds of frames per second, with the right configuration and region of
interest) could be captured without processing, and simply processed offline
using our code.

The dynamic range of our Stokes Camera is the limiting factor compared with
commercial polarimeters. Since an 8 bit image is used (i.e. [0; 255]), the
dynamic range is approximately 48 dB. Using a higher bit depth (for example,
10 bit) on the camera is unnecessary, as the extinction ratio is only 1:300,
which fundamentally only provides about a 1 dB improvement at the expense of
frame rate. In practice, we believe the measurements using this system are
likely to be superior to manually operated Stokes measurement systems often
used in experiments.

##  V Summary

This paper presents a self-contained tutorial aimed at simplifying the task of
real-time reconstruction of the transverse polarisation distribution of an
optical beam using a polarisation camera. One aim of this tutorial is to
popularise Stokes polarimetry amongst a broader audience, encompassing both
researchers and laboratory classes. To this end, we start by clarifying basic
concepts from a theoretical standpoint, including the polarisation ellipse,
the Poincaré sphere, and the Stokes parameters. Furthermore, we outline how a
specific polarisation state might change after passing through single or
multiple optical polarisation components, such as linear polarisers, wave
retarders, and q-plates, using straightforward matrix multiplication. We also
provide a brief overview of vector beams—general light states characterised by
a non-homogeneous polarisation distribution. Our focus on vector beams is
twofold: they have recently gained considerable attention in the field, and
they exemplify the effectiveness of polarisation cameras in reconstructing the
transverse polarisation state of any optical beam, regardless of its
polarisation complexity. Crucially, real-time reconstruction of the transverse
polarisation distribution of vector beams has potential applications in
optical metrology, especially for characterising dynamic systems sensitive to
polarisation. This highlights the importance of timely and accurate methods
for such beam polarisation distribution reconstructions. The early sections of
the tutorial are supported with experimental results, where we detail the use
of the polarisation camera and describe the necessary experimental setup to
simultaneously capture all intensity images needed to compute the Stokes
parameters. For hands-on application, we have included a set of MATLAB scripts
to streamline the polarisation reconstruction process.

## Funding

NRF TTK2204011621; NRF NLCRPP221109667491; Wits Friedel Sellschop Award

## Disclosures

The author declares no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available
at this time but may be obtained from the authors upon reasonable request.

## References

  * Pinnell  _et al._ (2020) J. Pinnell, I. Nape, B. Sephton, M. Cox, V. Rodriguez-Fajardo, and A. Forbes, Journal of the Optical Society of America A 37, 146 (2020).
  * Tidwell  _et al._ (1990) S. C. Tidwell, D. H. Ford, and W. D. Kimura, Applied Optics 29, 2234 (1990).
  * Niziev  _et al._ (2006) V. G. Niziev, R. S. Chang, and A. V. Nesterov, Appl. Opt. 45, 8393 (2006).
  * Passilly  _et al._ (2005) N. Passilly, F. Treussart, R. Hierle, R. de Saint Denis, K. Aït-Ameur, and J.-F. Roch, Journal of the Optical Society of America A 22, 984 (2005).
  * Mendoza-Hernández  _et al._ (2019) J. Mendoza-Hernández, M. F. Ferrer-Garcia, J. A. Rojas-Santana, and D. Lopez-Mago, Opt. Express 27, 31810 (2019).
  * Marrucci  _et al._ (2006) L. Marrucci, C. Manzo, and D. Paparo, Phys. Rev. Lett. 96, 163905 (2006).
  * Naidoo  _et al._ (2016) D. Naidoo, F. S. Roux, A. Dudley, I. Litvin, B. Piccirillo, L. Marrucci, and A. Forbes, Nature Photonics 10, 327 (2016).
  * Devlin  _et al._ (2017) R. C. Devlin, A. Ambrosio, N. A. Rubin, J. B. Mueller, and F. Capasso, Science 358, 896 (2017).
  * Cox  _et al._ (2022) M. A. Cox, T. Celik, Y. Genga, and A. V. Drozdov, Applied Optics 61, D1 (2022).
  * Rosales-Guzmán and Forbes (2017) C. Rosales-Guzmán and A. Forbes,  _How to shape light with spatial light modulators_ , SPIE.SPOTLIGHT (SPIE Press, 2017).
  * Davis  _et al._ (2000) J. A. Davis, D. E. McNamara, D. M. Cottrell, and T. Sonehara, Appl. Opt. 39, 1549 (2000).
  * Maurer  _et al._ (2007) C. Maurer, A. Jesacher, S. Fürhapter, S. Bernet, and M. Ritsch-Marte, New J. Phys. 9, 78 (2007).
  * Moreno  _et al._ (2012) I. Moreno, J. A. Davis, T. M. Hernandez, D. M. Cottrell, and D. Sand, Opt. Express 20, 364 (2012).
  * Mitchell  _et al._ (2017) K. J. Mitchell, N. Radwell, S. Franke-Arnold, M. J. Padgett, and D. B. Phillips, Opt. Express 25, 25079 (2017).
  * Rosales-Guzmán  _et al._ (2017) C. Rosales-Guzmán, N. Bhebhe, and A. Forbes, Opt. Express 25, 25697 (2017).
  * Rong  _et al._ (2014) Z.-Y. Rong, Y.-J. Han, S.-Z. Wang, and C.-S. Guo, Opt. Express 22, 1636 (2014).
  * Liu  _et al._ (2018a) S. Liu, S. Qi, Y. Zhang, P. Li, D. Wu, L. Han, and J. Zhao, Photon. Res. 6, 228 (2018a).
  * Zhan (2009) Q. Zhan, Adv. Opt. Photonics 1, 1 (2009).
  * Rosales-Guzmán  _et al._ (2018) C. Rosales-Guzmán, B. Ndagano, and A. Forbes, J. Opt. 20, 123001 (2018).
  * Yao-Li  _et al._ (2020) Yao-Li, X.-B. Hu, B. Perez-Garcia, Bo-Zhao, W. Gao, Z.-H. Zhu, and C. Rosales-Guzmán, Appl. Phys. Lett. 116, 221105 (2020).
  * Bo  _et al._ (2022) Z. Bo, R.-F. Valeria, H. Xiao-Bo, I. H.-A. Raul, P.-G. Benjamin, and R.-G. Carmelo, Nanophotonics 11, 681 (2022).
  * Rosales-Guzmán  _et al._ (2022) C. Rosales-Guzmán, X. Hu, V. Rodríguez-Fajardo, R. I. Hernandez-Aranda, A. Forbes, and B. Perez-Garcia, J. Opt. 23, 034004 (2022).
  * Hu  _et al._ (2021a) X. Hu, B. Perez-Garcia, V. Rodríguez-Fajardo, R. Hernandez-Aranda, A. Forbes, and C. Rosales-Guzmán, Photon. Res. 9 (2021a).
  * Yang  _et al._ (2021) Y. Yang, Y. Ren, M. Chen, Y. Arita, and C. Rosales-Guzmán, Adv. Photonics 3 (2021).
  * Michihata  _et al._ (2009) M. Michihata, T. Hayashi, and Y. Takaya, Applied Optics 48, 6143 (2009).
  * Kozawa and Sato (2010) Y. Kozawa and S. Sato, Optics Express 18, 10828 (2010).
  * Donato  _et al._ (2012) M. G. Donato, S. Vasi, R. Sayed, P. H. Jones, F. Bonaccorso, A. C. Ferrari, P. G. Gucciardi, and O. M. Maragó, Optics Letters 37, 3381 (2012).
  * Zhao  _et al._ (2005) Y. Zhao, Q. Zhan, Y. Zhang, and Y.-P. Li, Optics Letters 30, 848 (2005).
  * Donato  _et al._ (2021) M. G. Donato, F. Patti, R. Saija, M. A. Iatì, P. G. Gucciardi, F. Pedaci, G. Strangi, and O. M. Maragò, J. Quant. Spectrosc. Radiat. Transfer 258, 107381 (2021).
  * Hao  _et al._ (2010) X. Hao, C. Kuang, T. Wang, and X. Liu, Journal of Optics 12, 115707 (2010).
  * Török and Munro (2004) P. Török and P. Munro, Optics Express 12, 3605 (2004).
  * Segawa  _et al._ (2014a) S. Segawa, Y. Kozawa, and S. Sato, Optics Letters 39, 3118 (2014a).
  * Segawa  _et al._ (2014b) S. Segawa, Y. Kozawa, and S. Sato, Optics Letters 39, 4529 (2014b).
  * Ndagano  _et al._ (2017) B. Ndagano, B. Perez-Garcia, F. S. Roux, M. McLaren, C. Rosales-Guzmán, Y. Zhang, O. Mouane, R. I. Hernandez-Aranda, T. Konrad, and A. Forbes, Nature Phys. 13, 397 (2017).
  * Hu  _et al._ (2019) X.-B. Hu, B. Zhao, Z.-H. Zhu, W. Gao, and C. Rosales-Guzmán, Opt. Lett. 44, 3070 (2019).
  * Töppel  _et al._ (2014) F. Töppel, A. Aiello, C. Marquardt, E. Giacobino, and G. Leuchs, New J. Phys. 16, 073019 (2014).
  * Berg-Johansen  _et al._ (2015) S. Berg-Johansen, F. Töppel, B. Stiller, P. Banzer, M. Ornigotti, E. Giacobino, G. Leuchs, A. Aiello, and C. Marquardt, Optica 2, 864 (2015).
  * Trichili  _et al._ (2020) A. Trichili, M. Cox, B. Ooi, and M.-S. Alouini, Journal of the Optical Society of America B 37 (2020), 10.1364/JOSAB.399168.
  * Cox  _et al._ (2021) M. A. Cox _et al._ , IEEE Journal of Selected Topics in Quantum Electronics 27, 1 (2021), arXiv:2005.14586 .
  * Barreiro  _et al._ (2008) J. T. Barreiro, T.-C. Wei, and P. G. Kwiat, Nature Physics 4, 282 (2008).
  * D’Ambrosio  _et al._ (2016) V. D’Ambrosio, G. Carvacho, F. Graffitti, C. Vitelli, B. Piccirillo, L. Marrucci, and F. Sciarrino, Physical Review A 94, 030304 (2016).
  * Otte  _et al._ (2020) E. Otte, I. Nape, C. Rosales-Guzmán, C. Denz, A. Forbes, and B. Ndagano, J. Opt. Soc. Am. B 37, A309 (2020).
  * Ndagano  _et al._ (2018) B. Ndagano, I. Nape, M. A. Cox, C. Rosales-Guzmán, and A. Forbes, J. Light. Technol. 36, 292 (2018).
  * Milione  _et al._ (2015a) G. Milione, M. P. Lavery, H. Huang, Y. Ren, G. Xie, T. A. Nguyen, E. Karimi, L. Marrucci, D. A. Nolan, R. R. Alfano, _et al._ , Optics Letters 40, 1980 (2015a).
  * Milione  _et al._ (2015b) G. Milione, T. A. Nguyen, J. Leach, D. A. Nolan, and R. R. Alfano, Optics Letters 40, 4887 (2015b).
  * Wang and Song (2017) X. Wang and Y. Song, Optics Express 25, 29342 (2017).
  * Liu  _et al._ (2018b) J. Liu, S.-M. Li, L. Zhu, A.-D. Wang, S. Chen, C. Klitis, C. Du, Q. Mo, M. Sorel, S.-Y. Yu, X.-L. Cai, and J. Wang, Light: Science & Applications 7, 17148 (2018b).
  * Willner (2018) A. E. Willner, Light: Science and Applications 7, 18002 (2018).
  * Lerner  _et al._ (2012) V. Lerner, D. Shwa, Y. Drori, and N. Katz, Opt. Lett. 37, 4826 (2012).
  * Ren  _et al._ (2015) Y.-X. Ren, R.-D. Lu, and L. Gong, Annalen der Physik 527, 447 (2015).
  * Mitchell  _et al._ (2016) K. J. Mitchell, S. Turtaev, M. J. Padgett, T. Čižmár, and D. B. Phillips, Opt. Express 24, 29269 (2016).
  * Gong  _et al._ (2014) L. Gong, Y. Ren, W. Liu, M. Wang, M. Zhong, Z. Wang, and Y. Li, Journal of Applied Physics 116, 183105 (2014).
  * Manthalkar  _et al._ (2020) A. Manthalkar, I. Nape, N. T. Bordbar, C. Rosales-Guzmán, S. Bhattacharya, A. Forbes, and A. Dudley, Opt. Lett. 45, 2319 (2020).
  * Zhao  _et al._ (2020) B. Zhao, X.-B. Hu, V. Rodríguez-Fajardo, A. Forbes, W. Gao, Z.-H. Zhu, and C. Rosales-Guzmán, Appl. Phys. Lett. 116, 091101 (2020).
  * Cox and Drozdov (2021) M. A. Cox and A. V. Drozdov, Applied Optics 60, 465 (2021).
  * Rosales-Guzmán  _et al._ (2020) C. Rosales-Guzmán, X.-B. Hu, A. Selyem, P. Moreno-Acosta, S. Franke-Arnold, R. Ramos-Garcia, and A. Forbes, Sci. Rep. 10, 10434 (2020).
  * Hu  _et al._ (2021b) X.-B. Hu, S.-Y. Ma, and C. Rosales-Guzmán, J. Opt. 23, 044002 (2021b).
  * Scholes  _et al._ (2019) S. Scholes, R. Kara, J. Pinnell, V. Rodríguez-Fajardo, and A. Forbes, Opt. Eng. 59, 1 (2019).
  * Zhao  _et al._ (2019) B. Zhao, X.-B. Hu, V. Rodríguez-Fajardo, Z.-H. Zhu, W. Gao, A. Forbes, and C. Rosales-Guzmán, Opt. Express 27, 31087 (2019).
  * Hu and Rosales-Guzmán (2022) X.-B. Hu and C. Rosales-Guzmán, Journal of Optics 24, 034001 (2022).
  * Toninelli  _et al._ (2019) E. Toninelli, B. Ndagano, A. Vallés, B. Sephton, I. Nape, A. Ambrosio, F. Capasso, M. J. Padgett, and A. Forbes, Adv. Opt. Photon. 11, 67 (2019).
  * Cox and Rosales-Guzmán (2023) M. A. Cox and C. Rosales-Guzmán, “Real-time stokes polarimetry using a polarisation camera,”  (2023), https://github.com/WitsOCLab/StokesCamera.
  * Poincaré’ (1892) H. Poincaré’,  _Theorie Mathematique de la Lumiere_ , Vol. 2 (Gauthiers-Villars, Paris, 1892).
  * Goldstein (2011) D. H. Goldstein,  _Polarized light_ (CRC Press, 2011).
  * Cox  _et al._ (2016) M. A. Cox, C. Rosales-Guzmán, M. P. J. Lavery, D. J. Versfeld, and A. Forbes, Optics Express 24, 18105 (2016).
  * Marrucci  _et al._ (2011) L. Marrucci, E. Karimi, S. Slussarenko, B. Piccirillo, E. Santamato, E. Nagali, and F. Sciarrino, Journal of Optics 13, 064001 (2011).
  * Yamazaki  _et al._ (2016) T. Yamazaki, Y. Maruyama, Y. Uesaka, M. Nakamura, Y. Matoba, T. Terada, K. Komori, Y. Ohba, S. Arakawa, Y. Hirasawa, Y. Kondo, J. Murayama, K. Akiyama, Y. Oike, S. Sato, and T. Ezaki, in  _2016 IEEE International Electron Devices Meeting (IEDM)_ (IEEE, 2016) pp. 8.7.1–8.7.4.
  * (68) Sony Semiconductor Solutions Group, “Polarization Image Sensor,” https://www.sony-semicon.com/en/products/is/industry/polarization.html.


Face recognition (FR) has seen significant advancements due to the utilization of large-scale datasets. Training deep FR models on large-scale datasets with multiple GPUs is now a common practice. In fact, computing power has evolved into a foundational and indispensable resource in the area of deep learning. It is nearly impossible to train a deep FR model without holding adequate hardware resources. Recognizing this challenge, some FR approaches have started exploring ways to reduce the time complexity of the fully-connected layer in FR models. Unlike other approaches, this paper introduces a simple yet highly effective approach, Moving Haar Learning Rate (MHLR) scheduler, for scheduling the learning rate promptly and accurately in the training process. MHLR supports large-scale FR training with only one GPU, which is able to accelerate the model to $1/4$ of its original training time without sacrificing more than $1\%$ accuracy. More specifically, MHLR only needs $30$ hours to train the model ResNet100 on the dataset WebFace12M containing more than $12$M face images with $0.6$M identities. Extensive experiments validate the efficiency and effectiveness of MHLR.

§ INTRODUCTION

Face recognition (FR) plays an important role in real-life applications. In recent years, the scale of FR training datasets have grown from $0.5$M (Million) images with $0.01$M IDs (Identities) in CASIA Webface [28] to $42$M images with $2$M IDs in WebFace42M [29], which indicates the training of deep FR models have entered the large-scale datasets era, called large-scale FR training.

Deep FR models have been greatly advanced because of large-scale training datasets, such as MegaFace [22], MS1M [9], and WebFace42M [29]. Yet, the basic requirement for training a deep FR model gradually grows from $1$ node (Personal computer or server) with $4\times$GPUs to $1$ node with $8\times$GPUs. The state-of-the-art FR approach, Partial FC (PFC) [1], even employs at most $8$ nodes with $8\times$GPUs for training a FR model, which is $64\times$GPUs in total. That shuts many researchers out the door of large-scale FR training. In fact, computing power has become a fundamental and critical resource in the deep learning area. It is nearly impossible to train a deep FR model without holding enough hardware resources, which can provide massive computing power.

Indeed, many cutting-edge FR approaches, e.g. PFC [1], Faster Face Classification (F$^{2}$C) [26], Dynamic Class Queue (DCQ) [14], and Virtual Fully-connected (Virtual FC) [15], have already started to consider the problem of large-scale FR training. Nevertheless, those researches concentrate on the space and time complexity of the fully-connected layer in the model, which means large-scale FR training with multiple GPUs is still a basic requirement, in order to guarantee the feasible training time. Interestingly, the majority of researchers have rarely considered large-scale FR training on $1$ node with only $1\times$GPU. Therefore, we would like to find out if it is possible to do large-scale FR training with $1\times$GPU while assuring its training time and accuracy.

[The loss curve]

[The loss curve with MHLR]

The loss curve of training ResNet100 on MS1MV3. All curves are normalized to $[0,1]$ for putting in one figure. With MHLR, the loss after training $2e5$ steps is close to that after training $8e5$ steps without MHLR. Thus, only $1/4$ of its original training time is required and we can train models with $1\times$GPU.

Till now, the common practice of large-scale FR training is to train models for more than $20$ epochs with at least $4\times$GPUs. In addition, a learning rate (LR) scheduler is employed for stabilizing the convergence of models. However, we identified that the loss curve would be shaped like a staircase while doing multiple step LR scheduling, as shown in Fig. <ref>(a). This phenomenon indicates that manual LR scheduling is vulnerable to training time waste, since many subsequences in the loss curve keep stationary. It is safe to schedule the LR early once the loss curve goes stationary. Thus, the intuition is to eliminate those stationary subsequences of the loss curve efficiently and effectively, as shown in Fig. <ref>(b). In this way, the deep FR model can converge faster without sacrificing much accuracy.

In this paper, we propose a simple yet effective approach, Moving Haar Learning Rate (MHLR) scheduler, which schedules the learning rate promptly and accurately. MHLR is able to guarantee the loss curve does not contain long stationary subsequences, in order to reduce the training time of the model. Thus, the required time for large-scale FR training with $1\times$GPU can be similar to $8\times$GPUs. The contribution of this paper can be summarized as follows:

* Efficient. MHLR is able to train the model for $5$ epochs on $1\times$GPU, which costs $1/4$ of its original training time. To be specific, on WebFace12M, MHLR requires $30$H (hours) compared to the model trained for $20$ epochs on $8\times$GPUs, which requires around $24$H.

* Effective. MHLR sacrifices negligible accuracy for reducing training epochs from $20$ to $5$. Specifically, for ResNet100 trained on WebFace12M, the accuracy of MHLR on the dataset IJB-C is $97.20$ while PFC-0.3 is $97.58$. For ResNet100 trained on WebFace42M, the accuracy of PFC-0.3 on IJB-C is $97.82$. Thus, we conclude the accuracy loss of MHLR is less than $1\%$.

* Impactful. MHLR proves that large-scale FR training with $1\times$GPU is not a mission impossible. It gives the possibility for many researchers to enter the door of large-scale FR training without massive hardware resources.

§ RELATED WORK

In this section, we mainly reviewed the FR approaches relating to PFC [1], since it is one of the state-of-the-art approaches aiming at large-scale FR training.

Back in 2017, SphereFace [17] had already started to train FR models with $4\times$GPUs. It first introduced mapping the face image onto a hyper-sphere in the feature space. CosFace [25] proposed to normalize both the feature vector and the class vector in order to stabilize the training process. It adopted $8\times$GPUs to train FR models. UniformFace [8] introduced to learn equidistributed feature vectors in order to maximize the utilization of the feature space. It employed $4\times$GPUs to train FR models. ArcFace [3] proposed additive angular margin in order to maximize the inter-class separation and the intra-class compactness. It trained FR models with $4\times$GPUs. CurricularFace [12] embedded the idea of curriculum learning into the loss function, in order to address easy samples in the early training stage and hard ones in the later stage. It utilized $4\times$GPUs in the FR model training. MagFace [19] suggested building a connection between the quality of images and the magnitude of their corresponding feature vectors, in order to improve the generalization ability of the model. It leveraged $8\times$GPUs in the FR model training. ElasticFace [2] proposed to allow flexibility in the push for inter-class separability, in order to give the decision boundary chances to extract and retract to allow space for flexible class separability learning. It used $4\times$GPUs to train FR models. Also, there are some papers in which we failed to find the hardware information, such as Large-Margin Softmax (L-Softmax) loss [16], Additive Margin Softmax (AM-Softmax) loss [24], Circle loss [23], and AdaFace [13]. To sum up, all FR approaches employ multiple GPUs for training FR models, except the hardware information of some approaches is not given. Therefore, we can conclude that training FR models with multiple GPUs is a common practice in this area.

Recently, some FR approaches have already started to consider large-scale FR training. Dynamic Class Queue (DCQ) [14] proposed to dynamically select and generate class vectors instead of using them all in the training process. It adopted $8\times$GPUs to train FR models. Similar to DCQ, Faster Face Classification (F$^{2}$C) [26] introduced Dynamic Class Pool (DCP) to store and update class vectors dynamically. It employed $8\times$GPUs for the FR model training. Instead of maintaining a queue or a pool, Partial Fully-Connected (PFC) [1] suggested randomly selecting negative class vectors. It leveraged at least $8\times$GPUs for the FR model training. Virtual Fully-connected (Virtual FC) [15] has a similar goal to our approach, which is doing large-scale FR training with $1\times$GPU. Virtual FC focuses on reducing the space and time complexity of the FC layer, yet we focus on training epochs, which are different and have no conflict. In summary, all the mentioned approaches above concentrate on the FC layer. Yet in this paper, we aim at training epochs, which means our approach is compatible with the other approaches.

§ PROBLEM STATEMENT

In this section, we first revisit the large-scale FR training. After that, the problem in this paper is discussed.

A face image is defined as $\mathbf{x}\in\mathbb{R}^{N\times N\times C}$, where $N$ is the size of the image and $C$ is the channel of it. Thus, a deep FR model is defined as $\mathbf{f}=\mathcal{F}(\mathbf{x})$, where $\mathbf{f}\in\mathbb{R}^{D}$ is the feature vector extracted from $\mathbf{x}$. In order to train the model, the equation of combined margin loss function used in ArcFace and Partial FC is reviewed as follows:

\begin{equation}
\label{eq:cml}
\begin{aligned}
&\mathcal{L}=-\frac{1}{B}\sum_{i=1}^{B}\log\frac{e^{s\mathcal{M}(\theta_{y_{i},i})}}{e^{s\mathcal{M}(\theta_{y_{i},i})}+\sum_{j=1,j\notin y_{i}}^{C}e^{s\cos\theta_{j,i}}}\\
\end{aligned}
\end{equation}

where $B$ denotes the batch size, $C$ is the number of classes, $\mathcal{M}$ means the marginal function, $m_{1}$ stands for the multiplicative angular margin, $m_{2}$ represents the additive angular margin, $m_{3}$ denotes the additive cosine margin, $s$ is the scale factor controlling the radius of the hyper-sphere in feature space, $\theta_{j,i}$ represents the angle between the class vector $\mathbf{w}_{j}$ and the feature vector $\mathbf{f}_{i}$, $\mathbf{w}_{j}\in\mathbb{R}^{D}$ stands for the $j$-th class vector, and $\mathbf{f}_{i}$ denotes the feature vector of the $i$-th image belonging to the $y_{i}$-th class. Thus, $\theta_{y_{i},i}$ can be calculated by the equation $\theta_{y_{i},i}=\arccos(\mathbf{w}_{y_{i}}\mathbf{f}_{i}/\Vert\mathbf{w}_{y_{i}}\Vert\Vert\mathbf{f}_{i}\Vert)$, and $\mathcal{M}(\theta_{y_{i},i})$ can be computed afterwards. Yet, since the time complexity of $\arccos(\cdot)$ is pretty high, the variant in Eq. (<ref>) is a better option. Note that $m_{1}$ will change the frequency of the $\cos$ function and destroy its monotonicity. Thus, we set $m_{1}$ as $1$ in this paper and thus it is treated as a constant.

\begin{equation}
\label{eq:cos}
\begin{aligned}
&=\sin m_{2}\cos\theta_{y_{i},i}-\cos m_{2}\sqrt{1-\cos^2\theta_{y_{i},i}}-m_{3}
\end{aligned}
\end{equation}

where $m_{1}=1$, $\sin m_{2}$, $\cos m_{2}$ and $m_{3}$ are constants which can be pre-computed and stored, and $\cos\theta_{y_{i},i}=\mathbf{w}_{y_{i}}\mathbf{f}_{i}/\Vert\mathbf{w}_{y_{i}}\Vert\Vert\mathbf{f}_{i}\Vert$.

To this end, we are able to train the model $\mathcal{F}$. Specifically, $\mathcal{F}$ is trained with $E$ epochs, where each epoch contains $U$ training iterations. Thus, $\mathcal{F}$ is updated $T=E\cdot U$ times in total, where the model and the loss after $t$-th updating are denoted $\mathcal{F}_{t}$ and $\mathcal{L}_{t}$ respectively. Formally, all $\mathcal{F}_{t}$ forms a sequence (time-series) $\mathbf{F}=\{\mathcal{F}_{0}, \mathcal{F}_{1}, \dots, \mathcal{F}_{T}\}$. Similarly, all $\mathcal{L}_{t}$ forms another sequence $\mathbf{L}=\{\mathcal{L}_{0}, \mathcal{L}_{1}, \dots, \mathcal{L}_{T}\}$. Hence, we wish to identify a sequence of the learning rate $\Gamma=\{\gamma_{0}, \gamma_{1}, \dots, \gamma_{T}\}$, in order to minimize the final loss $\mathcal{L}_{T}$ with as less $T$ as possible.

In practice, the learning rate (LR) $\lambda$ is scheduled manually in the training process, which means the scheduling strategy is a hyper-parameter. The majority of FR approaches select Multiple Step LR (MultiStepLR) scheduler, where $\lambda_{t}$ is reduced by a decay factor $\delta$ after certain training iterations. The new LR $\lambda_{t+1}$ can be calculated by $\lambda_{t+1}=\lambda_{t}/\delta$. However, we discovered that MultiStepLR will always create many stationary subsequences in the loss curve, as already shown in Fig. <ref>(a). It represents that large-scale FR training suffers from training time waste, which indicates it is safe to schedule the LR early once the loss curve goes stationary. Therefore, we propose an easy yet effective LR scheduler, namely Moving Haar LR (MHLR) scheduler, in order to accelerate the model convergence.

§ THE PROPOSED APPROACH

In this section, Moving Haar Learning Rate (MHLR) scheduler is introduced from three subsections. The first subsection presents the moving average part of MHLR. The second subsection demonstrates the Haar kernel part of MHLR. The third subsection summarizes the whole algorithm of MHLR.

§.§ Exponential Moving Average

Given a loss sequence $\mathbf{L}=\{\mathcal{L}_{0}, \mathcal{L}_{1}, \dots, \mathcal{L}_{T}\}$, it is hard to detect the inside stationary subsequences as it contains too many noises, which are caused by the optimization algorithm, Mini-batch Gradient Descent. More specifically, Batch Gradient Descent can generate a more smoothed loss sequence $\mathbf{L}$, since it considers all samples (face images) in every iteration. Thus, the loss function $\mathcal{L}$ remains unchanged. On the contrary, Mini-batch Gradient Descent considers distinct $B$ samples in every iteration. Accordingly, the loss function $\mathcal{L}$ will always change in each iteration. This leads to the result that the generated loss sequence $\mathbf{L}$ is quite unstable, which contains too many noises, as shown in Fig. <ref>(a). Therefore, the first process in MHLR is smoothing the loss sequence $\mathbf{L}$, which filters noises in order to help the detection of stationary subsequences. The expected $\mathbf{L}$ after smoothing is shown in Fig. <ref>(b).

[The loss curve $\mathbf{L}$ before smoothing]

[The loss curve $\mathbf{L}^{EMA}$ after smoothing]

Smoothing the loss curve by EMA. As shown in (a), the original loss curve is hard to be analyzed since it contains too much noise. Thus, as shown in (b), we adopt EMA to smooth it. Note all curves are normalized

In this paper, we select Exponential Moving Average (EMA) to smooth $\mathbf{L}$. There are two reasons why we adopt EMA: 1) The loss sequence $\mathbf{L}$ is generated on-the-fly in the training process. EMA supports real-time updating, which is just compatible with our problem; 2) The time complexity of EMA in each iteration is only $\mathcal{O}(1)$, which indicates that it occupies negligible computations. The equation of EMA is given in Eq. (<ref>).

\begin{equation}
\label{eq:ema}
\mathcal{L}_{t}^{EMA}=
\begin{cases}
\mathcal{L}_{t},\hfill t=0\\
\alpha\mathcal{L}_{t}+(1-\alpha)\mathcal{L}_{t-1}^{EMA},\hfill t>0
\end{cases}
\end{equation}

where $\alpha\in[0,1]$ indicates the importance of the current loss $\mathcal{L}_{t}$ compared to the average of previous losses $\mathcal{L}_{t-1}^{EMA}$, and $t\in\{0,1,\dots,T\}$. Initially, when $t=0$, $\mathcal{L}_{0}^{EMA}=\mathcal{L}_{0}$ directly. After that, $\mathcal{L}_{t}^{EMA}$ is calculated by the current loss $\mathcal{L}_{t}$ and the average of previous losses $\mathcal{L}_{t-1}^{EMA}$. In this paper, we $\alpha$ is set as $0.001$.

Overall, EMA generates a smoothed loss sequence  $\mathbf{L}^{EMA}=\{\mathcal{L}_{0}^{EMA}, \mathcal{L}_{1}^{EMA}, \dots, \mathcal{L}_{T}^{EMA}\}$ on-the-fly while generating the loss sequence $\mathbf{L}$, in order to support the detection of stationary subsequences.

§.§ Haar Convolutional Kernel

With EMA calculating $\mathbf{L}^{EMA}$ on-the-fly, we are able to detect stationary subsequences on it in real-time. Formally, given a smoothed loss sequence $\mathbf{L}^{EMA}=\{\mathcal{L}_{0}^{EMA}, \mathcal{L}_{1}^{EMA}, \dots, \mathcal{L}_{T}^{EMA}\}$, we wish to identify if the subsequence of $\mathbf{L}^{EMA}$ starting from $S$ and ending at $T$, denoted $\mathbf{L}_{S,T}^{EMA}=\{\mathcal{L}_{S}^{EMA}, \mathcal{L}_{S+1}^{EMA}, \dots, \mathcal{L}_{T}^{EMA}\}$, keeps stationary statistically, where $0\leq S<T$. Note $\mathcal{L}_{T}^{EMA}$ is always the last value in $\mathbf{L}^{EMA}$ and it increases continuously, since $\mathbf{L}^{EMA}$ is updating in real-time.

In this paper, we adopt a simple Haar Convolutional Kernel (HCK) in the Haar-like Feature Descriptors to detect stationary subsequences. HCK is defined as $\mathcal{H}_{2s}$, where $s$ represents the half size of it, the first half values of $\mathcal{H}_{2s}$ from $1$ to $s$ are $-1$, the second half values of $\mathcal{H}_{2s}$ from $s+1$ to $2s$ are $1$. For example, $\mathcal{H}_{4}=\{-1,-1,1,1\}$. Thus, performing a convolutional operation on $\mathbf{L}^{EMA}$ using $\mathcal{H}_{2s}$, denoted $\mathcal{H}_{2s}*\mathbf{L}^{EMA}$, indicates calculating the difference between the sum of the first half values in $\mathbf{L}_{S,T}^{EMA}$ and the sum of the second half values in $\mathbf{L}_{S,T}^{EMA}$, where the size of $\mathbf{L}_{S,T}^{EMA}$ is $2s$. In this way, a stationary $\mathbf{L}_{S,T}^{EMA}$ will get a small difference value, while a $\mathbf{L}_{S,T}^{EMA}$ with a decreasing trend will get a relatively larger difference value. Hence, we are now able to distinguish the decreasing curve and the stationary curve. Nevertheless, a large $2s$ will cause two drawbacks:

* The time complex of computing a single step convolution between $\mathcal{H}_{2s}$ and $\mathbf{L}_{S,T}^{EMA}$ is $\mathcal{O}(2s)$. It is inefficient to compute a long $\mathcal{H}_{2s}$ in every training iteration. For example, if $s=5000$, a $10$K times multiplication and a $10$K times summation are needed in each iteration. In large-scale FR training, the total iterations $T$ is commonly more than $500$K, which results in an extra $10$B (Billion) times computations in total.

* A long $\mathcal{H}_{2s}$ requires a delay in detection. For instance, if $s=5000$, then we have to wait until $\mathbf{L}^{EMA}$ accumulates $10$K values before starting the detection, since the convolutional operation requires the size of $\mathbf{L}^{EMA}$ to be longer than $\mathcal{H}_{2s}$. However, we wish the stationary subsequence detection could start as soon as possible.

Thus, $s$ is set to $1$ in this paper, since it requires negligible computations and the shortest delay. Moreover, when $s=1$, this convolutional operation degrades into calculating the difference $\mathcal{D}_{t}$ of $\mathcal{L}_{t-1}^{EMA}$ and $\mathcal{L}_{t}^{EMA}$. We denote the sequence of $\mathcal{D}_{t}$ as $\mathbf{D}=\{\mathcal{D}_{1},\mathcal{D}_{2},\dots,\mathcal{D}_{T}\}$ which can be calculated on-the-fly. The equation of it is given below:

\begin{equation}
\label{eq:diff}
\mathcal{D}_{t} = \mathcal{L}_{t-1}^{EMA} - \mathcal{L}_{t}^{EMA}
\end{equation}

where $t\in\{1,2,\dots,T\}$ here. Nonetheless, considering only two values in each step leads to a result that $\mathbf{D}$ is sensitive to noises. As shown in Fig. <ref>(a), it looks like a meaningless sequence, which is hard to use.

[The curve of $\mathbf{D}$]

[The curve of $\mathbf{D}^{EMA}$]

Examples of $\mathbf{D}$ and $\mathbf{D}^{EMA}$. As shown in (a), it is hard to find a relationship between $\mathbf{L}$ and $\mathbf{D}$. By contrast, as shown in (b), $\mathbf{D}$ is smoothed by EMA and denoted as $\mathbf{D}^{EMA}$. It is clear to observe a spike on $\mathbf{D}^{EMA}$ every time when there is a decline on $\mathbf{L}$. Note all curves are normalized

Therefore, we leverage EMA the second time, which aims at smoothing $\mathbf{D}$. The smoothed $\mathbf{D}$ is denoted as $\mathbf{D}^{EMA}=\{\mathcal{D}_{1}^{EMA}, \mathcal{D}_{2}^{EMA}, \dots, \mathcal{D}_{T}^{EMA}\}$. The equation of calculating $\mathbf{D}^{EMA}$ is shown as follows:

\begin{equation}
\label{eq:dema}
\mathcal{D}_{t}^{EMA}=
\begin{cases}
\mathcal{D}_{t},\hfill t=1\\
\beta\mathcal{D}_{t}+(1-\beta)\mathcal{D}_{t-1}^{EMA},\hfill t>1
\end{cases}
\end{equation}

where $\beta\in[0,1]$ indicates the importance of the current difference $\mathcal{D}_{t}$ compared to the average of previous differences $\mathcal{D}_{t-1}^{EMA}$. In this paper, $\beta$ is set as $0.001$.

Consequently, $\mathbf{D}^{EMA}$ is able to help us identify stationary subsequences, as shown in Fig. <ref>(b). In detail, we are able to observe a spike on $\mathbf{D}^{EMA}$ every time when there is a decline on $\mathbf{L}$, which indicates a rapid drop in the loss curve. Thus, the subsequence with small values before the spike in $\mathbf{D}^{EMA}$ represents the stationary part in $\mathbf{L}$. We conclude that MHLR can schedule $\mathcal{L}_{t}$ once identifying a consecutive subsequence of small $\mathcal{D}_{t}^{EMA}$.

§.§ Moving Haar Learning Rate Scheduler

To this end, we propose to combine EMA (Eq. (<ref>)) and HCK (Eq. (<ref>)) together. Thus, the calculation and storage for $\mathbf{L}^{EMA}$ are omitted. The combined approach is called Moving Haar Learning Rate (MHLR) scheduler. MHLR is able to detect the stationary subsequences on-the-fly in the training process. Its time and space complexity is $\mathcal{O}(1)$ in each iteration.

First, with Eq. (<ref>) and Eq. (<ref>), we are able to get a variant of the equation for calculating $\mathcal{D}_{t}$, as given in Eq. (<ref>).

\begin{equation}
\label{eq:vdiff}
\begin{aligned}
\mathcal{D}_{t}&=\mathcal{L}_{t-1}^{EMA}-\mathcal{L}_{t}^{EMA}\\
\end{aligned}
\end{equation}

Next, it is easy to have $\mathcal{D}_{t-1}=\alpha\mathcal{L}_{t-2}^{EMA}-\alpha\mathcal{L}_{t-1}$ from Eq. (<ref>), where $\mathcal{D}_{t-1}$ just means one iteration before $\mathcal{D}_{t}$. After performing $\mathcal{D}_{t-1}-\mathcal{D}_{t}$, we have:

\begin{equation}
\label{eq:d2-d1}
\begin{aligned}
\mathcal{D}_{t-1}-\mathcal{D}_{t}&=\alpha\mathcal{L}_{t-2}^{EMA}-\alpha\mathcal{L}_{t-1}-\alpha\mathcal{L}_{t-1}^{EMA}+\alpha\mathcal{L}_{t}\\
\end{aligned}
\end{equation}

From Eq. (<ref>), we know $\mathcal{D}_{1}=\alpha\mathcal{L}_{0}^{EMA}-\alpha\mathcal{L}_{1}$. Then, since $\mathcal{L}_{0}^{EMA}=\mathcal{L}_{0}$ in Eq. (<ref>), we have $\mathcal{D}_{1}=\alpha(\mathcal{L}_{0}-\mathcal{L}_{1})$. At last, we are able to get the final equation of $\mathcal{D}_{t}$ from Eq. (<ref>), which is shown as follows:

\begin{equation}
\label{eq:mhlr}
\mathcal{D}_{t} =
\begin{cases}
\alpha(\mathcal{L}_{t-1}-\mathcal{L}_{t}),\hfill t=1\\
(1-\alpha)\mathcal{D}_{t-1}+\alpha(\mathcal{L}_{t-1}-\mathcal{L}_{t}),\hfill t>1
\end{cases}
\end{equation}

To this end, calculating $\mathcal{D}_{t}$ by Eq. (<ref>) does not include $\mathcal{L}_{t}^{EMA}$ anymore. In addition, it has three advantages: 1) $\mathcal{D}_{t}$ can be calculated on-the-fly. It requires only $1$ step delay, i.e. when $t=0$; 2) Its time complexity is $\mathcal{O}(1)$ in each iteration; 3) Its space complexity is also $\mathcal{O}(1)$, since it only needs to store the $T$-th and $(T-1)$-th values in sequences $\mathbf{D}$ and $\mathbf{L}$, instead of the whole sequences.

After calculating $\mathcal{D}_{t}$ by Eq. (<ref>) and then $\mathcal{D}_{t}^{EMA}$ by Eq. (<ref>), we set a threshold $\lambda$ in order to determine if $\mathcal{D}_{t}^{EMA}$ is small enough, which indicates the loss sequence $\mathbf{L}$ keeps stationary. In addition, a tolerance $\tau$ is defined in case that MHLR is too sensitive to $\lambda$. Specifically, MHLR will reduce the learning rate $\gamma_{t}$ if $\mathcal{D}_{t}^{EMA}$ keeps lower than $\lambda$, i.e. $\mathcal{D}_{t}^{EMA}<\lambda$, for $\tau$ iterations. After extensive experiments, we empirically set $\lambda=5e-5$ and $\tau=0.05T$. The algorithm of MHLR is shown in Alg. <ref>.

The algorithm of MHLR

Input: Previous learning rate $\gamma_{t-1}$, current loss $\mathcal{L}_{t}$, and previous loss $\mathcal{L}_{t-1}$

Parameter: The threshold $\lambda$, tolerance $\tau$,  $\alpha=0.001$, and $\beta=0.001$

Output: Current learning rate $\gamma_{t}$

Initialization: Iterations count $t=1$, tolerance count $c=0$, $\mathcal{D}_{t-1}=0$, and $\mathcal{D}_{t-1}^{EMA}=0$
[1]
$\mathcal{D}_{t}\gets\alpha(\mathcal{L}_{t-1}-\mathcal{L}_{t})$$\triangleright$ Eq. (<ref>)
$\mathcal{D}_{t}^{EMA}\gets\mathcal{D}_{t}$$\triangleright$ Eq. (<ref>)
$\mathcal{D}_{t}\gets(1-\alpha)\mathcal{D}_{t-1}+\alpha(\mathcal{L}_{t-1}-\mathcal{L}_{t})$$\triangleright$ Eq. (<ref>)
$\mathcal{D}_{t}^{EMA}\gets\beta\mathcal{D}_{t}+(1-\beta)\mathcal{D}_{t-1}^{EMA}$$\triangleright$ Eq. (<ref>)
$\gamma_{t}\gets \gamma_{t-1}$
$c\gets c+1$
$\gamma_{t}\gets \gamma_{t-1}/2$
$c\gets 0$
$t\gets t+1$
return $\gamma_{t}$

§ EXPERIMENTS

§.§ Implementation Details

Datasets. In this paper, we separately employ MS1MV2 [3], MS1MV3 [4], WebFace4M, WebFace8M, and WebFace12M [29] as the training sets. MS1MV2 is also called MS1M-ArcFace, which has more than $5.8$M images with $0.085$M identities. MS1MV3 is also called MS1M-RetinaFace, which contains more than $5$M images with $0.09$M IDs. WebFace42M [29] contains $42$M images with $2$M IDs. Afterwards, we randomly select $10\%$, $20\%$, and $30\%$ samples from WebFace42M as WebFace4M, WebFace8M, and WebFace12M respectively. Hence, WebFace12M can have more than $12$M images with $0.6$M IDs in total. Note that the reason we do not employ WebFace42M as our training set is that we have not applied PFC [1] yet. Large-scale FR training without PFC will lead to a huge fully-connected layer, which is infeasible to have enough memory and time to train it on a single GPU.

Many public test sets are adopted to evaluate the performance of MHLR, including LFW [11], CFP-FP [21], AgeDB-30 [20], IJB-B [27], and IJB-C [18].

Experimental Settings. The experiments in the paper are implemented on a computer equipped with an Intel Core i9-13900K 3.00 GHz, 32 GB memory, and an NVIDIA GeForce RTX 4090 GPU. The operating system is Windows 11. In addition, the development environment is Anaconda 23.7.2 with Python 3.11.4 and Pytorch 2.0.1. Following the settings in recent papers  [1, 3], input images are cropped and centered to the size of $112\times112\times3$. Each pixel in the image is normalized by subtracting $127.5$ and then being divided by $128$, which is mapped from $[0,255]$ to $[-1,1]$. Images are randomly flipped for data augmentation. Customized ResNet50 and ResNet100 [1, 10] are employed as the backbones. The ResNet model outputs the feature vector of dimension $512$. After that, feature vectors are input into the combined margin loss. For training set MS1MV3, $m_{1}=1$, $m_{2}=0.5$ and $m_{3}=0$. For training sets WebFace4M, WebFace8M, and WebFace12M, $m_{1}=1$, $m_{2}=0$ and $m_{3}=0.4$. The scale factor $s$ is set as $64$. In addition, the ResNet model is trained for $5$ epochs on different training sets with only one single GeForce RTX 4090 GPU. The batch size is set as $128$ samples in each updating. Stochastic Gradient Descent (SGD) with momentum $0.9$ and weight decay $5e-4$ is selected as the optimizer. The learning rate is set as $0.2$, the decay factor $\delta=2$, and it is automatically scheduled by MHLR. For EMA, $\alpha$ and $\beta$ are both set as $0.001$. The hyper-parameters in MHLR are set as $\lambda=1e-5$ and $\tau=0.05$.

For testing, two feature vectors extracted from the original image and its horizontal-flipped one are summed together as the final feature vector. In addition, True Accept Rate (TAR) @ False Positive Rate (FAR) $=1e-4$ is reported on IJB-B and IJB-C, denoted TAR@FAR$=1e-4$.

§.§ Ablation Study

Different FR approaches with/without MHLR. In Tab. <ref>, we re-implement CosFace and ArcFace. After that, we train ResNet100 by them on different datasets with/without MHLR. The CosFace and ArcFace with MHLR are denoted as Cos+MHLR and Arc+MHLR. All approaches are trained on $1\times$GPU. Without MHLR, ArcFace and CosFace need to train models for $20$ epochs in common practice, which takes $36$H (Hours) on MS1MV3 and $108$H on WebFace12M. That represents a high time cost for researchers training a FR model. On the contrary, MHLR trains models for only $5$ epochs, which takes $9$H on MS1MV3 and $30$H on WebFace12M. That is nearly $1/4$ of the training time required by ArcFace and CosFace. Moreover, on MS1MV3, the accuracy of ArcFace is $96.83$ while that of Arc+MHLR is $96.41$. On WebFace12M, the accuracy of CosFace is $97.55$ while that of Cos+MHLR is $97.20$. It means that the accuracy drop of MHLR is tiny, which is less than $1\%$. Here, we would like to propose a question if it is worth to chase the small accuracy increasing while shutting the door of large-scale FR training for the majority of researchers. In summary, MHLR gives a competitive result with only $1/4$ of the original training time. Note we have not re-implemented PFC yet, and thus it is infeasible to train models on WebFace42M with $1\times$GPU since the limitation of the GPU memory.

MHLR across different ResNet models In Tab. <ref>, we train different models on MS1MV3 and WebFace12M. All models are trained for $5$ epochs on $1\times$GPU. Note the training on MV1MV3 is ArcFace+MHLR while that on WebFace12M is CosFace+MHLR. We witness a time increasing from ResNet18 all the way to ResNet200. More specifically, on MS1MV3, the training time and throughput are from 4H and 1876 images/s to 16H and 450 images/s respectively. On WebFace12M, the training time and the throughput are from 18H and 1015 images/s to 30H and 590 images/s respectively. Note we are not training ResNet200 on WebFace12M, since it is more than a NVIDIA Geforce RTX 4090 can handle. In addition, we observe a more significant training time increasing for ResNet50 and ResNet100. It indicates the computational bottleneck is not on the model for ResNet18, ResNet34, and ResNet50. By contrast, the computational time required by ResNet50 and ResNet100 overwhelms the computational time required by others. For instance, on MS1MV3, the training time of ResNet18 and ResNet34 do not have much difference if we do not arrange more threads in reading images, since the bottleneck is the speed of data loading. We can also observe that the throughput on WebFace12M is lower than that on MS1MV3. The reason is that the number of identities in WebFace12M is more than that in MS1MV3. Thus, it requires more computational time for the FC layer in the model. For accuracy, it is observed an increasing trend from ResNet 18 to ResNet200.

Method     Dataset     Epoch     Time     IJB-B     IJB-C

ArcFace     MS1MV3         20     36H      95.03     96.73

ArcFace     WebFace4M      20     32H      94.91     96.75

ArcFace     WebFace12M     20     107H     95.22     97.31

CosFace     MS1MV3         20     36H      95.02     96.54

CosFace     WebFace4M      20     32H      94.60     96.68

CosFace     WebFace12M     20     108H     95.41     97.55

Arc+MHLR     MS1MV3         5     9H      94.98     96.21

Arc+MHLR     WebFace4M      5     8H      94.85     96.17

Arc+MHLR     WebFace12M     5     30H     95.37     96.99

Cos+MHLR     MS1MV3         5     9H      94.91     96.11

Cos+MHLR     WebFace4M      5     8H      94.16     96.04

Cos+MHLR     WebFace12M     5     30H     95.46     97.08

The performance of ArcFace and CosFace with/without MHLR. ResNet100 is selected as the model and trained with $1\times$GPU. TAR@FAR$=1e-4$ is reported as the accuracy ($\%$) on IJB-B and IJB-C

Model     Dataset     images/s     Time     IJB-B     IJB-C

ResNet18      MS1MV3         1876     4H      91.51     93.62

ResNet34      MS1MV3         1423     5H      93.63     95.31

ResNet50      MS1MV3         1203     6H      94.69     96.12

ResNet100     MS1MV3         786      9H      95.03     96.41

ResNet200     MS1MV3         450      16H     95.16     96.65

ResNet18      WebFace12M     1015     18H     92.13     94.43

ResNet34      WebFace12M     878      20H     94.12     95.87

ResNet50      WebFace12M     783      23H     95.25     96.85

ResNet100     WebFace12M     590      30H     95.68     97.20

The performance of MHLR training different models. All models are trained for $5$ epochs with $1\times$GPU. TAR@FAR$=1e-4$ is reported as the accuracy ($\%$) on IJB-B and IJB-C

§.§ Hyperparameter Study

Initial learning rate $\gamma_{0}$ and decay factor $\delta$  In Fig. <ref>, we grid search the initial learning rate $\gamma_{0}$ and decay factor $\delta$. All results are reported based on ResNet100 trained for $5$ epochs on MS1MV3 with $1\times$GPU, while the test set is IJB-C. We observe an accuracy drop while $\delta$ is increasing, where a large $\delta$ represents the learning rate $\gamma_{t}$ will decrease fast in the training process. That indicates a large $\delta$ is harmful for large-scale FR training, since the model will stop learning if $\gamma_{t}$ becomes tiny in a very short time, as shown in Fig. <ref>(b). We also witness that $\gamma_{0}=0.02$ achieves the highest accuracy, since a small $\gamma_{t}$ makes the learning stop while a big $\gamma_{t}$ makes the learning unstable. In summary, we set $\gamma_{0}=0.02$ and $\delta=2$ in this paper.

[The heatmap of $\gamma_{0}$ and $\delta$]

[The impact of decay factor $\delta$]

Grid searching for the initial learning rate $\gamma_{0}$ and decay factor $\delta$. In (a), it shows the heatmap of $\gamma_{0}$ and $\delta$, where TAR@FAR$=1e-4$ on IJB-C is reported as the accuracy ($\%$). In (b), it illustrates the impact of the decay factor $\delta$, where $\delta=8$ for $\mathbf{\Gamma}_{1}$ and $\delta=2$ for $\mathbf{\Gamma}_{2}$. Thus, $\mathbf{\Gamma}_{1}$ declines fast and $\mathbf{L}^{EMA}_{1}$ goes stationary early

[The heatmap of $\lambda$ and $\tau$]

[Relationship of $\lambda$ and $\mathcal{D}^{EMA}_{t}$]

Grid searching for the threshold $\lambda$ and tolerance $\tau$ in MHLR. In (a), it shows the heatmap of TAR@FAR$=1e-4$ on IJB-C is reported as the accuracy ($\%$). In (b), it illustrates the relationship of different $\lambda$ and $\mathcal{D}^{EMA}_{t}$. The curve of $\mathcal{D}^{EMA}_{t}$ under $\lambda$ will be treated as the signal to schedule the learning rate

2*Method     2*Dataset     2*Model     2*#GPUs     2*Epoch     LFW     CFP-FP     AgeDB-30     IJB-B     IJB-C

3c$1:1$ Verification Accuracy     2cTAR@FAR=1e-4

CosFace [25]            Private                   64-layer CNN     8              21     99.81     98.12      98.11     94.80     96.37

ArcFace [3]            MS1MV2                    ResNet100        4              16     99.83     98.27      98.28     94.25     96.03

CurricularFace [12]     MS1MV2                    ResNet100        4              24     99.80     98.37      98.32     94.80     96.10

Sub-center [5]          MS1MV3                    ResNet100        8              24     99.80     98.80      98.31     94.94     96.28

MagFace [19]             MS1MV2                    ResNet100        8              25     99.83     98.46      98.17     94.08     95.97

ElasticFace [2]      MS1MV2                    ResNet100        4              26     99.80     98.73      98.28     95.43     96.65

VPL [7]                 Cleaned MS1M              ResNet100        8              20     99.83     99.11      98.60     95.56     96.76

F$^{2}$C [26]            MS1MV2                    ResNet50         8              20     99.50     98.46      97.83     -         94.91

F$^{2}$C [26]            WebFace42M       ResNet100        8              20     99.83     99.33      98.33     -         97.31

Virtual FC [15]          Cleaned MS1M              ResNet100        1     16     99.38     95.55      -         61.44     71.47

PFC-0.3 [1]               WebFace4M                 ResNet100        8              20     99.85     99.23      98.01     95.64     96.80

PFC-0.3 [1]               WebFace12M                ResNet100        8              20     99.83     99.40      98.53     96.31     97.58

PFC-0.3 [1]               WebFace42M       ResNet100        32             20     99.85     99.40      98.60     96.47     97.82

MHLR     MS1MV2         ResNet100     1     5     99.80     98.39      98.11     94.98     96.33

MHLR     MS1MV3         ResNet100     1     5     99.80     98.53      98.12     95.03     96.41

MHLR     WebFace4M      ResNet100     1     5     99.65     98.63      97.57     94.36     96.18

MHLR     WebFace8M      ResNet100     1     5     99.85     99.24      97.90     95.58     97.09

MHLR     WebFace12M     ResNet50      1     5     99.80     99.09      97.92     95.25     96.85

MHLR     WebFace12M     ResNet100     1     5     99.83     99.17      98.02     95.68     97.20

Performance comparisons between MHLR and other state-of-the-art FR methods on various benchmarks. For LFW, CFP-FP, and AgeDB-30, 1:1 verification accuracy ($\%$) is reported. For IJB-B and IJB-C, TAR@FAR$=1e-4$ is reported as the accuracy ($\%$).

The threshold $\lambda$ and tolerance $\tau$ in MHLR In Fig. <ref>, we grid search the threshold $\lambda$ and tolerance $\tau$ in MHLR. All results are reported based on ResNet100 trained for $5$ epochs on MS1MV3 with $1\times$GPU, while the test set is IJB-C. We witness an accuracy decline while $\tau$ is decreasing, which indicates a small $\tau$ is too sensitive and hurried to schedule the learning rate. To the opposite, a big $\tau$ is slow in stationary subsequence detection. We also observe that $\lambda=5e-5$ has the highest accuracy. The reason is, as shown in Fig. <ref>(b), a small $\lambda$ will filter the majority of $\mathcal{D}^{EMA}_{t}$ while a big $\lambda$ considers too many $\mathcal{D}^{EMA}_{t}$ as the signal of scheduling the learning rate. In summary, we set $\lambda=5e-5$ and $\tau=0.05$ in this paper.

§.§ Comparison with SOTA methods

In Tab. <ref>, we compared MHLR with other state-of-the-art (SOTA) methods. Other methods employ at least $\times4$ GPUs and averagely $20$ epochs, while MHLR only uses $x1$ GPU and $5$ epochs. Thus, it requires only $1/4$ of the training time compared to other methods. For example, it only requires $30$H to train ResNet100 on WebFace12M with MHLR, as shown in Tab. <ref>. Comparing the accuracy of MHLR with other methods, its performance is competitive. For instance, on WebFace12M, the accuracy of MHLR is $97.20$ while that of PFC-0.3 is $97.58$. Yet, PFC-0.3 needs to be trained with $8\times$GPUs for $20$ epochs, which is around $24$ hours for training. Thus, we can conclude that it is meaningful doing large-scale training with $1\times$GPU.

In addition, PFC-0.3 and F$^{2}$C also train ResNet100 on WebFace42M, which contains $42$M images with $2$M identities. The accuracy of PFC-0.3 and F$^{2}$C on IJB-C are $97.82$ and $97.31$ respectively. In order to reach such a high performance, their training cost, including time and hardware resources, is huge. To be specific, it takes PFC nearly $25$H to train ResNet100 on WebFace42M with $32\times$GPUs, and it takes F$^{2}$C days to train it on WebFace42M with only $8\times$GPUs. Nonetheless, MHLR is able to train ResNet100 on WebFace12M, not Webface42M, to get $97.20$ on IJB-C by $30$H with only $1\times$GPU. The accuracy gap is in $1\%$. Therefore, we conclude that large-scale FR training now faces the law of diminishing marginal utility, which means the cost increase rapidly in order to improve a small amount of the performance for FR models.

Please note Virtual FC also tries to train ResNet100 on self-cleaned MS1M with $1\times$GPU. Yet, it trains the model for $16$ epochs, which requires more training time compared to MHLR for only $5$ epochs. Besides, its accuracy on IJB-B and IJB-C has a large gap compared to other methods as reported in [1]. By contrast, on MS1MV2 and MS1MV3, MHLR is competitive with other methods.

§ CONCLUSION

In this paper, we propose Moving Haar Learning Rate (MHLR) scheduler for large-scale face recognition training. MHLR schedules the learning rate promptly and accurately in order to reduce the converging time of FR models. As a result, MHLR is able to train the model with $1/4$ of its original training time on $1\times$GPU by sacrificing less than $1\%$ accuracy. We conclude that large-scale face recognition training now faces the law of diminishing marginal utility, which means the cost increase rapidly in order to improve a small amount of the performance for FR models.

Limitations. We have not implemented PFC yet, and thus training models on WebFace42M is still infeasible with $1\times$GPU. In the future, we will implement PFC+MHLR to validate the performance of trained models on WebFace42M with only $1\times$GPU. Afterwards, we will also conduct experiments on Masked Face Recognition (MFR) challenge [6].

Potential Societal Impacts. MHLR opens the gate of large-scale face recognition training to many researchers who do not have massive computing power, which is caused by lacking of expensive hardware sources. By attracting more researchers, this research direction will become thriving.

§ ACKNOWLEDGMENTS

This work was supported in part by the Guangdong Basic and Applied Basic Research Foundation under Grant 2022A515110020, the Guangzhou Basic and Applied Basic Research Foundation under Grant 202201010476, the Fundamental Research Funds for the Central Universities under Grant 21621017, and the National Natural Science Foundation of China under Grant 62271232.

[1]
Xiang An, Jiankang Deng, Jia Guo, Ziyong Feng, Xuhan Zhu, Jing Yang, and
Tongliang Liu.
Killing two birds with one stone: Efficient and robust training of
face recognition cnns by partial FC.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 4032–4041, 2022.

[2]
Fadi Boutros, Naser Damer, Florian Kirchbuchner, and Arjan Kuijper.
Elasticface: Elastic margin loss for deep face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition Workshops, CVPRW, pages 1577–1586, 2022.

[3]
Jiankang Deng, Jia Guo, Niannan Xue, and Stefanos Zafeiriou.
Arcface: Additive angular margin loss for deep face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 4690–4699, 2019.

[4]
Jiankang Deng, Jia Guo, Debing Zhang, Yafeng Deng, Xiangju Lu, and Song Shi.
Lightweight face recognition challenge.
In IEEE/CVF International Conference on Computer Vision
Workshops, ICCVW, pages 2638–2646, 2019.

[5]
Jiankang Deng, Jia Guo, Tongliang Liu, Mingming Gong, and Stefanos Zafeiriou.
Sub-center arcface: Boosting face recognition by large-scale noisy
web faces.
In Andrea Vedaldi, Horst Bischof, Thomas Brox, and Jan-Michael
Frahm, editors, European Conference on Computer Vision, ECCV, volume
12356, pages 741–757, 2020.

[6]
Jiankang Deng, Jia Guo, Xiang An, Zheng Zhu, and Stefanos Zafeiriou.
Masked face recognition challenge: The insightface track report.
In IEEE/CVF International Conference on Computer Vision
Workshops, ICCVW, pages 1437–1444, 2021.

[7]
Jiankang Deng, Jia Guo, Jing Yang, Alexandros Lattas, and Stefanos Zafeiriou.
Variational prototype learning for deep face recognition.
In IEEE Conference on Computer Vision and Pattern Recognition,
CVPR, pages 11906–11915, 2021.

[8]
Yueqi Duan, Jiwen Lu, and Jie Zhou.
Uniformface: Learning deep equidistributed representation for face
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 3415–3424, 2019.

[9]
Yandong Guo, Lei Zhang, Yuxiao Hu, Xiaodong He, and Jianfeng Gao.
Ms-celeb-1m: A dataset and benchmark for large-scale face
In European Conference on Computer Vision, ECCV, volume 9907
of Lecture Notes in Computer Science, pages 87–102, 2016.

[10]
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.
Deep residual learning for image recognition.
In IEEE Conference on Computer Vision and Pattern Recognition,
CVPR, pages 770–778, 2016.

[11]
Gary B Huang, Marwan Mattar, Tamara Berg, and Eric Learned-Miller.
Labeled faces in the wild: A database forstudying face recognition in
unconstrained environments.
In University of Massachusetts, Technical Report, 2007.

[12]
Yuge Huang, Yuhan Wang, Ying Tai, Xiaoming Liu, Pengcheng Shen, Shaoxin Li,
Jilin Li, and Feiyue Huang.
Curricularface: Adaptive curriculum learning loss for deep face
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 5900–5909, 2020.

[13]
Minchul Kim, Anil K. Jain, and Xiaoming Liu.
Adaface: Quality adaptive margin for face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 18729–18738, 2022.

[14]
Bi Li, Teng Xi, Gang Zhang, Haocheng Feng, Junyu Han, Jingtuo Liu, Errui Ding,
and Wenyu Liu.
Dynamic class queue for large scale face recognition in the wild.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 3763–3772, 2021.

[15]
Pengyu Li, Biao Wang, and Lei Zhang.
Virtual fully-connected layer: Training a large-scale face
recognition dataset with limited computational resources.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 13315–13324, 2021.

[16]
Weiyang Liu, Yandong Wen, Zhiding Yu, and Meng Yang.
Large-margin softmax loss for convolutional neural networks.
In Maria-Florina Balcan and Kilian Q. Weinberger, editors, International Conference on Machine Learning, ICML, volume 48, pages
507–516, 2016.

[17]
Weiyang Liu, Yandong Wen, Zhiding Yu, Ming Li, Bhiksha Raj, and Le Song.
Sphereface: Deep hypersphere embedding for face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 6738–6746, 2017.

[18]
Brianna Maze, Jocelyn C. Adams, James A. Duncan, Nathan D. Kalka, Tim Miller,
Charles Otto, Anil K. Jain, W. Tyler Niggel, Janet Anderson, Jordan Cheney,
and Patrick Grother.
IARPA janus benchmark - C: face dataset and protocol.
In International Conference on Biometrics, ICB, pages
158–165, 2018.

[19]
Qiang Meng, Shichao Zhao, Zhida Huang, and Feng Zhou.
Magface: A universal representation for face recognition and
quality assessment.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 14225–14234, 2021.

[20]
Stylianos Moschoglou, Athanasios Papaioannou, Christos Sagonas, Jiankang Deng,
Irene Kotsia, and Stefanos Zafeiriou.
Agedb: The first manually collected, in-the-wild age database.
In IEEE Conference on Computer Vision and Pattern Recognition
Workshops, CVPRW, pages 1997–2005, 2017.

[21]
Soumyadip Sengupta, Jun-Cheng Chen, Carlos Castillo, Vishal M Patel, Rama
Chellappa, and David W Jacobs.
Frontal to profile face verification in the wild.
In IEEE Winter Conference on Applications of Computer Vision
(WACV), pages 1–9, 2016.

[22]
Ira Kemelmacher Shlizerman, Steven M. Seitz, Daniel Miller, and Evan Brossard.
The megaface benchmark: 1 million faces for recognition at scale.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 4873–4882, 2016.

[23]
Yifan Sun, Changmao Cheng, Yuhan Zhang, Chi Zhang, Liang Zheng, Zhongdao Wang,
and Yichen Wei.
Circle loss: A unified perspective of pair similarity optimization.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 6397–6406, 2020.

[24]
Feng Wang, Weiyang Liu, Hanjun Dai, Haijun Liu, and Jian Cheng.
Additive margin softmax for face verification.
In International Conference on Learning Representations
Workshops, ICLRW, 2018.

[25]
Hao Wang, Yitong Wang, Zheng Zhou, Xing Ji, Dihong Gong, Jingchao Zhou, Zhifeng
Li, and Wei Liu.
Cosface: Large margin cosine loss for deep face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 5265–5274, 2018.

[26]
Kai Wang, Shuo Wang, Panpan Zhang, Zhipeng Zhou, Zheng Zhu, Xiaobo Wang,
Xiaojiang Peng, Baigui Sun, Hao Li, and Yang You.
An efficient training approach for very large scale face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 4073–4082, 2022.

[27]
Cameron Whitelam, Emma Taborsky, Austin Blanton, Brianna Maze, Jocelyn C.
Adams, Tim Miller, Nathan D. Kalka, Anil K. Jain, James A. Duncan, Kristen
Allen, Jordan Cheney, and Patrick Grother.
IARPA janus benchmark-b face dataset.
In IEEE Conference on Computer Vision and Pattern Recognition
Workshops, CVPRW, pages 592–600, 2017.

[28]
Dong Yi, Zhen Lei, Shengcai Liao, and Stan Z. Li.
Learning face representation from scratch.
CoRR, abs/1411.7923, 2014.

[29]
Zheng Zhu, Guan Huang, Jiankang Deng, Yun Ye, Junjie Huang, Xinze Chen, Jiagang
Zhu, Tian Yang, Jiwen Lu, Dalong Du, and Jie Zhou.
Webface260m: A benchmark unveiling the power of million-scale deep
face recognition.
In IEEE/CVF Conference on Computer Vision and Pattern
Recognition, CVPR, pages 10492–10502, 2021.