# MTSMAE: Masked Autoencoders for Multivariate Time-Series Forecasting

Peiwang Tang1,2, Xianchao Zhang3,4∗ 1 Institute of Advanced Technology,
University of Science and Technology of China, China  2 G60 STI Valley
Industry & Innovation Institute, Jiaxing University, China  3 Key Laboratory
of Medical Electronics and Digital Health of Zhejiang Province,
Jiaxing University, China  4 Engineering Research Center of Intelligent Human
Health Situation Awareness of Zhejiang Province,
Jiaxing University, China<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Large-scale self-supervised pre-training Transformer architecture have
significantly boosted the performance for various tasks in natural language
processing (NLP) and computer vision (CV). However, there is a lack of
researches on processing multivariate time-series by pre-trained Transformer,
and especially, current study on masking time-series for self-supervised
learning is still a gap. Different from language and image processing, the
information density of time-series increases the difficulty of research. The
challenge goes further with the invalidity of the previous patch embedding and
mask methods. In this paper, according to the data characteristics of
multivariate time-series, a patch embedding method is proposed, and we present
an self-supervised pre-training approach based on Masked Autoencoders (MAE),
called MTSMAE, which can improve the performance significantly over supervised
learning without pre-training. Evaluating our method on several common
multivariate time-series datasets from different fields and with different
characteristics, experiment results demonstrate that the performance of our
method is significantly better than the best method currently available.

###### Index Terms:

Autoencoder, Pre-Training, Time-Series, Forecasting

(a) Pre-Training

(b) Fine-Tuning

Figure 1: Our MTSMAE architecture. In the pre-training, our model consists of
encoder and Rdeocder (the decoder responsible for recovering the original
input); in the fine-tuning, our model consists of encoder and Pdeocder (the
decoder responsible for predicting future data).

##  I Introduction

With the fast evolvement of deep learning in recent years [1, 2, 3], training
a model is anticipated to accommodate hundreds of millions of labeled data
[4]. The demand for large scale data processing has been solved by self-
supervised pre-training in natural language processing (NLP) and computer
vision (CV) fields [6, 5]. Most of these solutions are based on masked
modeling, such as masked language modeling in NLP [7, 8] or masked image
modeling in CV [9, 10, 11]. Their ideas are conceptually simple: firstly mask
parts of the data based on the original data and then enable to recovery these
parts by learning [12, 13]. Masked modeling encourages the model to infer the
deleted parts according to the context information, so that the model can
learn the deep semantics, which has become the benchmark of self-supervised
pre-training in NLP and CV fields [6, 5]. These pre-trained masked modeling
has been proved to be well applied to various downstream tasks, among which a
simpler and more effective way is masked autoencoders (MAE) [5]. However,
despite a widely interests in this idea from academia and industry following
the success of MAE, the progress of autoencoder methods in the field of
multivariate time-series data (MTSD) lags behind other fields.

One of the main reasons is that the information density of MTSD is different
from that of CV and NLP. The local information of MTSD seems to be heavy
spatial redundancy, but the multivariate information within each time point
has high specificity. Missing information can be easily learned from
information at adjacent time points with little high-level understanding. In
order to overcome this difference and encourage learning more useful features,
we use the idea of Vision Transformer (ViT) [4], patch MTSD, and mask more
random patches than the original MAE, e.g. 85$\%$. This simple strategy has a
ideal performance in MSTD, which can reduce redundancy effectively and futher
increase the overall understanding beyond low-level information of the model.

Another reason is the design of the decoder. The decoder of the autoencoder
maps the latent representation back to the input. In CV, the decoder can
reconstruct the pixel level representation in patch. In NLP, the decoder
predicts missing words. In MTSD, the decoder recovers completely different
data with information specificity and a dimension that even highly to 321. For
different types of MTSD, the dimensions of data at each time point are
different, which may have 321 dimensions or only 7 dimensions. We find that
the design of decoder plays a key role in the latent representation of
learning for MTSD.

Based on the above analysis, we propose a very simple and effective MTSMAE for
MTSD representation learning. Our MTSMAE idea is simple: In the pre-training,
patch MTSD, masks random patches from the input and recover the missing
patches; in the fine-tuning, take out the encoder trained in the previous
step, and the input of the decoder is redesigned. In addition, our encoder
only calculates visible patches. Unlike the decoder with only one layer of MLP
in Bert [6], we design different levels of decoder according to different
MTSD, but compared with MAE, our decoder are all lightweight. We have
conducted extensive experiments on four different datasets of three types. The
final experimental results show that our proposed MTSMAE can significantly
improve the accuracy of prediction, and is superior to other state-of-the-art
models.

##  II Related work

###  II-A Masked modeling for Self-supervised learning

Self-supervised learning approaches have aroused great interest in natural
language processing and computer vision, often focusing on different pretext
tasks for pre-training [14, 15, 16]. BERT [6] and GPT [7, 8, 17] are very
successful masked modeling for pre-training in NLP. They learn representations
from the original input corrupted by masking, which have been proved by a
large amount of evidence to be highly extensible [17], and these pre-trained
representations can be well extended to various downstream tasks. As the
original method applied BERT to CV, BEIT [10] firstly “tokenize” the original
image as a visual token and then randomly mask some image patches and input
them into the backbone Transformer [2]. The goal of pre-training is to restore
the original visual token from the damaged image patch. Based masked modeling
as all autoencoder [18], MAE [5] uses the encoder to map the observed signal
to the potential representation, and the decoder to reconstruct the original
signal from the latent representation. In turn, different from the classic
autoencoder, MAE adopts an asymmetric design to allow the encoder to operate
on only part of the observed signal (without mask token) and use a lightweight
decoder to reconstruct the complete signal from the latent representation
along with mask tokens.

###  II-B Models for Time-Series Forecasting

Forecasting is one of the most important applications of time-series. LSTNet
[19] uses convolutional neural network (CNN) [1] and recurrent neural network
(RNN) to extract short-term local dependence patterns between variables and
find long-term patterns of time-series trends. In addition, the traditional
autoregressive model is used to solve the scale insensitivity problem of
neural network model. The temporal convolutional network (TCN) [20] is
proposed to make the CNN have time-series characteristics, as a variant of CNN
that deals with sequence modeling tasks. It mixes RNN and CNN architecture to
use the causal convolution to simulate temporal causality. Reformer [21] is
proposed to replace the original dot-product attention with a new one using
locality-sensitive hashing. It decreases the complexity from
$\mathcal{O}(L_{2})$ to $\mathcal{O}(L\log{L})$ and makes the storage
activated only once in the training process rather than $n$ times (here $n$
refers to the number of layers) by using reversible residual layers to replace
standard residuals. LogTrans [22] proposes a convolutional self-attention
mechanism, which uses causal convolution to process the local context of the
sequence and calculate the query / key of self-attention. Further more, it
uses the logsparse Transformer architecture to ensure that each cell can
receive signals from other cells in the sequence data, while reducing the time
complexity of the architecture. Informer [23] is proposed ProbSparse self-
attention and encoder with self-attention distilling. The former is based on
query and key similarity sampling dot-product pairs, which reduces the
computational complexity of Transformer and allows it to accept longer input.
The latter adopts the concept of distillation to design encoder, so that the
model can continuously extract the feature vectors that focus on, while
reducing the memory occupation.

Figure 2: Patch embedding method.

##  III Methodology

The problem of multivariate time-series forecasting is to input the past
sequence
$\mathcal{X}^{t}=\left\\{x^{t}_{1},\cdots,x^{t}_{L_{x}}|x^{t}_{i}\in\mathbb{R}^{d_{x}}\right\\}$
at time $t$, and output the predict the corresponding future sequence
$\mathcal{Y}^{t}=\left\\{y^{t}_{1},\cdots,y^{t}_{L_{y}}|y^{t}_{i}\in\mathbb{R}^{d_{y}}\right\\}$,
where $L_{x}$ and $L_{y}$ are the lengths of input and output sequences
respectively, and $d_{x}$ and $d_{y}$ are the feature dimensions of input
$\mathcal{X}$ and output $\mathcal{Y}$ respectively.

Our masked autoencoders (MTSMAE) is a simple autoencoding method and the
training process is divided into two stages, as shown in the Fig. 1. As all
autoencoders, there is an encoder and a decoder in our method. The encoder
maps the observed signal to a latent representation and the decoder
reconstructs the original signal from the latent representation in the pre-
training, or output $\mathcal{Y}$ in the fine-tuning. In the pre-training, the
encoder only operates the partial signal (without mask token) observed.

###  III-A Patch embedding

As shown in the Fig. 2, the input embedding consists of three parts, a scalar
projection SP (we project the scalar context $x^{t}_{i}$ into $d_{model}$-dim
vector with 1-D convolutional filters. The kernel width of one-dimensional
convolutional filter is 3, and stride is 1), a local position embedding PE [2]
and a global time stamp embedding SE [23]:

$\mathcal{X}=SP+PE+SE$ (1)

which

$SP=Conv1d(x^{t}_{i})$ (2)
$PE_{\left(i,2j\right)}=sin\left(i/10000^{2j/d_{model}}\right)$ (3)
$PE_{\left(i,2j+1\right)}=cos\left(i/10000^{2j/d_{model}}\right)$ (4)
$SE=E(month)+E(day)+E(hour)+E(minute)$ (5)

where $i\in\left\\{1,\dots,L_{x}\right\\}$, $j\in\left\\{1,\dots,\left\lfloor
d_{model}/2\right\rfloor\right\\}$, $d_{model}$ is the feature dimension after
embedding, $E$ is a learnable stamp embeddings with limited vocab size (up to
60, namely taking minutes as the finest granularity).

The original MAE continues the idea of ViT, processing image data
$\mathcal{X}\in\mathbb{R}^{H\times W\times C}$, where $(H,W)$ is the
resolution of the original image, $C$ is the number of channels. For MTSD,
$\mathcal{X}\in\mathbb{R}^{L_{x}\times d_{x}}$, the original patch embedding
method is no longer applicable. Therefore, unlike the method of patch image
data in ViT, we patch MTSD in the direction of time after embedding:

$\mathcal{X}_{h}=Conv1d(\mathcal{X})\in\mathbb{R}^{L_{x}/P\times d_{model}}$
(6) $\mathcal{X}_{pte}=Conv1d(\mathcal{X}_{h})\in\mathbb{R}^{L_{x}/P^{2}\times
d_{model}}$ (7)

where the kernel width of one-dimensional convolutional filter and $stride=P$,
$\mathcal{X}_{pte}$ is the final result of the patch embedding. We use two-
step one-dimensional convolution to control the resolution of each MSTD patch
to $(P,d_{model})$, so the length of the last input sequence $L=L_{x}/P^{2}$.

###  III-B Model inputs

In the pre-training, the input of the encoder is the output
$\mathcal{X}_{pte}$ of the past sequence $\mathcal{X}^{t}$ after patch
embedding. However, only unmasked patches are input. A random masking method
is adopt, that is, the patches are randomly sampled without replacement, and
follow the uniform distribution. The random sampling can tremendously remove
the information redundancy of MSTD by deleting a large number of patches
(i.e., high masking rate). Thus, a task that can not easily recover
information from the visible nearby patches is created. Finally, the highly
sparse input provides an opportunity to train an efficient encoder. In the
pre-training, the decoder only inputs the output of the encoder.

In the fine-tuning, the input of the encoder is all the output
$\mathcal{X}_{pte}$ of the past sequence $\mathcal{X}^{t}$ after patch
embedding. As shown in the Fig. 1(b), the input of the decoder includes not
only the output of the decoder, but also the label part
$\mathcal{X}^{t}_{label}$ and the forecast part $\mathcal{X}_{0}$, $nonpatch$
$embedding$ is a common embedding without convolution kernel processing, such
as the embedding in Informer. It can be formulated as follows:

