Introduction

Reservoir computing1,2,3 is one of the successful derivatives of recurrent neural networks (RNN), which are composed of an input layer, a reservoir, and an output layer (Fig. 1a). The reservoir is a dynamical system which maps the input data onto a higher-dimensional information space nonlinearly to mimic the roles of hidden layers in a RNN. In reservoir computing, only the weight matrix \(\varvec{W}_{\textrm out}\) connecting the reservoir and the output nodes is required to be trained, instead of training all the weight matrices linking the layers in a RNN. In this way, the reservoir computing efficiently resolves the difficulties of RNNs regarding the time-consuming and potentially unstable trainings. Aside from purely dynamical reservoir models, many physical reservoirs have been proposed to date4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33. Among them, magnetic reservoirs20,21,22,23,24,25,26,27, including the systems involving skyrmions28,29,30,31,32,33,34, have several advantages such as nonvolatility35 (namely, the magnetic system should retain its initial texture after the removal of inputs to fulfill the reproducibility), durability36,37, low energy consumption as compared to CMOS architectures38,39,40, and quick responses21.

Figure 1
figure 1

(a) Architecture of the reservoir computing. (b) Magnetization configuration of the skyrmion reservoir and the layout of input (red circles) and readout (yellow circles) nodes. In-plane magnetization vectors \((m_{ix}, m_{iy})\) are presented by arrows at sites (\(i_x\), \(i_y\)) when \(\text {mod}(i_x, 2)=\text {mod}(i_y, 2)=0\), while out-of-plane components \(m_{iz}\) are presented by colors. Inset shows an enlarged view of one readout node.

Full size image

Magnetic skyrmions are nanometric topological spin textures, which form a hexagonal lattice under the application of a static magnetic field in chiral magnets with broken inversion symmetry41,42. Skyrmions are robust against environmental agitations36 because of their topological protection43,44. A magnetic skyrmion crystal exhibits specific spin-wave modes at microwave frequencies45,46,47. In this sense, a skyrmion crystal works as a set of series-connected spin-torque oscillators (STOs)20,21,22,23,24 and is expected to possess the characteristics required by a physical reservoir. In contrast to STOs, one advantage of the skyrmion reservoirs is that advanced nanofabrication and complicated manufacturing processes are not required in their production. In our previous work48, we have demonstrated that a skyrmion lattice in a thin-plate chiral magnet possesses the fundamental requirements of a reservoir, including the generalization ability, short-term memory, and nonlinearity, inherently carried by the spin wave dynamics in the skyrmion lattice excited by the locally applied magnetic-field pulses as inputs1,2,3.

Skyrmion-based reservoirs have been proposed in literature mostly by random skyrmion configurations both numerically28,29,30,31 and experimentally32. The major purpose of using a random texture is to enhance the nonlinearity of the physical dynamics to increase the potential of separating the linearly inseparable input classes in tasks such as pattern classifications. In this paper, alternatively, we propose to use a (slightly distorted) skyrmion lattice as a reservoir, without introduction of random pinning sites or anisotropies, and we demonstrate that by utilizing the nonlinear interferences of spin waves in such a skyrmion crystal, it suffices to generate high performances in the recognition task of handwritten digits extracted from the Modified National Institute of Standards and Technology (MNIST) database.

For two unrepeated sets of randomly chosen 6600 and 3300 digits as the training and testing datasets, respectively, with 3530 or even smaller amounts of optimized components in the weight matrix, a recognition rate of 88.2% for the testing dataset can be reached by the skyrmion lattice. This performance is higher than that of 79.3% by using a dynamical echo-state network model instead of the skyrmion crystal to transform the input data into the reservoir state49,50,51 and 50% by directly using the greyscale data of handwritten digits as state vectors without the operation of any reservoir32. All training algorithms for these three cases are the same, providing an unbiased comparison of their performances. Importantly, we reveal that this great performance is attributable to the highly nonlinear transformation of input data carried by spin-wave interferences in the skyrmion lattice, instead of merely linear transformations of the input into the amplitudes of spin wave dynamics.

Although the CMOS architecture has been maturely applied to practical handwritten digit recognitions, semiconductor devices are vulnerable to environmental stimuli and cause considerable energy consumption, therefore there is a strong urge to find spintronics alternatives39,40,52. Our result shows that the skyrmion lattice is a promising material for spintronics reservoirs in machine learning applications.

