\section{Introduction} \label{sec:intro}

The majority of existing methods for EEG data processing and source analysis falls into the following categories: parametric inversions (the dipole fitting), current-reconstruction methods (\cite{PM1994,PM2002,Hamalainen1994,Lin2006,Baillet2001,Grech2008}) and, somewhat less frequently, beamforming methods (\cite{VanVeen1997,Mosher1998,Hillebrand2005, Sekihara2008}. The approaches are very different (sometimes, fundamentally). All of them, however, use the EEG measurements as enter data and, thus, highly depend on its quality.

The measured electrical potential is strongly distorted by outer (with respect to the brain) head tissues.
The current research is aimed to improve the EEG data via partial excluding the influence of the outer tissues on the EEG signal.

Since the outer compartments (tissues) do not contain any sources producing the objective signal, the electric potential related to the brain sources satisfies the Laplace's equation there.
Measuring the potential on the part of the head surface makes it possible to state the Cauchy problem for the Laplace's equation in order to map the data from the scalp to the brain surface. Such mapping is an auxiliary problem, allowing to sufficiently increase the accuracy of source localization due to better quality and spatial resolution of enter data.

The Cauchy problem for the data mapping in the context of EEG was considered in \cite{clerc2007}, where authors propose to use the boundary elements method (BEM) in order to construct the cost functional for further minimization of it. The method, however,  probably may have some drawbacks common for the BEM techniques. It is known that the conventional BEM (i.e., double layer formulation) accuracy degrade when the distance between the source and one of the surfaces becomes smaller \cite{DeMunck}. Beyond this, when two surfaces approach each other, the resulting system of linear equations tends to be singular. Further, conventional BEM cannot incorporate conductivity anisotropies. Finally, numerical comparison presented in \cite{Hyde} indicated much higher modelling errors in the BEM versus finite element method (FEM).

In the current paper, we use the technique based on the mixed quasi-reversibility (MQR) method for linear finite elements and proposed by L.Bourgeois in \cite{Bourgeois}.
The method allows to reduce the Cauchy problem to the system of linear equations per compartment, built using two regularization parameters, the choice of which is also considered in \cite{Bourgeois, Bourgeois_2005}.
As we show later, the application of our algorithm enables us to 'focus' a measured signal, sufficiently improving its spatial resolution even in complicated cases of in-brain source currents.

Regarding practical use of the proposed method, it can be applied either directly to detect active parts of the cortex or as an intermediate procedure to propagate data from scalp to cortex before employing other processing techniques. Though in this study we focused only on mathematical and numerical aspects, we believe that the proposed method may serve as a basis for fast and accurate algorithms for analysis of real EEG data.

The paper is organized as follows. In Section \ref{sec:equations}, we describe the conductivity model of the human head, discuss the approximations being used and write out the governing equation of the EEG. In Section \ref{sec:Cauchy}, we discuss the Laplace's equation for the sourceless compartments of the head and present some important notes about our implementation of the MQR method. The Section \ref{sec:numerical_results} presents some numerical results. Concluding remarks are given in Section 5.

\section{Mathematical description of EEG} \label{sec:equations}
The current section is aimed to introduce the notations, which will be used in paper, and remind the reader about some aspects of EEG together with its assumptions and equations.

\subsection{The computational domain and sources} \label{sec:domains}

The head volume in general and the brain in particular consist of many parts and organs: scalp, skull bone, cerebrospinal fluid, the brain itself and its components.
Using the macroscopic approximation, we assume the conductivity to be identical across the particular compartment of the head, i.e., we expect the conductivity to be the piecewise-constant distribution over the computational domain.

Consider the domain $\Omega\subset \mathbb{R}^3$ with a piecewise-smooth
boundary $\partial\Omega$. The domain under consideration represents the volume of the head and is assumed to consist of several subdomains
\begin{equation*}
	\Omega_i, i=1...N_d: \Omega = \bigcup_{i=1}^{N_d} \Omega_i,
\end{equation*}
representing corresponding $N_d$ organs inside the head (see Fig.\ref{fig:domain_structure}). Each organ expected to have the constant electrical conductivity $\sigma_i$ (see \cite{tissue_conductivity1}).

\textbf{We further assume a nested domain topology of the introduced partitioning.}
We denote the outer boundary of each volume $\Omega_m$ as $\partial\Omega_m$. For convenience, indexing can be done starting from outer areas (scalp, skull) to inner areas. In such way, the outer subdomain $\Omega_1$, for example, is bounded with the surfaces $\partial\Omega_1$ (outer surface) and $\partial\Omega_2$ (inner surface). The outer head surface in such indexing is $\partial\Omega \equiv \partial\Omega_1$.

\begin{figure}[h!]
	\begin{center}
		\includegraphics[width=0.5\linewidth]{head-crop}
	\end{center}
	\caption{Computational domains and boundaries.
	\label{fig:domain_structure}
	}
\end{figure}

In the EEG technique, the measurements are being performed in a bounded area of the surface located commonly at the top of a head. In further consideration, this part of the boundary, which we call \emph{accessible part of the boundary}, is represented by the subdomain of the outer surface $\Gamma_1\subset\partial\Omega_1$ (see Fig. \ref{fig:domain_structure}). Generally, if the domain $\Omega$ consists of more than two sourceless domains, we define the area $\Gamma_i$ for each of them as follows:
\begin{equation*}
	\Gamma_i = \begin{cases}
		\Gamma_1, & i=1; \\
		\partial\Omega_i, & i\ge2.
	\end{cases}
\end{equation*}

For further discussion we also need to introduce the rest of the boundary $\Pi_i$ for the subdomain $\Omega_i$, which we call \emph{inaccessible part of the boundary}. This part of the boundary contains unknown data, which we will need to restore.
\begin{equation}
	\Pi_i = \begin{cases}
		\big(\partial\Omega_1 \setminus \Gamma_1\big) \cup \partial\Omega_2, & i=1; \\
		\partial\Omega_i, & i\ge2.
	\end{cases}
\end{equation}

EEG method measures the electrical potential evoked by the biochemical currents inside a head. There are two main kinds of such currents: currents generated in the brain and the muscle currents. Since the EEG technique aimed to study the electrical activity of the brain, the signal of interest is generated by brain currents. Muscle currents have much higher frequencies and thus can be easily filtered. Therefore, we can assume the current sources are located only inside the brain (more accurately - on the brain surface), while all other areas do not contain electrical sources: $supp\big(\mathbf{J}(\mathbf{x})\big) \subset \Omega_{N_d}$.

The structure of the calculation domain for simplified head models is depicted in Fig.\ref{fig:domain_structure}.

\begin{remark}\label{remark:domains}
	For simplicity, in the current article, we take into account only simple head models consisting of two domains: the brain is represented by the inner domain $\Omega_2$, and the outer tissues (scalp, skull, cerebrospinal fluid, etc.) are described by the only one outer domain $\Omega_1$. Despite that, all other reasoning is correct with respect to numerous sourceless domains, and, thus, the algorithm described in further consideration can be easily applied to a real head model containing a massive amount of sourceless domains.
\end{remark}

\subsection{Governing Equation}

As mentioned above, the conductivity $\sigma(\mathbf{x}), x\in \Omega$ assumed to be a piecewise-constant distribution over the volume under consideration (head). The conductivity of the media outside the head expected to be zero:
		\begin{equation}
			\sigma(\mathbf{x}) = \begin{cases}
				\sigma_i, & \mathbf{x}\in \Omega_i \subset \Omega; \\
				0, & \mathbf{x}\notin \Omega.
			\end{cases}
		\end{equation}

Given the frequencies of the neural currents within the interval of 1-100 Hz, the quasi-static approximation is justified. This was proved in numerous well-known works on EEG, e.g. \cite{Sarvas, DeMunck, Plonsey_1969, Hamalainen, Zakharova}.

Let the volumetric distribution of a primary neuronal current density (the source) be denoted as $\mathbf{J}^p$. Under the assumptions, the electrical potential satisfies the equation
\begin{eqnarray}
	\nabla\cdot(\sigma\nabla U) = \nabla\cdot \mathbf{J}^p, \quad \mathbf{x}\in\Omega. \label{U_equation} \
\end{eqnarray}
The potential further satisfies Neumann boundary condition on $\partial\Omega$:
\begin{equation}\label{V_boundary_cond}
	\mathbf{n}\cdot\nabla U = 0, \quad x\in \partial\Omega.
\end{equation}
This boundary condition follows from the assumption of zero conductivity of the media outside the head \cite{DeMunck}.

We assume the existence of the described in subsection \ref{sec:domains} \emph{accessible part of the boundary} $\Gamma_1 \subset \Omega$, on which the electrical potential $U(\mathbf{x})$ is known:
\begin{equation}\label{boundary_cond_V_B}
	U|_{\Gamma_1} = u(\mathbf{x}),  \quad \Gamma_1 \subset
		\partial\Omega \subset \mathbb{R}^3.
\end{equation}

\begin{remark}\label{remark:conductivity}
Equation (\ref{U_equation}) contains the coefficient $\sigma(\mathbf{x})$ - distribution of the tissues conductivity over a head volume. This distribution can be obtained using methods such as magnetic resonance imaging (MRI) or computed tomography (CT). There are also other techniques for conductivity distribution acquisition. Additionally, some standard models of such distribution are often employed. Regardless of the approach, in the current work we assume the distribution $\sigma(\mathbf{x})$ to be known at every point $\mathbf{x}\in\Omega$.
In biomedical applications, the discussed conductivity varies within a considerable interval. For example, the conductivity of the brain is about one order greater than the conductivity of the tissues located outside the brain (see, e.g., \cite{tissue_conductivity1, tissue_conductivity2}). This fact causes distortion of the electrical potential measured on the outer surface of the domain (head). This situation will be discussed later, in Section \ref{sec:numerical_results}. In the current paper, the conductivity of the inner part of the model, representing the brain, was taken as a constant equal to $2.2$ S/m, and the conductivity of the outer sourceless part was taken as $0.2$ S/m.
\end{remark}

\section{The Cauchy problem for sourceless domains} \label{sec:Cauchy}

Areas outside a brain volume consist of skin, skull bones, muscles. Muscles in a tensed state can also produce rather strong currents and, consequently, electric field. However, they generate currents with much higher frequencies, which allows to filter it out. Accordingly, we consider the "useful" signal to be produced only by brain currents.
Thus, in our consideration we can expect the tissues outside a brain do not contain any sources, which assumption is important for further reasoning.

In sourceless homogeneous volumes the equation (\ref{U_equation}) will take a form:
\begin{eqnarray}
	\Delta U = 0, & \mathbf{x}\in \partial\Omega_i, & i=1,...,N_{sl},\label{laplace_U}\
\end{eqnarray}
where $N_{sl}$ - the number of sourceless domains. Indeed, the equation above can be obtained from (\ref{U_equation}) by nulling the right-hand side due to absence of sources with the assumption that the conductivity $\sigma_i = const$.

We can also write the Dirichlet and Neumann boundary conditions for sourceless domains:
\begin{eqnarray}
	U(\mathbf{x})|_{\Gamma_i} = u_i(\mathbf{x}), \\
	\mathbf{n}\nabla U(\mathbf{x})|_{\Gamma_i} = \frac{\sigma_{i-1}}{\sigma_{i}} \mathbf{n}\nabla \widetilde{U}|_{\Gamma_{i}} \equiv g_i(\mathbf{x});
\end{eqnarray}
where $\widetilde{U}$ denotes the potential on the outer side of the interface $\partial\Omega_i$. The boundary conditions are obtained from requirements for potential and its normal derivative to be continuous inside the computational domain $\Omega$. Assuming the conductivity $\sigma(\mathbf{x}) = 0, \mathbf{x}\notin\Omega$, we can easily obtain zero Neumann condition, common for EEG problem.

\subsection{Statement of the problem} \label{sec:statement}

The equation (\ref{laplace_U}) enable us to state an auxiliary problem, which we named the \emph{propagation problem}.

\textbf{Propagation problem for the electrical potential}. \emph{Assume the function $u(\mathbf{x}), \quad \mathbf{x}\in\Gamma_1$ to be known. Find the function $U(\mathbf{x})$, such that:}
\begin{eqnarray}\label{Cauchy_U}
	\Delta U = 0, & \mathbf{x} \in \Omega_i; \label{extra1}\\
	U|_{\Gamma_1} = u_i(\mathbf{x}), & \mathbf{x} \in \Gamma_i; \label{extra2}\\
	\mathbf{n}\cdot\nabla U|_{\Gamma_i} = g_i(\mathbf{x}), & \mathbf{x}\in \Gamma_i \label{extra3}.
\end{eqnarray}

\begin{remark}\label{remark:medical applications}
In medical applications, the interface between the brain and cerebrospinal fluid contains sources. Moreover, the mentioned interface is the only area, where the electrical currents are possible except muscle currents. Due to this fact, from the mathematics' point of view, the Cauchy problem for the Laplace's equation can be stated in this case only for open domain, which does not contain the brain surface itself. Despite this, our numerical experiments showed that its solution (i.e., backpropagation of the potential to the brain surface) is accurate enough and can be used in practice. This caused by the fact we do not compute the potential only on boundary, but compute it in the whole sourceless domain. Thus, we can estimate the electrical potential at the points arbitrary closed to the desirable surface. The respective numerical experiments will be shown in Section \ref{sec:numerical_results}.
\end{remark}

The Cauchy problem, in general, is a well-known problem, and the Cauchy problem for the Laplace's equation with incomplete boundary data is one of the classical examples of an ill-posed problem. During the last years, a number of validated approaches to its solution were presented. The approaches can be divided into three main categories: approaches based on the optimization (least square, Tiknonov's regularization etc), methods based on construction via quasi-reversebility of the linear systems of equations and solving it and iterative approaches (see, for example, \cite{Engl}).