$\mathcal{X}_{eni}=Patch(\mathcal{X}^{t})$ (8)
$\mathcal{X}_{dei}=(\mathcal{X}_{eno},Concat(Patch(\mathcal{X}^{t}_{label}),NonPatch(\mathcal{X}_{0})))$
(9)

where $\mathcal{X}_{0}$ denotes the placeholder filled with zero,
$\mathcal{X}_{eno}$ indicates the output of the encoder, $Patch$ refers to
patch embedding and $NonPatch$ refers to $nonpatch$ $embedding$.

###  III-C Encoder and Decoder

Our encoder is the encoder of transformer. In the pre-training, our encoder
embeds only visible, unmasked patches through patch embedding, and then
processes the output data through a series of transformer encoder blocks. Our
encoder only operates on a small part of the whole set, e.g., only 15$\%$,
which can greatly reduce the redundancy of information and increase the
overall understanding of the model beyond low-level information. In the fine-
tuning, our encoder can see all the patches. Except for this, it is no
different from the encoder in the pre-training. The Transformer encoder layers
are composed of two sub-blocks. The first is a multi-head self-attention
mechanism (MSA), and the second is a simple, position-wise fully connected
feed-forward network (MLP). Residual connections [30] are used around each of
the two sub-blocks , and layer normalization (LN) [31] is then performed.

$\textbf{z}^{\prime}_{n}=LN(\textbf{z}_{n-1}+MSA(\textbf{z}_{n-1}))$ (10)
$\textbf{z}_{n}=LN(\textbf{z}^{\prime}_{n}+MLP(\textbf{z}^{\prime}_{n}))$ (11)

where $n=1\dots N$, $z_{0}=\mathcal{X}_{eni}$, $\mathcal{X}_{eno}=z_{n}$, $N$
is the number of layers of the encoder, and MLP consists of two linear
transformations with a ReLU activation in between.

$MLP(x)=max(0,xW_{1}+b_{1})W_{2}+b_{2}$ (12)

$W_{1}$, $W_{2}$, $b_{1}$, $b_{2}$ are all parameters to be learned, and the
parameters of each MLP are different.

In the pre-training, we set up the decoder as MAE. The input to the decoder is
a complete set including the visible patches output by encoder and mask
tokens, where as the vector of learning, masked tokens are the data to be
recovered. Except different inputs, there is no difference between the decoder
and encoder in the pre-training. In the fine-tuning, our decoder is also the
decoder of the transformer. In addition to the two sub-blocks in each encoder
layer, the decoder also inserts a third sub-block in the two sub-blocks, which
performs multi-head attention on the output of the encoder layers. In order to
prevent positions from attending to subsequent positions, MSA is modified to
masked multi-head self-attention (MMSA). This masking ensures that the
prediction of position $i$ can only rely on the known outputs of positions
less than $i$.

$\textbf{z}^{\prime}_{n}=LN(\textbf{z}_{n-1}+MMSA(\textbf{z}_{n-1}))$ (13)
$\textbf{z}^{\prime\prime}_{n}=LN(\textbf{z}^{\prime}_{n}+MSA(\textbf{z}^{\prime}_{n},\mathcal{X}_{eno}))$
(14)
$\textbf{z}_{n}=LN(\textbf{z}^{\prime\prime}_{n}+MLP(\textbf{z}^{\prime\prime}_{n}))$
(15)

In the pre-training, our MTSMAE reconstructs the input by recovering the
specific value of each masking patch. Each element output by the decoder is a
vector that can represent a patch. The last layer of the decoder is a linear
projection, whose output channel is $P\times D$, $P$ is the length of the
patch, and $D$ is the dimension of the time-series. In the fine-tuning, each
element output by the decoder represents the data $y^{t}_{i}$, and the output
channel of the last layer of linear projection is $D$.

Our loss function is calculated by the mean square error (MSE) between the
model output data $y_{o}$ (recovery, prediction) and the real data $y$. But in
the pre-training, we only calculate the loss on the masking patch, similar to
Bert, MAE.

$\textbf{Loss}(y_{o},y)=\frac{1}{L}\sum_{1\leq i\leq
L}(y\widehat{}_{i}-y_{i})^{2}$ (16)

Where $L$ represents the number of masking patch or the length of prediction
sequence.

##  IV Experiment

###  IV-A Datasets

We conducted extensive experiments on four datasets of three types to evaluate
the progressiveness of the proposed MTSMAE.

####  IV-A1 ETT(Electricity Transformer Temperature)

The ETT is a key indicator for long-term deployment of the electric power,
which collects data from two counties in China from July 2016 to July 2018 for
a total of two years. ETT is divided into two types, and each of which has two
different datasets including 1-hour-level {ETTh1, ETTh2} and 15-minute-level
{ETTm1, ETTm2}. The data at each time point is composed of the target value
“oil temperature” and six power load features. The size of training set,
validation set and test set are 12/4/4 months respectively. We choose ETTm1
and ETTh2 as the experimental dataset.

####  IV-A2 ECL (Electricity Consuming Load)

It collects the electricity consumption (Kwh) of 321 customers, that is, each
time point includes the respective electricity consumption of them. Due to the
lack of data [22], we follow the settings in Informer [23] to convert the
dataset to hourly consumption for 2 years. The size of training set,
validation set and test set are 15/3/4 months respectively.

####  IV-A3 WTH(Weather)

This dataset contains the local climatologicale data of nearly 1600 locations
in the United States. The data are collected by once an hour from 2010 to
2013. The data at each time point includes the target value “wet bulb” and 11
climate features. The size of training set, validation set and test set are
28/10/10 months respectively.

###  IV-B Experimental Details

(a) MTSMAE (b) Informer (c) LogTrans (d) LSTM

(e) Transformer (f) Reformer

Figure 3: Prediction cases from the ECL dataset under the $L_{y}=96$ setting.
Blue lines are the ground truth and orange lines are the model prediction.

####  IV-B1 Baselines

We selected five methods as comparison, including Transformer [2], three
latest state-of-the-art transformer-based models: Reformer [21], LogTrans
[22], Informer [23], and one RNN-based models: LSTM (Long Short-Term Memory
networks) [24].

####  IV-B2 Experiment setting

We fixed the random seed to ensure the repeatability of the experiment. Our
experiment was implemented in Pytoch [29], and all the experiments are
conducted on a single Nvidia RTX 3090 GPU (24GB memory). The input of each
dataset is zero-mean normalized. We use two evaluation metrics, including mean
square error (MSE) and mean absolute error (MAE):

$MSE=\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{d}\frac{(y-\hat{y})^{2}}{d}$ (17)
$MAE=\frac{1}{n}\sum_{i=1}^{n}\sum_{j=1}^{d}\frac{\left|y-\hat{y}\right|}{d}$
(18)

Where $n$ is the length of the sequence and $d$ is the dimension of data at
each time point of MTSD. We use these two evaluation metrics on each
prediction window to calculate the average of multivariate forecasts and roll
the whole set with $stride=1$.

TABLE I: Pre-training setting.

config | value
---|---
optimizer | AdamW [25]
base learning rate | 1e-3
actual learning rate | 2.5e-4
weight decay | 0.05
optimizer momentum |  $\beta_{1}$, $\beta_{2}$=0.9, 0.95 [9]
batch size | 64
learning rate schedule | cosine decay [27]
loss function | MSE
input length | 784
patch size | 4
encoder layers | 3
decoder layers | 1
masking ratio | 85%

The pre-training default setting is in Table I. The input sequence length in
the model is 784, and the patch size is 4, that is, we have a total of 196
patches. Due to the limitation of GPU memory, we changed the batch size to 64,
and then we chose the Adam optimizer for optimization with the same strategy
as MAE, that is, the configuration of the optimizer has not been modified. We
follow the setting of the original MAE and use the linear learning rate
scaling rule [28]:$lr=base\\_lr\times batchsize/256$. The number of encoder
layers is 3 and the number of decoder layers is 1. In pre-training, we choose
a high proportion of masking ratio by about 85%. These configurations are the
same on all datasets without making corresponding changes according to
different datasets.

TABLE II: Fine-tuning setting

config | value
---|---
optimizer | Adam [26]
initial learning rate | 1e-4
optimizer momentum |  $\beta_{1}$, $\beta_{2}$=0.9, 0.999
batch size | 32
early stopping | 3
learning rate schedule | exponential decay
loss function | MSE

TABLE III: MSTD forecasting results on four datasets Methods | MTSMAE | PatchTrans | Informer | Reformer | LogTrans | Transformer | LSTM
---|---|---|---|---|---|---|---
Metric | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE | MSE | MAE
ECL | 24 | 0.237 | 0.344 | 0.243 | 0.347 | 0.264 | 0.362 | 0.294 | 0.393 | 0.223 | 0.332 | 0.242 | 0.350 | 0.341 | 0.417
48 | 0.257 | 0.357 | 0.278 | 0.369 | 0.299 | 0.389 | 0.304 | 0.397 | 0.291 | 0.374 | 0.261 | 0.359 | 0.335 | 0.412
96 | 0.263 | 0.362 | 0.272 | 0.370 | 0.302 | 0.388 | 0.293 | 0.382 | 0.299 | 0.379 | 0.286 | 0.380 | 0.329 | 0.408
192 | 0.278 | 0.371 | 0.328 | 0.415 | 0.322 | 0.379 | 0.349 | 0.426 | 0.279 | 0.370 | - | 0.326 | 0.407
384 | 0.289 | 0.382 | 0.368 | 0.441 | 0.338 | 0.411 | - | 0.331 | 0.404 | - | 0.319 | 0.402
ETTm1 | 24 | 0.282 | 0.356 | 0.311 | 0.379 | 0.362 | 0.407 | 0.397 | 0.423 | 0.905 | 0.656 | 0.339 | 0.392 | 1.253 | 0.837
48 | 0.412 | 0.445 | 0.455 | 0.477 | 0.461 | 0.458 | 0.558 | 0.535 | 0.880 | 0.662 | 0.582 | 0.579 | 1.262 | 0.838
96 | 0.558 | 0.544 | 0.657 | 0.604 | 0.618 | 0.576 | 0.703 | 0.611 | 0.562 | 0.553 | 0.855 | 0.752 | 1.197 | 0.829
192 | 0.603 | 0.597 | 0.531 | 0.546 | 0.765 | 0.659 | 0.867 | 0.678 | 0.564 | 0.527 | - | 1.200 | 0.824
384 | 0.672 | 0.634 | 0.818 | 0.708 | 0.965 | 0.763 | - | 0.694 | 0.612 | - | 1.199 | 0.822
ETTh2 | 24 | 0.648 | 0.661 | 0.689 | 0.686 | 1.567 | 1.022 | 0.893 | 0.764 | 0.851 | 0.725 | 0.968 | 0.802 | 3.239 | 1.371
48 | 1.126 | 0.874 | 1.266 | 0.948 | 3.326 | 1.500 | 1.240 | 0.940 | 1.165 | 0.904 | 2.123 | 1.218 | 3.359 | 1.403
96 | 2.561 | 1.304 | 3.167 | 1.555 | 7.756 | 2.321 | 4.161 | 1.637 | 2.795 | 1.431 | 6.597 | 2.163 | 3.309 | 1.391
192 | 3.341 | 1.588 | 4.687 | 1.800 | 5.527 | 1.993 | 4.425 | 1.791 | 5.135 | 1.946 | - | 3.417 | 1.415
384 | 3.107 | 1.403 | 3.375 | 1.626 | 3.581 | 1.587 | - | 3.156 | 1.561 | - | 3.456 | 1.422
WTH | 24 | 0.301 | 0.360 | 0.305 | 0.367 | 0.366 | 0.421 | 0.306 | 0.360 | 0.387 | 0.436 | 0.348 | 0.390 | 0.941 | 0.771
48 | 0.377 | 0.414 | 0.389 | 0.423 | 0.478 | 0.490 | 0.391 | 0.428 | 0.501 | 0.519 | 0.449 | 0.459 | 0.946 | 0.773
96 | 0.496 | 0.501 | 0.513 | 0.514 | 0.613 | 0.583 | 0.469 | 0.480 | 0.551 | 0.541 | 0.536 | 0.517 | 0.943 | 0.772
192 | 0.546 | 0.534 | 0.589 | 0.557 | 0.665 | 0.607 | 0.564 | 0.531 | 0.577 | 0.558 | - | 0.945 | 0.772
384 | 0.573 | 0.563 | 0.577 | 0.564 | 0.660 | 0.612 | - | 0.635 | 0.591 | - | 0.939 | 0.769
Count | 31 | 1 | 0 | 2 | 5 | 0 | 1