Results

Skyrmion spin-wave reservoir

Figure 1b shows the design of our skyrmion spin-wave reservoir. Sixteen input nodes (red circles labeled from 1 to 16) are installed near the top and bottom edges of a rectangular skyrmion lattice. Another sixteen nodes (yellow circles labeled from 1 to 16) are installed in between as readout nodes (detectors). We propose in experiments to install circular current loops underlying the nodes to apply and detect local fields via electromagnetic induction (as in reference32). Every node has a radius of 5 lattice constants and contains 80 sites inside. Taking the lattice constant as the length unit, the system size is \(128\times 64\). The center of the first input node is located at (8.5, 55.5), and the distances between centers of neighboring nodes are 16 in both of the x and y directions. The underlying skyrmion lattice is described by the classical Heisenberg model, and we solve the magnetization dynamics by numerically solving the Landau-Lifshitz-Gilbert (LLG) equation (see Methods).

Figure 2
figure 2

Schematics of the input procedure. (a) One example of the original 28 \(\times\) 28 greyscale image extracted from MNIST database. (b) This image is trimmed and rotated by \(\pi /2\) to become a 22 \(\times\) 16 image. The greyscale data in each of the 16 rows is simultaneously injected into the 16 input areas through magnetic-field pulses, respectively.

Full size image

Input injection and reservoir state detection

The procedure of preprocessing and injection of input data begins with the original 28 \(\times\) 28 pixel images randomly extracted from the MNIST database with one example as shown in Fig. 2a, each gray-scale pixel representing one of the integers in the range of [0, 255]. The image is reduced in size to 16 \(\times\) 22 by trimming, i,e, three, three, six, and six lines of mostly white pixels are removed from the top, bottom, left, and right edges, respectively. Then each of the trimmed images is rotated by \(\pi /2\) as shown in Fig. 2b, and the sixteen rows of greyscales are injected simultaneously into the sixteen input nodes, respectively. Each of the rotated trimmed images is described by a 16 \(\times\) 22 input matrix \(\varvec{S}_{\textrm in}\) whose components \(S_{{\textrm in},mp}\) are numbers representing the grayscales. In this study, we consider a system with multiple input nodes to reduce the computational cost for the numerical simulations by simultaneously injecting the sixteen sets of sequential pixel data into the skyrmion-lattice reservoir. We note that systems with a fewer input nodes or even with a single input node may also show good performance in the reservoir computing. Such systems might be advantageous in terms of the easy fabrication. Optimization of the number and locations of input and detector nodes to maximize the performance is left for future study.

The 22 components in each of the 16 rows of \(\varvec{S}_{\textrm in}\) are sequentially entered into the skyrmion reservoir via each of the 16 input nodes as time series of local magnetic-field pulses \(H_{\textrm in}(\varvec{r}_i,t)\hat{\varvec{z}}\). For the input \(S_{{\textrm in},mp}\) (\(1\le m \le 16, 1\le p\le 22\)), the magnetic-field pulse is applied to the sites i within an area of the mth input node for \((22-p)\Delta t \le t < (23-p)\Delta t\) to excite the magnetization dynamics. The magnitude of the pulse is set to be \(H_{\textrm in}(\varvec{r}_i,t)= 5\times 10^{-4}S_{{\textrm in},mp}\), while the duration of the pulse is fixed at \(\Delta t=2.5\). During and after the input injection with these field magnitude and duration, we find the skyrmion lattice retains its overall topological number; namely, the initial number of skyrmions as fifteen in Fig. 1b remains the same. Therefore, the magnetization excitations belong to the spin wave regime, and the spin waves eventually reach the detectors after experiencing complex reflections and interferences.

