11institutetext:  Université de Paris, Institut de Physique du globe de Paris,
CNRS UMR 7154, F-75005 Paris, France

# Triskels and Symmetries of Mean Global Sea-Level Pressure

F. Lopes 11 V. Courtillot 11 J-L. Le Mouël 11

The evolution of mean sea-level atmospheric pressure since 1850 is analyzed
using singular spectrum analysis. Maps of the main components (the trends)
reveal striking symmetries of order 3 and 4. The northern hemisphere (NH)
displays a set of three positive features, forming an almost perfect
equilateral triangle. The southern hemisphere (SH) displays a set of three
positive features arranged as an isosceles triangle, with a possible fourth
(weaker) feature. This geometry can be modeled as Taylor-Couette flow of mode
3 (NH) or 4 (SH). The remarkable regularity and order three symmetry of the NH
triskel occurs despite the lack of cylindrical symmetry of the northern
continents. The stronger intensity and larger size of features in the SH is
linked to the presence of the annular AAO. In addition to the dominant trends,
quasi-periodic components of $\sim$130, 90, 50, 22, 15, 4, 1.8, 1, 0.5, 0.33,
and 0.25 years, i.e. the Jose, Gleissberg, Hale and Schwabe cycles, the annual
cycle and its first three harmonics are identified.

###### Key Words.:

Global Sea Level Pressure, Taylow-Couette flow, Triskel pattern

##  1 Introduction

Understanding the mean circulation of air masses is one of the most ancient
problems in meteorology (cf. Lorentz (1967); Lindzen et Hou (1988)). Solar
insolation drives the first order structure and motions of atmospheric masses.
Early views thought them to be organized in planetary scale circulation cells,
arranged roughly symmetrically with respect to the equator. In the
troposphere, warm winds blow from the equator, forming the top of the Hadley
cells (cf. Hadley (1735)); air cools, becomes denser and sinks near 30° (N and
S), generating two belts of high pressures. The circuit is closed by return
winds blowing towards the (low pressure) equator. Cold winds in the lower
atmosphere blow from the poles to warmer regions, their density decreases as
they approach 60° (N and S) latitude, where they rise. The circuit is closed
by return winds in the troposphere, back to the polar high pressures, forming
the polar cells. The Ferrel cells (cf. Ferrel (1856)) located between the
polar and Hadley cells, extend between 30° and 60° (N and S).

Lindzen et Hou (1988) summarize the situation at the end of the first half of
the 20th century: “By the early part of this century (Jeffreys 1926), the idea
was being put forth that the zonally averaged circulation might, in large
measure, be forced by eddies” …“Starr (1948) was going so far as to suggest
that the symmetric circulation was inconsequential”. Schneider and Lindzen
(1977) and Schneider (1977) were probably the first to propose an explicit
calculation of the symmetric circulation; they showed that purely symmetric
circulations could maintain strong subtropical jets and contribute to the
maintenance of surface winds. The symmetry of the mean heat received by Earth
was sufficient to force the symmetry of the cell structure, the symmetry of
the subtropical currents and the general wind pattern.

Based on a decade of observations (1963-1973), Lindzen et Hou (1988) found
that as soon as the peak heating is a few degrees in latitude off the equator,
profound asymmetries in the Hadley circulation result, with the summer cell
becoming negligible. The annually averaged Hadley circulation is much larger
than that forced by the annually averaged heating.

We know that the heat coming from the Sun’s activity varies with time (e.g.
Usoskin et al. (2003); Solanki et al. (2004); Lockwood et al. (2011); Vieira
et al. (2011); Abreu et al. (2012); Kutiev et al. (2013); Thuillier et al.
(2014); Le Mouël et al. (2020a, b); Courtillot et al. (2021)) and so does the
inclination of our planet’s rotation axis (e.g. Barnes et al. (1983);
Rochester (1977); Lambeck (2005); Schindelegger et al. (2013); Lopes et al.
(2017); Barkin et al. (2019); Le Mouël et al. (2019a); Krylov et al. (2020);
Le Mouël et al. (2021a); Lopes et al. (2021)). Consequences of these two
variabilities are expected to affect the atmosphere, and in particular the
climate (e.g. Mörth et Schlamminger (1979); Barnes et al. (1983); Lindzen
(1994); Solanki et al. (2004); Gray et al. (2010, 2013); Roy et Haigh (2010);
Abreu et al. (2012); Kutiev et al. (2013); Schindelegger et al. (2013);
Johnstone and Mantua (2014); Le Mouël et al. (2019b); Gruzdev and Bezverkhnii
(2020); Cionco et al. (2021); Connolly et al. (2021); Drews et al. (2021);
Sonechkin et Vakulenko (2021)). To these already complex interactions, one
must add those due to the ocean, the largest heat exchanger that warms and
cools on its own (longer) time scales. Lindzen (1978) was the first to show
how (to first order) atmospheric circulation was forced by the physics of the
ocean surface.