Our fine-tuning follows common practice of Informer training, and the default
setting is in Table II. In the fine-tuning, we use Adam optimizer for
optimization with a learning rate starts from $1e^{-4}$, decaying two times
smaller every epoch. There is no limit to the total number of epochs, with
appropriate early stopping, that is, when the loss of the validation set does
not decrease on the three epochs, the training will be stopped. We set the
comparison methods according to the suggestion, and the batch size is 32.

###  IV-C Results and Analysis

####  IV-C1 Multivariate Time-series Forecasting

In order to be consistent with the setting of pre-training, we choose the same
input length, i.e. $L_{x}=784$, and we gradually extend the size of the
prediction window $L_{y}$ as a higher requirement of prediction capacity,
i.e., {24, 48, 96, 192, 384}, representing {6h, 12h, 24h, 48h, 96h} in ETTm1,
{1d, 2d, 4d, 8d, 16d} in {ETTh2, ECL, WTH}, and we set the length of the label
to double $L_{y}$, i.e., {48, 96, 192, 384, 768}. The parameter sensitivity of
input length and output length will be discussed in subsequent chapters. In
the pre-training, the default configuration selected by our MTSMAE is shown in
the Table I,

Table III summarizes the multivariable evaluation results of all metrics on
{ECL, ETTm1, ETTh2, WTH}. The best results are highlighted in bold, where ‘-’
represents out-of-memory, and PatchTrans indicates our model without pre-
training. From Table III, we observe that the proposed model MTSMAE has better
performance than other methods. Compared with PatchTrans without pre-training,
the MTSMAE improves the prediction performance in almost all datasets with
different prediction lengths (except in the case of $L_{y}=192$ in the ETTm1
dataset).

In order to better evaluate the prediction ability of different models, we
drew one dimension of the prediction results of the $test$ $set$ from the ECL
dataset for qualitative comparison. Fig. 3 intuitively shows the prediction
results of different models. All models can roughly predict the development
trend of MTSD, but our model shows the best performance in different models.
Since prediction results are more consistent with the ground truth, our model
can accurately predict the periodicity and long-term changes. Compared with
other Transformer based models, prediction results of LSTM are smoother, so it
can’t well predict the trend of some inflection points, which is why LSTM
can’t be used in long time-series prediction.

TABLE IV: Different input lengths for two prediction lengths Predicition length | 336 | 480
---|---|---
Encoder’s input | 336 | 480 | 720 | 960 | 1200 | 1440 | 480 | 720 | 960 | 1200 | 1440
Informer | MSE | 3.339 | 3.256 | 3.423 | 3.460 | 3.371 | 3.099 | 3.374 | 3.504 | 3.441 | 3.441 | 3.400
MAE | 1.393 | 1.379 | 1.421 | 1.418 | 1.401 | 1.337 | 1.403 | 1.434 | 1.421 | 1.417 | 1.408
Reformer | MSE | 2.471 | 3.692 | 4.561 | 3.479 | 3.677 | - | 3.078 | 3.052 | 2.873 | 3.271 | -
MAE | 1.229 | 1.482 | 1.692 | 1.489 | 1.664 | 1.381 | 1.473 | 1.451 | 1.324
LogTrans | MSE | 3.217 | 3.917 | 3.391 | - | - | - | 3.431 | - | - | - | -
MAE | 1.461 | 1.655 | 1.554 | 1.579
PatchTrans | MSE | 3.588 | 3.996 | 3.919 | 2.963 | 2.937 | 3.588 | 3.453 | 3.069 | 2.900 | 2.789 | 3.087
MAE | 1.613 | 1.681 | 1.673 | 1.465 | 1.408 | 1.563 | 1.615 | 1.423 | 1.438 | 1.370 | 1.451
MTSMAE | MSE | 3.507 | 3.025 | 3.386 | 2.946 | 2.691 | 2.883 | 2.974 | 2.692 | 2.769 | 2.746 | 2.575
MAE | 1.572 | 1.492 | 1.569 | 1.397 | 1.377 | 1.365 | 1.506 | 1.373 | 1.365 | 1.313 | 1.362

####  IV-C2 Input Length

We perform the sensitivity analysis (input length) of the proposed MTSMAE
model on ETTh2 dataset. Compared with the previous experiment, $L_{x}=784$, in
this experiment, we gradually extended the length of the input sequence, i.e.,
{336, 480, 720, 960, 1200, 1440}, while keeping the prediction sequence
unchanged. We choose two prediction lengths, 336 and 480, and the length of
the label sequence is consistent with the length of the prediction sequence.

The experiment results are shown in the Table IV. It can be shown that
increasing the input length will lead to the decline of MSE and MAE, because
it brings repeated short-term patterns, and longer encoder inputs may contain
more dependencies. However, with the increasing of input sequence, although
the input information will increase, the influence of noise will also
increase. Some models can not effectively remove the impact of these noises,
so MSE and MAE may increase. The performance of our proposed MTSMAE in short
sequence input is not the best, but with the increase of input sequence, the
prediction performance of MTSMAE is getting better and better. There is no
problem of out-of-memory due to the long input sequence and sudden decline in
prediction performance, which demonstrates that MTSMAE has good denoising
ability.

####  IV-C3 Masking ratio

(a) MSE

(b) MAE

Figure 4: Results of MSE and MAE at different prediction lengths with
different masking ratios.

Fig. 4 shows the influence of different masking ratio. The dataset of this
experiment is ETTh2. The default encoder and decoder layers are selected
during pre-training, i.e., the number of decoder layers is 3 and decoder
layers is 1. Different broken lines in the figure represent different
prediction lengths. The optimal ratio is beyond our expectation, which is
higher than the original MAE’s optimal ratios (75%), and this masking ratio is
also opposite to the self-supervised method BERT in the NLP field, which
selects a masking ratio of 15% (a very low masking ratio).

High ratio masking is conducive to the fine-tuning of subsequent downstream
tasks, which related to the data characteristics of multivariate time-series.
High ratio masking can lose more information, greatly reduce the redundancy of
data, and increase the overall understanding of the model beyond low-level
information. However, a too high masking ratio, such as 95%, means that a
large amount of data is lost. The data that the model can learn is finite,
which is not conducive to the model’s understanding of the data.

####  IV-C4 Decoder depth

TABLE V: Different decoder depths

Blocks | 1 | 2 | 3 | 4
---|---|---|---|---
24 | MSE | 0.648 | 0.619 | 0.609 | 0.584
MAE | 0.661 | 0.644 | 0.649 | 0.628
48 | MAE | 1.126 | 1.058 | 0.975 | 1.111
MAE | 0.874 | 0.824 | 0.807 | 0.856
96 | MSE | 2.561 | 2.680 | 2.461 | 2.806
MAE | 1.304 | 1.380 | 1.306 | 1.346
192 | MSE | 3.341 | 3.262 | 3.501 | 3.752
MAE | 1.588 | 1.571 | 1.541 | 1.566
384 | MSE | 3.107 | 2.931 | 2.787 | 2.932
MAE | 1.403 | 1.933 | 1.344 | 1.398

As shown in the Table V, our MTSMAE decoder can be designed flexibly. We
changed the depth of the decoder (the number of Transformer decoder blocks),
i.e., {1, 2, 3, 4 }, keep the fixed depth of the encoder as 3. The last few
layers of the autoencoder are dedicated to reconstruction and have little
correlation with acquiring semantics. During pre-training, our encoder is also
used to recover the original data. A reasonable depth decoder can explain the
specialization of reconstruction, leaving the latent representation at a more
abstract level.

The dataset selected in this experiment is ETTh2, and we choose 85% masking
ratio. As shown in Table V, the deeper the decoder has the correspondingly
better performance of downstream prediction tasks. However, when the number of
layers of the decoder reaches 4 (more than the number of layers of the
encoder), the prediction performance of MTSMAE begins to decline. A deep
enough decoder is extremely important for the downstream prediction task, but
too deep decoder is easy to affect the performance of the encoder to obtain
semantics.

##  V Conclusion

In this paper, we studied the problem of multivariable time-series
forecasting, and proposed a model named MTSMAE for multivariable time-series
data. Specifically, firstly according to the data characteristics of
multivariable time-series, we proposed a new patch embedding method, which can
reduce the use of memory, so that the model can process longer sequences.
Then, we further proposed a self-supervised pre-training method based on
masked autoencoders (MAE). A large number of experiments on real datasets have
proved the effectiveness of MTSMAE in improving multivariable time-series
forecasting.