To construct the reservoir state of the skyrmion reservoir, during and after the input procedure for each of the kth digit, the out-of-plane magnetizations \(m_{iz}\) averaged over the sites within each of the mth detector area are measured and denoted by \(M_{mq}^{(k)}\). Specifically, we record \(m_{iz}\) data at instants with a constant interval of \(2\Delta t\) for sequentially 22 times (namely, \(1\le q \le 22\)) starting from the initial time when the first input is injected. Note that the input sequence has a total period of \(22\Delta t\) as described in the last paragraph. Therefore, half of the measured \(m_{iz}\) data lies in the input procedure, while the other half lies in a period of \(22\Delta t\) after the last input is fully injected. In such a way, we expect that the magnetization dynamics dominated by complex spin wave interferences post the input procedure can enhance the nonlinearity and hence the digit recognition rate. (In the end of this paper we design and perform a numerical simulation to consolidate this point.) We then define the reservoir state vector \(\varvec{\psi }^{(k)}\) for each of the kth input using all elements in the measured signals \(\varvec{M}^{(k)}\) and an additional constant bias element. The dimension of \(\varvec{\psi }^{(k)}\) is therefore 353(\(=16 \times 22+1\)) with 16 detectors each with 22 temporal nodes, plus one constant bias.

Comparison of skyrmion reservoir and echo state network

In Fig. 3, we compare the recognition rates for the testing set performed by the skyrmion reservoir and the echo state network (ESN) as a famous dynamical reservoir model, plotted as functions of the number of iterations in the training procedure at intervals of 100 iterations (see Methods). Although the recognition rate by ESN increases faster than that by the skyrmion reservoir, it stops increasing in very initial few hundreds of iterations, and eventually the skyrmion reservoir reaches a recognition rate of 88.18%, higher by about 10% than that by ESN (79.3%). Note that the dimensions of \(\varvec{W}_{\textrm out}\) for both cases are the same, 3530 (= \(10\times 353\)), for a fair comparison of both methods. For a more detailed visualization, in Fig. 4a,b the confusion matrices for both methods are shown, in which the x and y axes represent the predicted and true digits, respectively, and the colormap shows the number of predicted digits in each class. It is apparent that the skyrmion reservoir leads to a better result compared with ESN.

To test the ability of the skyrmion reservoir, we have reduced the temporal sampling of \(m_{iz}\) from 22 to 10 for each detectors, leading to a number of 1610(\(=10\times (16\times 10+1)\)) weight matrix elements. The recognition rate in this case only slightly drops to \(87\%\), still higher than that by ESN with 3530 weight matrix elements. This achievement is comparable to a recent work33 in which the curren-driven motion of a single skyrmion was utilized for the digit recognition. With 1970 weight parameters, the authors obtained a 87.6% recognition rate. These facts, together with our previous work48 showing that spin waves in skyrmion lattices are endowed with the input-estimation ability, short-term memory, and nonlinearity, comprehensively reveal that the skyrmion lattice is a suitable candidate for reservoir computing applications.

Figure 3
figure 3

Recognition rate during the training procedure performed by the skyrmion reservoir (red curve) and the echo state network (blue curve).

Full size image
Figure 4
figure 4

Confusion matrices for (a) skyrmion lattice and (b) echo state network. The diagonal represents the correctly recognized digits. The color intensity reflects the number of recognized digits for each digit class.

Full size image
Figure 5
figure 5

(a) The layout of input (red circles) and detector (yellow circles) nodes for the comparison of linearity and nonlinearity. (b) Schematics of the procedure of input processing and injection. (c) Time profiles of space-averaged out-of-plane magnetizations excited at input nodes 1–7. Profiles of input pulses (gray boxes) are also presented. (d) Same as (c) for magnetizations measured at detectors 1–7. (e) Three different arrangements of temporal-node sets examined in the comparison. The magnetizations are measured at these moments at each detectors.

Full size image
Figure 6
figure 6

Recognition rate of the testset during the training procedure performed for three temporal-node sets A–C in Fig. 5e.

Full size image

Comparison of linearity and nonlinearity

In general, the reservoir is required to contain nonlinearity encoded in its states with respect to input data, in order to separate the linearly inseparable classes (digits 0–9) of input data. However, the reservoir is often treated as a black box with complex dynamics in which the actual key factors that alter the content of nonlinearity are unclear. In the skyrmion lattice, we find it is possible to compare the performances by two kinds of reservoir states (the \(m_{iz}\) data) that are dominated by either a linear or a nonlinear transformation of the input greyscales signals.