Interest in recent climate warming has led to question how the convection
cells, and in particular the Hadley cells, being the main ones, react to a
temperature increase (e.g. Chang (1995); Dima et Wallace (2003); Frierson et
al. (2007); Hu et Fu (2007); Kharin et al. (2007); Lu et al. (2007); Tandon et
al. (2013); Shepherd (2014); Tao et al. (2016); Grise and Davis (2020)). Also,
how do these variations affect extreme meteorological events (e.g. Schaeffer
et al. (2005); Stott et al. (2010); Rahmstorf and Coumou (2011); Rummukainen
(2012); Trenberth et al. (2015))? Most authors who have studied the space-time
evolution of convection cells call upon the analysis in terms of components of
sea-level pressure (SLP). SLP is directly related to climate indices (e.g.
Vaideanu et al. (2020); Le Mouël et al. (2019b)). Decomposition methods often
use “Empirical Orthogonal Functions” (EOF). In this paper, we select an
alternate method: the time decomposition of SLP using “Singular Spectrum
Analysis” (SSA, e.g. Golyandina and Zhigljavsky (2013), a method that we have
used extensively with success in some recent works (e.g. Lopes et al. (2017);
Le Mouël et al. (2019a, 2020a, 2020b); Courtillot et al. (2021); Le Mouël et
al. (2021a); Lopes et al. (2021)

The SLP data we use are described in section 2, their analysis using SSA in
section 3 and the results are discussed in section 4.

##  2 The SLP Data

The SLP data are maintained by the Met Office Hadley
Centre111https://www.metoffice.gov.uk/hadobs/hadslp2/data/download.html. They
are available in map form for global pressure every month since 1850 to the
present, with a spatial sampling of 5°x5°. As indicated by the Met Office, the
series is not homogeneous in time (two intervals from 1850 to 2004 and 2005 to
2020 in which the means are homogeneous but the variances differ). Allan and
Ansell (2006) discuss the locations and number of ground observations they use
to build the series. We show the map of 1850 to 2020 mean pressures in Figure
1.

Figure 1: 1850 to 2020 mean pressure at sea level (SLP) (Hammer-Aitoff
projection).

The map of Figure 1 is very similar to that obtained from simply averaging
monthly SLP maps, demonstrating the very stable geometry of the global
atmospheric structure. Well-known patterns are clearly rendered by the map: in
the Southern hemisphere, the structure is dominated by three large positive
features in the southern parts of the main oceans (Pacific, Atlantic and
Indian) separated by the southern ends of the three southern continents (South
America, South Africa and Australia). This approximate three-fold symmetry
actually seems to “leak” into a four-fold symmetry: the positive feature
extending from Australia to western South America can be described as
exhibiting two weaker features, the one over Australia and East of it being
the weakest.

The pattern in the Northern hemisphere is similar to that in the Southern
hemisphere, with two features extending over the northern Pacific and Atlantic
oceans, similar to the southern structure, but with the third positive anomaly
lying over the Asian continents (central Asia and Tibet). This overall
structure of global positive pressure features could be roughly described as
the intersection of a series of three strong and one weaker cylinder, parallel
to the Earth’s axis of rotation, with the earth’s surface. South of 40°S
latitude, the features turn negative and form the quasi-circular Antarctic
oscillation (AAO), also known as the Southern annular mode (SAM), a belt of
low pressures surrounding the frozen continent (e.g. Marshal (2003); Gillett
et al. (2006)). A zonal average of Figure 1 is not really consistent with the
classical latitudinal description of the Hadley, Ferrel and polar cells, with
relatively moderate low pressures over the equatorial belt, high pressures
around 30°N and lower values at higher latitudes, particularly over Antarctica
(see also Figure 4). The overall averaged zonal structure is not symmetrical
with respect to the equator (Lindzen et Hou 1988).

Winds are proportional to the space derivative of pressure (Figure 2, which is
the map most people are familiar with), and thus in space phase quadrature
with pressure features. The space derivative is a high-pass operator which
exacerbates discontinuities and promotes oscillations one order higher than
the causative features. Figure 2 shows how winds (the gradients of pressure)
are slowed or suppressed by continental masses; also, the red and blue ”
anomalies ” are often paired, which reminds one of Laplace’s statement that
the vector sum of winds must be close to zero at any instant.

Figure 2: 1850 to 2020 time average of the spatial pressure gradients (i.e.
winds up to a constant) at sea level (Hammer-Aitoff projection).

##  3 The SSA of SLP Data

Rather than building a map of spatial likelihood of pressure structures at a
given time, we build the time series of pressure for each couple of
(longitude, latitude) coordinates. These are analyzed using the method of
Singular spectrum analysis (iSSA , e.g. Golyandina and Zhigljavsky (2013)):
each time series is decomposed into a sum of a trend and periodic or quasi-
periodic components.

Note: the trend is the first component “extracted” by SSA, it is a sort of
mean with respect to time, a non-oscillating function of time; it should not
be confused with the spatial mean of maps in time (see for example Figure 2).

The trend ($\sim$1008 hP), that is the first and largest SSA component,
represents more than 70% of the total variance (sv) of the original series.
The sequence of the next quasi-periodic components is, in decreasing order of
periods, $\sim$130 years ($\sim$0.7 hP, $\sim$0.06% of the sv; Jose (1965)),
$\sim$90 yr ($\sim$21 hP, $\sim$1.9% of the sv; Gleissberg (1939); Le Mouël et
al. (2017)), $\sim$50 yr ($\sim$0.2 hP, $\sim$0.02% of the sv), $\sim$22 yr
($\sim$0.50 hP, $\sim$0.04% of the sv; Hale cycle, Usoskin (2017)), $\sim$15
yr ( $\sim$0.2 hP, $\sim$0.02% of the sv, upper bound of the Schwabe cycle;
Schwabe (1844), $\sim$4 yr ($\sim$0.3 hP, $\sim$0.03% of the sv) , $\sim$1.8
yr ($\sim$0.3 hP, $\sim$0.03% of the sv), then 1 yr ( $\sim$93 hP, $\sim$8.3%
of the sv), 0.5 yr ($\sim$65 hP, $\sim$5.8% of the sv), 0.33 yr ($\sim$44 hP,
$\sim$3.9% of the sv) and 0.25 yr ($\sim$21 hP, $\sim$ 1.9% of the sv).

These are readily recognized as the Jose, Gleissberg, Hale and Schwabe cycles,
the quasi biennal oscillation and the annual cycle followed by its first three
harmonics. Figure 3 represents the mean (from 1850 to 2020) of the SSA trends
of SLP in Hammer-Aitoff and North and South stereographic polar projections.
This mean is representative, since the observed (decreasing) overall variation
of pressure is only 0.1% in the 170 years of the study. We illustrate the
stability of the SLP trend both in space and time with Figure 4, which shows
the variation in the trend of SLP since 1850 as a function of latitude.
Indeed, the latitudinal structure of mean pressure (SSA trend) is remarkably
stable over 170 years (a rather rare feature in geophysics and geodynamics,
both internal and external to the solid Earth), with its typical asymmetry
between the two hemispheres.

Figure 3: 1850 to 2020 mean of trends of pressure at sea level (SLP) from,
determined using Singular spectrum analysis (SSA) (top: Hammer-Aitoff
projection; lower left: stereographic projection of the NH; lower right:
stereographic projection of the SH).

The stereographic projections of Figure 3 give a clear image of the large
scale atmospheric circulation. They are close to the original mean of Figure
1, and the same large scale features appear, with sharper contours. The polar
projections are particularly revealing. There is a strong negative feature
south of Greenland. In the northern hemisphere, three sharp positive features
lie on a circle at 40°N latitude and are located at the three apices of an
equilateral triangle. The center of this 3-fold symmetry (that we will call a
“triskel”, a Celtic symbol) is located near 77°N, 90°W (13° away from the
North pole). In the southern hemisphere, the pattern is closer to 4-fold
symmetry, with one of the four features much weaker than the other three.

Figure 4: Evolution of the SSA trend of SLP as a function of time (170 yr
since 1850) and latitude, illustrating the stability of the zonal pressure
structure in time.

##  4 Discussion

The physics that is appropriate to analyze the motions of fluid masses in the
atmosphere, ocean (and mantle) is that of stationary turbulent flow (e.g.
Chandrasekhar (1961); Frisch et Kolmogorov, (1995)). The large Hadley, Ferrel
and polar cells are generally interpreted in terms of Taylor-Couette flow.
Analytical solutions are available in the case of a cylindrical geometry (e.g.
Taylor (1923)), but raise problems in the spherical case (e.g. Schrauf (1986);
Mamun et Tuckerman (1995); Nakabayashi et Tsuchida (1995); Hollerbach et al.
(2006); Malhoul2016; Garcia et al. (2019); Mannix et Mestel (2021)).

Forbes and Bassom (2018) show that a simple solution for Taylor-Couette flow
of a viscous fluid between two concentric cylinders with radii $a$ and $b$
rotating in opposite directions and submitted to an azimuthal perturbation
with angle $\theta$, can be written as:

$r(\xi)=r_{n}+\varepsilon a\cos(q\xi)$ (1a)
$\theta(\xi,t)=\xi-(b^{2}+a^{2})*\omega*t*[1-\dfrac{r_{n}^{2}}{r^{2}(\xi)}]/(b^{2}-a^{2})$
(1b)

in which $r_{n}$ is the neutral radius and $q$ is an integer corresponding to
the flow mode. One has:

$r_{n}=\sqrt{\dfrac{2a^{2}*b^{2}}{a^{2}+b^{2}}}$

The solutions of system (1a,1b) are displayed in Figure 5 for stationary flows
with integers q = 1 to 4, with a = 0.5, b = 1, $\omega*t$ = 1 and
$\varepsilon$ = 0.1. Modes corresponding to parameter $q$ evolve and $q$ is
fixed by the nature and symmetry of the flow.

Figure 5: Solutions of system of equations (1a,1b), for integer mode values q
= 1 to 4 (see text).

Figure 6(a) shows the superimposition of a pattern of Taylor-Couette flow with
mode $q$ = 3 and the SLP trend map of Figure 3 (lower left) for the northern
hemisphere. We invert equations (1a,1b) by simulated annealing (Kirkpatrick et
al. 1983) for $\omega*t$ = 1 (flow stationary and constant in time). The
amplitude of the flow perturbation $\varepsilon$ is on the order of 0.100
$\pm$ 0.002, and the triskel pattern is centered on 91.3 $\pm$ 0.1°W - 77.1 ±
0.2°N. The fit is excellent and compatible with the shearing due to the
direction of winds (East to West) at mid latitudes, opposite to the direction
of Earth’s rotation from West to East.

Figure 6(b) shows the superimposition of a pattern of Taylor-Couette flow with
mode $q$ = 4 and the SLP trend map of Figure 3 (lower right) for the southern
hemisphere. The amplitude of the flow perturbation $\varepsilon$ is reduced to
0.05 $\pm$ 0.01, and the symmetry 4 is centered on 88.1 $\pm$ 0.2°W - 7.5
$\pm$ 0.1°N.

(a) Superimposition of the triskel Taylor-Couette pattern of flow for mode $q$
= 3 and observed mean pattern of sea-level pressure in the North hemisphere
since 1850

(b) Superimposition of the Taylor-Couette pattern of flow for mode $q$ = 4 and
observed mean pattern of sea-level pressure in the South hemisphere since
1850.

Figure 6: Superimposition of the Taylor-Couette patterns both North and South
hemispheres.

The Hadley, Ferrell or polar type cells change with season. This is why the
zonal structure of the atmosphere is shown in seasonal (monthly, bi-monthly)
maps that can be averaged over several years. These oscillations cannot be
analyzed through spatial decomposition methods such as EOF or spherical
harmonics. On the other hand, SSA is a scalar time analysis and can extract
from the data temporal features. We have shown that atmospheric surface
pressures display simple geometric/geographic features that are stable both in
space and time. The oscillatory components belong to the harmonic families of
the Schwabe (33, 22, 11, 5.5, 3.3 and 1.6 yr) and annual (1, 1/2, 1/3, 1/4 and
1/5 yr) cycles, plus the dominant 90 yr Gleissberg cycle.

Finally, one can only be struck by the kinship between the axial cylinders
that form the main components of the sea-level pressure structure (with the
quasi-symmetry of the outline of these cylinders, yet the failure of perfect
symmetry, with 3 cylinders in the northern hemisphere and more like 4 in the
southern hemisphere). These are geometrically reminiscent of the axial
cylinders found by Gubbins and Bloxham (1987) and Gubbins and Kelly (1993) in
their analyses of flow in the Earth’s molten outer core that generates the
geomagnetic field. The physics in both cases appears quite distinct but the
analogy may indicate a path to finding the mechanism that generates these
features.

##  5 Conclusion

In this note, we have analyzed the time evolution of mean sea-level pressure
SLP since 1850, using the iSSA method of spectral analysis. The method
decomposes the original time series into a set of time dependent components
and a trend. The quasi-periodic components are in decreasing order of periods,
cycles of approximately 130, 90, 50, 22, 15, 4, 1.8, 1, 0.5, 0.33, and 0.25
years. These are recognized as the Jose, Gleissberg, Hale and Schwabe “long”
(decadal and multi-decadal) cycles (with some harmonics of the well-known
Schwabe “solar” cycle222The Schwabe quasi-cycle is traditionally given as
being 11 yr long but actually spans the 9 to 14 yr period range. ) and the
“short” annual cycle and its first three harmonics. These periods are
encountered in many solar and terrestrial phenomena and are now attributed to
commensurable periods of the Jovian planets (Mörth et Schlamminger 1979; Lopes
et al. 2021).

As far as the general trends of sea-level pressure extracted with SSA are
concerned, they have fluctuated by less than 0.1% in 170 years (or 6.10 -3
hP/yr). The mid latitudes of each hemisphere of the Earth’s surface harbor a
set of positive features (on the order of 20 hP), all but one (over central
Asia and Tibet) lying over the main ocean basins (Figure 3, top). This is in
full agreement with Schneider and Lindzen (1977), who see the oceans as the
main engine driving the 194 general circulation of the atmosphere. In
stereographic polar projection, the northern hemisphere displays a set of
three positive features, forming an almost perfect equilateral triangle or
triskel (Figure 03, lower left). The southern hemisphere also features a set
of three main positive features but arranged as an isosceles triangle, with a
possible fourth (but much fainter) feature (Figure 3, lower right). A
preliminary analysis suggests that the atmosphere in the two hemispheres could
be the site of Taylor-Couette (sheared) differential flow of mode 3 (N
hemisphere) or mode 4 (S hemisphere). The mode is imposed by the geometry of
the oceanic vs continental nature of the Earth’s surface. The remarkable
regularity and order three symmetry of the northern hemisphere triskel occurs
despite the lack of cylindrical symmetry of the northern continents. The
stronger intensity and larger size of features in the southern hemisphere is
due to the presence of the annular flow of cold air over the southern ocean
and the obstacles due to the southern parts or tips of South America, Southern
Africa and Australia.

Singular spectrum analysis is a scalar analysis as a function of time, which
is useful in detecting time varying patterns and features. Spatial components
that SSA extracts are not artefacts, contrary to what could happen when using
spherical harmonics or EOF. A component is present only if there is spatial
coherency in the observational data. In the case of sea-level pressure, we
find stable (non oscillatory) components/structures/features (as we have
called them), The stability of pressure trends over the 170 years of the
analysis implies that the idea of eddies maintaining the cells and forcing
their oscillations is no more warranted. Taylor-Couette flow between rotating
cylinders provides a good fit to the observations. The geometrical analogy
with the cylinders generating the geodynamo may be a promising path to solving
the dynamics of the large scale circulation in the Earth’s atmosphere.

## References

  * Abreu et al. (2012) Abreu, J. A., Beer, J., Ferriz-Mas, A., McCracken, K. G., and Steinhilber, F., “Is there a planetary influence on solar activity?”, Astronomy & Astrophysics,548, A88.
  * Allan and Ansell (2006) Allan, R., and Ansell, T., “A new globally complete monthly historical gridded mean sea level pressure dataset (HadSLP2): 1850–2004”, Journal of Climate, 19(22), 5816-5842.
  * Barnes et al. (1983) Barnes, R. T. H., Hide, R., White, A. A., and Wilson, C. A., “Atmospheric angular momentum fluctuations, length-of-day changes and polar motion”, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences, 387(1792), 31-73.
  * Barkin et al. (2019) Barkin, M. Y., Krylov, S. S. and Perepelkin, V. V., “Modeling and analysis of the Earth pole motion with nonstationary perturbations”, In Journal of Physics: Conference Series, (Vol. 1301, No. 1, p. 012005). IOP Publishing.
  * Chandrasekhar (1961) Chandrasekhar, S., “hydrodynamic and hydromagnetic stability”, Oxford University Press.
  * Chang (1995) Chang, E. K., “The influence of Hadley circulation intensity changes on extratropical climate in an idealized model”, Journal of Atmospheric Sciences, 52(11), 2006-2024.
  * Cionco et al. (2021) Cionco, R. G., Kudryavtsev, S. M. and Soon, W. H., “Possible Origin of Some Periodicities Detected in Solar‐Terrestrial Studies: Earth’s Orbital Movements”, Earth and Space Science, 8(8), e2021EA001805.
  * Connolly et al. (2021) Connolly, R., Soon, W., Connolly, M., Baliunas, S., Berglund, J., Butler, C. J., … and Zhang, W, “How much has the Sun influenced Northern Hemisphere temperature trends? An ongoing debate”, Research in Astronomy and Astrophysics, 21(6), 131, 2021.
  * Courtillot et al. (2021) Courtillot, V., Lopes, F., and Le Mouël, J. L, “On the prediction of solar cycles”, Solar Physics, 296(1), 1-23.
  * Dima et Wallace (2003) Dima, I. M. and Wallace, J. M., “On the seasonality of the Hadley cell”, Journal of the atmospheric sciences, 60(12), 1522-1527.
  * Drews et al. (2021) Drews, A., Huo, W., Matthes, K., Kodera, K. and Kruschke, T., “The Sun’s Role for Decadal Climate Predictability in the North Atlantic”, Atmospheric Chemistry and Physics Discussions, 1-17.
  * Ferrel (1856) Ferrel, W., “Essay on the winds and ocean currents”, Nashville J. of Medicine and Surgery, 11, 287-301.
  * Frisch et Kolmogorov, (1995) Frisch, U. and Kolmogorov, A. N., “Turbulence: the legacy of AN Kolmogorov”, Cambridge university press.
  * Frierson et al. (2007) Frierson, D. M., Lu, J. and Chen, G., “Width of the Hadley cell in simple and comprehensive general circulation models”, Geophysical Research Letters, 34(18).
  * Forbes and Bassom (2018) Forbes, L. K. and Bassom, A. P., “Interfacial behaviour in two-fluid Taylor–Couette flow”, The Quarterly Journal of Mechanics and Applied Mathematics, 71(1), 79-97.
  * Garcia et al. (2019) Garcia F, Seilmayer M, Giesecke A, Stefani F., ”Modulated rotating waves in the magnetised spherical Couette system”, Journal of Nonlinear Science. 29(6), 2735-59.
  * Gillett et al. (2006) Gillett, N. P., Kell, T. D., and Jones, P. D., “Regional climate impacts of the Southern Annular Mode”, Geophysical Research Letters, 33(23).
  * Gleissberg (1939) Gleissberg, W., “A long-periodic fluctuation of the sunspot numbers”, Observatory 62, 158.
  * Golyandina and Zhigljavsky (2013) Golyandina, N. and Zhigljavsky, A., “Singular Spectrum Analysis for time series”, (Vol. 120). Berlin: Springer, 2013
  * Gray et al. (2010) Gray, L. J., Beer, J., Geller, M., Haigh, J. D., Lockwood, M., Matthes, K., … and White, W., “Solar influences on climate”, Reviews of Geophysics, 48(4), 2010.
  * Gray et al. (2013) Gray, L. J., Scaife, A. A., Mitchell, D. M., Osprey, S., Ineson, S., Hardiman, S., … and Kodera, K., “A lagged response to the 11 year solar cycle in observed winter Atlantic/European weather patterns”, Journal of Geophysical Research: Atmospheres, 118(24), 13-405.
  * Grise and Davis (2020) Grise, K. M. and Davis, S. M., “Hadley cell expansion in CMIP6 models”, Atmospheric Chemistry and Physics, 20(9), 5249-5268.
  * Gruzdev and Bezverkhnii (2020) Gruzdev, A. N. and Bezverkhnii, V. A., “Manifestation of the 11-year solar cycle in the North Atlantic climate”, In IOP Conference Series: Earth and Environmental Science, (Vol. 606, No. 1, p. 012018). IOP Publishing.
  * Gubbins and Bloxham (1987) Gubbins, D., and Bloxham, J., “Morphology of the geomagnetic field and implications for the geodynamo“. Nature, 325(6104), 509-511.
  * Gubbins and Kelly (1993) Gubbins, D., and Kelly, P., ”Persistent patterns in the geomagnetic field over the past 2.5 Myr ”, Gubbins, 365(6449), 829-832.
  * Hadley (1735) Hadley, G., “VI. Concerning the cause of the general trade-winds”, Philosophical Transactions of the Royal Society of London, 39(437), 58-62, 1735.
  * Hollerbach et al. (2006) Hollerbach, R., Junk, M. and Egbers, C., “Non-axisymmetric instabilities in basic state spherical Couette flow”, Fluid dynamics research, 38(4), 257.
  * Hu et Fu (2007) Hu, Y. and Fu, Q., “Observed poleward expansion of the Hadley circulation since 1979”, Atmospheric Chemistry and Physics, 7(19), 5229-5236.
  * Jeffreys (1926) Jeffreys, H., “On the dynamics of geostrophic winds”, Quart. J. Roy. Meteor. Soc, 52(217), 85-104, 1926.
  * Johnstone and Mantua (2014) Johnstone, J. A., and Mantua, N. J., “Atmospheric controls on northeast Pacific temperature variability and change, 1900–2012”, Proceedings of the National Academy of Sciences, 111(40), 14360-14365.
  * Jose (1965) Jose, P. D., “Sun’s motion and sunspots”, The Astronomical Journal, 70(3), 193-200.
  * Kharin et al. (2007) Kharin, V. V., Zwiers, F. W., Zhang, X. and Hegerl, G. C., “Changes in temperature and precipitation extremes in the IPCC ensemble of global coupled model simulations”, Journal of Climate, 20(8), 1419-1444.
  * Kirkpatrick et al. (1983) Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P., “Optimization by simulated annealing”, Science, 220(4598), 671-680.
  * Krylov et al. (2020) Krylov, S. S., Perepelkin, V. V. and Filippova, A. S., “Long-Period Lunar Perturbations in Earth Pole Oscillatory Process: Theory”, In Advances in Theory and Practice of Computational Mechanics: Proceedings of the 21st International Conference on Computational Mechanics and Modern Applied Software Systems, (Vol. 173, p. 315). Springer Nature.
  * Kutiev et al. (2013) Kutiev, I., Tsagouri, I., Perrone, L., Pancheva, D., Mukhtarov, P., Mikhailov, A., … and Torta, J. M., “Solar activity impact on the Earth’s upper atmosphere”, Journal of Space Weather and Space Climate, 3, A06.
  * Lambeck (2005) Lambeck, K., “The Earth’s variable rotation: geophysical causes and consequences”, Cambridge University Press.
  * Le Mouël et al. (2017) Le Mouël, J. L., Lopes, F. and Courtillot, V., “Identification of Gleissberg cycles and a rising trend in a 315-year-long series of sunspot numbers”, Solar Physics, 292(3), 43.
  * Le Mouël et al. (2019a) Le Mouël, J. L., Lopes, F., Courtillot, V. and Gibert, D., “On forcings of length of day changes: From 9-day to 18.6-year oscillations”, Physics of the Earth and Planetary Interiors, 292, 1-11, 2019a.
  * Le Mouël et al. (2019b) Le Mouël, J. L., Lopes, F. and Courtillot, V., “A solar signature in many climate indices”, Journal of Geophysical Research: Atmospheres, 124(5), 2600-2619.
  * Le Mouël et al. (2020a) Le Mouël, J. L., Lopes, F., and Courtillot, V., “Characteristic time scales of decadal to centennial changes in global surface temperatures over the past 150 years”, Earth and Space Science, 7(4), e2019EA000671.
  * Le Mouël et al. (2020b) Le Mouël, J. L., Lopes, F., and Courtillot, V., “Solar turbulence from sunspot records”, Monthly Notices of the Royal Astronomical Society, 492(1), 1416-1420.
  * Le Mouël et al. (2021a) Le Mouël, J. L., Lopes, F. and Courtillot, V., “Sea-Level Change at the Brest (France) Tide Gauge and the Markowitz Component of Earth’s Rotation”, Journal of Coastal Research.
  * Lindzen (1978) Lindzen, R. S., “On a calculation of the symmetric circulation and its implications for the role of eddies”, The General Circulation-Theory, Modeling, and Observations, 257-281.
  * Lindzen et Hou (1988) Lindzen, R. S. and Hou, A. V., “Hadley circulations for zonally averaged heating centered off the equator”, Journal of Atmospheric Sciences, 45(17), 2416-2427.
  * Lindzen (1994) Lindzen, R. S., “Climate dynamics and global change”, Annual Review of Fluid Mechanics, 26(1), 353-378.
  * Lorentz (1967) Lorentz, E. N., ”The nature and theory of the general circulation of the atmosphere”, WMO 218, TP-115, 161 pp., World Meteorol., Organ., Geneva, Switzerland.
  * Lockwood et al. (2011) Lockwood, M., Owens, M. J., Barnard, L., Davis, C. J., and Steinhilber, F., “The persistence of solar activity indicators and the descent of the Sun into Maunder Minimum conditions”, Geophysical Research Letters, 38(22).
  * Lopes et al. (2017) Lopes, F., Le Mouël, J. L. and Gibert, D., “The mantle rotation pole position. A solar component”, Comptes Rendus Geoscience, 349(4), 159-164.
  * Lopes et al. (2021) Lopes, F., Le Mouël, J. L., Courtillot, V., and Gibert, D., “On the shoulders of Laplace”, Physics of the Earth and Planetary Interiors, 316, 106693.
  * Lu et al. (2007) Lu, J., Vecchi, G. A. and Reichler, T., “Expansion of the Hadley cell under global warming”, Geophysical Research Letters, 34(6), 2007.
  * Mahloul et al. (2016) Mahloul, M., A. Mahamdia, and Kristiawan, M., ”The spherical Taylor–Couette flow”, European Journal of MechanicsB/Fluids, 59 (2016): 1-6.
  * Mamun et Tuckerman (1995) Mamun, C. K., and Tuckerman, L.S., ”Asymmetry and Hopf bifurcation in spherical Couette flow.”, Physics of Fluids, 7.1 (1995): 80-91.
  * Mannix et Mestel (2021) Mannix, P. M., and A. J.,Mestel, “Bistability and hysteresis of axisymmetric thermal convection between differentially rotating spheres”, Journal of Fluid Mechanics, 911.
  * Marshal (2003) Marshall, G. J., “Trends in the Southern Annular Mode from observations and reanalyses”, Journal of climate,16(24), 4134-4143.
  * Mörth et Schlamminger (1979) Mörth, H. T. et Schlamminger, L., “Planetary motion, sunspots and climate”, In Solar-terrestrial influences on weather and climate, (pp. 193-207). Springer, Dordrecht.
  * Nakabayashi et Tsuchida (1995) Nakabayashi, K. and Tsuchida, Y., “Flow-history effect on higher modes in the spherical Couette system”, Journal of Fluid Mechanics, 295, 43-60.
  * Rahmstorf and Coumou (2011) Rahmstorf, S., and Coumou, D., “Increase of extreme events in a warming world”, Proceedings of the National Academy of Sciences, 108(44), 17905-17909.
  * Rochester (1977) Rochester, M. G., “Causes of fluctuations in the rotation of the Earth”, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences,313(1524), 95-105.
  * Roy et Haigh (2010) Roy, I., and Haigh, J. D., “Solar cycle signals in sea level pressure and sea surface temperature”, Atmospheric Chemistry and Physics, 10(6), 3147-3153.
  * Rummukainen (2012) Rummukainen, M., “Changes in climate and weather extremes in the 21st century”, Wiley Interdisciplinary Reviews: Climate Change, 3(2), 115-129.
  * Schaeffer et al. (2005) Schaeffer, M., Selten, F. M., and Opsteegh, J. D., “Shifts of means are not a proxy for changes in extreme winter temperatures in climate projections”, Climate Dynamics, 25(1), 51-63.
  * Schindelegger et al. (2013) Schindelegger, M., Böhm, S., Böhm, J., and Schuh, H., “Atmospheric effects on Earth rotation”, In Atmospheric Effects in Space Geodesy (pp. 181-231), Springer, Berlin, Heidelberg.
  * Schneider (1977) Schneider, E. K., “Axially symmetric steady-state models of the basic state for instability and climate studies. Part II, Nonlinear calculations”, Journal of Atmospheric Sciences, 34(2), 280-296.
  * Schneider and Lindzen (1977) Schneider, E. K. and Lindzen, R. S., “Axially symmetric steady-state models of the basic state for instability and climate studies. Part I. Linearized calculations”, Journal of Atmospheric Sciences, 34(2), 263-279.
  * Schrauf (1986) Schrauf, G., “The first instability in spherical Taylor-Couette flow”, Journal of Fluid Mechanics, 166, 287-303.
  * Schwabe (1844) Schwabe, H., “Sonnen-Beobachtungen im Jahre 1843”, Astron. Nachr., 21(495), 233–236.
  * Shepherd (2014) Shepherd, T. G., “Atmospheric circulation as a source of uncertainty in climate change projections”, Nature Geoscience, 7(10), 703-708.
  * Solanki et al. (2004) Solanki, S. K., Usoskin, I. G., Kromer, B., Schüssler, M., and Beer, J., “Unusual activity of the Sun during recent decades compared to the previous 11,000 years”, Nature, 431(7012), 1084-1087.
  * Sonechkin et Vakulenko (2021) Sonechkin, D. M., and Vakulenko, N. V., “Polyphony of Short-Term Climatic Variations”, Atmosphere, 12(9), 1145.
  * Starr (1948) Starr, V. P., ”An essay on the general circulation of the earth atmosphere”, Journal of Atmospheric Sciences, 5.2: 39-43, 1948.
  * Stott et al. (2010) Stott, P. A., Gillett, N. P., Hegerl, G. C., Karoly, D. J., Stone, D. A., Zhang, X., and Zwiers, F., “Detection and attribution of climate change: a regional perspective”, Wiley Interdisciplinary Reviews: Climate Change, 1(2), 192-211, 2010.
  * Tandon et al. (2013) Tandon, N. F., Gerber, E. P., Sobel, A. H., and Polvani, L. M., “Understanding Hadley cell expansion versus contraction: Insights from simplified models and implications for recent observations”, Journal of climate, 26(12), 4304-4321.
  * Tao et al. (2016) Tao, L., Hu, Y., and Liu, J., “Anthropogenic forcing on the Hadley circulation in CMIP5 simulations”, Climate dynamics, 46(9-10), 3337-3350.
  * Taylor (1923) Taylor, G. I., VIII. “Stability of a viscous liquid contained between two rotating cylinders”, Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 223(605-615), 289-343.
  * Trenberth et al. (2015) Trenberth, K. E., Fasullo, J. T., and Shepherd, T. G., “Attribution of climate extreme events”, Nature Climate Change, 5(8), 725-730.
  * Thuillier et al. (2014) Thuillier, G., Bolsée, D., Schmidtke, G., Foujols, T., Nikutowski, B., Shapiro, A. I., … and Schmutz, W., “The solar irradiance spectrum at solar activity minimum between solar cycles 23 and 24”, Solar Physics, 289(6), 1931-1958.
  * Usoskin et al. (2003) Usoskin, I. G., Solanki, S. K., Schüssler, M., Mursula, K. and Alanko, K. “Millennium-scale sunspot number reconstruction: Evidence for an unusually active Sun since the 1940s”, Physical Review Letters, 91(21), 211101.
  * Usoskin (2017) Usoskin, I. G., “A history of solar activity over millennia”, Living Reviews in Solar Physics, 14(1), 1-97.
  * Vaideanu et al. (2020) Vaideanu, P., Dima, M., Pirloaga, R. and Ionita, M., “Disentangling and quantifying contributions of distinct forcing factors to the observed global sea level pressure field”, Climate Dynamics, 54(3), 1453-1467.
  * Vieira et al. (2011) Vieira, L. E. A., Solanki, S. K., Krivova, N. A., and Usoskin, I. “Evolution of the solar irradiance during the Holocene”, Astronomy & Astrophysics, 531, A6, 2011.


# Factor Graph Molecule Network for Structure Elucidation

Le Trung Hieu, 1 Xu Yiqing, 1 Lee Wee Sun, 1

###### Abstract

Designing a network to learn a molecule structure given its physical/chemical
properties is a hard problem, but is useful for drug discovery tasks. In this
paper, we incorporate higher-order relational learning of Factor Graphs with
strong approximation power of Neural Networks to create a molecule-structure
learning network that has strong generalization power and can enforce higher-
order relationship and valence constraints. We further propose methods to
tackle problems such as the efficient design of factor nodes, conditional
parameter sharing among factors, and symmetry problems in molecule structure
prediction. Our experiment evaluation shows that the factor learning is
effective and outperforms related methods.

## Introduction

The identification of a Chemistry Molecule Structure given its measured
physical/chemical properties is an important task for drug discovery.
Discovering the structure will aid scientists to have better information about
how the drugs will behave in various scenarios.

In this paper, we will elucidate chemistry structure by leveraging the
measurement of the molecule mass spectrum, which is a plot of the ion signal
as a function of the mass-to-charge ratio, an indication of detected ions
abundance. Traditional methods require a chemist to match components of
spectrum against a library of known spectra, then subsequently do result
combination. Such process is tedious and time-consuming. Figure 1 contains an
example mass spectrum and the corresponding molecule structures. The task is,
given the mass spectrum and the molecule formula (e.g C4H5O6), we need to
predict the chemical structure of the compound.

Figure 1: Mass Spectrum

(Kitson, Larsen, and McEwen 1996)

Works have been done for the above task using traditional Machine Learning
methods (Lim et al. 2018), but results are still not satisfactory. Deep
Learning has led to breakthroughs in Computer Vision (Krizhevsky, Sutskever,
and Hinton 2012) and Natural Language Processing (Vaswani et al. 2017), so it
is natural to ask whether it would perform well for other areas such as Drug
Discovery. However, many of the promising works in this field (Maziarka et al.
2020) mainly focus on predicting the molecule properties given its structure,
instead of structure elucidation given some observed measurements. Indeed,
latter task is much harder than the former, because of the large search space
given the large number of potential atom edge combinations, and the lack of
effective learning methods to capture the higher-order relationship between
the nodes as well as the valence constraints.

To tackle the challenge, in this paper we propose the Molecule Factor Graph
Network (MFGN), which incorporates Factor Graph into a Neural Network to
better leverage the higher-order relationship learning between the nodes. We
adapt the Low Rank Belief Propagation Neural Network (Dupty and Lee 2020) to
the chemical structures prediction task by proposing various effective design
for different types of graph nodes and factor nodes tailored for the task.
These architecture elements capture the higher-order relationships and
enforces hard constraints such as the valence constraint of each atom. We also
remove the need for substructure generated by external softwares, relying
mainly on a learning approach. The experiment result shows performance
improvement using the proposed model. We also conduct various ablation study
and analysis to understand how each of our components behave.

## Related Work

### Chemical Structure Elucidation (CSE)

One notable work on molecular structure prediction is computer-aided CSE from
mass spectrometry (MS) by matching substructures (Lim et al. 2018), because
substructures might be responsible for certain m/z peaks in a mass spectra. In
early work (Varmuza and Werther 1996), the authors developed a set of mass
spectral classifiers to recognize the presence or absence of substructures or
more general structural properties in a molecule, then they used the predicted
substructures to reduce the number of candidates. Later work (Kerber et al.
2001)(Schymanski et al. 2008)(Schymanski, Meringer, and Brack 2011) combined a
software named MOLGEN-MS (Kerber et al. 2001) with additional structural
information and made the number of candidates much smaller. The latest state-
of-the-art work (Lim et al. 2018) combines the substructure classification
results to the final list of candidate structures through graph isomorphism
methods.

Such approaches require domain knowledge to generate the possible
substructures to feed the classifiers, usually through third-party softwares.
Our paper contributes a new approach that eliminates the need to rely on
external software for the task.

### Deep Learning for CSE

Deep Learning has been successfully applied to the field of Chemistry and Drug
Discovery, but limited to molecule property prediction task given its known
structures. Two current state-of-the-art methods in this field are variants of
Message Passing Neural Network (MPNN) (Gilmer et al. 2017) and Transformer
(Maziarka et al. 2020), in which their neural architectures allow efficient
interaction between atoms based on the bond and distance of pairwise edge
connection. We will extend the MPNN and Transformer application above to our
task, as the baselines.

In general, the extension of deep learning for structure prediction is
difficult (Lim et al. 2018) with current state-of-the-art models that use
Graph Neural Network (GNN) style architecture. The reason is that, GNNs mostly
capture first order approximations in aggregating information from the
neighbouring nodes, but fails to capture the more complex higher-order
constraints in a molecule structure prediction task (e.g valence constraints
of each atom, or the joint-information from each mass spectrum peak to a set
of other atoms that represent a likely substructure)

### Factor Graph Neuralization

Node interactions have also been modeled with Probabilistic Graphical Models
(PGM). A common graphical representation of PGM that is able to captures
higher-order constraints of nodes is the Factor Graph, which include factor
nodes that represents a higher-order relationship between a subset of graph
nodes. There is rich background of theoretically-grounded PGM algorithms to
infer a state of a node, given its statistical relation with the factors.
Knowledge of these inference algorithm produce good inductive bias that can
lead to better generalization of neural networks, if incorporated correctly
(Battaglia et al. 2018).

A successful attempt at Factor Graph Neuralization comes from (Dupty and Lee
2020), where the factors are approximated using tensor decomposition, and in
turn can be abstracted into a neural network module that can be learned end-
to-end within a bigger architecture. We will cover more details in Preliminary
section. Inspired by their idea, we will design a better network that enforce
the necessary special higher-order constraints in CSE task.

## Preliminary

### Basic Chemistry

We will briefly introduce about two important chemistry concepts here. Mass
spectrum requires a soft higher-order relationship between each spectrum peak
to all of the molecule atoms, while valence constraint enforce a hard high-
order constraint between the sum of edges bonds for each atom. These two
higher-order relationships will be captured in our architecture.

  * •

Mass spectrum: The mass spectrum represents the distribution of ions by mass
in a molecule to reflect the relative abundance of detected ions (Figure 1).
The x-axis is the mass-to-charge scale (that represents a relationship between
the mass of a given ion and the number of elementary charges that it carries),
and the y-axis is the signal intensity scale of the ions. The mass spectrum
will contain peaks that represent fragment ions as well as the molecular ion
(Kitson, Larsen, and McEwen 1996). In Figure 1, the x-axis value of 186
represents the molecule mass of decanoic acid substructure C11H22O2 (atom mass
sum equal 186), but if we analyze it manually, it could be some other
fragments with the same total mass like C10H18O3. Therefore, it is important
for the network to capture the higher-order relationship between each mass
spectrum peak to all the atoms of the molecules.

  * •

Valence Constraint: From the structural view, a molecule consists of the
aggregation of different atoms held together by valence forces, hence the
structure must also obey the valence rules. For example, the valence of Carbon
is 4 and valence of H is 1, therefore 1 Carbon cannot connect with more than 4
Hs. So the sum of edge bonds connecting to each atom must be equal to its
valence, and enforcing this hard constraint is difficult . Inspired from error
correcting code, we will design a factor node enforcing valence constraint in
our network as explained later.

### Factor Graph Neuralization

#### Graphical Models and Loopy Belief Propagation (LBP)

Consider a factor graph $G=(V,E)$ where vertices $V={1,2,...,n}$ and
hyperedges $E\in 2^{V}$. Each hyperedge represent a factor node that connect
to a subset of vertices. Every vertex $v_{i}\in V$ is associated with discrete
random variable $x_{i}$, while every hyperedge $e\in E$ is associated with a
potential function (of a factor node $a$) $f_{a}\in F$. In this paper, we will
index vertices as integers ${0,1,2,...,i,j}$ and factor nodes as ${a,b,c,...}$

Loopy Belief Propagation (LBP) is the process of computing and passing
messages between factor nodes and vertices to compute approximate marginal
probability $p(x_{i})$ for $x_{i}\in X$. LBP initializes two kinds of
messages, factor-to-node $m_{a\rightarrow i}(x_{i})$ and node-to-factor
$m_{i\rightarrow a}(x_{i})$. The recursive update rule for the message is as
follows:

$m_{i\rightarrow a}(x_{i})=\prod_{c\in N(i)\setminus a}m_{c\rightarrow
i}(x_{i})$ (1)

$m_{a\rightarrow i}(x_{i})=\sum_{X_{a}\setminus{x_{i}}}f_{a}(X_{a})\prod_{j\in
N(a)\setminus{i}}m_{j\rightarrow a}(x_{j})$ (2)

whereas $N(i)$ is the set of neighbours of $i$. After sufficient number of
recursive iterations, the belief of variable is computed as :

$b_{i}(x_{i})=f_{i}(x_{i})\prod_{a\in N(i)}m_{a\rightarrow i}(x_{i})$ (3)

#### Low-rank LBP

In (Dupty and Lee 2020), each factor is approximated as a low-rank tensor. By
applying low-rank tensor decomposition, the message update rule becomes:

$m_{a\rightarrow i}=W_{a_{i}}[\bigodot_{j\in
N(a)\setminus{i}}(W_{a_{j}}^{T}m_{j\rightarrow a})]$ (4)

$m_{i\rightarrow a}=\bigodot_{c\in N(i)\setminus{a}}m_{c\rightarrow i}$ (5)

where $\bigodot$ denotes component-wise multiplication.

With this derivation, the computational complexity of messages grow linearly
with addition of number of variables to factors. Therefore, this inference
algorithm is efficient enough to be used in an end-to-end neural network.

#### Neuralizing low-rank LBP

Figure 2: (a) Factor graph with $x_{1}$ and its neighbours. (b) Message
passing update for node $x_{1}$(Dupty and Lee 2020)

With the efficient derivations in (4) and (5), we can derive a neuralization
that incorporates structural inductive bias to the neural network using the
factor graph inference learning (Dupty and Lee 2020). The LBP algorithm
neuralization in (4) and (5) contains chained multiplication of numerous terms
leading to numerical instability. Empirically, Dupty and Lee (2020) found that
replacing the multiplication operation in the message update function with a
sum operation, followed by a MLP works well in learning useful
representations. Let graph $G=(V,E)$ with nodes and edges, and $F$ be the
factor graph defined on top of $G$. Let $h_{x_{i}}^{t}$ be hidden state of
node $x_{i}$ at iteration t. Then, the network will have the following message
passing update:

$h_{x_{i}}^{t+1}=h_{x_{i}}^{t}+MLP(\sum_{a\in N(i)}W_{a_{i}}[\bigodot_{j\in
N(a)\setminus{i}}W_{a_{j}}^{T}h_{x_{j}}^{t}])$ (6)

Here we are using node embedding $h_{x_{i}}$ as approximation for message
embedding, which is simpler and works well empirically (Dupty and Lee 2020).
The component of this message update neuralization is depicted in Figure 3

## Molecule Factor Graph Network

Figure 3: FGMN Architecture

In this section, we will describe our proposed model architecture together
with the pitfalls and design considerations. Following (Dupty and Lee 2020)
(our Preliminary Section), we will design a Graph Neural Network with
different types of factors. There will be 3 types of graph vertices here: atom
nodes, edge nodes, mass spectrum nodes. We first recap the problem
formulation: Given k mass spectrum nodes representing the peaks
$z_{1},z_{2},\dots,z_{k}$, predict the adjacency matrix $\mathcal{A}$
representing the molecule structures, where vertices in the matrix are list of
$n$ atom nodes $x_{1},x_{2},\dots,x_{n}$, and entries in the matrix are list
of $m$ edge nodes $y_{1},y_{2},\dots,y_{m}$.

With these three types of graph nodes, we need to define the factors that
capture the higher-order relationship among them. In Preliminary section
(Basic Chemistry), we mentioned two most important constraint here, first is
that a mass spectrum peak indicates the presence of a certain substructure,
and the second is the valence constraint of sum of edges bonds connecting to
the same atom.

### Graph Factors

  * •

MSP-Atom factor (Type C): The $k_{th}$ factor C $f_{C}^{k}$ will connect with
mass spectrum peak node $z_{k}$ and all the $n$ atoms of this current molecule
$x_{1},x_{2},\dots,x_{n}$. This factor node design is based on the fact that
each peak index represents the total mass of a substructure fragment.
Therefore, when connecting a mass spectrum node to all the molecule atom nodes
(atom nodes embed the atom mass in its value as well), more information will
flow from this mass spectrum peak to the atoms in corresponding substructures
during message update, which aids the learning and inference.

  * •

Edge-Atom factor (Type B): Each edge node $y_{l}$ and the two atoms linked by
the edge are connected to a factor type B $f_{B}^{i,j}$. This factor is
responsible for the relationship of each edge and its two atoms.

  * •

Edge-Edge factor (Type A): Each atom will holds a valence constraint (e.g
Carbon valence is 4, Oxygen valence is 2). For each atom $x_{i}$, we find all
$t$ edges $y_{1},y_{2},\dots,y_{t}$ which have atom $x_{i}$ at one of two ends
of the edge, then factor type A $f_{A}^{i}$ will connect to all these $t$
edges, and enforce the sum of these edges value to equal the valence value of
the atom $x_{i}$.

### Graph Factors Implementation

  * •

Type B and Type C: The two factor types here connect to a low number of nodes.
In this work, we will approximate these factors as low-rank tensors and adapt
derivations from Equation (6).

  * •

Type A: This is a hard-constraint factor node that is enforcing the sum of the
nodes connected to a certain value. We will derive as following in the
following section.

#### Valence Factor Type A Formula

We will derive the message from factor to each variable through a Dynamic
Programming formula. The idea is inspired by Parity Check in Error-Correcting
code formula derivation (MacKay 1999), in which each of its factor tries to
enforce the XOR operations of all nodes connected to zero.

Let $\mu^{valence}_{m}$ be the message from a valence factor type A to a
specific atom variable with valence constraint equal to $valence$, using only
the first $m$ variables connected to the factor ($m\in{1,2,\dots}$. Let
$g_{m}(x)$ be the message from variable $m$ to the factor, with value of
variable as $x$. Denote $b$ as maximum value of $x$. Then the message from
factor to variable is equivalent to the following derivation (Note that
$g_{i}(x)$ is the message from variable $i$ to this current factor):

$\begin{split}\mu^{v}_{m}&=\sum_{{x_{i}}:\sum_{i=1}^{m}x_{i}=v}\prod_{i=1}^{m}(g_{i}(x_{i}))\\\
&=g_{m}(0)*\sum_{{x_{i}}:\sum_{i=1}^{m-1}x_{i}=v}\prod_{i=1}^{m-1}(g_{i}(x_{i}))\\\
&\qquad+g_{m}(1)*\sum_{{x_{i}}:\sum_{i=1}^{m-1}x_{i}=v-1}\prod_{i=1}^{m-1}(g_{i}(x_{i}))\\\
&\qquad+\dots\\\
&\qquad+g_{m}(b)*\sum_{{x_{i}}:\sum_{i=1}^{m-1}x_{i}=v-b}\prod_{i=1}^{m-1}(g_{i}(x_{i}))\\\
&=g_{m}(0)*\mu^{v}_{m-1}+g_{m}(1)*\mu^{v-1}_{m-1}+...+g_{m}(b)*\mu^{v-b}_{m-1}\end{split}$
(7)

Therefore, from a brute force formula that is hard-to-implement and
computationally expensive $\mathcal{O}(max\\_bond^{max\\_atoms})$, we have
reduced to a Dynamic Programming formula efficient in performance and
implementation with time complexity
$\mathcal{O}(max\\_bond\\_type*max\\_atoms*max\\_value)$

$\begin{split}\mu^{v}_{m}=&g_{m}(0)*\mu^{v}_{m-1}+g_{m}(1)*\mu^{v-1}_{m-1}\\\
&+\dots+g_{m}(b)*\mu^{v-b}_{m-1}\end{split}$ (8)

We also observe that this Dynamic Programming (DP) formula will contain a long
chain of multiplication of the terms leading to numerical instability due to
overflow and underflow issues. Because normalization does not alter final
beliefs in Loopy Belief Propagation (Pearl 2014), we are normalizing the
hidden state of $g$ at each DP iteration.

Table 1: Results performance of each methods Model | Train | Test
---|---|---
|  loss |  acc |  loss |  acc
MASE with Encoder (non-positional encoding) |  0 |  1 |  0.403 |  0.8823
MASE with Encoder (positional encoding) |  0 |  1 |  0.53 |  0.8901
MASE with Encoder & Decoder |  1.212 |  0.787 |  0.86 |  0.804
MMPNN |  0.278 |  0.908 |  0.713 |  0.841
FGMN-Low |  0.166 |  0.932 |  0.481 |  0.877
FGMN-Medium |  0.059 |  0.989 |  0.275 |  0.914
FGMN-High |  0.202 |  0.928 |  0.563 |  0.871

### Factor weights conditional sharing

The factor weights are shared by conditioning on the nodes connected to the
factor. We will define amount of shareable parameters as low, medium and high
for each factor type B and type C (type A has no learnable parameter):

#### Factor Type B Parameter Sharing

  * •

Low: Regardless of nodes connected, set of factor weights remains the same.

  * •

Medium: Factor weight set is conditioned on the atom type of two atoms
connected to the factor, therefore number of different factor weight sets is
$\text{max\\_num\\_atom\\_types}^{2}$,

  * •

High: Factor weight set is conditioned on the atom index in the molecule,
therefore number of different factor weight sets is
$\text{max\\_num\\_atom\\_per\\_molecule}^{2}$

Table 2: Valence Satisfaction after Valence Factor Decoding Algorithm (Factor Type A), performed on various noise levels for ground truth Beta Noise Level | Initial | Sum Factor | Multiply Factor
---|---|---|---
0.1 | (Acc:1.000,valence:0.995) | (Acc:1.000,valence:0.997) | (Acc:1.000,valence:1.000)
0.2 | (Acc:0.991,valence:0.881) | (Acc:0.981,valence:0.768) | (Acc:0.987,valence:0.847)
0.5 | (Acc:0.811,valence:0.073) | (Acc:0.889,valence:0.267) | (Acc:0.893,valence:0.318)
1 | (Acc:0.565,valence:0) | (Acc:0.743,valence:0.128) | (Acc:0.775,valence:0.168)
2 | (Acc:0.398,valence:0.001) | (Acc:0.580,valence:0.027) | (Acc:0.650,valence:0.090)

#### Factor Type C Parameter Sharing

  * •

Low: Regardless of nodes connected, set of factor weights remains the same.

  * •

Medium: We will do a $K$-means clustering on all the peaks we have from
training dataset, then select top $k$ peaks. Then for each factor, the set of
parameter is chosen based on the closest center from the clustering to the
current mass spectrum node connected. Therefore, number of sharable factor
weight set is $K$.

  * •

High: Dependent on the mass spectrum node value, which has around 1000
distinct values in our case.

### Symmetry Problem

It is not easy to distinguish between first and subsequent atoms of the same
type in the same molecule. First carbon connecting to first oxygen in the list
of molecule atoms, or second carbon connecting to second oxygen in the list
might actually represent the same molecule structure. Therefore, given the
same inputs configurations, the training dataset might have different label
structure while in fact just permutation matrix of each other. So, we need to
devise good ordering strategy to approximate a canonical representation for
isomorphic label structure matrices for during training in order to make
learning easier.

We adopt the following ordering ways for the input data and label data. For
the input data of edges, we use the atom mass based ordering, e.g. with
lighter atom C at the front and heavier O at the back. Take the molecule
Methyl decanoate C11H22O2 for instance. The atoms list will be [C, C, C,
C,C,C,C,C,C,C,C,O,O] and the input edges will be constructed by each pair of
atoms. For the label data of edges, we make use of the SMILES ordering since
during training we know the molecular SMILES representation, which contains
the information of molecular graph. The SMILES representation of this molecule
is CCCCCCCCCC(=O)OC. The atoms list is [C,C,C,C,C,C,C,C,C,C,O,O,C]. If we
index the position, the list now is
[C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,O11,O12,C13]. The structure matrix is also
ordered based on this atoms list with the $i$th row and column being $i$th
atom in the atoms list. We can find that for the neighboring atoms in the
atoms list, there is always a link. This is also how SMILES is constructed.
Since the input data of edges is in the order of C first then O. We change the
atoms list of SMILES odering from [C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,O11,O12,C13]
to [C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,O1,O2]. Even though the SMILES ordering
and mass based ordering are two ordering, they are very similar in our
experiment due to the reason that we currently focus on atom C, H, and O. The
number of C is usually more than O, so usually the previous several atoms will
not change.

With this ordering of the edges in the training label structure matrix ground
truth, the network can be guided to learn in a more systematic learning
direction. We also tried Singular-Value Decomposition of the label structure
matrix and rank the atoms based on corresponding eigenvalues as an alternative
ordering for SMILES but did not get good result, so we will omit the details
on it here.

### Other details

  * •

Due to the large number of hydrogens per molecule which makes the whole graph
very computational and memory expensive, we will omit them in the SMILES
ordering and the graph design as well. Instead, we will add one fake hydrogen
atom for each molecule. For example, if the first Carbon is connected to 3
different hydrogen atoms, we will replace it with the Fake Hydrogen Atom
connecting to this Carbon with a bond value of three.

  * •

For each atom nodes and mass spectrum nodes, there are multiple messages
coming in so we combine the messages through a summation followed by MLP
layer, an approximation similar to (Dupty and Lee 2020). For each edge node,
number of messages are limited and we can use a multiplication of messages as
the original LBP formula.

  * •

We add the previous hidden state with the new message as residual connection
to improve network learning.

### Overall Architecture

The architecture of the module is depicted in Figure 3. Note that initial
state of the nodes are initialized through a Message Passing Neural Network as
well.

## Experiments

### Experiment Settings

The dataset is from the NIST database that contains mass spectrum files and
molecular structure information. In our experiment, we focus on the ester
functional group with only C H and O atoms. Since we know the molecule is in
ester group, we know the substructure of C(=O)O, rendering roughly 1800
molecules in our dataset. So we will extract this substructure and put C(=O)O
in the first three position in atoms list.

### Results

\- Evaluation: For single output, the loss may be 0 and 1 for a classification
task, with accuracy or recall as evaluation. Since we have multiple outputs
for a structure matrix instead of one output, we will use the accuracy to
evaluate performance of a single molecule, and average the accuracy of all the
molecules to evaluate the model. For the triangle matrix of 13 atoms, the
average accuracy of all-0 output is around 0.84. For loss function we are
using is cross-entropy, and we will find minimum loss among the permutation of
predicted output matrix.

We use Molecule Attention Structure Elucidation (MASE) as first baseline, an
adoption of Molecule Attention Transformer (MAT) (Maziarka et al. 2020) in
which we also include the Mass Spectrum Node and Edge Node in. MASE/MAT has no
recurrence or convolution operator, so we need to embed the relative or
absolute position of the atoms inside nodes as well. We will vary MASE with
options of having positional encoding and decoder or not. Besides, we also use
Molecule- Message-Passing-Neural-Network (MMPNN) (Gilmer et al. 2017) as
baselines to compare with our FGMN-Low, Medium, High Network, as illustrated
in Table 1.

Note: FGMN-Low means that the conditional parameter sharing for type B and
type C is low, likewise for FGMN-Medium and FGMN-High. Also, SMILES ordering
is used instead of Singular-Value-Decomposition to order the atoms since it is
observed to boost better performance.

The result is depicted in Table 1, and our FGMN-Medium is the best performance
model, where accuracy beats the second-best MASE model by a margin of 2.6%

### Ablation studies/ Analysis

#### Valence Factor Type A Effect

We will analyze the effect of our valence factor type A. We first design a
synthetic experimental settings: We will add random exponential distribution
noise (represented by beta value) to the ground truth one-hot vector of each
edge node, then try if the inference algorithm can decode the result. It is
analogous to experiments to show to what extent of noise can Error-Correcting
code Inference algorithm can reconstruct the original values.

Table 2 depicts the noise level and the percentage of atoms that have valence
constraint perfectly satisfied. Each atom node will have messages coming from
multiple factor type A, so Sum Factor refers to messages being summed up while
Multiply Factor refers to messages being multiplied. We can see that our
decoding algorithm is eliminating the noise to some extent and boosting the
performance.

By removing the valence factor type A node from our FGMN-Medium model, the
accuracy and percentage of atoms satisfying valence constraints also reduce
from 0.914 and 0.387 to 0.892 and 0.295 respectively.

#### Factor Type B and C effect

By training MPNN directly with Valence Factor without factor type B and C, we
obtain a testing accuracy as 0.862, with (1.7%, 5.7%, 1.0%) performance margin
reduction for FGMN-Low, Medium, High models respectively. Therefore, these
factors are important since through iterations of message update, it will flow
the corresponding information from each mass spectrum peak to the atom nodes
of certain substructures, and subsequently flowing to the edges nodes
corresponding to these atoms.

#### Parameter Sharing Conditioning

As explained in (Dupty and Lee 2020), designing good parameter sharing among
factors is important. If set of shareable parameter weights is too low, model
will lose approximation power to fit the training dataset. Otherwise, if set
of shareable parameter weights is too high, we may overfit the training set.

In our experiment (Table 1), FGMN-Medium beats FGMN-Low and FGMN-High by a
margin of 4.2% and 4.9% respectively.

## Conclusion & Future Work

In this paper, we have proposed a new architecture that uses various high-
order relationships and results in a good network for learning the molecule
structure, beating other baselines. Our model also removes the need for using
third-party software to get candidate substructures to rank and combine,
unlike existing works.

There are various aspects that can be explored to improve this work. For
example, we know that permutation of the result structure matrix is a problem,
so we can explore ways to embed and learn this knowledge during the training
process as well. We might also explore using Reinforcement Learning to
sequentially label the edges, further improving the enforcement of the valence
constraint.

## References

  * Battaglia et al. (2018) Battaglia, P. W.; Hamrick, J. B.; Bapst, V.; Sanchez-Gonzalez, A.; Zambaldi, V.; Malinowski, M.; Tacchetti, A.; Raposo, D.; Santoro, A.; Faulkner, R.; et al. 2018.  Relational inductive biases, deep learning, and graph networks.  _arXiv preprint arXiv:1806.01261_ .
  * Dupty and Lee (2020) Dupty, M. H.; and Lee, W. S. 2020.  Neuralizing Efficient Higher-order Belief Propagation.
  * Gilmer et al. (2017) Gilmer, J.; Schoenholz, S. S.; Riley, P. F.; Vinyals, O.; and Dahl, G. E. 2017.  Neural message passing for quantum chemistry.  _arXiv preprint arXiv:1704.01212_ .
  * Kerber et al. (2001) Kerber, A.; Laue, R.; Meringer, M.; and Varmuza, K. 2001.  MOLGEN-MS: Evaluation of low resolution electron impact mass spectra with MS classification and exhaustive structure generation.  _Adv Mass Spectrom_ 15(939-940): 22.
  * Kitson, Larsen, and McEwen (1996) Kitson, F. G.; Larsen, B. S.; and McEwen, C. N. 1996.  _Gas chromatography and mass spectrometry: a practical guide_.  Academic Press.
  * Krizhevsky, Sutskever, and Hinton (2012) Krizhevsky, A.; Sutskever, I.; and Hinton, G. E. 2012.  Imagenet classification with deep convolutional neural networks.  In _Advances in neural information processing systems_ , 1097–1105.
  * Lim et al. (2018) Lim, J.; Wong, J.; Wong, M. X.; Tan, L. H. E.; Chieu, H. L.; Choo, D.; and Neo, N. K. N. 2018.  Chemical structure elucidation from mass spectrometry by matching substructures.  _arXiv preprint arXiv:1811.07886_ .
  * MacKay (1999) MacKay, D. J. C. 1999.  Good Error Correcting Codes based on Very Sparse Matrices.  _IEEE Transactions on Information Theory_ 45(2): 399–431.  ISSN 0018-9448.
  * Maziarka et al. (2020) Maziarka, Ł.; Danel, T.; Mucha, S.; Rataj, K.; Tabor, J.; and Jastrzębski, S. 2020.  Molecule Attention Transformer.  _arXiv preprint arXiv:2002.08264_ .
  * Pearl (2014) Pearl, J. 2014.  _Probabilistic reasoning in intelligent systems: networks of plausible inference_.  Elsevier.
  * Schymanski et al. (2008) Schymanski, E. L.; Meinert, C.; Meringer, M.; and Brack, W. 2008.  The use of MS classifiers and structure generation to assist in the identification of unknowns in effect-directed analysis.  _Analytica chimica acta_ 615(2): 136–147.
  * Schymanski, Meringer, and Brack (2011) Schymanski, E. L.; Meringer, M.; and Brack, W. 2011.  Automated strategies to identify compounds on the basis of GC/EI-MS and calculated properties.  _Analytical Chemistry_ 83(3): 903–912.
  * Varmuza and Werther (1996) Varmuza, K.; and Werther, W. 1996.  Mass spectral classifiers for supporting systematic structure elucidation.  _Journal of Chemical Information and Computer Sciences_ 36(2): 323–333.
  * Vaswani et al. (2017) Vaswani, A.; Shazeer, N.; Parmar, N.; Uszkoreit, J.; Jones, L.; Gomez, A. N.; Kaiser, Ł.; and Polosukhin, I. 2017.  Attention is all you need.  In _Advances in neural information processing systems_ , 5998–6008.


# Linear programming word problems formulation using EnsembleCRF NER labeler
and T5 text generator with data augmentations

JiangLong He, Mamatha N, Shiv Vignesh, Deepak Kumar, Akshay Uppal
Infrrd AI Lab
Infrrd.ai
<EMAIL_ADDRESS>

###### Abstract

We propose an ensemble approach to predict the labels in linear programming
word problems. The entity identification and the meaning representation are
two types of tasks to be solved in the NL4Opt competition. We propose the
ensembleCRF method to identify the named entities for the first task. We found
that single models didn’t improve for the given task in our analysis. A set of
prediction models predict the entities. The generated results are combined to
form a consensus result in the ensembleCRF method. We present an ensemble text
generator to produce the representation sentences for the second task. We
thought of dividing the problem into multiple small tasks due to the overflow
in the output. A single model generates different representations based on the
prompt. All the generated text is combined to form an ensemble and produce a
mathematical meaning of a linear programming problem.

##  1 Introduction

Linear programming problems are used in business applications to optimize the
solution found for the problem. We define the linear programming problem as a
maximization or minimization problem. We may have to obtain maximum profit or
reduce operational costs based on the requirements. People conceived the idea
of formulating a mathematical problem to run a business and drive society to
the next level. The problem arises when there are multiple ingredients or
resources which scale only in a linear manner. We like to blend two or more
materials to maximize profit or minimize cost, an example of linearity in the
materialistic view.

The formulation of linear programming problems is a powerful concept in
solving a set of business problems consisting of monetary benefits. The
solution for the linear programming problems is still a critical concept to
understand by non-experts. Dantzig developed a method known as a simplex
method to solve linear programming problems Nash (2000). Here a set of
constraints defines the feasible region to optimize the linear programming
problem. We locate the corners of the feasible region and evaluate them in
this method to find the optimal solution. The solutions to the linear
programming problems got better with interior point methods. In an interior
point-based approach, the valuation begins from inside the feasible region,
which reduces the computational time compared to the simplex method when the
number of variables is large Karmarkar (1984). Now, we have readily available
software to solve the problem CPLEX (2022). The software helps to solve
mathematically framed linear programming problems. It still requires an expert
to formulate a linear programming problem from the word problem. The NL4Opt
competition NL4Opt (2022) creates an opportunity to develop an automated
solution to convert a word problem into a mathematically framed problem which
is then passed on to readily available software to find the optimal solution
Ramamonjison et al. (2022). The competition conceptualizes the hurdle faced by
the non-experts while formulating a linear programming problem. There are two
tasks in this competition, and they are named entity recognition and
generating the precise meaning representation. A method has to identify all
the entities present in the word problems in the entity recognition task.
Another solution uses the located entities to generate a mathematical
representation from the previous method. However, the text generation method
depends on the ground truth of the first task and acts as a reference for
evaluation.

Our contributions and observations are summarized as follows:

  * •

The number of samples in the training set is less to train a large language
model. So, we introduce the data augmentation strategies: LwTR, SR, MR, and
SiS Dai and Adel (2020) as a preprocessing step. Thus, we increase the number
of samples in the data to provide a more plausible output.

  * •

We introduce an ensemble of models to learn the weights for each label
predicted from the single models.

  * •

We introduce multi-task learning and train a large language model to generate
text through different prompts. An entity wrapper is placed around an entity
to enhance the data. We combine the task-specific text from the same model to
form a meaningful representation for linear programming solvers.

##  2 Sub-task 1: named entity recognition

Our approach follows the flair framework Akbik et al. (2019) and uses PyTorch
and hugging face Face (2022) transformers to build and experiment. Our base
model is a ‘Text Embedding + BiLSTM + CRF’ transformer model as shown in
Figure 1. We had used the base model to train a multilingual transformer for
MultiCONER competition Malmasi et al. (2022a, b). The text embedding layer was
XLM-Roberta embeddings for MultiCONER competition He et al. (2022). The F1
score of the model is equivalent to the baseline results when the text
embedding layer is Roberta-base Liu et al. (2019). Since the number of
training samples available to train the model is less, we performed data
augmentation techniques on the training samples.

Figure 1: The main architecture of the model where the text embedding layer is
replaced with different types of embeddings available in hugging face.

###  2.1 Data augmentation

A lot of labeled data is required by transformer-based language models to
produce a good performance. Here, we are working on linear programming word
problems, and the availability of the word problems is scarce. Apart from
this, we have a time-consuming annotation task that requires expert knowledge.
In this case, both of them are very critical. We need to overcome this data
deficiency problem. So, we made use of data augmentation techniques to replace
token-level embeddings Dai and Adel (2020). The method expands the number of
training samples by adding additional transformed samples without changing the
labels. There are four techniques in this method: Label-wise token replacement
(LwTR), Synonym replacement (SR), Mention replacement (MR), and Shuffle within
segments (SiS). Here, we have used all the techniques to perform data
augmentation in the training samples. However, the development set was not
utilized and kept as it was for evaluation purposes. Table 1 illustrates an
example of generation of augmented data.

Table 1: The different types of data augmentation techniques: Label-wise token replacement (LwTR), Synonym replacement (SR), Mention replacement (MR), and Shuffle within segments (SiS). Method | Instance |  |  |  |  |  |  |
---|---|---|---|---|---|---|---|---
Original | A | serving | of | chicken |  | costs | $ | 10
| O | O | O | B-VAR |  | B-OBJ_NAME | O | B-PARAM
LwTR | A | flour | of | cereal |  | profit | $ | 25
| O | O | O | B-VAR |  | B-OBJ_NAME | O | B-PARAM
SR | A | serving | of | volaille |  | cost | $ | X
| O | O | O | B-VAR |  | B-OBJ_NAME | O | B-PARAM
MR | A | serving | of | Durian | TV | savings | $ | 0.32
| O | O | O | B-VAR | I-VAR | B-OBJ_NAME | O | B-PARAM
SiS | serving | of | a | chicken |  | costs | $ | 25
| O | O | O | B-VAR |  | B-OBJ_NAME | O | B-PARAM

####  2.1.1 Single model accuracy after data augmentation

We started evaluating the models with different types of embeddings and a
combination of embeddings. The results are tabulated in Table 2. The XLM-
Roberta-base model Conneau et al. (2019) was the initial choice as
experimented with in the baseline results. Quickly, we moved to the next set
of embeddings due to an increase in the F1 score for different types of
embeddings. The Roberta-base embeddings provided better results compared to
the XLM-Roberta-base. We performed ablation study by changing the text
embeddings and their combination. The transformer model has complex contextual
information between the tokens. Whether the contextual information from the
transformer model is sufficient or not is evaluated in the single model by
removing a few of the components in the pipeline. In the seventh and eighth
rows, nobilstm and nobilstmcrf refer to without the BiLSTM module and
BiLSTM+CRF module from the model. The F1 score was better without the BiLSTM
module giving a glimpse of contextual information gained from the transformer
rather than the separate BiLSTM module. However, we couldn’t distinguish a
single model that is better than other.

Table 2: The performance of a single model with different text embeddings on the development set. Model | Text embeddings | Development
---|---|---
number |  | F1 score
1 | roberta-base + no augmentation | 0.8117
2 | roberta-base + LwTR & SR augmentation | 0.8641
3 | xlm-roberta-base + 4 augmentations | 0.8898
4 | roberta-base + 4 augmentations | 0.9127
5 | roberta-base + glove + 4 augmentations | 0.9151
6 | only glove + 4 augmentations | 0.787
7 | roberta-base + glove + nobilstm + 4 augmentations | 0.9154
8 | roberta-base + glove + nobilstmcrf + 4 augmentations | 0.9034

###  2.2 Ensemble architecture

A simple ensemble strategy of Majority Voting is developed by He et al.
(2022). Given a set of M sequence labeling models denoted as
$C=\\{c_{1},c_{2},...,c_{M}\\}$ and an input sentence denoted as
$S=\\{w_{1},w_{2},...,w_{n}\\}$, where each $w$ is a word from S. Each model
from $C$ will output a sequence of predictions for each word $w$ in sentence
$S$. Let
$O^{s}_{c_{i}}=\\{O^{c_{i}}_{w_{1}},O^{c_{i}}_{w_{2}},...,O^{c_{i}}_{w_{n}}\\}$
denote the prediction output of model $C_{i}$ on sentence $S$.
$O^{c_{i}}_{w_{j}}$ denotes the prediction of model $C_{i}$ on word $w_{j}$ in
IOB format. The set of outputs for all models in $C$ on sentence $S$ will be
denoted as

$O_{S}=\\{O^{s}_{c_{1}},O^{s}_{c_{2}},...,O^{s}_{c_{M}}\\}$ (1)

The Majority Voting strategy takes all model’s predictions of word $w_{j}$ and
outputs the most frequent prediction as to the final prediction for $w_{j}$.
An obvious issue with Majority Voting is the IOB scheme constraint. The final
ensemble result is not guaranteed to be passing all the constraints where (I)
tag must follow and (B) tag and the entity of neighbor (B) and (I) tag must be
the same. We introduce an ensemble learning approach via sequence labeling
called ‘EnsembleCRF’ as shown in Figure 2.

Figure 2: The architecture of EnsembleCRF model. Given a sentence as input,
each of the sequence labeling models will output the name entity prediction in
IOB format. One hot encoder combines them to generate ensemble output from a
Conditional Random Field model.

The model outputs are stacked together and passed through a one-hot encoder,
three linear layers, and CRF. The CRF layer is trained to optimally combine
the model predictions to form a new set of predictions. The addition of the
three linear layers helped in the performance improvement. The EnsembleCRF
model is of the form

$\begin{array}[]{c}C_{en}=EnsembleCRF(C=\\{c_{1},c_{2},...,c_{M}\\},\\\
D_{en}=\\{X,Y\\})\end{array}$ (2)

$D_{en}$ is the ensemble learning dataset composed of $X$ and $Y$.
$X=\\{O_{s_{1}},O_{s_{2}},...,O_{s_{k}}\\}$ is created by using model set $C$
and set of input sentences $\\{S\\}=\\{S_{1},S_{2},...,S_{k}\\}$ with size K.
Each element of $X$ is defined as equation (1). $Y$ is the ground truth entity
label in IOB format. During the training phase for some models in set $C$, we
included only the training datasets. Thus, we choose to perform the ensemble
learning using an augmented train set created using the data augmentation
strategies explained in Section 2.1. Since the second layer classifier is a
CRF layer, we solved the problem of breaking the IOB constraints. By learning
to optimally combine the model predictions, EnsembleCRF also learns to avoid
mistakes made by single sequence labeling models.

We experimented with creating $D_{en}$ with the augmented training dataset.
However, we found that there is no positive correlation between the number of
models in $C$ and the macro-averaged F1 score on the development dataset.
Treating every possible combination of set $C$ as a hyperparameter to optimize
will yield the optimal result.

####  2.2.1 Ensemble model accuracy

In the MultiCONER competition Malmasi et al. (2022a, b), an ensemble model
provided a boost in the F1 score. We followed a similar approach and used the
same architecture. The single models, which were evaluated in the earlier
section, are combined to form an ensemble. The results are tabulated in Table
3. We found all six models with 4 augmentations gave an F1 score of 0.9173 on
the development set, but the best model evaluated on the test set had a
different combination, which is represented, in bold format. The results
varied at the third decimal digit position indicating that the combination is
providing improvement but is not sufficient to bring in the change in the f1
score. However, the results are consistent with the ensemble model. One of the
main concerns with the ensemble model is the number of single models to use to
gain the maximum F1 score. Since it has to be evaluated in a trial-and-error
approach.

Table 3: The performance of ensemble with different text embeddings on the development set. Model | Text embeddings | Development
---|---|---
number |  | F1 score
1 | roberta-base + 4 augmentations | 0.9154
(3 models) | roberta-base + glove + 4 augmentations |
| xlm-roberta-base + 4 augmentations |
2 | roberta-base + 4 augmentations | 0.9164
(3 models) | roberta-base + glove + 4 augmentations |
| roberta-base + glove + nobilstm + 4 augmentations |
3 | roberta-base + 4 augmentations | 0.9161
(4 models) | roberta-base + glove + 4 augmentations |
| roberta-base + glove + nobilstm + 4 augmentations |
| roberta-base + glove + nobilstmcrf + 4 augmentations |
4 | roberta-base + 4 augmentations | 0.9167
(5 models) | roberta-base + glove + 4 augmentations |
| roberta-base + glove + nobilstm + 4 augmentations |
| roberta-base + glove + nobilstmcrf + 4 augmentations |
| xlm-roberta-base + 4 augmentations |
5 | roberta-base + 4 augmentations | 0.9173
(6 models) | roberta-base + glove + 4 augmentations |
| roberta-base + glove + nobilstm + 4 augmentations |
| roberta-base + glove + nobilstmcrf + 4 augmentations |
| xlm-roberta-base + 4 augmentations |
| glove + 4 augmentations |

###  2.3 Results of sub-task 1

The test results are tabulated in Table 4. Our ensemble approach is placed at
the top with an F1 score of 0.939. The results are tabulated after performing
a reproducibility test by the competition organizers which helps us confirm
that our programs are reproducible and consistent with the F1 score.

Table 4: The performance of different submissions on the test set NL4Opt (2022). Rank | Team Name | Affiliation(s) | F1 score
---|---|---|---
1 | Infrrd AI Lab | Infrrd | 0.939
2 | mcmc | OPD | 0.933
3 | PingAn-zhiniao | PingAn Technology | 0.932
4 | Long | BDAA-BASE | 0.931
5 | VTCC-NLP | Viettel | 0.929
6 | Sjang | POSTECH | 0.927
7 | DeepBlueAI | DeepBlueAI | 0.921
8 | TeamFid | Fidelity | 0.920
9 | KKKKKi | Netease | 0.917
10 | holajoa | Imperial College London | 0.910
11 | Dream |  | 0.884
| Baseline (xlm-roberta-base) | Nl4Opt | 0.906

##  3 Sub-task 2: text generation

The generation of declarations was conceptualized as entity relationship
mapping. The relationship between the ground truth entities was missing, and
our refined objective was to map the relationship to capture the declaration
that exists in the natural language. To do so, we tried to implement the
entity relationship way of mapping the entities. We couldn’t progress further
and moved to the next step, the text generation approach. We observed that the
text generation process was much easier than entity relationship mapping in
our experiments.

We use a text generator to generate the declarations. We have used text to
text transfer transformer (T5) in our experiments Raffel et al. (2019). The
architecture has visible tokens for the input text, and the output text is
observed based on the past predicted outputs. We fine-tuned the regularly
available T5 model from the hugging face library Face (2022). In the first
stage, we provided the raw input text as input to the transformer, and the
expected output was as defined in sub-task 2 NL4Opt (2022). Our output is a
generation of declaration, and the results generated were on par with the
baseline results shared by the competition organizers.

Figure 3: The generation of text using T5 model for different types of
prompts.

We observed that the score was par but didn’t have much-anticipated
improvements. We performed an initial analysis and found the two points
identified during our discussion. They were the number of samples available to
train the model and the repeated occurrence of declarations in terms of
constraints. To address these points, we decompose each sample into pieces for
multi-task learning and, thus, implicitly increase the number of training
samples in our dataset.

We believe the output was low due to the small number of samples. We split the
single task into multiple tasks as shown in Figure 3. In this way, we
invariably introduced a lot of examples in the training set. This approach
helped the model to map the declaration with the probable entities in the
process. First, we trained the model to predict only the objective of the
problem, and it was then scaled to generate the constraint declarations. There
were seven different types of constraints in the competition. They are sum
constraints, upper bound constraints, lower bound constraints, linear
constraints, ratio constraints, xby constraints, and xy constraints. Thus, we
combine all the declarations by repeating the training seven more times.

Every linear programming problem has a definite objective statement, but the
constraints declarations may be missing in a few. In our multi-task learning
setting, the missing constraint statements in a problem act as negative
samples in the training process. It helps the model in predicting the correct
constraints without much effort. Tables 5 and 7 show the original mapping and
the proposed mapping for sub-task 2 to add multiple samples to the training
dataset.

Table 5: The original mapping for text generation in sub-task 2 NL4Opt (2022).
Original mapping:
---
`<s>`
`<DECLARATION>`
`<OBJ_DIR> maximize </OBJ_DIR>`
`<OBJ_NAME> number of coconuts </OBJ_NAME> [is]`
`<VAR> rickshaws </VAR> [TIMES] <PARAM> 50 </PARAM>`
`<VAR> ox carts </VAR> [TIMES] <PARAM> 30 </PARAM>`
`</DECLARATION>`
`<DECLARATION>`
`<CONST_DIR> at most </CONST_DIR>`
`<OPERATOR> LESS_OR_EQUAL </OPERATOR>`
`<LIMIT> 200 </LIMIT>`
`<CONST_TYPE> [LINEAR_CONSTRAINT] </CONST_TYPE> [is]`
`<VAR> rickshaws </VAR> [TIMES] <PARAM> 10 </PARAM>`
`<VAR> ox carts </VAR> [TIMES] <PARAM> 8 </PARAM>`
`</DECLARATION>`
`<DECLARATION>`
`<CONST_DIR> must not exceed </CONST_DIR>`
`<OPERATOR> LESS_OR_EQUAL </OPERATOR>`
`<CONST_TYPE> [XY_CONSTRAINT] </CONST_TYPE>`
`<VAR> ox carts </VAR> [is] <VAR> rickshaws </VAR>`
`</DECLARATION>`
`</s>`
Table 6: The exploratory data analysis of number of tokens for different training strategies. Type of training | Prefix + input length (max) | Output mapping length (max)
---|---|---
Original: prefix=generate | Train: 203 | Train: 716
linear program mapping | Dev: 214 | Dev: 601
augmented input | Train: 545 | Train: 716
| Dev: 529 | Dev: 601
multi-task | Train: 248 | Train: 564
| Dev: 258 | Dev: 438
augmented input + multi-task | Train: 591 | Train: 564
| Dev: 570 | Dev: 438
Table 7: The proposed multi-task mapping for text generation. Multi-task
mapping:
---
`prompt <OBJ_DECLARATION> </OBJ_DECLARATION>:`
`<OBJ_DECLARATION>`
`<OBJ_DIR> maximize </OBJ_DIR>`
`<OBJ_NAME> number of coconuts </OBJ_NAME> [is]`
`<VAR> rickshaws </VAR> [TIMES] <PARAM> 50 </PARAM>`
`<VAR> ox carts </VAR> [TIMES] <PARAM> 30 </PARAM>`
`</OBJ_DECLARATION>`
`prompt <CONST_DECLARATION> [LINEAR_CONSTRAINT] </CONST_DECLARATION>:`
`<CONST_DECLARATION>`
`<CONST_DIR> at most </CONST_DIR>`
`<OPERATOR> LESS_OR_EQUAL </OPERATOR>`
`<LIMIT> 200 </LIMIT>`
`<CONST_TYPE> [LINEAR_CONSTRAINT] </CONST_TYPE> [is]`
`<VAR> rickshaws </VAR> [TIMES] <PARAM> 10 </PARAM>`
`<VAR> ox carts </VAR> [TIMES] <PARAM> 8 </PARAM>`
`</CONST_DECLARATION>`
`prompt <CONST_DECLARATION> [XY_CONSTRAINT] </CONST_DECLARATION>:`
`<DECLARATION>`
`<CONST_DIR> must not exceed </CONST_DIR>`
`<OPERATOR> LESS_OR_EQUAL </OPERATOR>`
`<CONST_TYPE> [XY_CONSTRAINT] </CONST_TYPE>`
`<VAR> ox carts </VAR> [is] <VAR> rickshaws </VAR>`
`</DECLARATION>`
negative samples
`prompt <CONST_DECLARATION> [SUM_CONSTRAINT] </CONST_DECLARATION>:`
`prompt <CONST_DECLARATION> [XY_CONSTRAINT] </CONST_DECLARATION>:`
`prompt <CONST_DECLARATION> [RATIO_CONSTRAINT] </CONST_DECLARATION>:`
`prompt <CONST_DECLARATION> [UPPER_BOUND] </CONST_DECLARATION>:`
`prompt <CONST_DECLARATION> [LOWER_BOUND] </CONST_DECLARATION>:`

We fine-tune the model through prompting. We prefix the input text to generate
the appropriate output depending on the prefix task. We performed exploratory
data analysis to compute the maximum number of tokens after prefix at the
input and the output. The numbers are tabulated in Table 6 and exceeded the
limit of 512 tokens in the original training set. After splitting the main
task into multiple tasks and increasing the data, we observed that the results
started moving up and away from the baseline results. The results were not as
expected after this step, and another round of analysis was carried out to
identify the reasons for not picking up the correct declarations from the
statements.

Figure 4: An example word problem with the wrapped entities used in the training. Table 8: The declaration-level accuracy on the development set for different types of training with the hyperparameters. Type of training | Hyperparameters | Declaration-level accuracy
---|---|---
Original: prefix=generate | max_seq_length: 600 | 62.8%
linear program mapping | max_output_length: 750 |
| training_epochs: 20 |
augmented input | max_seq_length: 600 | 75.9%
| max_output_length: 750 |
| training_epochs: 24 |
multi-task | max_seq_length: 600 | 65.6%
| max_output_length: 750 |
| training_epochs: 17 |
augmented input + multi-task | max_seq_length: 600 | 83.5%
| max_output_length: 750 |
| training_epochs: 50 |

We thought of teaching the model by explicitly mentioning the actual locations
of the entities in the input text. Thus, we prepared a wrapper for all the
entities with a notation shown in Figure 4. The entities present in the input
text get wrapped before training the model. The number of tokens doubled with
the addition of the wrapper tokens. Table 6 shows the increase in the number
of tokens processed by the model. We observed the maximum number of tokens in
the training set before and after the addition of the wrapper are 240 and 590,
respectively. With the addition of new tokens, the sequence length crossed the
actual limit of 512, as suggested in the paper Raffel et al. (2019). We have
to fix the size of the transformer to 600 and perform training on the model.
Since the size is not a power of 2, the performance deteriorates to some
extent. Table 8 shows the ablation study for different training strategies.

###  3.1 Results of sub-task 2

The results of sub-task 2 on the test set is tabulated in Table 9. We stand in
the fifth position after evaluating the program independently by the
competition organizers and thus results are updated by reproducing the model
by the organizers. There is a slight drop in the accuracy due to the retrieval
of multiple instances for the same prompt question.

Table 9: The performance of different submissions on the test set NL4Opt (2022). Rank | Team Name | Affiliation(s) | F1 score
---|---|---|---
1 | UIUC-NLP | UIUC | 0.899
2 | Sjang | POSTECH | 0.878
3 | Long | BDAA-BASE | 0.867
4 | PingAn-zhiniao | PingAn Technology | 0.866
5 | Infrrd AI Lab | Infrrd | 0.780
6 | KKKKKi | Netease | 0.634
| Baseline (BART) | NL4Opt | 0.608

###  3.2 Discussion

The additional wrapper tokens not only increase the length of the input
sequence and redundant for many repeated items. We have to replace each
wrapper token with a unique token to limit the size of tokens doesn’t cross
the prescribed limit of 512. We have performed data augmentation in sub-task
1. We may need to use this approach to increase the number of samples in the
training set. The objective and constraints declarations consist of a lot of
numerical values. It is difficult to identify whether to perform the
replacement of a numerical value or word reflecting the entity in a problem.
The retrieval of the model through prompting has one serious concern to be
addressed. When the entity appears at multiple locations, the model may pick
one of them and predict it as an output. If we have different values at
several locations for the same entity, the model picks one and doesn’t
retrieve all of them. The values present in the input text differ for the same
entity. It has to be addressed in the future to make the model adapt to
picking up all of them belonging to the same entity in the training procedure.

##  4 Conclusion

We proposed an ensemble of models to predict the labels and generate the
texts. The approaches are consistent with better results compared to a single
model. We would like to explore the decoders in the ensemble since the text
generation had the decoder part. We would like to use the large models in
place of the base models to know whether the models are better at predicting
the labels. The OBJ NAME entity is difficult to capture in the present set of
models which has to be addressed in a different way. We need to overcome the
drawbacks of prompting like two instances in the same question to improve the
score.

## References

  * Akbik et al. [2019] Alan Akbik, Tanja Bergmann, Duncan Blythe, Kashif Rasul, Stefan Schweter, and Roland Vollgraf.  Flair: An easy-to-use framework for state-of-the-art nlp.  In _NAACL 2019, 2019 Annual Conference of the North American Chapter of the Association for Computational Linguistics (Demonstrations)_ , pages 54–59, 2019.
  * Conneau et al. [2019] Alexis Conneau, Kartikay Khandelwal, Naman Goyal, Vishrav Chaudhary, Guillaume Wenzek, Francisco Guzmán, Edouard Grave, Myle Ott, Luke Zettlemoyer, and Veselin Stoyanov.  Unsupervised cross-lingual representation learning at scale.  _arXiv preprint arXiv:1911.02116_ , 2019.
  * CPLEX [2022] CPLEX.  IBM ILOG CPLEX Optimization Studio, 2022.  URL www.ibm.com/products/ilog-cplex-optimization-studio.
  * Dai and Adel [2020] Xiang Dai and Heike Adel.  An analysis of simple data augmentation for named entity recognition.  _arXiv preprint arXiv:2010.11683_ , 2020.
  * Face [2022] Hugging Face.  Hugging face.  https://huggingface.co/docs/transformers/index, 2022.  Accessed: 2022-12-23.
  * He et al. [2022] Jianglong He, Akshay Uppal, Mamatha N, Shiv Vignesh, Deepak Kumar, and Aditya Kumar Sarda.  Infrrd.ai at semeval-2022 task 11: A system for named entity recognition using data augmentation, transformer-based sequence labeling model, and ensemblecrf.  In _Proceedings of the 16th International Workshop on Semantic Evaluation (SemEval-2022)_ , pages 1501–1510, Seattle, United States, July 2022\. Association for Computational Linguistics.  doi: 10.18653/v1/2022.semeval-1.206.  URL https://aclanthology.org/2022.semeval-1.206.
  * Karmarkar [1984] N Karmarkar.  A new polynomial-time algorithm for linear programming.  _Combinatorica_ , 4:373–395, 1984.  doi: https://doi.org/10.1007/BF02579150.
  * Liu et al. [2019] Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov.  Roberta: A robustly optimized bert pretraining approach.  _arXiv preprint arXiv:1907.11692_ , 2019.
  * Malmasi et al. [2022a] Shervin Malmasi, Anjie Fang, Besnik Fetahu, Sudipta Kar, and Oleg Rokhlenko.  MultiCoNER: a Large-scale Multilingual dataset for Complex Named Entity Recognition.  2022a.
  * Malmasi et al. [2022b] Shervin Malmasi, Anjie Fang, Besnik Fetahu, Sudipta Kar, and Oleg Rokhlenko.  SemEval-2022 Task 11: Multilingual Complex Named Entity Recognition (MultiCoNER).  In _Proceedings of the 16th International Workshop on Semantic Evaluation (SemEval-2022)_. Association for Computational Linguistics, 2022b.
  * Nash [2000] John C Nash.  The (dantzig) simplex method for linear programming.  _Computing in Science & Engineering_, 2(1):29–31, 2000.
  * NL4Opt [2022] NL4Opt, 2022.  URL https://nl4opt.github.io/.  Accessed: 2022-12-23.
  * Raffel et al. [2019] Colin Raffel, Noam Shazeer, Adam Roberts, Katherine Lee, Sharan Narang, Michael Matena, Yanqi Zhou, Wei Li, and Peter J. Liu.  Exploring the limits of transfer learning with a unified text-to-text transformer, 2019.  URL https://arxiv.org/abs/1910.10683.
  * Ramamonjison et al. [2022] Rindranirina Ramamonjison, Haley Li, Timothy T. Yu, Shiqi He, Vishnu Rengan, Amin Banitalebi-Dehkordi, Zirui Zhou, and Yong Zhang.  Augmenting operations research with auto-formulation of optimization models from problem descriptions, 2022.  URL https://arxiv.org/abs/2209.15565.