rescaled such that the maximum singular value of the weight matrix, whfined as
W = (w1 w2 · · · wMx ) ,“Forecasting the geomagnetic secular variation by machine learning”
8485
8687
8889

becomes 0.99. This rescaling is applied to satisfy the so-called “echo sta
to theEarthESNand SpaceScience
which guaranteesmanuscriptthat thesubmitted
state ofis notaffected by distant pa
The output of the ESN at time tk , y k , is then obtained from xk as follo

H. Toh, S. Sato [Division of Earth & Planetary Sciences, Kyoto University] &S. Nakano [Department of Statistical Modelling, The Institute of Statistical Mathematics; Center for Data Assimilation Research and Applications, ROIS-DS]
136
90
137

91
138
139
92
140
141
142

94

143
144

Results145

The objective function of this study:
95
146
147

ESN was found so effective for short-term forecast of GSV that it didn't need long training data.

67

Data covariance does affect fits of the machine learning method to the observation.
148

Echo State Network is a variant of reservoir computing applicable to machine learning for time-series prediction.

We conduct hindcast experiments to reproduce the temporal evolution of thewhere Γ denotes
Thetheoutputy k corresponds
geomagnetic
main fieldthe
afterweighttraining matrix.
the ESN using
COV-OBS.x2
model. The to a pre
COV-OBS.x2
model gives
the temporal
evolution of the Gauss coefficients of thethe observation
at timetk .
scalar potential V as an ensemble of field models sampled from a continuous stochastic model.
Here we obtain
the Gauss coefficients
temporalDenoting
the observation
at timefortkevery
as yeardk , andthe thematrixΓ is determimm
evolution of ∆gn (tk ) and ∆hn (tk ) for each year is modelled with the ESN. Whenmmmizing
the following
objective function:∆gn (tk ) and ∆hn (tk ) are obtained as the temporal difference for a 1-year interval, their typical scale is of the order of 10 nT. The
for a suf&inputs& ESN requires
Earth and Space K
Science&&2 the#Γ#2!
ficient number of timemanuscript
steps submitted
beforeto its
outputcan beobservadkcompared− ΓT xkwith
F2J ∆g= nm (tk ) and ∆hm,tions. Hence, we use the observations of(tk ) from +1840 to 1860
σnk21861 to 2005. We
λ2 then2 Method
for spin-up
and train the ESN using the observations
fromk=1
Following
many modelsof the)Earth’smagneticthe IGRF
model, the hindcast of g m (t )predict ∆gnm(tk ) and∆hm
2006field,to including
2015 andobtaink fromkn (tnwe represent the magnetic field
B with a scalar potential V asandhmk ) accordingly.
n (tthewheresecond term on the right-hand side of this equation is a regula

ESNによるΔg10とΔg21の出力
Discussion I:
93

Echo State Network (ESN) Method

y k = Γ T xk ,3 Hindcast experiments

149
96

97
150
151
98

Table 1. Hindcast errors of ESN as a function of training data length and time resolution of the COV-OBS geomagnetic field model.

152

99

153

100
154
155

101

156

B = −∇V.

(1)

termToto
avoid
and#Γ#Frobenius
norm
determineq i using(12),
theparameters
σk,i theand λmust be given
in of the mThe
potential
V isoverfitting
expanded
into Eq.
spherical
harmonics:
F denotes# λ correspondadvance.
The parameterto the uncertaintyof theof
observationdk,ik,i" acorresponds
The parameters
σ! σ!
andto the scales
uncertainties
for thV (r, θ, φ, t) = a k[g (t) cos mφ + h (t) sin mφ] P (cos θ)(2)and it was determined by ther variance for each coefficient gnm (tk ) and hm
n (tk ) given used intions and constraints,respectively.
Thedeterminedvalues ofbytheparameters$by the COV-OBS.x2
model.
The Theparameter
λ was field
the maximizawhere a denotes the Earth’s
mean radius.
SV of the geomagnetic
is represented as the firstin
time thederivativesof the Gauss
coefficients gDecomposing(t) and h (t).
are ofdescribed
nextsection.
dk beandΓ as dtion
themarginal
likelihood
the minimizationof Ji canregarded
ask = dk,1"#because
Nnn+1

mnmn

mnn=1 m=068

mn

69mn

We model their temporal variations by the ESN model. The state of the sysa Bayesianestimationdk,iand λi canberegarded
as the variance of theat time
represented
by state vector
xk . The numberof statevariables
Mx can be rewritten as:and Γtem
=γtk1inis,this. .study.
. ,problem
γAtM,qofk,respectively,Eq.
(6)
y stepis setto 1000distribution
time
Gaussianpriorfori . the i-th element of xk , xk,i , is updated
7071
7273

as follows:$+ ξ tanh wT x%k,ik−1,ik−1
k
i'i K
(
$among the toT be %2005
First, we set thestart timeof thei M
hindcast
experiments
2 and the2yfor connecting
where z k denotes the input vector, wi is a weight vector!
!variables,d2005.
γ i parameter
xk
#γ
158
ESN wastrained
the vector
COV-OBS.x2
model
λ wask,istatevariables, uwith
for connecting with
the inputuntiland −ξ The
i is a weighti #2102
J the= marginal+ used.
is the leakage rate (Jaeger et al., 2007; Lukoševičius,2012). We fixed the value of ξ159
determinedby the maximization oflikelihood, which
is often
2at 0.5 in this study. The weights wi and ui are given in advance and are fixed. We σ
λ2 in
90% of the weights(e.g.,{wi } and{ui } (randomlychosen)
to zero.
Thevalues ofbecause
the ki=1k=1160
Bayesianset estimation
Morris,1983;
Casella,1985),
the minimization
remaining non-zero elements of ui are drawn randomly from a normal distribution
161
of Ji canwithbemean
regarded
as deviationa Bayesianestimation
with a Gaussian prior for
0 and standard
σu . The standard
deviation problem
σu is set to adjust
the rangeof the input
z as described
The valuesforof theele162
q i . Figure
1 shows
thevariables
resultsof the later.
hindcast
g10non-zero,} gare11 ,thenh11 , g20 , g21 , h12 , g22 , h22 , andments of wi are also drawn from a normal distribution. The weights {wi163
g30 . In eachblue line
indicates
results
ofwhichtheishindcast
conducted withrescaledpanel,such thatthe
the maximum
singularvalue of the weightmatrix,defined as164
the ESN,the red line indicates
the COV-OBS.x2, and the–3–
greenline indicates theW = (w1 w2 · · · wM ) ,(4)165
predictionof 0.99.
the This10th-generation
IGRF
model
(Maus et al., 2005). The
becomes
rescaling is applied to satisfy
the (IGRF-10)
so-called “echo state
property”
which guarantees
thatistheregarded
state of the ESN
not affectedby distant
inputs.
166
COV-OBS.x2
model
asistheactualSV. pastSince
the IGRF-10 was released inThe output of the ESN at time tk , y k , is then obtained from xk as follows:
167
2005, the prediction by the IGRF-10T is also shown as a benchmark of the predictiony =Γ x ,(5)168from 2005. The ESN provided a kbetterk predictions for g10 , g11 , g20 , g21 , h22 . Althoughwhere Γ denotes the weight matrix. The output y k corresponds to a prediction of169
the prediction
with the ESN deviated from the actual SV as described by the COVthe observation at time tk .
170
OBS.x2 forDenoting
g22 andg30 of which
the trends changed around 2005, the performance ofthe observation
at time tk as dk , the matrix Γ is determined by minifollowing objectiveto
function:
171
the ESNmizing
wasthecomparable
or better than that of the IGRF-10 even for those co&2&K &dT!
172
efficients. This result suggests
thatxk &2 ESN#Γ#2F has significant potential for predictingk − Γthe+,(6)The objective function, J, of this study to determine the matrix G containing
Wout is given by J = k=1
sk is theso-called
λ2 , where173
SV with satisfactory accuracy.
Weσk2also examined
the performance for the case when‘data covariance’. One interest is how the data covariance affects theresult.
The
figure
above
shows
terms
of
fits totwo specificwherethe secondtermon the right-handof it
thisequation
is a regularization
174
the start
timeof thehindcast
wassidesettoin
be2010.
Figure2 shows
the result of theterm to avoid overfitting and #Γ#F denotes the Frobenius norm of the matrix Γ.
175hindcast
initiated2010
with tothe
Thetomodel and(thethe 11thTheparameters
σincorrespondthe scales
uncertainties
forCOV-OBS.x2
theobserva- COV-OBS
k and λconstanttrained differential Gauss coefficients, Dg10 and Dg21, from 1860 through2005.
Althoughsblue
k‘sof line.
seem
followtions
andconstraints,respectively.
Thevaluesoftheparameters
usedin this study
176
generation IGRF (IGRF-11) model (Finlay et al.,are
$ 2010)% also indicated with thedescribedin the next
section. Decomposing dk and Γ as dk = dk,1 , . . . , dk,M"# Because
blue solid line) much better, it does not necessarily mean that the variablek‘sareareworse.
older
data
arelikely
tohavelarger
177
redsandgreenlines,respectively,forreference.
AlthoughtheIGRF-11
model welland Γ = γ 1 , . . . , γ M , respectively, Eq. (6) can be rewritten as:
178
predicted
theSVfrom2010,
theESNpredicted
theSV
betterthanIGRF-11'($%
observation errors, it is rather desired to minimize the effect of older data with large errors. ESNseems
indifferent
to data olderthe
thatnM2K!
!
dk,i − γ T#γ #2i xk
179
model.
J=+ i2 2 .
(7)
274

157

= (1 − ξ)xx

75

+ uT z + η

(3)

76

manuscript submitted to Earth and Space
77
78
79
80
81
82
83
84
85
86

x

87
88
89

Table 1 clearly shows that ESN doesn't require longer training data for prediction of short-term GSVs of ourplanet's main magnetic field. Rather, the shorter training data of 105 yrs yielded better accuracy. The dependence
202of time resolution, i.e., the total number of data, is weak enough to give the best estimates by a combination ofshorter training data and half a year resolution. The improvement in estimation error sums up to 20 to 30 %. 203

9091years, while higher coefficients change increasingly over shortthat the prediction errors by ESN!are dependent on harmonshows the temporal evolution of !dPn,k for each n from 200
figuretheprediction
dP, are smaller for hig
1960,
whereas it thatneatly follows
the differential
COV-OBS modelerrors,
for the latest years (1995andn,konwards).
errors themselves increase as time passes from the starting y
Discussion II:
929394

95

96
97
98

Two hindcast results

20499

100
101

i=1180
181

206

182

Data Showcase(Old Magnetic Storms)

Conclusions

l We conducted the hindcasts of GSV for Epochs 2005 and 2010 using ESN, one of the simplest
reservoir computing methods in machine learning. In doing so, we adopted the COV-OBS (Gilletet al., 2013) geomagnetic field model since 1840 as ‘observation’.
l ESN turned out so effective that it was able to reproduce the short-term GSV without longer
training data with a typical 10-y mean prediction error of less than 100 nT.
l It was found that the data covariance may have a large impact in forecasting GSVs.
l It, however, may not be possible to forecast occurrence of geomagnetic jerks even by themachine learning method used here. This needs further elaboration such as inclusion of Length ofDay in input data.

k=1

σk

λWe also evaluate the accuracy of the prediction with the root-mean-square–3–(RMS) differences (Whaler & Beggan,2017; Minami et al., 2020):"# n&'#%!
('( )m (t ) − g m (t ) 2 + hm (t ) − hm (t ) 2 ,dPn,k = $(n + 1) gn,pkn,o kn,p kn,o k

The degree-dependent forecast error:

The so-called geomagnetic secular variation (GSV) is, at a glance, dominated by smooth and gradual changes
in the vector geomagnetic field. However, it still contains highly non-linear temporal variations such as thegeomagnetic jerk characterized by abrupt changes in the first time derivatives of, especially, the eastward
geomagnetic components. Presence of those non-linear variations makes the prediction of our planet's
magnetic field difficult enough to force the geomagnetic research communities revise their global geomagnetic
field models unexpectedly urgently (e.g., Witze, 2019).
Conventional prediction of GSV has been made empirically taking advantage of the linearity of its principal
element. However, this doesn't work especially when the aforementioned geomagnetic jerks occur. To
circumvent this, the number of GSV predictions using a combination of geodynamo simulation and data
assimilation is increasing recently (e.g., Minami et al., 2020). The combination turned out effective fornon-linear predictions of GSV except for its computer-intensive nature in both geodynamo simulation and dataassimilation. Namely, it is not until working with an enough number of ensemble members (each of which is a
result of a fairly long geodynamo run) that data assimilation for GSV forecast becomes effective.
On the other hand, machine learning is an emerging technique, which has a promising potential ofapplication to highly non-linear processes such as the geomagnetic jerks. For example, Koopman Mode
Decomposition is a variant of machine learning methods based on a decomposition of observed time-series into
multiple modes, while the key procedure of Recurrent Neural Network is feedback of various internal states tothe present status. In this study, we adopted ESN (e.g., Nakano and Kataoka, 2022) to perform hindcasts ofshort-term GSVs from 2005 and 2010 using a geomagnetic field model (COV-OBS; Gillet et al., 2013) stemmingfrom 1840 because ESN is known as cost-effective reservoir computing and works even with relatively few
training data.

yyy205102m=0"# N#%

dP = $
dP ,Smaller forecast errors for higher-degree Gauss coefficients!

183184The panel above compares the hindcast results for Epoch 2005 (Left) and 2010 (Right).
Although overall forecasts are fairly good for both epochs, ESN seems to have failed to follow
abrupt changes in the internal Gauss coefficients (e.g., g22 at around 2005). The hindcast,
therefore, does not necessarily guarantee the flexibility of ESN to reproduce non-linear temporalvariations such as the ‘geomagnetic jerks’.

Acknowledgements
This study is supported byROIS-DS-JOINT Program #003RP2023 of
Research Organization of Information and207Systems, Japan.
208

REFERENCES

Figure 5.

!

dPkn,k(15)

(16)

n=1

It’s quite–5– 2
10 natural thatthe forecasterrors
increase withtime.
101 However, they
do possess
degreedependence
because the100 higher degree
Gausscoefficientsvary more
rapidly, which
10-1 does appear
even in shortteaching data.

of each degree n for the predn,kü Alken, P., Thébault, E., Beggan, C.D. et al. International Geomagnetic
Reference Field: the thirteenth generation.
Earth Planets Space 73, 49 (2021). https://doi.org/10.1186/s40623-020-01288-x.
ü Gillet, N., Jault, D., Finlay, C. C., and Olsen, N. (2013), Stochastic modeling of the Earth's magnetic field:
Inversion for covariances over the observatory era, Geochem. Geophys. Geosyst., 14, 766- 786,
doi:10.1002/ggge.20041.
ü Minami, T., Nakano, S., Lesur, V. et al. A candidate secular variation model for IGRF-13 based on MHD dynamo
simulation and 4DEnVar data assimilation. Earth Planets Space 72, 136 (2020). https://doi.org/10.1186/s40623020-01253-8.
ü Takahashi, F. (2012) Implementation of a high-order combined compact difference scheme in problems ofthermally driven convection and dynamo in rotating spherical shells, Geophysical & Astrophysical Fluid
Dynamics, 106:3, 231-249, DOI: 10.1080/03091929.2011.565337
ü Witze A. Earth's magnetic field is acting up and geologists don't know why. Nature. 2019 Jan;565(7738):143-144.
doi: 10.1038/d41586-019-00007-1.

4 Discussion

To predict the evolution of a dynamical system for gen