The first category includes classical Tikhonov's regularization \cite{Berntsson, Li}, which is well-established and mostly used, but have some issues with the solution accuracy \cite{JIIP-2016}. Additionally, it includes a set of powerful methods proposed by the group of M. Klibanov together with one of the authors of the current paper. The approach based on the construction of strictly convex Tikhonov-like functional weighted with a problem-depending Carleman weight functions \cite{JIIP-2016, NONRWA-2016}. Despite a great accuracy of the solution, the method, however, is rather hard to implement and is better suitable for nonlinear problems, for which it was designed. Other solutions based on the optimization approaches are presented, for example, in \cite{clerc2007}.

The second category is probably the most known category of the approaches. It based on usage of quasi-reversibility method, which is being used in order to state a weak formulation of a problem, which, after discretization, is turned into linear system of algebraic equations. Some methods are presented, for example, in \cite{Bourgeois, Bourgeois_2005, Cao, Cao_Klib_Per_2009, Darde-2013}.

\subsection{The method} \label{sec:method}

The equations (\ref{extra1}-\ref{extra3}) are the Cauchy problems for the Laplace's equation. Uniqueness of this problem was proved in many ways. The most effective and beautiful proof can be found in \cite{LavRomShish}. Solution of the problem considered in a set of works (see, for example, \cite{Bourgeois, Engl, Berntsson, Li, NONRWA-2016, Cao, Cao_Klib_Per_2009}).

In the current paper, we use the MQR method provided for linear finite-element approximations in \cite{Bourgeois} and \cite{Bourgeois2006}. As one can see below, the method depends on two regularization parameters. The most efficient way of choosing the regularization parameters is the balancing principle, which is described in \cite{Cao_Klib_Per_2009}.

For simplicity, we will next consider only the case of two compartments: 'inactive' sourceless compartment $\Omega_1$ and 'active' subdomain $\Omega_2$ containing the electrical sources. Consider the Cauchy problem (\ref{Cauchy_U}). In our consideration the part of the boundary $\Gamma_1$ is shown, e.g., in Fig.\ref{fig:domain_structure}.

We introduce the following notations:
\begin{eqnarray}\label{MQR_subspaces}
	V_0 = \{h\in H^1(\Omega_1) : h|_{\Gamma_1} = 0 \}, \label{V0}\\
	V_1 = \{h\in H^1(\Omega_1) : h|_{\Pi_1} = 0 \}, \label{V1}\\
	\tilde{V}_0 = \{h\in H^2(\Omega_1) : h|_{\Gamma_1} = u(x)\}. \label{V00}
\end{eqnarray}

Let $0<\epsilon\ll1, 0<\delta\ll1$ be two fixed small numbers (regularization parameters).

\textbf{Weak formulation of MQR method \cite{Bourgeois}}. \emph{Find a pair of functions $(U, \lambda) \in \tilde{V}_0 \times V_1$, such that:}
\begin{eqnarray}\label{weak_QR}
	\epsilon \int\limits_{\Omega_1} \nabla U\cdot \nabla h dx + \epsilon \int\limits_{\Omega_1} uhdx +
		\int\limits_{\Omega_1} \nabla h \cdot \nabla \lambda dx = 0, & \forall h\in V_0; \label{weak_qr1}\\
	\int\limits_{\Omega_1} \nabla U \cdot\nabla \mu dx - \delta\int\limits_{\Omega_1}
		\nabla\lambda\cdot\nabla\mu dx - (1+\delta)\int\limits_{\Omega_1} \lambda\mu dx =
		\int\limits_{\Gamma_1} g_1 \mu d\Gamma_1, & \forall \mu \in V_1. \label{weak_qr2}
\end{eqnarray}

In \cite{Bourgeois} one can find proofs of existence and uniqueness of the solution, its stability and convergence of the solution of (\ref{weak_QR}) and the exact solution of the problem (\ref{Cauchy_U}) with $(\epsilon, \delta) \to 0$, when the error in the enter data $u(\mathbf{x}) \to 0$.

\paragraph{The finite element approximation}
In order to make further consideration shorter we start with the fact the second term of the equation (\ref{weak_qr1}) and the third term in (\ref{weak_qr2}) can be omitted, the justification of which one can find in \cite{Bourgeois}. We also note that since we consider the Cauchy problem for the outer domain $\Omega_1$, the function $g(\mathbf{x}) \equiv 0$.

Omitting the index under the function $u_1(\mathbf{x})$, we can rewrite the system under consideration as follows:

\begin{eqnarray}\label{final_system}
	\epsilon \int\limits_{\Omega_1} \nabla U\cdot \nabla h dx +
		\int\limits_{\Omega_1} \nabla h \cdot \nabla \lambda dx = 0, & \forall h\in V_0, \label{weak_qr1_2}\\
	\int\limits_{\Omega_1} \nabla U \cdot\nabla \mu dx - \delta\int\limits_{\Omega_1}
		\nabla\lambda\cdot\nabla\mu dx = 0, & \forall \mu \in V_1. \label{weak_qr2_2}
\end{eqnarray}
where the function $u(\mathbf{x})$ is known on $\Gamma_1$. In our consideration the part of the boundary $\Gamma_1$ is shown, e.g., in Fig.\ref{fig:domain_structure}.

Among many different methods available to numerical solution of the diffusion equation (finite-difference, finite-volume, see e.g. \cite{Yavich-2007}, finite-element etc.), we used the finite-element method with linear basis functions on tetrahedrons as it combines simplicity and flexibility.

Let $X_h$ be the linear finite element space with the basis $h_i(\mathbf{x}), i=1,...,N$ (here $N$ is the number of nodes). Denote the following subspaces, which are the finite-dimensional analogs of the corresponding subspaces (\ref{MQR_subspaces}):
\begin{eqnarray}
	X_{0} = \{h\in X_h: h|_{\Gamma_1} = 0\}, \\
	X_1 = \{h\in X_h: h|_{\Pi_1} = 0\}, \\
	\tilde{X}_0 = \{h\in X_h : h - U_0 \in X_0\}.