For this purpose, we design a numerical test based on a simplified layout with seven input and seven detector nodes as shown in Fig. 5a. The center of the first input node is located at (16.5, 37.5), and distance between neighboring nodes is 16 (20) in the \(x\ (y)\) direction. To illustrate our idea with a reduced simulation time, we manually choose 400 (200) less distorted digits for the training (testing) dataset, and adopt the preprosseing and input procedure as shown in Fig. 5b. The digits are trimmed to 15 \(\times\) 21 pixels by removing four, three, seven, and six lines of pixels from the top, bottom, left, and right edges, respectively. Then we coarse grain them into 5 \(\times\) 7 pixels as the input matrices \(\varvec{S}^{(k)}_{\textrm in}\) by taking averages of the greyscales in each 3 \(\times\) 3 squares. Each \(\varvec{S}^{(k)}_{\textrm in}\) is injected into the skyrmion lattice via the seven input nodes as magnetic-field pulses \(H_{\textrm in}(\varvec{r}_i,t)\hat{\varvec{z}}\) with magnitude being \(S^{(k)}_{{\textrm in},mp} \times 10^{-4}\), while the duration of the pulse is fixed at \(\Delta t=10\).

Figure 5c shows time profiles of space-averaged out-of-plane magnetization \(M_p(t)\) at seven input nodes (p = 1–7) excited by the applied field pulses for an example digit image of “6”. The figure also shows sequences of the applied five field pulses with gray boxes. The width of box corresponds to the duration of one pulse. We find that the profiles of \(M_p(t)\) exhibit the same trend as the instantaneous pulse magnitudes. This is because \(M_p(t)\) is directly excited by the applied field pulses at each input node. On the contrary, as seen in Fig. 5d, \(M_q(t)\) at seven detectors (q = 1–7) exhibit delayed responses to sequences of the input pulses because the detectors are located away from the input nodes. For detectors 1 and 7, the trends are less correlated, possibly due to edge effects on the spin wave propagations.

Since \(M_q(t)\) roughly follows the trends of input streams with a small delay in the early time domain of \(0 \le t \le 5\Delta t\), it indicates that the magnetization dynamics in this early time domain depends almost linearly on the input signals, because they mostly reach the detectors directly without experiencing significant interferences. On the contrary, \(M_q(t)\) in the later time domain of \(5\Delta t \le t \le 10\Delta t\) is dominantly affected by nonlinear spin wave interferences since the input signals have already stopped. These observations motivate us to compare the performances dominated by linearity and nonlinearity inherent in the reservoir states.

We examine the recognition task for different temporal-node sets at which instants we measure the \(m_{iz}\) data for the reservoir state. The sets B and C shown in Fig. 5e contain ten temporal nodes in the early and later time domains, respectively, whose recognition results are plotted in Fig. 6 with blue and orange curves. We clearly find that the nonlinearity leads to a better recognition rate (orange curves) by about 10% than the linearity (blue curves). We note the recognition rate by nonlinearity is even slightly larger than that by temporal-node set A carried out by both linearity and nonlinearity (green curves) with the same number of nodes. This result clearly reveals the power of nonlinearity and, more importantly, shows that its content is potentially adjustable in the spin wave dynamics in skyrmion reservoirs by specialized designs of the input-detector setup, which could be a valuable guideline for future experimental realizations.

For practical applications to general pattern recognitions, we expect that a magnet with a smaller Gilbert damping constant can realize a better performance because the spin waves will have substantially larger amplitudes to enhance the nonlinearity encoded in the magnetization dynamics. Another aspect concerns randomness of the skyrmion lattice. A random skyrmion texture has been proposed as a suitable reservoir by virtue of the large content of nonlinearity28,30. It is generally suggested that when the system is located on the verge of the phase boundary between ordered and chaotic states (so-called edge of chaos), the highest nonlinearity manifests itself53,54. These scopes will be our future studies.

Comparison with other works