## References

  * [1] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” _Advances in neural information processing systems_ , vol. 25, 2012.
  * [2] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” _Advances in neural information processing systems_ , vol. 30, 2017.
  * [3] P. Tang and X. Zhang, “Features fusion framework for multimodal irregular time-series events,” _arXiv preprint arXiv:2209.01728_ , 2022.
  * [4] A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai, T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, and S. a. Gelly, “An image is worth 16x16 words: Transformers for image recognition at scale,” in _International Conference on Learning Representations_ , 2021.
  * [5] K. He, X. Chen, S. Xie, Y. Li, P. Dollár, and R. Girshick, “Masked autoencoders are scalable vision learners,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 16 000–16 009.
  * [6] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” _arXiv preprint arXiv:1810.04805_ , 2018.
  * [7] A. Radford, K. Narasimhan, T. Salimans, I. Sutskever _et al._ , “Improving language understanding by generative pre-training,” 2018.
  * [8] A. Radford, J. Wu, R. Child, D. Luan, D. Amodei, I. Sutskever _et al._ , “Language models are unsupervised multitask learners,” _OpenAI blog_ , vol. 1, no. 8, p. 9, 2019.
  * [9] M. Chen, A. Radford, R. Child, J. Wu, H. Jun, D. Luan, and I. Sutskever, “Generative pretraining from pixels,” in _International conference on machine learning_. PMLR, 2020, pp. 1691–1703.
  * [10] H. Bao, L. Dong, and F. Wei, “Beit: Bert pre-training of image transformers,” _arXiv preprint arXiv:2106.08254_ , 2021.
  * [11] Z. Xie, Z. Zhang, Y. Cao, Y. Lin, J. Bao, Z. Yao, Q. Dai, and H. Hu, “Simmim: A simple framework for masked image modeling,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2022, pp. 9653–9663.
  * [12] X. Geng, H. Liu, L. Lee, D. Schuurams, S. Levine, and P. Abbeel, “Multimodal masked autoencoders learn transferable representations,” _arXiv preprint arXiv:2205.14204_ , 2022.
  * [13] X. Li, W. Wang, L. Yang, and J. Yang, “Uniform masking: Enabling mae pre-training for pyramid-based vision transformers with locality,” _arXiv preprint arXiv:2205.10063_ , 2022.
  * [14] C. Doersch, A. Gupta, and A. A. Efros, “Unsupervised visual representation learning by context prediction,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 1422–1430.
  * [15] X. Wang and A. Gupta, “Unsupervised learning of visual representations using videos,” in _Proceedings of the IEEE international conference on computer vision_ , 2015, pp. 2794–2802.
  * [16] M. Noroozi and P. Favaro, “Unsupervised learning of visual representations by solving jigsaw puzzles,” in _European conference on computer vision_. Springer, 2016, pp. 69–84.
  * [17] T. Brown, B. Mann, N. Ryder, M. Subbiah, J. D. Kaplan, P. Dhariwal, A. Neelakantan, P. Shyam, G. Sastry, A. Askell _et al._ , “Language models are few-shot learners,” _Advances in neural information processing systems_ , vol. 33, pp. 1877–1901, 2020.
  * [18] G. E. Hinton and R. Zemel, “Autoencoders, minimum description length and helmholtz free energy,” _Advances in neural information processing systems_ , vol. 6, 1993.
  * [19] G. Lai, W.-C. Chang, Y. Yang, and H. Liu, “Modeling long-and short-term temporal patterns with deep neural networks,” in _The 41st international ACM SIGIR conference on research & development in information retrieval_, 2018, pp. 95–104.
  * [20] S. Bai, J. Z. Kolter, and V. Koltun, “An empirical evaluation of generic convolutional and recurrent networks for sequence modeling,” _arXiv preprint arXiv:1803.01271_ , 2018.
  * [21] N. Kitaev, L. Kaiser, and A. Levskaya, “Reformer: The efficient transformer,” in _ICLR_ , 2019.
  * [22] S. Li, X. Jin, Y. Xuan, X. Zhou, W. Chen, Y.-X. Wang, and X. Yan, “Enhancing the locality and breaking the memory bottleneck of transformer on time series forecasting,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [23] H. Zhou, S. Zhang, J. Peng, S. Zhang, J. Li, H. Xiong, and W. Zhang, “Informer: Beyond efficient transformer for long sequence time-series forecasting,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 35, no. 12, 2021, pp. 11 106–11 115.
  * [24] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” _Neural computation_ , vol. 9, no. 8, pp. 1735–1780, 1997.
  * [25] I. Loshchilov and F. Hutter, “Decoupled weight decay regularization,” in _ICLR_ , 2019.
  * [26] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _International Conference on Learning Representations_ , 2015. [Online]. Available: http://arxiv.org/abs/1412.6980
  * [27] I. Loshchilov and F. Hutter, “Sgdr: Stochastic gradient descent with warm restarts,” in _ICLR_ , 2017.
  * [28] P. Goyal, P. Dollár, R. Girshick, P. Noordhuis, L. Wesolowski, A. Kyrola, A. Tulloch, Y. Jia, and K. He, “Accurate, large minibatch sgd: Training imagenet in 1 hour,” _arXiv preprint arXiv:1706.02677_ , 2017.
  * [29] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga _et al._ , “Pytorch: An imperative style, high-performance deep learning library,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [30] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [31] J. L. Ba, J. R. Kiros, and G. E. Hinton, “Layer normalization,” _arXiv preprint arXiv:1607.06450_ , 2016.


Optimal designs for experiments for scalar-on-function linear models

D. Michaelides111Damianos Michaelides<EMAIL_ADDRESS>School of
Mathematical Sciences, University of Southampton, Southampton, SO17 1BJ, UK. ,
M. Adamou, D. C. Woods & A. M. Overstall

Southampton Statistical Sciences Research Institute

University of Southampton, United Kingdom

The aim of this work is to extend the usual optimal experimental design
paradigm to experiments where the settings of one or more factors are
functions. For these new experiments, a design consists of combinations of
functions for each run of the experiment along with settings for non-
functional variables. After briefly introducing the class of functional
variables, basis function systems are described. Basis function expansion is
applied to a functional linear model consisting of both functional and scalar
factors, reducing the problem to an optimisation problem of a single design
matrix.

Keywords: Design of Experiments; Functional Linear Model; Profile Factors;
Basis functions

##  1 Introduction

In science and engineering, an increasing number of experiments involve the
investigation of the relationship between a response and functional or profile
variables, i.e. variables whose values can be varied as a function, usually of
time, within a single run of an experiment. A common example of a profile
variable is temperature, being varied monotonically or as a step function
through the run. The statistical design problem then becomes choosing suitable
functions that determine how each profile variable varies during each run. For
ease of exposition, throughout this paper we shall assume time
$t\in[0,\mathcal{T}]$ is the continuous single input to functional variables.
The methods extend naturally to situations where there are functional
variables with multiple inputs, e.g. spatio-temporal studies.

Statistical modelling with functional data is well established in the
statistics literature, see Ramsay and Silverman (1997) and Ramsay and
Silverman (2005) for an introduction. However, the design of experiments for
such models has received much less attention, with two main approaches being
proposed: response surface methods using dimension-reduction techniques
(Georgakis 2013, Roche 2015, Roche 2018 and Klebanov and Georgakis 2016) and
optimal design for dynamic models, typically derived from differential
equations (Balsa-Canto and Banga 2007 and Uciński and Bogacka 2007). We
develop methodology related to the response surface approach to find optimal
functions for profile variables assuming a scalar-on-function linear model.
Our approach leverages the power of standard linear model optimal design
methodology and it’s flexibility allows designs to be obtained for various
different scenarios, including common experiments in the pharmaceutical
industry, and using different optimality criteria.

In Section 2, with a view to model and design experiments involving profile
factors, a functional linear model along with an optimal design strategy is
developed. Models involving both profile and scalar factors are tackled and a
basis function expansion is used for both the profile factors and functional
parameters. Functional models turn out to be extension of the traditional
linear model and the model developed is capable of incorporating an unlimited
number of profile and scalar factors. In Section 3, examples with only profile
and both profile and scalar variables are presented, for which it is assumed
that control of the profile variables is represented by step functions. In
Section 4, the set of possible functions for the profile variables is extended
by the use of spline basis expansions and the performance of the designs is
compared to the designs found under the step function basis. Section 5
provides a conclusion of the findings and future research goals.

##  2 The scalar-on-function linear model and optimal design

###  2.1 The scalar-on-function linear model

Let $y_{i}$ denote a scalar response observed at time $\mathcal{T}$, the end
of the experiment, for the $i$th run of the experiment ($i=1,\ldots,n$). An
additive scalar-on-function linear model of the following form is assumed for
$y_{i}$:

$y_{i}=\beta_{0}+\sum_{j=1}^{f_{1}}\int_{0}^{\mathcal{T}}\beta_{j}(t)x_{ij}(t)\;dt+\sum_{k=f_{1}+1}^{f_{1}+f_{2}}\beta_{k}x_{ik}+\varepsilon_{i}\,,$
(1)

where $f_{1}$ is the number of functional variables involved in the
experiment, $f_{2}$ is the number of scalar variables involved in the
experiment, $\beta_{j}(t):[0,\mathcal{T}]\rightarrow\mathbb{R}$ is an unknown
function of time $0\leq t\leq\mathcal{T}$,
$x_{ij}(t):[0,\mathcal{T}]\rightarrow[u,v]$ the function controlling the $j$th
functional variable in the interval $[u,v]$,
$\beta_{0},\beta_{f_{1}+1},...,\beta_{f_{1}+f_{2}}$ are unknown constant
parameters, $x_{ik}\in[u_{s},v_{s}]$ is the value for the $k$th scalar factor
lying in the interval $[u_{s},v_{s}]$ and
$\varepsilon_{i}\sim\mathcal{N}(0,\sigma^{2})$, with errors from different
experiment runs assumed independent. The scalars $u,v,u_{s},v_{s}$ are
constants representing the bounds for the functional and scalar variables,
respectively.

Viewed via an increasingly fine discretisation of $t$, $t_{1},\dots,t_{m}$
with $m\rightarrow\infty$, model (1) has essentially an infinite number of
unknowns $\beta_{j}(t_{k})$. Hence, as the experiment only returns a finite
amount of data, the system is underdetermined and there will be an infinite
number of solutions
$\hat{\bm{\beta}}^{\rm{T}}(t)=(\hat{\beta}_{0},\hat{\beta}_{1}(t),\ldots,\hat{\beta}_{f_{1}}(t),\hat{\beta}_{f_{1}+1},\ldots,\hat{\beta}_{f_{1}+f_{2}})$
that each give a perfect fit to the data. Hence, inference from a given
experiment requires some restriction on the function space. We implicitly
define spaces of functions $\beta_{j}(t)\in\mathcal{B}_{j}[0,\mathcal{T}]$ and
$x_{ij}(t)\in\mathcal{X}_{j}[0,\mathcal{T}]$ via basis function expansions,
with $\mathcal{B}_{j}[0,\mathcal{T}]$ and $\mathcal{X}_{j}[0,\mathcal{T}]$ the
function spaces of the functional parameters and profile variables
respectively, from time 0 to time $\mathcal{T}$. For the functional parameters

$\displaystyle\beta_{j}(t)$ $\displaystyle=$
$\displaystyle\sum_{l=1}^{n_{\beta,j}}\theta_{jl}b_{jl}(t)$ (2)
$\displaystyle=$ $\displaystyle\bm{b}_{j}(t)^{\rm{T}}\bm{\theta}_{j}\,,$

where $\bm{b}_{j}(t)^{\rm{T}}=[b_{j1}(t),\ldots,b_{jn_{\beta,j}}(t)]$ are
known basis functions with corresponding (unknown) coefficients
$\bm{\theta}_{j}^{\rm{T}}=(\theta_{j1},\ldots,\theta_{jn_{\beta,j}})$. Hence,
we have reduced the problem of estimating an unknown function to the problem
of estimating $n_{\beta,j}$ coefficients. We can define the functions
$x_{ij}(t)$ via a similar expansion,

$x_{ij}(t)=\sum_{l=1}^{n_{x,j}}\gamma_{ijl}c_{l}(t)\,,$ (3)

and for $\bm{x}_{j}(t)^{\rm{T}}=[x_{j1}(t),\ldots,x_{jn}(t)]$ we have that
$\bm{x}_{j}(t)=\Gamma_{j}\bm{c}_{j}(t)$, where
$\bm{c}_{j}(t)^{\rm{T}}=[c_{j1}(t),\ldots,c_{jn_{x,j}}(t)]$ are known basis
functions, $\Gamma_{j}$ is a $n\times n_{x,j}$ coefficient matrix with $il$th
entry $\gamma_{ijl}$. Where necessary, linear transformations of the basis
expansions are used to ensure the desired image of the function $x_{ij}(t)$ is
obtained.

A substitution of the basis expansions into model (1) gives a familiar model
form for the response vector $\bm{y}^{\rm{T}}=(y_{1},\ldots,y_{n})$:

$\begin{split}\bm{y}&=\beta_{0}\bm{1}_{n}+\sum_{j=1}^{f_{1}}\int_{0}^{\mathcal{T}}\beta_{j}(t)\bm{x}_{j}(t)\;dt+\sum_{k=f_{1}+1}^{f_{1}+f_{2}}\beta_{k}\bm{x}_{k}+\bm{\varepsilon}\\\
&=\beta_{0}\bm{1}_{n}+\sum_{j=1}^{f_{1}}\int_{0}^{\mathcal{T}}\bm{\Gamma}_{j}\bm{c}_{j}(t)\bm{b}_{j}(t)^{T}\bm{\theta}_{j}\;dt+\sum_{k=f_{1}+1}^{f_{1}+f_{2}}\beta_{k}\bm{x}_{k}+\bm{\varepsilon}\\\
&=\beta_{0}\bm{1}_{n}+\bm{J}_{\gamma
cb}\bm{\theta}+\bm{X\beta}+\bm{\epsilon}\\\ \end{split}$ (4)