\end{eqnarray}

Here $U_0$ is the approximation of the function with the following properties:
\begin{eqnarray}
	U_0|_{\Gamma_1} = u(\mathbf{x}), \quad \partial_n U_0|_{\Gamma_1} = 0, \quad U_0 \in H^1(\Delta, \Omega),
\end{eqnarray}
where
\begin{equation}
	H^1(\Delta, \Omega) = \{U\in H^1(\Omega): \Delta U \in L^2(\Omega)\}.
\end{equation}

The formulation of the MQR method is still the same as (\ref{weak_qr1_2}-\ref{weak_qr2_2}) with the difference that the pair of functions to be found now belongs to $(U, \lambda) \in \tilde{X}_0 \times X_1$.

Finite dimensional approximations of the functions $U$, $\lambda$, $g_1$ can be written as follows:

\begin{eqnarray}\label{femappr_func}
	U(\mathbf{x}) \approx \sum\limits_{i=0}^N U_i h_i(\mathbf{x}), & \mathbf{x} \in \Omega_1, \label{femappr1}\\
	\lambda(\mathbf{x}) \approx \sum\limits_{i=0}^N \lambda_i h_i(\mathbf{x}), &  \mathbf{x} \in \Omega_1,
	\label{femappr2}
\end{eqnarray}
where $N_{\Gamma_1}$ is the number of nodes belonging to $\Gamma_1$.

For convenience introduce the notation $i(D) = \{i: \mathbf{x}_i \in D\}$ to be indices of nodes of a finite element mesh, which belong to some domain $D\subset \mathbb{R}^3$. We also use the following notation:
\begin{equation*}
	A_{gh} = \int\limits_{\Omega} \nabla h_i(\mathbf{x}) \cdot \nabla h_m(\mathbf{x}) d\mathbf{x}, \quad g,h=1,...,N.
\end{equation*}
Substituting (\ref{femappr1}, \ref{femappr2}) into (\ref{weak_qr1}, \ref{weak_qr2}) and assuming $v = h_j \in X_0, \mu = h_k \in X_1$, we obtain:
\begin{eqnarray}
	& \epsilon \sum\limits_{i(\Omega_1)} U_i A_{im} +
		\sum\limits_{j(\Omega_1)} \lambda_j A_{jm} = 0,
		& m=m(\Omega_1\setminus\Gamma_1), \label{findim_system1} \\
	& \sum\limits_{l(\Omega_1)} U_l A_{lk} +
		\delta \sum\limits_{n(\Omega_1)} \lambda_n A_{nk} =
		b_k, & k=k(\Omega_1\setminus\Pi_1), \label{findim_system2}
\end{eqnarray}
where

\begin{equation*}
	b_k = \sum\limits_{i(\Gamma_1)} u_i \int\limits_{\Gamma_1} \nabla h_k(\mathbf{x}) \cdot \nabla h_k(\mathbf{x}) d\mathbf{x}, \quad k = k(\Omega_1\setminus\Pi_1).
\end{equation*}

Since $U(\mathbf{x}) \in \tilde{X_0}$, we can state that $U_i = u_i, i\in i(\Gamma_1)$, and the first term of (\ref{findim_system1}) takes the form:
\begin{equation}
	\epsilon \sum\limits_{i(\Omega_1)} U_i A_{im} = \epsilon \sum\limits_{i(\Omega_1\setminus \Gamma_1)} U_i A_{im} +
		\epsilon \sum\limits_{i(\Gamma_1)} u_i A_{im}, \quad m \in m(\Omega_1 \setminus \Gamma_1).
\end{equation}
Since the function $\lambda(\mathbf{x}) \in X_1$, the second term can be written as follows:
\begin{equation}
	\sum\limits_{j(\Omega_1)} \lambda_j A_{jm} = \sum\limits_{j(\Omega_1 \setminus \Pi_1)} \lambda_j A_{jm}.
\end{equation}
Reasoning in the same way, we obtain for both equations \ref{findim_system1} and \ref{findim_system2}, one can rewrite these equations as follows:
\begin{eqnarray}
	& \epsilon \sum\limits_{i(\Omega_1\setminus\Gamma_1)} U_i A_{im} +
		\sum\limits_{j(\Omega\setminus\Pi_1)} \lambda_j A_{jm} =-a_m,
		& m=m(\Omega\setminus\Gamma_1), \label{findim_system11} \\
	& \sum\limits_{l(\Omega\setminus\Gamma_1)} U_l A_{lk} +
		\delta \sum\limits_{n(\Omega\setminus\Pi_1)} \lambda_n A_{nk} =
		 - b_k, & k=k(\Omega\setminus\Pi_1); \label{findim_system22}
\end{eqnarray}
where
\begin{eqnarray*}
	a_m = \sum\limits_{i(\Gamma_1)} u_i A_{im}, & m\in m(\Omega_1\setminus\Gamma_1), \\
	b_k = \sum\limits_{i(\Gamma_1)} u_i A_{ik}, & k\in k(\Omega_1\setminus\Pi_1). \\
\end{eqnarray*}

Denote the number of nodes located in some domain $D$ as $N_D$. It is easy to see that in the representation above the discreet system contains $2N_{\Omega} - N_{\partial\Omega}$ equations and the same number of variables. It is thus can be solved by an iterative solver.

\section{Numerical results} \label{sec:numerical_results}
The current section is aimed to present the numerical results on the algorithm described above. We start with two cases of the simple model representing the brain with the inner sphere, and outer tissues with the outer spherical layer. Then, we present the numerical results for more realistic model of the head, based on simplified MRI data (\cite{iso2mesh}). In the current paper, we use simulated data. The simulations were implemented with the finite element numerical solution of the equation (\ref{U_equation}) for the sources/conductivity distributions presented below.

\begin{figure}
	\begin{tabular}{cc}
		\includegraphics[width=0.45\linewidth]{fig2_a} &
  		\includegraphics[width=0.45\linewidth]{fig2_b} \\
		(a) & (b) \\
		\includegraphics[width=0.45\linewidth]{fig2_c} &
		\includegraphics[width=0.45\linewidth]{fig2_d} \\
		(c) & (d)
	\end{tabular}
	\caption{Spherical model (sources located under the surface of the active domain): a) simulation of the electrical potential on the outer surface (scalp); b) simulation of the electrical potential on the "active" domain surface; c) simulation of the potential on the surface of the inner sphere; d) the potential on the inner sphere surface mapped from the data depicted in b) via the solution of the Cauchy problem ($\epsilon=0.0001$ and $\delta=0.001$).}
	\label{fig:results_1}
\end{figure}

\subsection{Spherical model with the sources located under the surface of the active domain}
The case considered in this subsection is unrealistic since the currents here are located under the surface of the "active" domain. In truth, the biochemical sources are placed on the interface between a brain and cerebrospinal fluid, i.e., on the brain surface. However, in this case, in the strict sense, the Cauchy problem for the Laplace's equation cannot be stated. Thus, we find it reasonable to start the presentation of our numerical results in the strict case of the sourceless boundary.

The head here is represented with a simple model, which consists of the spherical "active" domain $\Omega_2$ with the radius $r_2=0.8$ (see Fig. \ref{fig:domain_structure}a) and a spherical layer representing the outer sourceless area $\Omega_1$ with the radius $r_1=1$. The conductivity of the inner domain is $\sigma(\Omega_2) = \sigma_2 =2.2  S/m$, and the domain $\Omega_1$ can be characterized with the conductivity $\sigma(\Omega_1) = \sigma_1 = 0.1 S/m$. The electrical sources are depicted with the support (three spots) on Fig.\ref{fig:results_1} a). The support spots contain the constant currents inside it: $\mathbf{J}(\mathbf{x}) = (0, 10^{-7}, 0), \mathbf{x}\in supp(\mathbf{J})$ A.

Fig.\ref{fig:results_1}a) and b) depict the simulated electric potential on the outer surface and the surface of the "active" domain, respectively, while c) and d) show respectively the simulated potential on the 'brain' surface, and the result of its reconstruction using the algorithm presented above.

The finite element mesh for the sourceless domain $\Omega_1$ consists of $484946$ tetrahedra and $90435$ vertices. The calculation time for the Cauchy problem solver implemented with the Matlab is equal to $5.7$ seconds with a laptop equipped with $16$ Gb RAM and CPU Intel Core-i7 7820HQ.

\begin{figure}
	\begin{tabular}{cc}
		\includegraphics[width=0.45\linewidth]{fig3_a} &
		\includegraphics[width=0.45\linewidth]{fig3_b} \\
		(a) & (b) \\
		\includegraphics[width=0.45\linewidth]{fig3_c} &
		\includegraphics[width=0.45\linewidth]{fig3_d} \\
		(c) & (d)
	\end{tabular}
	\caption{Spherical model (sources located on the surface of the inner sphere): a) The current density (sources); b) simulation of the electric potential on the outer surface (scalp); b) simulation of the electric potential on the "active" domain surface; c) the potential on the inner sphere surface, obtain vial solution of the Cauchy problem with regularization parameters $\epsilon=0.00001$ and $\delta=0.001$
		}
	\label{fig:results_2}
\end{figure}

\subsection{Spherical model: more realistic case}
Geometrically, the model under consideration in the current subsection is identical to the one in the previous subsection. It also characterized with the same conductivity distribution. The difference is in the current. First of all, the currents here are located on the interface between active and inactive domains, which is closer to the real-life encephalographic situation. Strictly speaking, the domain $\Omega_1$, in this case, should not contain its $\partial\Omega_2$ boundary surface, and it should be the opened domain. However, since the potential can be calculated via the solution of the Cauchy problem for the surface, which is arbitrary close to the surface $\partial\Omega_2$, we were not surprised that the modelled and restored via the Cauchy problem solution potentials on an active domain surface are relatively close to each other. The following results can be seen in Fig.\ref{fig:results_2}.