The existing works on the handwritten digit recognition task using skyrmion-based reservoirs include references32,33. In reference32, a recognition rate of 94.7% is reached by using a skyrmion network. Although their recognition rate is higher than ours, 88.2%, the authors take nine separated skyrmion configurations, each being applied by distinct static magnetic field to generate their own skyrmion textures, and a total amount of 15,840 weight matrix components are required to be optimized. In our case, we use a single skyrmion crystal which is expected to significantly reduce the device size and fabrication difficulty. Moreover, it is also expected to enhance the recognition rate by increasing number of parallel skyrmion crystals and/or weight matrix parameters. Additionally, the authors in reference32 measure the Hall voltage as the reservoir state, which requires an electrical current applied into the crystals, that will cause Joule heating and consume more energy than the spin-wave skyrmion reservoir. This is another advantage of our reservoir compared to other skyrmion-based reservoirs28,29,30,31,32,33 which use voltage or current as inputs. In order to increase the recognition rate, we note that one of the major contributions of our work is to demonstrate the origin of good recognition, which is the nonlinear spin wave interferences that dominate the magnetization dynamics after the input is fully applied (Fig. 6). Therefore, our work has provided a guide to designing better reservoirs that can achieve more accurate recognitions. For instance, larger input magnetic field pulses and/or smaller-Gilbert damping materials support the larger amplitudes and longer lifetime of spin waves, and thus we expect the recognition can be improved in these situations.

The purpose of using large numbers of input areas and detectors in Fig. 1b is to reduce the simulation time; namely, in such a way we can simultaneously inject the input fields corresponding to all the 16 rows of greyscales into the crystal, then measure the magnetization response also simultaneously from the 16 detector areas. In the experiment32, the authors flatten out the two-dimensional digit greyscales to a one-dimensional input sequence to inject only one magnetic field to each of their skyrmion crystals. To reduce the fabrication complexity of our layout, we expect that a smaller number of local input and detector areas can also lead to comparable results, since the major origin of high recognition comes from the nonlinear spin wave interferences excited by input fields, and the input/detector layout could be further designed to reduce the fabrication complexity while maintaining the nonlinearity capacity contained in the magnetization dynamics. This will be our scope of future works.

In another work on the handwritten digit recognition in the reference33, they use 196 separate quasi-one-dimensional chiral magnets, inside each a single skyrmion is driven by the input current, and their positions are recorded as the reservoir state. This complex layout containing 196 parallel skyrmions can obviously be improved by using our single skyrmion lattice. Moreover, after one computation, all of the skyrmions in each magnets may need to be driven by current to move back to their initial positions, if the device is to be used for later computations. In our case, since the input magnetic fields do not drive significant skyrmion motions such that the lattice can relax much more easily to the initial configuration for the next task, the spin-wave skyrmion reservoir can overcome these difficulties. As compared to other spintronics reservoirs based on memristors or spin-torque oscillators with typical microsecond time scales, the skyrmion devices can operate in a time scale of nanoseconds55, thus higher computational speed is expected.

Discussion

We have theoretically studied the performance of a skyrmion spin-wave reservoir on the handwritten digit recognition task. A high recognition rate of 88.2% is achievable by the skyrmion reservoir with 3530 parameters in the weight matrix to be trained for a subset of image data randomly extracted from the MNIST database. Importantly, skyrmions emerge spontaneously in magnetic materials with broken inversion symmetry under application of static magnetic field. Therefore, the skyrmion reservoir requires no advanced nano-fabrications for production, in contrast to other spintronics reservoirs using, e.g., spin-torque oscillators and magnetic tunnel junctions20,21,22,23,24,25,26. Recently, even a zero-field skyrmion lattice is possible to be stabilized experimentally56,57,58, raising up more possibilities for its spintronics applications. Our work paves the way to realizing high-performance spintronics reservoirs in machine-learning applications.

Methods

The skyrmion lattice

The underlying skyrmion lattice is described by the classical Heisenberg model on a square lattice in the xy plane with open boundary conditions. The Hamiltonian for local magnetizations contains the nearest-neighbor ferromagnetic exchange interactions, the Zeeman interactions, and the Dzyaloshinskii-Moriya interactions (DMI) as,

$$\begin{aligned} \mathscr {H}= & {} -J\sum _{\langle i,j\rangle }\varvec{m}_i \cdot \varvec{m}_j -\sum _i [\varvec{H}_{\textrm ext}+\varvec{H}_{\textrm in}(\varvec{r}_i,t)] \cdot \varvec{m}_i \nonumber \\{} & {} +D\sum _i (\varvec{m}_i\times \varvec{m}_{i+\hat{\varvec{x}}} \cdot \hat{\varvec{x}} + \varvec{m}_i\times \varvec{m}_{i+\hat{\varvec{y}}} \cdot \hat{\varvec{y}}), \end{aligned}$$
(1)