where $\bm{1}_{n}$ is the $n$-vector with every entry equal to one,
$\bm{\theta}^{\rm{T}}=(\bm{\theta}_{1}^{\rm{T}},\ldots,\bm{\theta}_{f_{1}}^{\rm{T}})$,
$\bm{X}$ the $n\times f_{2}$ model matrix for the scalar factors,
$\bm{\beta}^{\rm{T}}=(\beta_{f_{1}+1},\ldots,\beta_{f_{1}+f_{2}})$ and
$\bm{\varepsilon}^{\rm{T}}=(\varepsilon_{1},\ldots,\varepsilon_{n})$. The
$n\times\sum n_{\beta,j}$ matrix $\bm{J}_{\gamma cb}$ is formed by column
binding the $f_{1}$ sub-matrices

$\begin{split}\bm{J}_{\gamma
cb_{j}}&=\int_{0}^{\mathcal{T}}\bm{\Gamma}_{j}\;\bm{c}_{j}(t)\;\bm{b}_{j}(t)^{T}\;dt\,\quad
j=1,\ldots,f_{1}\\\
&=\bm{\Gamma}_{j}\int_{0}^{\mathcal{T}}\bm{c}_{j}(t)\;\bm{b}_{j}(t)^{T}\;dt,\end{split}$
(5)

with the integral of the matrix defined as the matrix of integrals, i.e. the
integral of a vector $\bm{\alpha}(t)^{T}=[\alpha_{1}(t),...,\alpha_{n}(t)]$ is
defined as

$\int_{0}^{\mathcal{T}}\bm{\alpha}(t)\;dt=\Bigg{[}\int_{0}^{\mathcal{T}}\alpha_{1}(t)\;dt,...,\int_{0}^{\mathcal{T}}\alpha_{n}(t)\;dt\Bigg{]}^{T}.$

Hence model (4) takes the form of the familiar linear model,

$\bm{y}=\bm{Z\nu}+\bm{\varepsilon}\,,$

with $n\times(\sum n_{\beta,j}+f_{2}+1)$ model matrix

$\bm{Z}=[\bm{1\>\;J_{\gamma cb}\>\;X}]\,,$ (6)

and ($\sum n_{\beta,j}+f_{2}+1$)-vector $\bm{\nu}$ of coefficients

$\bm{\nu}^{T}=[\beta_{0}\>\;\bm{\theta}^{T}\>\;\bm{\beta}^{T}]\,.$ (7)

The least squares estimator for $\bm{\nu}$ hence takes the form

$\hat{\bm{\nu}}=(\bm{Z}^{T}\bm{Z})^{-1}\bm{Z}^{T}\bm{y}\,.$ (8)

with variance-covariance matrix

$\begin{split}\text{Var}(\hat{\bm{\nu}})&=\text{Var}\Big{[}(\bm{Z}^{T}\bm{Z})^{-1}\bm{Z}^{T}\bm{y}\Big{]}\\\
&=\sigma^{2}(\bm{Z}^{T}\bm{Z})^{-1}.\end{split}$ (9)

###  2.2 The reduced design problem

Representation of the functional variables and parameters via basis function
expansions in model (2.1) reduces the problem of finding a functional optimal
design to one of choosing coefficient matrices $\bm{\Gamma}_{j}$, together
with values for the scalar inputs, to minimise a functional of equation (9).
Here, we focus on point estimation of $\bm{\theta}$ using $A$-optimality, and
find designs that lead to information matrices $\bm{M}=\bm{Z}^{\rm{T}}\bm{Z}$
that minimise

$\begin{split}\Psi_{A}(\bm{M})&=\text{tr}\big{\\{}\text{Var}(\hat{\bm{\nu}})/\sigma^{2}\big{\\}}\\\
&=\text{tr}\big{\\{}\bm{M}^{-1}\big{\\}}\,.\end{split}$ (10)

Clearly, an $A$-optimal design will depend on the choice of bases for both the
functional parameters and inputs, as well as the number of functions in each
of these bases. In what follows, we investigate the impact of these choices
for a number different scenarios.

For two designs with the same size $n$ resulting in information matrices
$M_{1}$ and $M_{2}$, we define the relative $A$-efficiency of design 2 with
respect to design 1 as $\Psi_{A}(M_{1})/\Psi_{A}(M_{2})$.

##  3 Optimal designs using step functions

For the examples to be considered in this section, control of a functional
variable is obtained via choice of step functions defined as

$x_{ij}(t)=\sum_{l=1}^{n_{x,j}}\gamma_{ijl}\bm{1}_{\lambda_{l}}(t)\,,$

for breakpoints
$\lambda_{0}=0<\lambda_{1}<\cdots\lambda_{n_{x,j}}=\mathcal{T}$ and

$\bm{1}_{\lambda_{l}}(t)=\begin{cases}1&\mbox{if
}\lambda_{l-1}<t<\lambda_{l}\\\ 0&\mbox{otherwise}\,.\end{cases}$

An example of such a function with two interior break points, and hence
$n_{x,j}=3$, is given in Figure 1.

Figure 1: Step function for $x_{ij}(t)$ with two interior breakpoints,
$n_{x,j}=3$.

$A$-optimal designs are found for two choices of power series bases for
$\beta_{j}(t)$, a linear basis

$\beta_{j}(t)=\alpha_{0}+\alpha_{1}t\,,$ (11)

and a quadratic basis

$\beta_{j}(t)=\alpha_{0}+\alpha_{1}t+\alpha_{2}t^{2}\,.$ (12)

See Figure 2. Designs are sought for two experiments: (i) with a single
functional variable, and (ii) with one functional variable and three scalar
variables.

Figure 2: (a) Linear and (b) Quadratic basis functions for $\beta(t)$.

###  3.1 Example 1: Single profile factor

Under this example, $f_{1}=1$ and $f_{2}=0$ and model (4) is simplified such
that

$y_{i}=\beta_{0}+\int_{0}^{\mathcal{T}}\beta(t)x_{i}(t)\;dt+\varepsilon_{i}\,,$
(13)

for $i=1,...,n$. We enforce the constraint $-1\leq x_{i}(t)\leq 1$ via linear
transformation and construct step functions assuming $n_{x}+2$ equally-spaced
breakpoints on $[0,1]$. To ensure identifiability of model (13), we restrict
the number of basis functions in the expansion of $x_{i}(t)$ to be $n_{x}\
ge2$ when finding designs for the linear $\beta(t)$ basis and $n_{x}\geq 3$
when considering the quadratic $\beta(t)$ basis.

Designs for different combinations of $n$ and $n_{x}$ are found under both
basis expansions (11) and (12) using 1000 random starts of a coordinate
exchange algorithm (Goos and Jones 2011). The results are given in Table 1 for
the linear basis and Table 2 for the quadratic basis.

Some general findings can be observed from the results. For both functional
parameter bases, the $A$-optimality objective function values decrease with
experiment size $n$, as expected, with typically, the designs for $n=8$ and
$n=12$ containing replicates of the functions from the $n=4$ design. Also, the
$A$-efficiency of the designs (within each experiment size) generally decrease
with $n_{x}$. The exception is for the case $n_{x}=3$ for the linear
functional parameter basis, where the inability to make a step change in the
centre of the interval $[0,1]$ (as $\lambda=0.5$ is not included in the set of
breakpoints) leads to an increase in the average variance compared to the
design for $n_{x}=2$. Again, increasing design performance with $n_{x}$ is
generally expected, as increasing the breakpoints enlarges the set of possible
functions from which the design can be chosen. Changes in design performance
with $n_{x}$ are greater for the quadratic $\beta(t)$ basis, as more complex
functions are chosen in the $A$-optimal design. Representative example designs
are shown in Figures 3 and 4 for the linear and quadratic bases respectively.

|  $n=4$ |  $n=8$ |  $n=12$
---|---|---|---
$n_{x}$ |  A-opt |  A-eff |  A-opt |  A-eff |  A-opt |  A-eff
2 |  8.750 |  0.961 |  3.958 |  0.981 |  2.583 |  0.972
3 |  8.828 |  0.952 |  4.287 |  0.906 |  2.778 |  0.904
4 |  8.750 |  0.961 |  3.903 |  0.995 |  2.570 |  0.977
8 |  8.493 |  0.990 |  3.902 |  0.995 |  2.539 |  0.989
16 |  8.427 |  0.997 |  3.887 |  0.999 |  2.520 |  0.997
100 |  8.404 |  1.000 |  3.882 |  1.000 |  2.512 |  1.000
Table 1: Objective function values and A-efficiency values with $n=4,8,12$ for the linear basis for $\beta(t)$ and step basis for $x(t)$. |  $n=4$ |  $n=8$ |  $n=12$
---|---|---|---
$n_{x}$ |  A-opt |  A-eff |  A-opt |  A-eff |  A-opt |  A-eff
3 |  386.408 |  0.535 |  189.766 |  0.510 |  126.409 |  0.499
4 |  246.869 |  0.838 |  103.553 |  0.934 |  67.735 |  0.931
8 |  218.479 |  0.947 |  99.109 |  0.976 |  65.217 |  0.966
16 |  208.843 |  0.991 |  97.408 |  0.993 |  63.610 |  0.991
100 |  206.884 |  1.000 |  96.709 |  1.000 |  63.028 |  1.000
Table 2: Objective function values and A-efficiency values with $n=4,8,12$ for
the quadratic basis for $\beta(t)$ and step basis for $x(t)$. Figure 3: Four
run optimal design for $n_{x}=16$, linear basis for $\beta(t)$ and step basis
for $x(t)$. Figure 4: Four run optimal design for $n_{x}=8$, quadratic basis
for $\beta(t)$ and step basis for $x(t)$.

###  3.2 Example 2: Single profile and three scalar factors

For this example, the model is extended to include a single profile factor and
$f_{2}=3$ scalar factors. The functional linear model takes the form,

$y_{i}=\beta_{0}+\int_{0}^{\mathcal{T}}\beta_{1}(t)x_{i1}(t)\;dt+\bm{f}^{T}(\bm{x}_{i})\bm{\beta}+\epsilon_{i},$
(14)

for $i=1,...,n$, assuming $-1\leq x_{i}(t)\leq 1$ and,

$\bm{f}^{T}(\bm{x}_{i})\bm{\beta}=\begin{cases}\beta_{2}x_{i2}+\beta_{3}x_{i3}+\beta_{4}x_{i4}&\text{Case
}1\\\
\beta_{2}x_{i2}+\beta_{3}x_{i3}+\beta_{4}x_{i4}+\beta_{5}x_{i2}^{2}+\beta_{6}x_{i3}^{2}+\beta_{7}x_{i4}^{2}&\text{Case
}2\end{cases}$

being the model for scalar factors, i.e. main effects or main and quadratic
effects. Control of the single profile factor is assumed to be represented by
a step function basis, the functional parameter $\beta_{1}(t)$ is assumed to
be represented by both a linear and quadratic basis, similar to Section 3.1,
and $n=12$ runs is chosen. In this example, the sensitivity study raised is
partly similar to the study in Section 3.1. The number of basis functions for
the profile factor and functional parameters are varied as before and
additionally the scalar part with main and main and quadratic effects is
considered. The optimal designs for the single profile factor are exactly the
same as the designs found in Section 3.1, so no further elaboration is
required.