The finite element mesh for the sourceless domain $\Omega_1$ consists of $484946$ tetrahedra and $90435$ vertices. The calculation time for the Cauchy problem solver implemented with the Matlab is equal to $5.7$ seconds with a laptop equipped with $16$ Gb RAM and CPU Intel Core-i7 7820HQ.

\begin{figure}
	\begin{tabular}{cc}
		\includegraphics[trim=360 0 250 0,clip,width=0.45\linewidth]{fig4_a} &
		\includegraphics[trim=360 0 250 0   ,clip,width=0.45\linewidth]{fig4_b} \\
		(a) & (b) \\
		\includegraphics[trim=360 0 250 0	,clip,width=0.45\linewidth]{fig4_c} &
		\includegraphics[trim=360 0 250 0	,clip,width=0.45\linewidth]{fig4_d} \\
		(c) & (d)
	\end{tabular}
	\caption{MRI-based model of the head: case with piecewise-constant current distribution on the brain surface.
	(a) - the current, (b) - FEM-simulated potential on the outer surface of the head (enter data for the Cauchy problem,
	(c) - FEM-simulated potential; (d) - reconstruction of the potential on the brain surface on the base of the
	data, depicted with (b).}
	\label{fig:results_3}
\end{figure}

\subsection{Simplified head model based on real MRI data}
The model described in this subsection is a simplified but more realistic model of the head. The mesh was constructed using the iso2mesh software \cite{iso2mesh}. Here the "realistic" distribution of the current density over the brain surface was used.

In order to make the result more realistic, we employed two kinds of current density distributions. Fig \ref{fig:results_3} depicts the model with several 'active' spots (support areas) of the cortex, where the current, corresponding to the spot is constant inside it. The Fig \ref{fig:results_4} shows the situation, when the current is randomely defined in each point inside the support areas (areas itself are similar to the previous case).

The finite element mesh for the sourceless domain $\Omega_1$ consists of $444916$ tetrahedra and $82239$ vertices. The calculation time for the Cauchy problem solver implemented with the Matlab varies in the interval $6.1-11.2$ seconds in dependence of the current distribution used for modelling. All calculations have been provided with a laptop equipped with $16$ Gb RAM and CPU Intel Core-i7 7820HQ.

\begin{figure}
	\begin{center}
		\begin{tabular}{cc}
			\includegraphics[trim=360 0 250 0,clip,width=0.45\linewidth]{fig5_a} &
			\includegraphics[trim=360 0 250 0   ,clip,width=0.45\linewidth]{fig5_b} \\
			(a) & (b) \\
			\includegraphics[trim=360 0 250 0	,clip,width=0.45\linewidth]{fig5_c} &
			\includegraphics[trim=360 0 250 0	,clip,width=0.45\linewidth]{fig5_d} \\
			(c) & (d)
		\end{tabular}
	\end{center}
	\caption{MRI-based model of the head: case with random current, distributed within big areas of the brain surface.
	(a) - the current, (b) - FEM-simulated potential on the outer surface of the head (enter data for the Cauchy problem,
	(c) - FEM-simulated potential; (d) - reconstruction of the potential on the brain surface on the base of the
	data, depicted with (b).}
	\label{fig:results_4}
\end{figure}

\section{Conclusion}
In this paper we presented an approach for Scalp-to-Cortex data mapping, based on finite elements computational scheme. The application of the approach shows great spatial resolution of the result, which is comparable to intracranial EEG (iEEG, ECoG) techniques. The method based on solving the Cauchy problem for the Laplace's equation in homogeneous compartment and can be described as propagation of EEG data collected on the scalp through skull to the brain surface.
Our numerical implementation employs the approximation of the first order finite elements on tetrahedral grids, which enables to model sophisticated irregularly shaped conductors with high precision. The method is considered as a pre-processing procedure, after which classical inverse source localization problem solvers can be applied.
Notably, our algorithm successfully works with minimal soft- and hardware capabilities and does not require professional knowledge. All calculations were provided in Matlab, using only a quad-core laptop and without explicit parallelization; the computational algorithm thus has a big potential on optimization. This allows us to assume that with further improvements our method can be successfully implemented for usage with mobile platforms.
Speaking of future work, we believe the solution of the Cauchy problem described in this paper can be easily adopted for source localization in magnetoencephalography. The other line of our work is to consider the same problem but with more complex algorithms in order to achieve better accuracy and compare numerically the spatial resolution of conventional EEG, ECoG and our method. Both of these areas will be certainly reflected in our future publications.

\bibliographystyle{unsrt}


\section*{Introduction}
Forecasting a future system state is an important task in various fields such as neuroscience, ecology, finance, fluid dynamics, weather, and disaster prevention. If accurate system equations are unknown but a time series is observed, we can use model-free forecasting approaches.
Although there are various model-free methods, if the time series is assumed to be nonlinear and deterministic, one of the most common approaches is delay embedding. According to embedding theorems\cite{Takens1981, Sauer1991a}, we can reconstruct the underlying dynamics by delay embedding (see Methods for details). These theorems ensure that the map from the original attractor to the reconstructed attractor has a one-to-one correspondence. These embedding theorems have also been extended to multivariate data and nonuniform embeddings\cite{Deyle2011}.

When multivariate time series are observed, various embeddings can be obtained on the basis of the extended embedding theorem\cite{Deyle2011}. Ye and Sugihara\cite{Ye2016} exploited this property; they developed multiview embedding (MVE), which combines multiple forecasts based on the top-performing multiple embeddings scored by the in-sample error. MVE can yield a more accurate forecast than the single best embedding, especially for short time series\cite{Ye2016}. Although MVE works well with low-dimensional data, it cannot be simply applied to high-dimensional data. This is because MVE requires forecasts derived from all possible embeddings, whose total number combinatorially increases with the number of variables.
Although we can approximate MVE by randomly chosen embeddings if the total number of possible embeddings is beyond the computable one\cite{Okuno2019}, it is difficult to find high-performance embeddings for high-dimensional data.
Recently, Ma et al.\cite{Ma2018a} presented an outstanding framework, randomly distributed embedding (RDE), to tackle these high-dimensional data. Their key idea is to combine forecasts yielded by randomly generated ``nondelay embeddings'' of target variables. These ``nondelay embeddings'' can also reconstruct the original state space for some cases and drastically reduce the possible number of embeddings. Ma et al.\cite{Ma2018a} also showed that small embedding dimensions worked fine, even for high-dimensional dynamics, and successfully forecasted short-term high-dimensional data using the RDE framework.
Although RDE showed outstanding results, there is room for improvement for some specific tasks. For example, ``nondelay embedding'' may overlook important pieces of information if delayed observations carry such pieces of information, e.g., the upstream river height for flood forecasting. A random distribution may also be a problem for the case where only partial variables are valid system variables; in this case, most of the random embeddings are invalid and not suitable to forecast. In addition, several important hyperparameters---the embedding dimension and the number of embeddings to combine---must be manually or empirically selected.
Although ensembles tend to yield better results with significant diversity among its members\cite{Sollich96, Kuncheva2003}, RDE only aggregates forecasts for a fixed embedding dimension with a fixed number of embeddings to combine, regardless of their forecast performance.

In this study, we propose another forecasting framework that overcomes the disadvantages of MVE and RDE. Our key idea is to prepare diverse ``suboptimal embeddings'' via combinatorial optimization.
We combine the optimal number of these embeddings to maximize the performance of the combined forecast.
In contrast to the existing frameworks, the suitable embedding dimensions and the number of embeddings to combine are automatically determined through this procedure.
Our proposed framework achieves better forecast performance than the existing frameworks for high-dimensional toy models and a real-world flood dataset.

\section*{Results}
\subsection*{Proposed forecasting framework}
Our proposed framework yields a single forecast to combine multiple ones, which are obtained by multiple embeddings. These embeddings are obtained through two steps.
First, we obtain suboptimal embeddings, which are suboptimal solutions to minimize the forecast error.
The aim of this procedure is to yield diverse embeddings to improve the performance of the combined forecast because such diversity is important for the performance of the ensemble\cite{Sollich96, Kuncheva2003}.
Second, we make the in-sample forecasts based on all suboptimal embeddings following MVE and RDE but pick the optimal number of embeddings to minimize the error of the combined forecasts. A schematic of the procedure is illustrated in Figure~\ref{fig:procedure}. We describe each step in detail in the following sub-subsections.

\begin{figure}[ht]
    \centering
    \includegraphics[width=.8\linewidth]{fig/procedure}
    \caption{Schematic of the proposed forecasting procedure. We prepare a pool of suboptimal embeddings in the first step. We solve $K$ combinatorial optimization problems to obtain various embeddings in this step. Next, we pick $\hat{k}_p$ embeddings to minimize the error of the combined forecast in the second step. We combine the forecasts obtained by the $\hat{k}_p$ embeddings to test the time series.}
    \label{fig:procedure}
\end{figure}

\subsubsection*{First step: Preparing suboptimal embeddings}
In our framework, we first obtain diverse suboptimal embeddings to minimize the in-sample forecast error. With the observed time series $y(t) \in \mathbb{R}^n$, a set of possible embeddings $\mathcal{E}$, and a corresponding delay vector $v_e (t) \in \mathbb{R}^E$ with $e \in \mathcal{E}$, we can obtain the $p$-steps-ahead forecast at time $t$ by a map $\psi_e$ as follows:
\begin{equation}
    \hat{y}_f^e (t+p|t) = \psi_e (v_e (t)).
\end{equation}
We can obtain such maps by suitable fitting algorithms, e.g., the conventional method of analogues\cite{Lorenz1969}, as in MVE, and Gaussian process regression\cite{Rasmussen2006}, as in RDE. Throughout this paper, we apply the same algorithm, a variation of the method of analogues (see Supplementary Information), to all forecasting frameworks including the proposed one. This is because we need to compare the performance of the frameworks on a fair basis that is not affected by the performance of each fitting algorithm.

To obtain high-performance embeddings in an efficient manner, we solve combinatorial optimization problems instead of random selection or using brute force.
We split the set of training time indices $\mathcal{T}_{train}=\{t \mid t < 0\}$ into $K$ datasets. Then, we solve the following optimization problem for each $K' \in \{1,2,...,K\}$ to minimize the error:
\begin{equation} \label{eq:fitness}
    \minimize_{e \in \mathcal{E}} \ \sum_{t \in \mathcal{T}_{train}^{K'}} \sum_p \|\hat{y}_f^e (t+p|t) - y_f (t+p)\|,
\end{equation}
where $\mathcal{T}_{train}^{K'}$ is the $K'$th set of training time indices and $\|\cdot\|$ is an appropriate norm. We use the $L^2$ norm throughout this paper. Note that $\hat{y}_f^e$ is computed in an in-sample manner; we use all of $\mathcal{T}_{train}$, except for the current time, as its library.
When we consider embeddings up to $l$ lags to embed the latest observation of the target variable, the number of possible embeddings is $2^{nl-1}$. If $nl$ is sufficiently large, it is almost impossible to compute all possible embeddings because of combinatorial explosion.
Here, we take a straightforward approach to minimize the forecast error: an evolution strategy that interprets an embedding as a one-dimensional binary series\cite{Vitrano2001, Small2003}. See Methods for details. Throughout this paper, we applied the $(\mu+\lambda)$ evolution strategy\cite{Schwefel1981a}.

In the process of each optimization, we store not only the best embedding but also the ``Hall of Fame,'' which preserves the best individuals in each generation. From the Hall of Fame, we select the best $M$ solutions that satisfy the condition that the Hamming distance $d(\cdot,\cdot)$ is larger than or equal to a certain threshold $\theta$ as follows:
\begin{equation} \label{eq:embedding_condition}
    d(e_i, e_j) \geq \theta, \ \forall i,j; i \neq j,
\end{equation}
where $e_i$ is the $i$th embedding of the sorted Hall of Fame.
Namely, we sort the Hall of Fame by the fitness (the objective function of equation \ref{eq:fitness}) and pick the best embeddings in order from the sorted Hall of Fame to satisfy equation~(\ref{eq:embedding_condition}).
This condition prunes similar embeddings and keeps diverse embeddings for the next step.
To apply this procedure, we obtain $K \cdot M$ suboptimal embeddings in total.
See also the first part of Algorithm~\ref{alg:framework} in Methods.

The procedure of the $K$ times optimizations yields diverse but high-performance embeddings.
If a dataset contains a large number of variables, it is extremely difficult to obtain the exact solution by evolution strategies or any other metaheuristic algorithms.
These algorithms find suboptimal solutions for most cases, and these solutions can be easily changed by slight changes in the optimization conditions.
The procedure rather exploits this property; the $K$ times optimization can find $K$ sets of diverse suboptimal solutions.
Instead of minimizing the whole in-sample error with $K$ different conditions, we minimize each $K$ split dataset and utilize the history of solutions to significantly reduce the computational time.
We numerically show the effect of $K$ times optimization in Supplementary Information.

To create diverse embeddings and improve the performance of the combined forecast, we can optionally consider $L$ finite impulse response (FIR) filters for each suboptimal embedding; namely, $K \cdot L \cdot M$ embeddings are obtained in total. These filtered embeddings are justified by the Filtered Delay Embedding Prevalence Theorem\cite{Sauer1991a}, and a study has shown that the combination of various forecasts via linear filters improves the forecast result in some cases\cite{Okuno2017a}. Although we can compute suboptimal embeddings for each filter, we simply consider $L$ filters based on the obtained $K \cdot M$ embeddings. Empirically, the optimization of the embeddings for each filter does not realize a significant improvement in performance to justify the additional computational cost.
See Methods for the details of forecasting filtered time series.
In this paper, we set the filter coefficients as $h_i(0)=1, h_i(1)=\rho \ \forall i$ for all numerical examples.

\subsubsection*{Second step: Combining forecasts to minimize the in-sample error}
Next, we pick the optimal number of embeddings from the $K \cdot L \cdot M$ embeddings to minimize the error of the combined forecast.
We compute in-sample forecasts by the embeddings for all of $\mathcal{T}_{train}$, as done in the MVE or RDE aggregation schemes.
We determine the optimal number of embeddings to integrate forecasts, which is neither the square root of the possible number of embeddings (as MVE does) nor the embedding dimension (as the RDE aggregation scheme does).
Here, we define a tuple of the forecast indices sorted by the in-sample error for step $p$ as $\mathcal{I}_p$, and the $i$th best embedding is written as $\mathcal{I}_p(i)$. Note that $\mathcal{I}_p$ is separately computed for each forecast horizon because accurate embeddings may differ according to the forecast horizons, as discussed later. We obtain $\hat{k}_p$ to minimize the combined forecast as follows:
\begin{align}
    Y_k(t+p|t) := 1/k \sum_{i=1}^k \hat{y}_f^{\mathcal{I}_p(i)} (t+p|t), \\
    \hat{k}_p = \argmin_{k=1,2,...,P} \sum_{t \in \mathcal{T}_{train}} [Y_k(t+p|t) - y_f(t)]^2, \\
    \hat{y}_f(t+p|t) = Y_{\hat{k}_p}(t+p|t).
\end{align}
We can obtain $\hat{k}_p$ in a brute-force manner with a very small computational cost because the $k$th combined forecast $Y_k(t+p|t)$ can be recursively computed as follows:
\begin{equation}
    Y_{k+1}(t+p|t) = \left[k Y_k(t+p|t) + \hat{y}_f^{\mathcal{I}_p(k+1)}(t+p|t) \right] / (k+1),
\end{equation}
where $k=1,2,...,P-1,\  Y_1(t+p|t) = \hat{y}_f^{\mathcal{I}_p(1)}(t+p|t)$. We compute $Y_k$ and the corresponding error for all $k$ and pick the optimal value to minimize the squared error. Note that $\hat{k}_p$ can be different for each $p$; empirically, $\hat{k}_p$ tends to be larger for a longer step.

For unobserved samples $t \geq 0$, the forecast $\hat{y}_f (t+p|t)$ is computed as follows:
\begin{equation}
    \hat{y}_f (t+p|t) =  \sum_{k=1}^{\hat{k}_p} \hat{y}_f^{\mathcal{I}_p(k)}(t+p|t) / \hat{k}_p.
\end{equation}
Note that only $\hat{k}_p$ forecasts need to be computed rather than $K \cdot L \cdot M$ forecasts.
We summarize the whole algorithm in Algorithm~\ref{alg:framework} in Methods.

\subsection*{Numerical experiments}
\subsubsection*{Toy models}
We demonstrate our framework with toy models.
We first tested with the 10-dimensional Lorenz'96I model \cite{Lorenz1995a}. The model is described as 10-dimensional differential equations, and we used five of their system variables and five random walks as a dataset---namely, Lorenz'96I's $x_0, x_1, x_2, x_3, x_4$ and random walks $x_5, x_6, x_7, x_8, x_9$. We forecasted $x_0$ up to 10 steps.
We compared the proposed framework with existing frameworks---namely, randomly distributed embedding with an aggregation scheme (RDE), multiview embedding (MVE), a variation of MVE (state-dependent weighting\cite{Okuno2019} (SDW)), and the single best embedding via the $(\mu+\lambda)$-ES algorithm. See Supplementary Information for the detailed conditions.

We first tested the performance with the data length of 4000 for the database and 500 for the evaluation. Figure~\ref{fig:lorenz96_rmse}(a) shows that the proposed forecast achieved the best performance among the existing frameworks up to 10 steps ahead with the length of 4000.

We conducted a sensitivity analysis with the data length. We set the data length, i.e., the size of the training dataset, to \{100, 200, 500, 1000, 2000, 4000\} and tested the performance with 500 samples for all cases. Figure~\ref{fig:lorenz96_rmse}(b) shows that the proposed five-steps-ahead forecast yielded the minimum error among the existing frameworks for all the data lengths.

We also checked the robustness against observational noise. We added Gaussian observational noise scaled by the standard deviation of the target variable and tested the performance with a scale range of \{0.0, 0.1, 0.2, 0.3\}. We set the training data length as 500 and evaluated with 500 samples.
Figure~\ref{fig:lorenz96_rmse}(c) shows that a lower noise level results in better performance for the proposed framework compared with the other frameworks. Although the proposed framework achieved the best performance up to 10\% noise, the performance was degraded by high-level noise and is the second best for 30\% noise.

We profiled which variables were embedded for the combined forecasts for the case where the data length is 4000 without noise. We summed $\hat{k}_p$ embeddings for each variable, where 1 means embedded and 0 means not embedded. Then, we computed the proportion of embedding of each variable, normalizing the sum to one for each $p$. Figure~\ref{fig:lorenz96_profile}(a) shows that random walks $x_5, x_6, x_7, x_8$, and $x_9$ were seldom embedded with the proposed framework. This result suggests that our proposed framework can select the appropriate variables. The results also show that the proportion of embedding of the target variable $x_0$ is large for short-term prediction and small for longer-term prediction. This result is consistent with that of an existing study\cite{Chayama2016a}.

We also profiled the filter coefficient $\rho$ and the embedding dimension $E$ of the proposed combined forecasts.
We averaged $\rho$ and $E$ over the $\hat{k}_p$ combined members for each step. Figure~\ref{fig:lorenz96_profile}(b) shows that smaller dimensions are selected for shorter steps, and larger dimensions are selected for longer predictions. The figure also shows that smaller and larger values of $\rho$ are selected for shorter and longer steps, respectively.
Because we employ high-pass filtering $(\rho\leq0)$ with this model, a smaller $\rho$ focuses on a smaller time period and vice versa.

We investigated the number of combined forecasts and the forecast errors to check whether the proposed integration scheme works properly. Although the proposed framework determines $\hat{k}_p$ on the basis of the in-sample error, $\hat{k}_p$ also minimizes the test error in this example (see Figure~\ref{fig:lorenz96_profile}(c)).

We also demonstrate our framework with the Kuramoto--Sivashinsky equations\cite{Kuramoto1976, Sivashinsky1977}. We generated 20-dimensional time series: the values of the first 10 grids of the Kuramoto--Sivashinsky equations $x_0, x_1, ..., x_9$ and 10 random walks $x_{10}, x_{11}, ..., x_{19}$. We forecasted $x_0$ up to 10 steps using the same parameter values, data length, and noise scales as used in the Lorenz'96I example. See Supplementary Information for the detailed conditions.

Despite the increases in the variables, Figures~\ref{fig:kuramoto_rmse}(a) and (b) show that the proposed forecast yielded the minimum error among the existing frameworks, as in the Lorenz'96I example. Regarding the noise sensitivity, a lower noise level results in better performance for the proposed framework compared to the other frameworks. The proposed framework achieved the best performance for a wide data range and relatively small observational noise.

Finally, we tested the sensitivity to the number of variables. We changed the dimensions of Lorenz'96I to \{10, 20, 40, 80, 160, 320\}, and we substituted a half of each dataset with random walks. We set the training data length to 500 and evaluated with 500 samples.
For all ranges of variables, the proposed framework achieved the best performance (Figure~\ref{fig:lorenz96_variables}(a)).
We also tested the case where all variables are available and the datasets do not contain random walks; this is the case in which delay embedding is not needed because all system variables are observed.
In this case, the performance difference is negligible for all frameworks (except the single best embedding), and RDE achieved the best score for the case where the number of variables is 320.
Although this condition is too ideal because the observation function is an identity map, we can expect RDE to perform well for such cases that include many informative variables.

\begin{figure}[]
    \centering
    \includegraphics[width=1.\linewidth]{fig/lorenz96_rmse}
    \caption{Forecast performance with the Lorenz'96I model: comparisons of the performance (a) up to 10 steps ahead with the fixed data length (4000) without noise, (b) with different values of the data length, and (c) with different scales of observational noise. We computed the RMSE of the five-steps-ahead forecasts with randomly distributed embedding (RDE), multiview embedding (MVE), state-dependent weighting (SDW), single-best embedding based on the $(\mu+\lambda)$-ES algorithm (SBE), and the proposed framework. These tests were carried out with 20 datasets generated with different random initial conditions and noise. The median, upper quartile, and lower quartile are shown.}
    \label{fig:lorenz96_rmse}
\end{figure}

\begin{figure}[h]
    \centering
    \includegraphics[width=1.\linewidth]{fig/lorenz96_profile}
    \caption{Forecast profile of the Lorenz'96I model with the data length of 4000. Panel (a) shows the proportion of embedding of the proposed forecasts. The color indicates the proportion of embedding for each variable averaged over the number of combined forecasts for each step. Note that for each prediction step, the sum of the proportions of all variables is one. Variables $x_0,...,x_4$ are the variables of the 10-dimensional Lorenz'96I equations, and the other variables are random walks. Panel (b) shows the averaged filter coefficient $\rho$ and embedding dimension $E$ over the combined forecasts for each step. Note that the averaged dimensions can be a decimal fraction since they are averaged over the number of combined forecasts for each step. Panel (c) shows the relation between the number of combined forecasts and the errors. The solid line shows the in-sample error, and the dashed line shows the test error. The vertical solid line shows the selected number of combined forecasts $k_p$ to minimize the in-sample error.}
    \label{fig:lorenz96_profile}
\end{figure}

\begin{figure}[]
    \centering
    \includegraphics[width=1.\linewidth]{fig/kuramoto_rmse}
    \caption{Forecast performance with the Kuramoto--Sivashinsky equations: comparisons of the performance (a) up to 10 steps ahead with the fixed data length (4000) without noise, (b) with different values of the data length, and (c) with different scales of observational noise. We computed the RMSE of the five-steps-ahead forecasts with randomly distributed embedding (RDE), multiview embedding (MVE), state-dependent weighting (SDW), single-best embedding based on the $(\mu+\lambda)$-ES algorithm (SBE), and the proposed framework. These tests were carried out with 20 datasets generated with different random initial conditions and noise. The median, upper quartile, and lower quartile are shown.}
    \label{fig:kuramoto_rmse}
\end{figure}

\begin{figure}[]
    \centering
    \includegraphics[width=1.\linewidth]{fig/lorenz96_variables}
    \caption{Forecast performance with the Lorenz'96I model for various numbers of variables: cases where (a) a half of the variables are substituted with random walks and (b) all variables are available. We computed the RMSE of the five-steps-ahead forecasts with randomly distributed embedding (RDE), multiview embedding (MVE), state-dependent weighting (SDW), single-best embedding based on the $(\mu+\lambda)$-ES algorithm (SBE), and the proposed framework. These tests were carried out with 20 datasets generated with different random initial conditions and noise. The median, upper quartile, and lower quartile are shown.}
    \label{fig:lorenz96_variables}
\end{figure}

\subsubsection*{Flood dataset}
We tested the proposed framework with real-world data, namely the flood forecasting competition dataset at ``Artificial Neural Network Experiment (ANNEX 2005/2006)''\cite{Dawson2005a}. This dataset contains nine variables: the river stages of the target site and three upstream sites and five rain gauges for three periods. We forecasted 6, 12, 18, and 24 h ahead of the target river stage using all nine variables. See Supplementary Information for the detailed conditions.
We also compared our framework with four other conventional machine-learning methods used in an existing study\cite{Okuno2019}: a recurrent neural network with long short-term memory (LSTM)\cite{Hochreiter:1997:LSM:1246443.1246450}, support vector regression (SVR)\cite{Boser:1992:TAO:130385.130401}, and random forest regression\cite{Breiman2001}.

The results in Table~\ref{tab:flood} show that the proposed framework yielded the best root-mean-square error (RMSE) through 12--24 h ahead, and the differences compared to the other methods were expanded over the forecast horizon.
Although the 6 h (one step)-ahead forecast was not the best, the difference compared to the best score was negligible. The dataset contains only 1460 points, and long-term memory is not required for this task; hence, LSTM did not perform well.
Although the test data include the largest river stage on 1995-02-01, Figure~\ref{fig:flood}(a) shows that the proposed framework properly forecasted the inexperienced river stage.

We compared the combined forecast with its individual members obtained with the $(\mu+\lambda)$-ES algorithm. Regarding the individual members, the order of the in-sample score was different from the order of the test (see Figure~\ref{fig:flood}(b)). In contrast, the combined forecast achieved the minimum error by far for the both the in-sample and test errors. This forecast stability is one of the largest advantages of combining forecasts.

We also profiled $\rho$ and $E$ averaged over the $\hat{k}_p$ combined members for each step. The results in Figure~\ref{fig:flood}(c) suggest the same trends as the toy models, even for the real-world data. Specifically, smaller dimensions and smaller values of $\rho$ were selected for shorter steps.

\begin{figure}[htb]
    \begin{center}
        \includegraphics[width=.9\linewidth]{fig/flood}
        \caption{Forecast results for the flood dataset. Panel (a) shows a comparison of the ground truth and the proposed 24-h-ahead forecast. The proposed forecast did not underestimate the maximum river stage, which is the maximum value of the whole dataset. Panel (b) shows a comparison of the in-sample and test errors of the ensemble members. The results demonstrate the difficulty of selecting the best forecast because the best in-sample forecast does not always perform the best. In contrast, the combined forecast performed the best by far for both training and test data. Panel (c) shows the filter coefficient $\rho$ and embedding dimension $E$ averaged over the combined forecasts for each step.}
        \label{fig:flood}
    \end{center}
\end{figure}

\section*{Discussion}
Our numerical experiments suggest that our proposed framework works well for a wide range of the data length, including short time series such as 100. In contrast, the framework did not appear to work fine with very noisy data, and MVE achieved the best performance, as shown in the noise sensitivity analyses (Figures~\ref{fig:lorenz96_rmse}(c) and ~\ref{fig:kuramoto_rmse}(c)). This suggests that noise prevents the framework from selecting optimal embeddings on the basis of the squared errors, and it results in obtaining invalid embeddings. MVE can still be a good option for very noisy data.
Note that the performance of RDE was not satisfactory in some of our numerical experiments, but this is because the RDE framework was originally proposed for extremely high-dimensional data, e.g., hundreds of variables. As shown in Figure~\ref{fig:lorenz96_variables}(b), we can expect RDE to perform well for cases that include many informative variables, and delay embedding is not needed.

As shown in the numerical experiments, the proposed framework worked well for high-dimensional dynamics and datasets containing high-dimensional variables. On the other hand, it is not suitable to apply it to relatively low-dimensional data. We tested the framework with low-dimensional datasets, and the results suggest that the proposed framework does not always yield the best performance, especially for low-dimensional datasets such as the R\"ossler equations\cite{Rossler1976} (see Supplementary Information).
This is because it is easy to find suitable embeddings to forecast in a brute-force manner, especially for low-dimensional datasets, while it becomes combinatorially difficult as the number of variables increases.
In addition, the aim of $K$ times optimization is to yield diverse embeddings (see the effect of $K$ times optimization in Supplementary Information).
However, if a dataset contains a small number of variables, almost the exact solutions are obtained for every optimization.

The proposed method combines forecasts via the simple average, optimizing the number of forecasts. The simple average works fine for most cases because suboptimal embeddings can equally forecast well thanks to solving the combinatorial optimization problem, and the test performance order is not always the same as the corresponding training data, especially for real-world data (see Figure~\ref{fig:flood}(b)).
If suboptimal embeddings show quite different forecast performances, the weighted average based on the error distribution\cite{Ma2018a} may work better. We can also apply the expert advice framework\cite{Cesa-Bianchi2006a} if better forecasts are assumed to be changed over time.

We need to select the values of several parameters to apply the proposed framework, i.e., the number of split datasets $K$, the number of the suboptimal embeddings for each batch $M$, the criteria of the Hamming distance $\theta$, the number of filters $L$, and their components.
We can control the variations in forecasts with these parameters. In order to set $K$, we consider the number of samples in each divided dataset; empirically, the framework works well with more than 50 samples for each divided dataset. Samples that are too small prevent the framework from finding good embeddings, but samples that are too large reduce the diversity of the embeddings.
Although $M$ can be any large integer (up to the total number of individuals of the evolution strategy), $M\leq5$ is sufficient for most cases. Note that the value of $M$ that is too large is overkill and unnecessarily increases the computational cost.
The criteria for the Hamming distance $\theta$ directly controls the similarity of suboptimal embeddings. We set $\theta=3$ in this study, and it worked fine for all cases. Note that conventional MVE and RDE guarantee that the Hamming distance among all embeddings is larger than two (i.e., $\theta=2$) for a fixed embedding dimension.
In practice, appropriate filters depend on the prediction steps. If we predict shorter time periods, high-pass filters will work fine; if we predict longer time periods, low frequencies are needed to obtain fine forecasts, as the numerical experiments suggest. Although we treat high-pass filters in this study, low-pass filters may work for extremely noisy data.
We need to understand the effects of these parameters and set the appropriate values within reason, $P=KLM$, which is usually much smaller than the total number of embeddings evaluated in the first step.

Our numerical experiments suggest that the appropriate embeddings vary with the prediction steps. The results show that a smaller $E$ with a smaller $\rho$ is preferable for a short-term forecast, and a larger $E$ with a larger $\rho$ is preferable for a longer-term forecast. These results are consistent with those of an existing study\cite{Chayama2016a} and with the fact that the longer-step forecasts need more complex map from the current inputs to the forecasts. Therefore, more information is needed to express the complex map.

Although we applied the $(\mu+\lambda)$-ES algorithm to obtain the suboptimal embeddings, it is worth considering other methods to solve the combinatorial optimization problem. We can apply other evolution algorithms, e.g., simulated annealing\cite{Kirkpatrick83optimizationby} and ant colony optimization\cite{Dorigo:2004:ACO:975277}. Any algorithm may work as long as it can store the suboptimal embeddings through its optimization procedure. Testing and evaluating other algorithms are open problems.

Throughout this study, we used a variation of a local nonlinear prediction method to forecast the target variable. One of the most significant advantages of these local prediction methods is that we can compute in-sample forecasts very easily with a small computational cost; all we have to do is to exclude the current query from the nearest neighbors.
In addition, owing to recent progress in approximate neighboring search algorithms (e.g., Refs.~\citenum{Muja2014a, Fu2016a}), these algorithms are fast enough to apply the evolution algorithms. Moreover, it is also suitable to obtain suboptimal embeddings through the forecast errors because local nonlinear prediction methods are sensitive to the false nearest neighbors derived from invalid features.
However, we can apply other regression methods, e.g., Gaussian process regression\cite{Rasmussen2006}, as done in RDE\cite{Ma2018a}, for further improvements. Note that these methods incur a very large computational cost to yield in-sample forecasts, as mentioned in Ref.~\citenum{Ma2018a}. To reduce the computational cost, we can apply these methods to only the process in the second step and can approximate the in-sample error by the k-fold cross-validation error.

\section*{Methods}
\subsection*{Delay embedding}
Delay embedding is a method of reconstructing the original state space using the delay coordinates $v(t)$. Here, with the $n$-dimensional observed time series $y(t)=[y_1(t),y_2(t),...,y_n(t)]$, variable $\sigma_i \in \{1,2,...,n\}$, and time lag $\tau_i \in \{0,1,...,l-1\}$, $v(t) \in \mathbb{R}^E$ is defined as follows:
\begin{equation}
    v(t)=\left[ y_{\sigma_1}(t-\tau_1), y_{\sigma_2}(t-\tau_2),...,y_{\sigma_E}(t-\tau_E) \right],
\end{equation}
where $E$ is the embedding dimension, $\tau_j=0 \ \exists j \in \{1,2,...,n\}$, and no duplication is allowed for any element. According to embedding theorems\cite{Takens1981,Deyle2011}, an attractor reconstructed by $v(t)$ is an embedding with the appropriate $E$. Note that the number of possible reconstructions increase combinatorially.

\subsection*{Optimizing embedding by an evolution strategy}
One can optimize embedding to solve combinatorial optimization problems to treat embeddings as binary series\cite{Vitrano2001, Small2003}.
Here, we solve equation~(\ref{eq:fitness}) by a combinatorial optimization problem.
We treat embeddings as binary series, i.e., 1 meaning embedded and 0 meaning not embedded.
For example, with $n=3$ and $l=2$, the embedding $[y_1 (t),y_2 (t),y_3 (t-1)]$ can be expressed by the $nl=6$ binary code $e=[1,0,1,0,0,1]$: the first element 1 corresponds to $y_1 (t)$ to embed, the second element 0 corresponds to $y_1 (t-1)$ to embed, the third element 1 corresponds to $y_2 (t)$ to embed, and so on.
The optimization of a binary series is a typical combinatorial optimization problem.
In this paper, we apply the $(\mu+\lambda)$ evolution strategy\cite{Schwefel1981a}.
The parameter $\mu$ denotes the number of parents, $\lambda$ denotes the number of offspring, and $+$ denotes {\it plus selection}, which means that the next parents are selected from both the previous parents and offspring.
Precisely, the algorithm yields $\lambda$ new offspring from the top-performing $\mu$ individuals selected from the previous parents and offspring.
The algorithm is an elitist method since the algorithm retains the best individuals unless they are replaced by superior individuals. The algorithm tends to converge fast, but it is likely to become trapped at local optima. Therefore, the algorithm is suitable for our proposed framework.
However, it is worth considering the application of other existing methods\cite{Runge2015, Vlachos2009, Chen2010} to obtain suboptimal embeddings.

\subsection*{Forecasting through FIR filters}
We apply FIR filters to the original time series to obtain various embeddings and forecasts. We apply the following linear transformation for each variable:
\begin{equation}
    z_i(t) = \left[\sum_{k=0}^{N-1} h_i(k) y_i(t-k) - \mu_i \right] / \sigma_i,
\end{equation}
where $h_i$'s are filter coefficients with $N$ elements, $h_i(0) = 1 \ \forall i$, and $\mu_i$ and $\sigma_i$ are the coefficients for standardization. After forecasting $z_f(t)$, we restore $\hat{z}_f(t)$ for the original time series. For $p=1$, we restore $\hat{y}_f(t+p|t)$ using $\hat{z}_f(t+p|t)$ as follows:
\begin{equation}
    \hat{y}_f(t+p|t) = \sigma_f \hat{z}_f(t+p|t) + \mu_f - \sum_{k=1}^{N-1} h_f(k) y_f(t+p-k).
\end{equation}
For $p \geq 2$, we obtain the $p$-steps-ahead forecast $\hat{y}_f(t+p|t)$ by substituting forecasts for the unobserved $y_f(t \geq t+1)$ as follows:
\begin{equation} \label{eq:restore_filter2}
    \hat{y}_f(t+p|t) = \sigma_f \hat{z}_f(t+p|t) + \mu_f \left[-\sum_{k=1}^{p-1} (h_f(k) \hat{y}_f(t+p-k|t)) + \sum_{k=p}^{N-1} (h_f (k) y_f (t+p-k)) \right].
\end{equation}
We can restore $\hat{y}_f(t+p|t)$ for any $p$ to apply equation~(\ref{eq:restore_filter2}) recursively.

We can apply FIR filters such as the moving average filter to suppress noise. In contrast, an existing study\cite{Okuno2017a} focuses on the filtered delay coordinates constructed by blending the derivatives of time series and the original one to search the nearest neighbors. Although these are a type of high-pass FIR filters, such filters improve the prediction accuracy in less noisy time series.

\subsection*{Overall algorithm of the proposed framework}
We denote a multivariate time series as $y(t) \in \mathbb{R}^n$ and the target variable to be forecasted as $y_f$. The proposed forecasting framework is presented in Algorithm~\ref{alg:framework} in terms of the user-defined number of split training datasets $K$, the number of top-performing embeddings chosen for each split $M$, the linear filters, and the number of filters $L$.

\section*{Effect of multiple optimizations}
We solve multiple combinatorial optimization problems to divide the objective function into $K$ problems in the first step.
This procedure yields diverse embeddings to improve the performance of the combined forecast.
Here, we numerically show the effect of the multiple optimizations.

We profiled the effect of $K$ times optimization with an example of the 10-dimensional Lorenz'96I model with the data length of 4000.
We forecasted the data by two methods: one is the proposed $K$ optimizations, and the other is a single optimization to minimize the whole training in-sample error.
All other processes follow the proposed forecast framework.
Note that the computational time is almost equivalent with both methods.

As shown in Figures~\ref{fig:kdiv}(a) and (b), the procedure with $K$ optimizations enhanced the diversity of embeddings compared with the single optimization.
This result suggests that the $K$ times optimization found different but useful embeddings for each optimization process.
The $K$ times optimization also reduced the number of forecasts to combine (Figures~\ref{fig:kdiv}(c) and (d)).
This is because the single optimization process yielded similar embeddings, and the combination of similar forecasts yields a small performance improvement empirically\cite{Sollich96, Kuncheva2003}.
Reductions in the number of combined forecasts and the computational time are important, especially for real-time applications.

\begin{figure}[ht]
    \centering
    \includegraphics[width=1.\linewidth]{fig/effect_of_kdiv.pdf}
        \caption{Effect of multiple optimizations: the proportion of embedding of the proposed forecasts with (a) $K$ times optimization and (b) a single optimization to minimize the whole in-sample error. The color indicates the proportion of embedding for each variable averaged over the number of combined forecasts for each step. Panels (c) and (d) show the relation between the number of combined forecasts and the five-steps-ahead errors for (c) the $K$ times optimization and (d) the single optimization, respectively. The solid line shows the in-sample error, and the dashed line shows the test error. The vertical solid line shows the selected number of combined forecasts.}
        \label{fig:kdiv}
\end{figure}

\section*{Application to low-dimensional dynamics}
We applied the proposed forecast framework to low-dimensional datasets---namely, the Lorenz'63 equations\cite{Lorenz1963} (three variables), the R\"ossler equations\cite{Rossler1976} (three variables), and the six-dimensional Lorenz'96I equations\cite{Lorenz1995a} (six variables).
We set the data length to 4000 as the database and evaluated forecasts up to 10 steps ahead with 500 samples for all the cases. See the next section for the detailed conditions.

As shown in Figure~\ref{fig:lowdims}, the proposed framework did not always achieve the best performance. There is little difference between the proposed framework and state-dependent weighting (SDW) for the Lorenz'63 dataset (Figure~\ref{fig:lowdims}(a)), and SDW achieved the best performance for the R\"ossler dataset (Figure~\ref{fig:lowdims}(b)).
In contrast, when we increased the number of variables to six, the proposed framework yielded the best performance (Figure~\ref{fig:lowdims}(c)), and as stated in the main text, the proposed framework achieved much better results than the others for cases with 10 or 20 variables.

Interestingly, the performance of SDW was superior to that of MVE for all of the low-dimensional toy models mentioned above, although the performance of SDW for the high-dimensional dataset was not satisfactory, as stated in the main text. This result means that weighting based on the forecast performance and its corresponding state\cite{Okuno2019} worked fine with simple low-dimensional data.
On the other hand, it is difficult to improve the performance by weighting for high-dimensional data, and a simple average is sufficient for such cases.
In short, it is preferable to apply the proposed framework to complex high-dimensional and less noisy data, but it is worth considering the application of weighting methods such as SDW to simple low-dimensional data.

\begin{figure}[ht]
    \centering
    \includegraphics[width=1.\linewidth]{fig/lowdims}
    \caption{Forecast performance for low-dimensional datasets: RSMEs for (a) the Lorenz'63 dataset, (b) the R\"ossler dataset, and (c) the six-dimensional Lorenz'96I dataset. We compared the performance of multiview embedding (MVE), state-dependent weighting (SDW), and single-best embedding based on the $(\mu+\lambda)$-ES algorithm (SBE). These tests were carried out with 20 datasets generated with different random initial conditions and noise. The median, upper quartile, and lower quartile are shown.}
    \label{fig:lowdims}
\end{figure}

\section*{Detailed conditions of the numerical experiments}
\subsection*{Lorenz'96I equations}
The Lorenz'96I equations\cite{Lorenz1995a} are expressed as 10-dimensional differential equations as follows:
\begin{equation}
    \frac{d x_i}{d t} = x_{i-1} (x_{i+1}-x_{i-2}) - x_i + F,
    i=0,1,2,...,9,
\end{equation}
where $F$ is a forcing variable and $i$ is cyclic. We generated 10-dimensional time series: Lorenz'96I's $x_0, x_1, x_2, x_3$, and $x_4$ and random walks $x_5, x_6, x_7, x_8$, and $x_9$. We chose the initial condition from a normally distributed random value for each variable. We set the integration step to 0.001 and recorded every 50 points. Note that the initial transients were disregarded. We forecasted $x_0$ up to 10 steps. We set $K=10, M=3, \theta=3$, and $h_i(0)=1, h_i(1)=\rho \ \forall i$ for $\rho \in \{0.0, -0.2, -0.4, -0.6, -0.8, -1.0\}$. For the $(\mu+\lambda)$-ES algorithm, we set $\mu=50, \ \lambda=100$, the number of generations to 10, and the number of populations to 100.

We compared the proposed framework with existing frameworks---namely, randomly distributed embedding with an aggregation scheme (RDE), multiview embedding (MVE), SDW with $E=4$, and up to four lags for MVE and SDW. Because the possible number of embeddings combinatorially increases, we randomly generated 1000 embeddings instead of using brute-force calculation.
Note that the total number of embedding evaluations is almost the same as that for the proposed framework.
We also computed the single-best embedding via the $(\mu+\lambda)$-ES algorithm to minimize the total in-sample error within the whole training dataset.

\subsection*{Kuramoto--Sivashinsky equations}
The Kuramoto--Sivashinsky equations\cite{Kuramoto1976, Sivashinsky1977} are expressed as follows:
\begin{equation}
    \frac{\partial y}{\partial t} = -\frac{\partial^2 y}{\partial x^2} - \frac{\partial^4 y}{\partial x^4} - u \frac{\partial y}{\partial x}.
\end{equation}
We computed $y(x,t)$ with spatially periodic boundary conditions in the interval $[0, L]$, where $L=22$ and the number of uniform grids is 128. We generated 20-dimensional time series: the values of the first 10 grids of the Kuramoto--Sivashinsky equations $x_0, x_1, ..., x_9$ with a sampling time of 1.0 and 10 random walks $x_{10}, x_{11}, ..., x_{19}$.
We set $E=5$ and considered up to five lags for the existing schemes (RDE, MVE, and SDW).
We forecasted $x_0$ up to 10 steps using the same parameter values, data length, and noise scales as those in the Lorenz'96I example.

\subsection*{Flood dataset}
The flood forecasting competition dataset ``Artificial Neural Network Experiment (ANNEX 2005/2006)''\cite{Dawson2005a} contains river stage and rainfall data for three periods: 1993-10-01 to 1994-03-31, 1994-10-01 to 1995-03-31, and 1995-10-01 to 1996-03-31. We forecasted river stage $Q$ 6, 12, 18, and 24 h ahead using nine variables: the river stages of the target ($Q$) and three upstream sites ($US1$, $US2$, and $US3$) and five rain gauges ($RG1$, $RG2$, $RG3$, $RG4$, and $RG5$). All data were sampled by 6 h. We used the period of 1994-10-01 to 1995-03-31 for testing and the others for training. We set $K=6, M=3, \theta=3$, and $h_i(0)=1, h_i(1)=\rho \ \forall i$ for $\rho \in \{0.0, -1.0\}$. For the $(\mu+\lambda)$-ES algorithm, we set $\mu=50, \lambda=100$, the number of generations to 20, and the number of populations to 100.

\subsection*{Lorenz'63 equations}
The Lorenz'63 equations\cite{Lorenz1963} are expressed as three-dimensional differential equations as follows:
\begin{equation}
    \frac{dx}{dt} = p (y - x),
    \frac{dy}{dt} = x (r - z) - y,
    \frac{dz}{dt} = xy - bz.
\end{equation}
We generated three-dimensional time series: Lorenz'63's $x$ and $y$ and a random-walk series. We set the integration step to 0.001 and recorded every 100 points with $p=10,b=8/3,$ and $r=28$. Note that the initial transients were disregarded. We forecasted $x$ up to 10 steps. We set $K=10, M=3, \theta=3$, and $h_i(0)=1, h_i(1)=\rho \ \forall i$ for $\rho \in \{0.0, -0.2, -0.4, -0.6, -0.8, -1.0\}$. For the $(\mu+\lambda)$-ES algorithm, we set $\mu=50, \lambda=100$, the number of generations to 10, and the number of populations to 100.

We compared the proposed framework with MVE and SDW with $E=4$ up to five lags and the single-best embedding via the $(\mu+\lambda)$-ES algorithm to minimize the total in-sample error within the whole training dataset. Note that we did not carry out calculations with RDE because we were not able to prepare a sufficient number of ``nondelay embeddings'' for this type of low-dimensional data.

\subsection*{R\"ossler equations}
The R\"ossler equations\cite{Rossler1976} are expressed as three-dimensional differential equations as follows:
\begin{equation}
    \frac{dx}{dt} = -y - z,
    \frac{dy}{dt} = x + ay,
    \frac{dz}{dt} = b + z (x-c).
\end{equation}
We generated three-dimensional time series: R\"ossler $x$ and $y$ and a random-walk series. We set the integration step to 0.001 and recorded every 500 points with $a=0.36,b=0.4$, and $c=4.5$. Note that the initial transients were disregarded. We forecasted $x$ up to 10 steps with the same conditions as those of the Lorenz'63 equations.

\subsection*{Six-dimensional Lorenz'96I equations}
The six-dimensional Lorenz'96I dataset contains the six-dimensional Lorenz'96I series $x_0, x_1$, and $x_2$ and random walks $x_3, x_4$, and $x_5$. The conditions for numerical integration were the same as those of the 10-dimensional Lorenz'96I equations. We forecasted $x_0$ up to 10 steps with the same conditions as those of the Lorenz'63 equations.

\section*{Method of analogues}
We applied a variation of the method of analogues\cite{Lorenz1969} to obtain a $p$-steps-ahead forecast $\hat{y}_f (t+p|t)$ at time $t$. The method of analogues takes the forward paths of neighboring trajectories as the forecast. Here, we forecast $y_f (t+p|t)$ from the set of delay coordinates $\{v(t) \mid t \in \mathcal{T}_{train}\}$.
We first search the database to find the neighboring points of $v(t)$. Then, the method of analogues gives $\hat{v}(t+p|t)$ in terms of the set of neighboring time indices $\mathcal{I}(t)$ as follows:
\begin{equation}
    \hat{v}(t+p|t) = \sum_{t' \in \mathcal{I} (t)} \lambda(t') v(t'+p),
\end{equation}
where $\lambda(t') \in \mathbb{R}_{+}$ is a weight satisfying $\sum_{t' \in \mathcal{I}(t)} \lambda (t)$. We employ $\lambda (t') \propto \|v (t') - v (t)\|^2_2$ throughout this paper.

Note that we can apply more advanced methods such as that in Refs.~\citenum{Farmer1987, Hirata2014} or any other regression methods with our proposed forecasting framework.