where \(\varvec{m}_i\) denotes the classical unit-length magnetization vector at site i. For the Zeeman term, we consider time-dependent local fields \(\varvec{H}_{\textrm in}(\varvec{r}_i,t)=H_{\textrm in}(\varvec{r}_i,t)\hat{\varvec{z}}\), in addition to a static global field \(\varvec{H}_{\textrm ext}=H_{\textrm ext}\hat{\varvec{z}}\) which stabilizes the lattice. The time-dependent fields are applied locally to the sites within the sixteen input-node areas to inject the sequential input data into the skyrmion reservoir. The position vector \(\varvec{r}_i=(i_x, i_y)\) of site i represents integer coordinates in units of the lattice constant. We take \(J=1\) as the energy unit and set \(D=0.36\) and \(H_{\textrm ext}=0.06\) to produce a stable skyrmion lattice48. When \(J=1\) meV, the dimensionless time \(t=1\) and field strength \(H=1\) correspond to 0.66 ps and 8.64 T, respectively.

Note that the values of \(D=0.36\) and \(H_{\textrm ext}=0.06\) are derived from numerical calculation results with \(D = 0.09\) in reference45 via a scale transformation \(D\rightarrow aD, H_{\textrm ext}\rightarrow a^2H_{\textrm ext}\), and \(t\rightarrow t/a^2\) with \(a = 4\), which reduces the system size that requires us to simulate. Performing an inversed scale transformation, we are simulating an actual chiral magnet with the skyrmion size of about 50 nm and the pulse width of input field in the setup of Fig. 1b is \(a^2\Delta t=16\times 2.5\times 0.66\text { ps}\approx 26\) ps that corresponds to a frequency of 38 GHz in the microwave range, by which the local breathing modes of skyrmions can be excited.

Magnetization dynamics induced by the applied \(\varvec{H}_{\textrm in}\) fields are numerically simulated using the fourth-order Runge-Kutta method to solve the Landau-Lifshitz-Gilbert (LLG) equation,

$$\begin{aligned} \frac{d\varvec{m}_i}{dt} =\frac{-1}{1+\alpha _{\textrm G}^2}\left[ \varvec{m}_i\times \varvec{H}^{\textrm eff}_i +\alpha _{\textrm G}\varvec{m}_i \times (\varvec{m}_i \times \varvec{H}^{\textrm eff}_i) \right] , \end{aligned}$$
(2)

where \(\alpha _{\textrm G}(=0.001)\) is the Gilbert-damping constant, and \(\varvec{H}^{\textrm eff}_i \equiv -\partial \mathscr {H}/\partial \varvec{m}_i\) is the effective local field. The initial magnetization configuration (Fig. 1b) is obtained by the Monte Carlo thermalization with simulated annealing to low temperatures followed by a sufficient relaxation executed with the LLG equation.

Training and testing procedures

For the training and testing datasets respectively, we randomly choose unrepeated 6600 and 3300 digits from the MNIST database with equal numbers of digits from 0 to 9 in both sets. Note that we choose a rather small ratio between the numbers of training and testing datasets, namely 2:1 (in literature this ratio is typically 5:1 or higher33), as a verification to test the learning ability of the skyrmion lattice. The procedure for the handwritten digit recognition is described below.