What is interesting to investigate is how the designs for the scalar factors
behave from main to main and quadratic effects. For the Case 1 model, i.e. the
model with main effects, the optimal choices of the scalar factors are those
at the boundaries. This is similar to the profile factor behaviour with linear
basis for the parameters, where the function of the profile factor changes at
most once. For the Case 2 model, i.e. the model with both main and quadratic
effects, the optimal design for the scalar factors include boundary points and
centre points in order to be able to estimate the curvature. This is similar
to the profile factor behaviour with quadratic basis for the parameters, where
the function of the profile factor changes at most twice. For the Case 2
model, the columns for the scalar factors are orthogonal. The optimal designs
for the scalar factors with a linear basis function for the functional
parameter and $n_{x}=4$ are given in Table 3.

Case 1
---
$i$ | $x_{i2}$ | $x_{i3}$ | $x_{i4}$
1 | -1 | 1 | -1
2 | -1 | 1 | 1
3 | 1 | 1 | -1
4 | -1 | -1 | 1
5 | -1 | -1 | -1
6 | -1 | 1 | 1
7 | 1 | -1 | 1
8 | 1 | -1 | -1
9 | 1 | 1 | 1
10 | 1 | -1 | 1
11 | 1 | 1 | -1
12 | -1 | -1 | -1

Case 2
---
$i$ | $x_{i2}$ | $x_{i3}$ | $x_{i4}$
1 | 1 | 1 | 1
2 | 0 | 0 | -1
3 | 1 | 0 | 0
4 | 0 | -1 | 1
5 | -1 | 0 | 1
6 | 0 | 1 | 1
7 | -1 | -1 | 0
8 | -1 | 0 | 0
9 | 0 | 0 | 0
10 | 0 | 1 | 0
11 | 0 | 0 | 0
12 | -1 | 1 | -1

Table 3: Optimal designs for main effects (Case 1) and main & quadratic
effects (Case 2) for linear basis function for the functional linear model.

##  4 Spline basis functions

To assess whether design performance is improved by widening the set of
possible functions from which we choose the runs for the functional variable,
we now use a spline basis expansion for $x_{ij}(t)$. A spline is a piecewise
polynomial with continuity, and possibly continuity of derivatives, defined at
the breakpoints. We use a B-spline basis expansion (see De Boor et al. 1978)
of degree 1 on the interval $[0,1]$. B-splines have a number of desirable
properties as a spline basis, see Dierckx (1995, Chapter 1), including ease of
computation through recursive formulae, they sum to unity for any $t$ in the
domain, and they are more numerically stable than alternative bases, e.g. the
truncated power basis.

We repeat the example of Section 3.1 and find $A$-optimal designs for model
(13) with a single functional variable but assuming a linear B-spline for
$x_{i}(t)$. We find designs for a linear functional parameter, $n=4,8,12$ runs
and choose the breakpoint vector such that $n_{x}=3,4,8,16,100,201$, assuming
equally spaced breakpoints. Results are given in Table 4, and example designs
are displayed in Figures 5 and 6. As in the previous section, design
performance increases with $n_{x}$ for all values of $n$, but with much
smaller increases for $n>16$. The functions chosen for the design actually
become somewhat simpler for larger values of $n_{x}$, e.g. for $n_{x}=8$ the
chosen functions commonly have four changes in slope (see Figure 5) whereas
each function has only two changes when $n_{x}=16$ (Figure 6). This difference
demonstrates the importance of the location of the breakpoints, as with larger
values of $n_{x}$ there is more flexibility in the location of the changes in
slope. For larger $n_{x}$, each chosen function consists of two constant
sections, at the extremes of the design space, joined by a linear component
with finite slope. However, comparison of the objective function values for
the designs found assuming step functions (Table 1) shows that the added
complexity of the linear spline basis did not result in better designs.

|  $n=4$ |  $n=8$ |  $n=12$
---|---|---|---
$n_{x}$ |  A-opt |  A-eff |  A-opt |  A-eff |  A-opt |  A-eff
3 |  12.471 |  0.674 |  6.224 |  0.624 |  4.123 |  0.609
4 |  9.314 |  0.902 |  4.168 |  0.931 |  2.759 |  0.910
8 |  8.594 |  0.978 |  3.940 |  0.985 |  2.571 |  0.977
16 |  8.433 |  0.997 |  3.895 |  0.997 |  2.528 |  0.994
100 |  8.405 |  1.000 |  3.882 |  1.000 |  2.512 |  1.000
201 |  8.404 |  1.000 |  3.882 |  1.000 |  2.512 |  1.000
Table 4: Objective function values and A-efficiency values with $n=4,8,12$ for
the linear basis for $\beta(t)$ and B-spline basis for $x(t)$.

Figure 5: Four run optimal design for $n_{x}=8$, linear basis for $\beta(t)$
and B-spline basis for $x(t)$. Figure 6: Four run optimal design for
$n_{x}=16$, linear basis for $\beta(t)$ and B-spline basis for $x(t)$.

##  5 Conclusions

We have demonstrated a general methodology for finding optimal designs for
experiments involving functional variables and a scalar response. The
methodology uses basis function expansions of the functional variable and the
functional parameters in a scalar-on-function linear model. It is flexible and
can be applied assuming a variety of different bases for both functional
variables and parameters. For simple functional parameters, e.g. linear, only
simple forms for the functional variable are required.

Ongoing work is extending the results to different criteria tailored to
estimation of the functional parameters, finding designs for scalar-on-
function generalised linear models, and developing methods for choosing
optimal functions for inputs into dynamic models.

## References

  * (1)
  * Balsa-Canto and Banga (2007) Balsa-Canto, Eva, M. R.-F. and Banga, J. R. (2007), ‘Optimal design of dynamic experiments for improved estimation of kinetic parameters of thermal degradation’, Journal of Food Engineering 82(2), 178–188.
  * De Boor et al. (1978) De Boor, C., De Boor, C., Mathématicien, E.-U., De Boor, C. and De Boor, C. (1978), A practical guide to splines, Vol. 27, springer-verlag New York.
  * Dierckx (1995) Dierckx, P. (1995), Curve and surface fitting with splines, Oxford University Press.
  * Georgakis (2013) Georgakis, C. (2013), ‘Design of dynamic experiments: A data-driven methodology for the optimization of time-varying processes’, Industrial & Engineering Chemistry Research 52(35), 12369–12382.
  * Goos and Jones (2011) Goos, P. and Jones, B. (2011), Optimal design of experiments: a case study approach, John Wiley & Sons.
  * Klebanov and Georgakis (2016) Klebanov, N. and Georgakis, C. (2016), ‘Dynamic response surface models: a data-driven approach for the analysis of time-varying process outputs’, Industrial & Engineering Chemistry Research 55(14), 4022–4034.
  * Ramsay and Silverman (1997) Ramsay, J. O. and Silverman, B. W. (1997), Functional Data Analysis, 1 edn, Springer.
  * Ramsay and Silverman (2005) Ramsay, J. O. and Silverman, B. W. (2005), Functional Data Analysis, 2 edn, Springer.
  * Roche (2015) Roche, A. (2015), ‘Response surface methodologiy for functional data with application to nuclear safety’.
  * Roche (2018) Roche, A. (2018), ‘Local optimization of black-box functions with high or infinite-dimensional inputs: application to nuclear safety’, Computational Statistics 33(1), 467–485.
  * Uciński and Bogacka (2007) Uciński, D. and Bogacka, B. (2007), ‘A constrained optimum experimental design problem for model discrimination with a continuously varying factor’, Journal of statistical planning and inference .


Astronomy Letters, 2010 Vol. 36, No. 3

Searching for Stars Closely Encountering with the Solar System

V.V. Bobylev

Pulkovo Astronomical Observatory, Russian Academy of Sciences, St-Petersburg

Abstract–Based on a new version of the Hipparcos catalog and currently
available radial velocity data, we have searched for stars that either have
encountered or will encounter the solar neighborhood within less than 3 pc in
the time interval from $-2$ Myr to $+2$ Myr. Nine new candidates within 30 pc
of the Sun have been found. To construct the stellar orbits relative to the
solar orbit, we have used the epicyclic approximation. We show that, given the
errors in the observational data, the probability that the well-known star HIP
89 825 (GL 710) encountering with the Sun most closely falls into the Oort
cloud is 0.86 in the time interval $1.45\pm 0.06$ Myr. This star also has a
nonzero probability, $1\times 10^{-4},$ of falling into the region $d<1000$
AU, where its influence on Kuiper Belt objects becomes possible.

DOI:10.1134/S1063773710030060

## INTRODUCTION

Interest in the problem of close encounters of field stars with the Solar
system stems from the fact that the passage of a star can lead to various
kinds of perturbations of Solar system objects.

Thus, for example, the interaction of a star with the Oort comet cloud can
give rise to comet showers reaching the region of the major planets (Hills
1981; Weissman 1996). The formation time scale of a comet shower is $\approx$1
Myr after the passage of a star. According to Weissman (1996), the Oort cloud
is a spheroid with a semimajor axis of about $10^{5}$ AU ($\approx$0.5 pc,
since 1 pc=206 265 AU) elongated toward the Galactic center and a semiminor
axis of about $8\times 10^{4}$ AU. The question about close encounters of
stars with the Sun within $r<2-5$ pc was considered by Revina (1988), Matthews
(1994), and Mullari and Orlov (1996) using various ground-based observations
and by Garcia-Sanchez et al. (1999, 2001) based on Hipparcos (1997) data in
combination with stellar radial velocity data. As a result, 156 Hipparcos
stars from the solar neighborhood 50 pc in radius that either have encountered
or will encounter with the Solar system within $r<5$ pc in the time interval
$\pm 10$ Myr are known to date (Garcia-Sanchez et al. 2001). Having analyzed
these data, Garcia-Sanchez et al. (2001) estimated the frequency of close
(within 1 pc) encounters of stars with the Sun to be $2.3\pm 0.2$ encounters
per Myr and, after a correction for the Hipparcos incompleteness, this value
increases to $11.7\pm 1.3$ encounters per Myr.

Analyzing the possibility of even closer encounters is of current interest.
For example, Kenyon and Bromley (2004) argue that only the passage of a star
can explain the peculiarities of the orbit of the minor planet 2003 VB12
(Sedna). The various kinds of influences of such a star on Kuiper Belt objects
were shown to manifest themselves at characteristic distances $r<1000$ AU
($\approx$0.005 pc). The Kuiper Belt proper extends from the Sun to a boundary
of $\approx$50 AU.

The goal of this study is to search for candidate stars closely encountering
with the Sun based on a new version of the Hipparcos catalog (van Leeuwen,
2007) and currently available stellar radial velocity data. We solve the
problem of statistical simulations by taking into account the random errors in
the input data and estimate the probability of a star penetrating into the
Oort cloud region and into the region of a possible influence of the star on
Kuiper Belt objects.

##  1 THE DATA

We use stars from the Hipparcos catalog (ESA 1997) while taking the new proper
motions and parallaxes from a revised version of this catalog (van Leeuwen,
2007); the stellar radial velocities are taken from the Pulkovo Compilation of
Radial Velocities (PCRV) (Gontcharov 2006) created as a result of implementing
the OSACA project (Bobylev et al., 2006) and containing radial velocity data
for about 35 000 stars. Note that, in contrast to the larger CRVAD-2 catalog
of radial velocities (Kharchenko et al., 2007), the PCRV catalog contains only
the stars with random errors in their radial velocities within 10 km s-1.

##  2 THE METHODS

###  2.1 Orbit Construction

We use a rectangular Galactic coordinate system with the axes directed away
from the observer toward the Galactic center ($l=0^{\circ},b=0^{\circ},$ the
$X$ axis), in the direction of Galactic rotation ($l=90^{\circ},b=0^{\circ},$
the $Y$ axis), and toward the North Pole ($b=90^{\circ},$ the $Z$ axis). The
corresponding space velocity components of an object $U,V,W$ are also directed
along the $X,Y,Z$ axes. The epicyclic approximation (Lindblad, 1927) allows
the stellar orbits to be constructed in a coordinate system rotating around
the Galactic center. We apply this method in the form given in Fuchs et al.
(2006):