First, we adopt the one-hot representation for the targets of input digits. Specifically, the kth target, i.e., the expected output for the kth input, is the integer \(N^{(k)}\) corresponding to the digit, which is represented by a ten-dimensional vector \(\varvec{x}^{(k)}\) where \(x^{(k)}_\ell =\delta _{\ell ,N^{(k)}+1}\; (\ell =1,2,\ldots ,10)\). Here \(\delta _{a,b}\) is the Kronecker’s delta variable, and \(N^{(k)}(=0,1,\ldots ,9)\) is the correct number for the kth input image. Second, since one of the fundamental requirements of the reservoir is its nonlinear mapping of the input data into a higher-dimensional space, we expect all the nonlinear transformations are encoded in the magnetization dynamics in the skyrmion lattice. Therefore, for the output function, only a linear transformation of the reservoir state, \(y^{(k)}_\ell = \sum _j W_{\text {out},\ell j} \psi ^{(k)}_j\), is executed, with \(\varvec{W}_{\textrm out}\) being a 10 \(\times\) 353 output weight matrix. This is distinct from the use of a postprocessing nonlinear softmax output function (such as further defining the output as \(y'^{(k)}_\ell = \exp [y^{(k)}_\ell ]/\sum _{j=1}^{10} \exp [y^{(k)}_j]\)) usually done in RNNs.

Finally, in the training procedure, the components of \(\varvec{W}_{\textrm out}\) are optimized by the gradient descent method with the Adam algorithm59 to minimize the loss function taken as a mean-square error between the output and target data defined by,

$$\begin{aligned} L=\frac{1}{N}\sum _{k=1}^{N} \sum _{\ell =1}^{10} (x^{(k)}_\ell - y^{(k)}_\ell )^2. \end{aligned}$$
(3)

Here \(N=6600\) is the number of training inputs. In this training procedure, we reset the first and second moments to zero after every 100 iterations in the Adam algorithm defined in59 to achieve an empirically faster convergence of the loss function. After the training procedure, we calculate \(\varvec{y}^{(k)}\) and the recognition rate using the optimized \(\varvec{W}_{\textrm out}\) for both training and testing datasets. The digit of the kth image is recognized as the number \(\ell -1\) when the \(\ell\)th component of \(\varvec{y}^{(k)}\) is the largest among all of its ten components. The recognition rate is defined as the ratio of the number of correctly recognized images to the total number of images.

Echo state network

As a reference for the performance of the skyrmion-lattice reservoir, we use the famous echo state network (ESN) model49,50 to solve the recognition task for the same set of the handwritten digits. For the kth digit, the reservoir state \(\varvec{\psi }^{(k)}_{\textrm ESN}\) and output function \(\varvec{y}^{(k)}_{\textrm ESN}\) in an ESN are defined via a nonlinear hyperbolic tangent as

$$\begin{aligned} \phi ^{(k)}_{n,p+1}= & {} \text {tanh}\Big [W_{\text {in},nm}S^{(k)}_{\text {in},mp}+W_{\text {res},nm}\phi ^{(k)}_{m,p}\Big ],\nonumber \\ \varvec{\psi }^{(k)}_{\textrm ESN}\equiv & {} (\phi ^{(k)}_{1,2},\ldots , \phi ^{(k)}_{1,23},\phi ^{(k)}_{2,2},\ldots , \phi ^{(k)}_{16,2},\ldots ,\phi ^{(k)}_{16,23}, 1)^{\textrm T},\nonumber \\ \varvec{y}^{(k)}_{\textrm ESN}= & {} \varvec{W}_{\textrm out}\varvec{\psi }^{(k)}_{\textrm ESN}. \end{aligned}$$
(4)

Here in the first line the repeated index is summed over, and \(\varvec{\phi }^{(k)}\) is a 16 \(\times\) 23 matrix with the first column being set as \(\phi ^{(k)}_{m,1}=0\) for all k, m51. This zero first column is disregarded when constructing the reservoir state as in the second line. Both \(\varvec{W}_{\textrm in}\) and \(\varvec{W}_{\textrm res}\) are fixed 16 \(\times\) 16 matrices with elements randomly distributed in \((-1,1)\), while \(\varvec{W}_{\textrm out}\) is a 10 \(\times\) 353 output weight matrix to be optimized, whose dimension is the same as that in the skyrmion reservoir described above. To ensure the echo state property, we normalize \(\varvec{W}_{\textrm res}\) by dividing all its elements by 1.01r, with r being the spectral radius (absolute value of the largest eigenvalue) of the initial unnormalized random matrix24,49,50,51, in such a way to make the spectral radius of the normalized matrix become less than one.

As another reference, we consider the digit recognition performed without the reservoir as well. In this case the reservoir state is simply formed by flattening out the 16 \(\times\) 22 input matrix \(\varvec{S}^{(k)}_{\textrm in}\) into a 353-component vector with an additional constant bias element. This is purely a linear transformation of the input data. Following the same training procedure described above, the recognition rate is about 50%. In literature the recognition rates by this linear model range from about 10% for the handwritten digit recognition32 to 70% for spoken digit recognition31. It may depend on the dimensions of weight matrix, ratio between training and testing data numbers, and/or the specific gradient descent algorithms. Here the same training algorithm as described above is applied for all three performances (i) by the skyrmion lattice, (ii) by the echo state network, and (iii) without reservoir, for an unbiased comparison.