$\halign{\hbox to\displaywidth{$\hfil\displaystyle#\hfil$}\cr
0.0pt{\hfil$\displaystyle\hfill X(t)=X(0)+{U(0)\over\kappa}\sin(\kappa
t)+\hfill\cr 0.0pt{\hfil$\displaystyle\hfill+{V(0)\over 2B}(1-\cos(\kappa
t)),\hfill\hbox to0.0pt{\hss(}1)\cr 0.0pt{\hfil$\displaystyle\hfill
Y(t)=Y(0)+2A\biggl{(}X(0)+{V(0)\over 2B}\biggr{)}t-\hfill\cr
0.0pt{\hfil$\displaystyle\hfill-{\Omega_{0}\over B\kappa}V(0)\sin(\kappa
t)+{2\Omega_{0}\over\kappa^{2}}U(0)(1-\cos(\kappa t)),\hfill\cr
0.0pt{\hfil$\displaystyle\hfill Z(t)={W(0)\over\nu}\sin(\nu t)+Z(0)\cos(\nu
t),\hfill\crcr}}}}}}$

where $t$ is the time in Myr (we proceed from the ratio pc/Myr = 0.978 km
s-1); $A$ and $B$ are the Oort constants; $\kappa=\sqrt{-4\Omega_{0}B}$ is the
epicyclic frequency; $\Omega_{0}$ is the angular velocity of Galactic rotation
for the local standard of rest, $\Omega_{0}=A-B$; $\nu=\sqrt{4\pi G\rho_{0}}$
is the vertical oscillation frequency, where $G$ is the gravitational constant
and $\rho_{0}$ is the star density in the solar neighborhood.

The parameters $X(0),Y(0),Z(0)$ and $U(0),V(0),W(0)$ in the system of
equations (1) denote the current stellar positions and velocities. The
displacement of the Sun from the Galactic plane is taken to be $Z(0)=17$ pc
(Joshi 2007). We calculate the velocities $U,V,W$ relative to the local
standard of rest (LSR) with $(U,V,W)_{LSR}=(10.00,5.25,7.17)$ km s-1 (Dehnen
and Binney 1998).

We took $\rho_{0}=0.1M_{\odot}$/pc3 (Holmberg and Flinn 2004), which gives
$\nu=74$ km s-1 kpc-1. We used the following Oort constants found by Bobylev
et al. (2008): $A=15.5\pm 0.3$ km s-1 kpc-1 and $B=-12.2\pm 0.7$ km s-1 kpc-1;
$\kappa=37$ km s-1 kpc-1 corresponds to these values.

Note that we neglect the gravitational interaction between the star and the
Sun.

As was shown by Mullari and Orlov (1996), even the simplest linear
approximation provides a sufficient accuracy in the time interval from $-2$
Myr to $+2$ Myr compared to the method of constructing the relative orbits
(the star relative to the Sun) in the Galactic potential. We choose the
epicyclic approximation, because we are planning to perform statistical
simulations with the computation of hundreds of thousands of orbits, while the
application of more complex methods requires a huge computational time.

###  2.2 Statistical Simulations

In accordance with the method of Monte Carlo statistical simulations, we
compute a set of orbits for each object by taking into account the random
errors in the input data. For each star, we compute the encounter parameter,
the shortest distance between the stellar and solar orbits
$d_{min}=\sqrt{\Delta X^{2}(t)+\Delta Y^{2}(t)+\Delta Z^{2}(t)}.$ The stellar
parameters are assumed to be distributed normally with a dispersion $\sigma$.
We add the errors to the equatorial coordinates, proper motion components,
parallax, and radial velocity of the star using the nominal errors.

##  3 RESULTS

First, we considered the solar neighborhood 30 pc in radius using $\approx$35
000 stars with known space velocities. For each star, we constructed the orbit
relative to the Sun in the time interval from $-2$ Myr to $+2$ Myr. In
addition to the stars from the list by Garcia-Sanchez et al. (2001), we found
several other Hipparcos stars. Data on the nine stars found with an encounter
parameter $d<3$ pc are presented in Table 1; Fig. 1a shows the trajectories
for six of these stars.

Figure 1b shows the model trajectories of four stars computed by taking into
account the random errors in the observational data; 300 realizations were
obtained for each of the three stars indicated in the figure.

Table 2 gives characteristics of the stars common to the list by Garcia-
Sanchez et al. (2001) that have the closest ($d<2$ pc) encounters with the
Sun.

Table 3 gives characteristics of the common list of stars that have the
closest ($d<2$ pc) encounters with the Sun. Note that this is currently the
most complete list with the restrictions that we used. In contrast to Tables 1
and 2, Table 3 gives the encounter parameters $d_{min}$ and $t_{min}$ with
their errors estimated through statistical simulations.

Table 3 also contains two stars from Table 1. The first of them is the well-
known white dwarf WD 0310.688 (HIP 14754). Note that the error in its radial
velocity (Pauli et al. 2006) is largest among all of the stars listed in
Tables 1 and 2. The data for the second star, HIP 27887, are very reliable. It
is contained in the catalog of F and G dwarfs of the Geneva-Copenhagen survey
(Nordstrom et al. 2004; Holmberg et al. 2007).

Note that the compilations by Gontcharov (2006) and Kharchenko et al. (2007)
give an obsolete radial velocity for HIP 3829, $V_{r}=263\pm 5$ km s-1, taken
from the GCRV (Wilson 1953). A very close encounter with the Sun is obtained
with this value: $d_{min}=0.96$ pc at $t_{min}=-16$ thousand years. However,
this star is the well-known white dwarf WD 0046+051 (van Maanen s star, vMa
2). According to more recent observations, its radial velocity is $V_{r}=54$
km s-1 (Greenstein and Trimble 1967) or $V_{r}(LSR)=-41$ km s-1 (Aannestad et
al. 1993). In this case, the encounter is considerably less close, $d_{min}>3$
pc (see also Garcia-Sanchez et al. 1999, 2001).

Figure 2 shows the trajectories of HIP 89825 (GL 710) relative to the Sun
computed by taking into account the random errors in the observational data.

We made only 300 realizations for a clearer graphical presentation. It is
obvious from the figure that very close encounters with the Solar system are
possible for this star. More detailed simulations show that, for example, out
of the one million model orbits found for GL 710, $d_{min}\leq 0.48$ pc in 855
902 cases and, hence, the probability of its falling into the Oort cloud is
$P_{1}=0.86$ and $d_{min}\leq 0.005$ pc in 94 cases and, hence, the
probability of the star approaching the boundary of its possible influence on
Kuiper Belt objects is $P_{2}=0.94\times 10^{-4}$.

Figure 3 show the distribution of encounter parameters of HIP 89825 (GL 710)
with the Sun obtained using one million model orbits.

##  4 DISCUSSION

(1) Among the Hipparcos stars we consider, there are nine stars common to the
list by Mullari and Orlov (1996). The encounter parameters for some of the
stars found by Mullari and Orlov (1996) are: $d_{min}=1.35$ pc, $t_{min}=1050$
thousand years for GL 710; $d_{min}=1.64$ pc, $t_{min}=-529$ thousand years
for GL 208; and $d_{min}=1.89$ pc, $t_{min}=88.7$ thousand years for GL 860 A.
The encounter parameters agree well with our results (Table2). There is only
one exception, GL 710, for which the new data give a closer encounter, which
was also pointed out by Garcia-Sanchez et al. (1999).

Note that the star HIP 99461 (GL 783) we selected is present in the list by
Mullari and Orlov (1996) but is absent in the lists by Garcia-Sanchez et al.
(1999, 2001).

(2) For an overwhelming majority of the stars, the encounter parameters
derived here agree well with those from Garcia-Sanchez et al. (2001).

The differences are most noticeable for the star GL 217.1, more specifically,
$|\Delta d_{min}|=0.38$ pc and $|\Delta t_{min}|=182$ thousand years.
According to Garcia-Sanchez et al. (2001), who adopted the radial velocity
$V_{r}=20.0\pm 3.7$ km s-1 that differs markedly from our value (Table 2), the
encounter parameters for GL 217.1 are: $d_{min}=1.65\pm 0.27$ pc and
$t_{min}=-(1046\pm 163)$ thousand years.

Computing the model epicyclic orbits for the star GL 217.1 using the initial
data from Garcia-Sanchez et al. (2001) with the adopted radial velocity error
$e_{V_{r}}=2$ km s-1 yields $d_{min}=1.61\pm 0.20$ pc and $t_{min}=-(1066\pm
116)$ thousand years, and these parameters almost coincide with those from
Garcia-Sanchez et al. (2001).

This leads us to conclude that the existing discrepancy is undoubtedly related
to the difference in input data.

As we see from Table 3, GL 217.1 is the most massive, $M=2M_{\odot}$, star
among our stars. Its improved encounter parameters, $d_{min}=1.28\pm 0.06$ pc
and $t_{min}=-(861\pm 40)$ thousand years, make it attractive for studying the
close passages that could take place in the past.

(3) The star GL 710 is of greatest interest to us, because it has a record
close encounter with the Sun among all of the candidates known to date.
According to Garcia-Sanchez et al. (2001), the encounter parameters of GL 710
are: $d_{min}=0.34\pm 0.18$ pc and $t_{min}=(1358\pm 41)$ thousand years, in
good agreement with those we found (Table 3).

GL 710 is the only star with a significant difference between the data of
Tables 2 and 3. This suggests that the conditions for its encounter with the
Sun make the stellar orbit very sensitive to small changes in such initial
data as the parallax and radial velocity. We clearly see from Fig. 3a that the
distribution of minimum encounter distance $d_{min}$ differs from a Gaussian
one, the mode of this distribution is $d_{min}=0.27$ pc, which is lower than
the median value of $d_{min}$ given in Table 3, while the distribution of
encounter times is nearly Gaussian (Fig. 3b).

(4) Our simulations show that among the candidates considered, only the star
GL 710 has a high probability of penetrating into the Oort cloud region.

##  5 CONCLUSIONS

Based on currently available space velocity data for about 35 000 Hipparcos
stars, we searched for stars closely encountering with the Solar system. For
this purpose, we took stars within 30 pc of the Sun and determined their
orbits relative to the Sun based of the epicyclic approximation in the time
interval from $-2$ Myr to $+2$ Myr. We found nine new candidates in addition
to the well-known list of such stars (Garcia-Sanchez et al. 2001). The use of
the PCRV (Gontcharov 2006) and improved Hipparcos stellar parallaxes (Leeuwen
2007) made this possible.

As a result, we compiled the currently most complete (given the restrictions)
list of Hipparcos stars that have close ($d<2$ pc) encounters with the Sun.

For the star GL 217.1, a well-known candidate for a passage close to the Sun,
the new observational data were shown to change noticeably its previously
known encounter parameters with the Sun. The encounter parameters found here
are: $d_{min}=1.28\pm 0.06$ pc and $t_{min}=-(861\pm 40)$ thousand years.
Improving the radial velocity for the white dwarf WD 0310–688 (HIP 14754)
whose orbit passed at a distance $d_{min}=1.61\pm 0.19$ pc from the solar
orbit about 300 thousand years ago is of current interest. Our statistical
simulations showed that the star GL 710 has not only a high probability of
penetrating into the Oort cloud, $P_{1}=0.86,$ but also a nonzero probability,
$P_{2}=1\times 10^{-4},$ of penetrating into the region where the influence of
the passing star on Kuiper Belt objects is significant.

ACKNOWLEDGMENTS

I wish to thank Prof. V.V. Orlov for a careful reading of the manuscript and
his remarks and A.T. Bajkova for a discussion of the results. The SIMBAD
searchable database was very helpful in the work. This study was supported by
the Russian Foundation for Basic Research (project no. 08–02–00400) and in
part by the “Origin and Evolution of Stars and Galaxies” Program of the
Presidium of the Russian Academy of Sciences.

REFERENCES

Aannestad P.A., S.J. Kenyon, G.L. Hammon, et al., Astron. J. 105, 1033 (1993).

Anosova J., V.V. Orlov, and N.A. Pavlova, Astron. Astrophys. 292, 115 (1994).

Bobylev V.V., A.T. Bajkova, and A.S. Stepanishchev, Astron. Lett. 34, 515
(2008).

Bobylev V.V., G.A. Gontcharov, and A.T. Bajkova, Astron. Rep. 50, 733 (2006).

Dehnen W., and J.J. Binney, MNRAS 298, 387 (1998).

Fuchs B., D. Breitschwerdt, M.A. Avilez, et al., MNRAS 373, 993 (2006).

Garcia-Sanchez J., R.A. Preston, D.L.Jones, et al., Astron. J. 117, 1042
(1999).

Garcia-Sanchez J., P.R.Weissman, R.A.Preston, et al., Astron. Astrophys. 379,
634 (2001).

Giampapa M.S., R. Rosner, V. Kashyap, et al., Astrophys. J. 463, 707 (1996).

Gontcharov G.A., Astron. Lett. 32, 759 (2006).

Greenstein J.L., and V.L. Trimble, Astrophys. J. 149, 283 (1967).

Hills J.G., Astron. J. 86, 1730 (1981).

The Hipparcos and Tycho Catalogues, ESA SP-1200 (1997).

Holmberg J., and C. Flinn, MNRAS 352, 440 (2004).

Holmberg J., B. Nordstrom, and J. Andersen, Astron. Astrophys. 475, 519
(2007).

Joshi Y.C., MNRAS 378, 768 (2007).

Kenyon S.J., and B. Bromley, Nature 432, 598 (2004).

Kharchenko N.V., R.-D. Scholz, A.E. Piskunov, et al., Astron. Nachr. 328,
(2007).

van Leeuwen F., Astron. Astrophys. 474, 653 (2007).

Lindblad B., Arkiv Mat., Astron., Fysik A 20, (17), (1927).

Matthews R.A.J., R. Astron. Soc. Quart. J. 35,1 (1994).

Mullari A.A., and V.V. Orlov, Earth, Moon, and Planets (Kluwer Netherlands,
1996), Vol. 72, p. 19.

Nordstrom B., M. Mayor, J. Andersen, et al., Astron. Astrophys. 419, 989
(2004).

Pauli E.M., R. Napiwotzki, U. Heber, et al., Astron. Astrophys. 447, 173
(2006).

Revina I.A., Analysis of the Motion of Celestial Bodies and Estimation of the
Accuracy of Their Observations (Latvian Univ., Riga, 1988), p. 121.

Weissman P.R., Earth, Moon, and Planets (Kluwer, Netherlands, 1996), Vol. 72,
p. 25.

Wilson R.E., General Catalogue of Stellar Radial Velocities (Carnegie Inst.,
Washington, DC, 1953), p. 601.

Translated by N. Samus’

Table 1: Hipparcos stars within 30 pc of the Sun found here

HIP | SP | $d_{min},$ | $t_{min},$ | $\pi\pm e_{\pi},$ | $V_{r}\pm e_{V_{r}},$ | n | $\varepsilon_{V_{r}},$
---|---|---|---|---|---|---|---
|  | pc |  $10^{3}$ yr | mas | km s-1 |  | km s-1
1 | 2 |  3 | 4 | 5 | 6 | 7 | 8
14754 | DA | $1.60$ | $-296$ | $97.66\pm 1.85$ | $33.8\pm 3.2$ | 1 |
27887 | K3V | $1.97$ | $-413$ | $76.95\pm 0.37$ | $30.7\pm 0.3$ | 1 |
99461* | K2V | $2.05$ | $41.1$ | $166.25\pm 0.27$ | $-129.0\pm 0.2$ | 2 | 0.8
24186 | M0V | $2.15$ | $-10.8$ | $255.66\pm 0.91$ | $244.4\pm 0.2$ | 3 | 1.1
105766 | G5 | $2.20$ | $341$ | $37.87\pm 0.42$ | $-76.6\pm 0.3$ | 2 | 0.2
30344 | K0V | $2.30$ | $-2035$ | $34.10\pm 0.62$ | $14.6\pm 0.3$ | 1 |
26373 | K0V | $2.33$ | $-762$ | $39.82\pm 1.36$ | $32.5\pm 0.3$ | 1 |
104214 | K5V | $2.79$ | $18.6$ | $286.82\pm 6.78$ | $-65.9\pm 0.1$ | 5 | 0.7
104217 | K7V | $2.80$ | $19.6$ | $285.88\pm 0.54$ | $-64.3\pm 0.1$ | 3 | 0.3

Note. $n$ is the number of catalogs used to calculate the mean $V_{r};$
$\varepsilon_{V_{r}}$ is the radial velocity error calculated from the
external convergence (Gontcharov 2006); (*) HIP 99461 (GL 783) is present in
the list by Mullari and Orlov (1996) but is absent in the lists by Garcia-
Sanchez et al. (1999, 2001).

Table 2: Data on the known Hipparcos stars from the solar neighborhood 30 pc
in radius encountering with the solar orbit within $d<2$ pc

HIP |  | $d_{min},$ | $t_{min},$ | $\pi\pm e_{\pi},$ | $V_{r}\pm e_{V_{r}},$ | n | $\varepsilon_{V_{r}},$
---|---|---|---|---|---|---|---
|  | pc |  $10^{3}$ yr | mas | km s-1 |  | km s-1
1 |  2 | 3 | 4 | 5 | 6 | 7 | 8
89825 | GL 710 | $0.21$ | $1444$ | $51.12\pm 1.63$ | $-13.8\pm 0.3$ | 4 | 8.1
70890 | Proxima Cen | $0.89$ | $27.4$ | $771.64\pm 2.60$ | $-25.1\pm 0.9$ |  |
71683 |  $\alpha$ Cen A/B * | $0.91$ | $28.4$ | $754.81\pm 4.11$ | $-24.7\pm 0.4$ | 1 |
57544 | AC+79 3888 | $1.06$ | $46.0$ | $186.86\pm 1.70$ | $-111.6\pm 0.2$ | 5 | 3.9
87937 | Barnard’s star | $1.15$ | $9.8$ | $548.31\pm 1.51$ | $-110.6\pm 0.2$ | 5 | 2.2
27288 | GL 217.1 | $1.27$ | $-862$ | $46.28\pm 0.16$ | $24.7\pm 1.1$ | 2 | 5.9
54035 | Lalande 21185 | $1.43$ | $20.5$ | $392.64\pm 0.67$ | $-85.8\pm 0.1$ | 7 | 3.1
26335 | GL 208 | $1.54$ | $-500$ | $88.97\pm 1.02$ | $22.0\pm 0.2$ | 5 | 2.1
57548 | Ross 128 | $1.92$ | $72.4$ | $298.04\pm 2.30$ | $-31.0\pm 0.2$ | 4 | 1.1
110893 | GL 860 A | $1.93$ | $90.0$ | $249.94\pm 1.87$ | $-34.0\pm 0.1$ | 6 | 2.3
92403 | Ross 154 | $1.96$ | $157$ | $336.72\pm 2.03$ | $-10.7\pm 0.2$ | 5 | 1.9
38228 | HD 63433 | $2.02$ | $1366$ | $45.45\pm 0.53$ | $-16.1\pm 0.1$ | 6 | 1.2

Note. $n$ is the number of catalogs used to calculate the mean $V_{r};$
$\varepsilon_{V_{r}}$ is the radial velocity error calculated from the
external convergence (Gontcharov 2006); (*) the encounter parameters were
calculated from component A.

Table 3: Hipparcos stars from the solar neighborhood 30 pc in radius
encountering with the solar orbit within $d<2$ pc

HIP |  | SP | $M/M_{\odot}$ | Ref |  ${\overline{d}}_{min},$ $10^{3}$ yr |  ${\overline{t}}_{min},$ pc
---|---|---|---|---|---|---
1 |  2 | 3 | 4 | 5 | 6 | 7
89825 | GL 710 | K7V | $0.4-0.6$ | (1) | $0.311\pm 0.167$ | $1447\pm 60$
70890 | Proxima Cen | M5Ve | $0.15\pm 0.02$ | (2) | $0.890\pm 0.019$ | $27.4\pm 0.1$
71683 |  $\alpha$ Cen A/B * | G2V/K1V | $2.0\pm 0.1$ | (2) | $0.910\pm 0.012$ | $28.4\pm 0.1$
57544 | AC+79 3888 | M4 | $\approx 0.15$ |  | $1.059\pm 0.013$ | $46.0\pm 0.3$
87937 | Barnard’s star | sdM4 | $0.144$ | (3) | $1.148\pm 0.006$ | $9.8\pm 0.1$
27288 | GL 217.1 | A2Vann | $2.0$ | (1) | $1.275\pm 0.064$ | $-861\pm 40$
54035 | Lalande 21185 | M2V | $0.39$ | (4) | $1.426\pm 0.005$ | $20.5\pm 0.1$
26335 | GL 208 | K7 | $0.47$ | (4) | $1.537\pm 0.050$ | $-500\pm 8$
14754 | WD 0310-688 | DA | $\approx 0.5$ |  | $1.607\pm 0.190$ | $-297\pm 29$
57548 | Ross 128 | M4.5V | $0.24$ | (4) | $1.920\pm 0.028$ | $72.8\pm 0.7$
110893 | GL 860 A | M2V | $0.56$ | (4) | $1.929\pm 0.028$ | $90.3\pm 0.5$
92403 | Ross 154 | M3.5Ve | $0.23$ | (4) | $1.959\pm 0.030$ | $157\pm 1$
27887 | HD 40307 | K3V | $0.75^{+0.78}_{-0.71}$ | (5) | $1.974\pm 0.035$ | $-413\pm 5$
38228 | HD 63433 | G5IV | $0.88^{+0.98}_{-0.84}$ | (5) | $2.038\pm 0.099$ | $1366\pm 20$

Note. The stellar mass estimates were taken from the following papers: (1)
Garcia-Sanchez et al. (1999); Anosova et al. (1994); (3) Giampapa et al.
(1996); (4) Mullari and Orlov (1996); and (5) Holmberg et al. (2007). If there
is no reference, then a typical mass is given; (*) the sum of the masses is
given.

Fig. 1. (a) Stellar trajectories relative to the Sun. The trajectories of the
stars found here are highlighted by the thick lines: HIP 27887 (1), HIP 14754
(2), HIP 24186 (3), the pair of HIP 104214 and HIP 104217 (4), HIP 99461 (5),
and HIP 105766 (6). (b) The model trajectories of four stars computed by
taking into account the random errors in the observational data (300
realizations). The trajectories hatch the $3\sigma$ confidence regions; the
Oort cloud region is shaded.

Fig. 2. Model trajectories of the star HIP 89825 (GL 710) relative to the Sun
computed by taking into account the random errors in the observational data
(300 realizations). The trajectories hatch the $3\sigma$ confidence region;
the Oort cloud region is shaded.

Fig. 3. Model encounter parameters of the star HIP 89825 (GL 710) with the
Sun. (a) The distribution of minimum distance $d_{min}$; (b) the histogram of
encounter times $t_{min}$.