Léo Gautheron

Pascal Germain
email: pascal.germain@ift.ulaval.ca

Amaury Habrard

Guillaume Metzler

Emilie Morvant

Marc Sebban

Valentina Zantedeschi
email: vzantedeschi@gmail.com

Landmark-based Ensemble Learning with Random Fourier Features and Gradient Boosting

Keywords: Gradient boosting, Random Fourier features, Kernel learning

This paper jointly leverages two state-of-the-art learning strategies-gradient boosting (GB) and kernel Random Fourier Features (RFF)-to address the problem of kernel learning. Our study builds on a recent result showing that one can learn a distribution over the RFF to produce a new kernel suited for the task at hand. For learning this distribution, we exploit a GB scheme expressed as ensembles of RFF weak learners, each of them being a kernel function designed to fit the residual. Unlike Multiple Kernel Learning techniques that make use of a pre-computed dictionary of kernel functions to select from, at each iteration we fit a kernel by approximating it from the training data as a weighted sum of RFF. This strategy allows one to build a classifier based on a small ensemble of learned kernel "landmarks" better suited for the underlying application. We conduct a thorough experimental analysis to highlight the advantages of our method compared to both boosting-based and kernel-learning state-of-the-art methods.

Introduction

Kernel methods are among the most popular approaches in machine learning due to their capability to address non-linear problems, their robustness and their simplicity. However, they exhibit two main flaws in terms of memory usage and time complexity. Landmark-based kernel approaches [START_REF] Balcan | Improved guarantees for learning via similarity functions[END_REF] can be used to drastically reduce the number of instances involved in the comparisons, but they heavily depend on the choice and the parameterization of the kernel. Multiple Kernel Learning [START_REF] Wu | Boosting based multiple kernel learning and transfer regression for electricity load forecasting[END_REF] and Matching Pursuit methods [START_REF] Vincent | Kernel matching pursuit[END_REF] can provide alternative solutions to this problem but these require the use of a pre-defined dictionary of base functions. Another strategy to improve the scalability of kernel methods is to use approximation techniques such as the Nyström [START_REF] Drineas | On the Nyström method for approximating a gram matrix for improved kernel-based learning[END_REF] or Random Fourier Features (RFF) [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. The latter is probably the most used thanks to its simplicity and ease of computation. It allows the approximation of any shift-invariant kernel based on the Fourier transform of the kernel. Several works have extended this technique by allowing one to adapt the RFF approximation directly from the training data [START_REF] Agrawal | Data-dependent compression of random features for large-scale kernel approximation[END_REF][START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF][START_REF] Sinha | Learning kernels with random features[END_REF]. Among them, the recent work of Letarte et al. [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF] introduces a method to obtain a weighting distribution over the random features by a single pass over them. This strategy is derived from a statistical learning analysis, starting from the observation that each random feature can be interpreted as a weak hypothesis in the form of trigonometric functions obtained by the Fourier decomposition. However, in practice, this method requires the use of a fixed set of landmarks selected beforehand and independently from the task before being able to learn the representation in a second step. This leads to three important limitations: (i) the need for a heuristic strategy for selecting relevant landmarks, (ii) these latter and the associated representation might not be adapted for the underlying task, and (iii) the number of landmarks might not be minimal w.r.t. that task, inducing higher computational and memory costs.

We propose in this paper to tackle these issues with a gradient boosting approach [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF]. Our aim is to learn iteratively the classifier and a compact and efficient representation at the same time. Our greedy optimization method is similar to Oglic & Gärtner's one [START_REF] Oglic | Greedy feature construction[END_REF], which at each iteration of the functional gradient descent [START_REF] Mason | Functional gradient techniques for combining hypotheses[END_REF] refines the representation by adding the base function minimizing a residual-based loss function. But unlike our approach, their method does not allow to learn a classifier at the same time. Instead, we propose to jointly optimize the classifier and the base functions in the form of kernels by leveraging both gradient boosting and RFF. Interestingly, we further show that we can benefit from a significant performance boost by (i) considering each weak learner as a single trigonometric feature, and (ii) learning the random part of the RFF. Organization of the paper. Section 2 describes the notations and the necessary background knowledge. We present our method in Section 3 as well as two efficient refinements before presenting an extensive experimental study in Section 4, comparing our strategy with boosting-based and kernel learning methods.

Notations and Related Work

We consider binary classification tasks from a d-dimensional input space R d to a label set Y ={-1, +1}. Let S= (x i , y i ) n i=1 be a training set of n points. We focus on kernel-based algorithms that rely on pre-defined kernel functions k : R d ×R d → R assessing the similarity between any two points of the input space. These methods present a good performance when the parameters of the kernels are learned and the chosen kernels are able to fit the distribution of the data. However, selecting the right kernel and tuning its parameters is computationally expensive, in general. To reduce this overhead, one can resort to Multiple Kernel Learning techniques [START_REF] Wu | Boosting based multiple kernel learning and transfer regression for electricity load forecasting[END_REF] which boil down to selecting the combination of kernels that fits the best the training data: a dictionary of T base functions {k t } T is composed of various kernels associated with some fixed parameters, and a combination is learned, defined as

H(x, x ) = T t=1 α t k t (x, x ), (1)
with α t ∈R the weight of the kernel k t (x, x ). As shown in Section 3, our main contribution is to address this issue of optimizing a linear combination of kernels by leveraging RFF and gradient boosting (we recall basics on it in Section 3.1).

To avoid the dictionary of kernel functions in Equation ( 1) from being precomputed, we propose a method inspired from Letarte et al. [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF] to learn a set of approximations of kernels tailored to the underlying classification task. Unlike Letarte et al., we learn such functions so that the representation and the classifier are jointly optimized. We consider landmark-based shift-invariant kernels relying on the value δ = x t -x ∈ R d and usually denoted by abuse of notation by

k(δ) = k(x t -x) = k(x t , x)
, where x t ∈ R d is a point-called landmark-lying on the input space which all the instances are compared to, and that strongly characterizes the kernel. At each iteration of our gradient boosting procedure, we optimize the kernel function itself, exploiting the flexibility of the framework of Letarte et al., where a kernel is a weighted sum of RFF [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] defined as

k q t (x t -x) = K j=1 q t j cos ω j • (x t -x) , (2)
where the ω j are drawn from the Fourier transform of a shift-invariant kernel k denoted by p(ω) and defined as

p(ω) = 1 (2π) d R d k(δ)e -iω•δ dδ. (3)
When q t is uniform, we retrieve the setting of RFF and we have k(δ) k q t (δ) where larger number of random features K give better approximations [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. Letarte et al. [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF] aim to learn the weights of the random Fourier features q t . To do so, they consider a loss function that measures the quality of the similarities computed using the kernel k q t . Their theoretical study on leads to a closed-form solution for q t computed as ∀j ∈ {1, . . . , K}, q t j =

1 Z t exp -β √ n n n i=1 (h t ω j (x i )) , (4)
with β ≥ 0 a parameter to tune, h t ω (x) = cos (ω•(x t -x)), and Z t a normalization constant such that K j=1 q t j =1. They learn a representation of the input space of n L features where each of them is computed using k q t with the landmark (x t , y t ) selected randomly from the training set. Once the new representation is computed, a (linear) predictor is learned from it, in a second step. It is worth noticing that this kind of procedure exhibits two limitations. First, the model can be optimized only after having learned the representation. Second, the landmarks have to be fixed before learning the representation. Thus, the constructed representation is not guaranteed to be compact and relevant for the learning algorithm considered. To tackle these issues, we propose in the following a strategy that performs both steps at the same time through a gradient boosting process that allows to jointly learn the set of landmarks and the final predictor.

Gradient Boosting Random Fourier Features

The approach we propose follows the widely used gradient boosting framework first introduced by Friedman [START_REF] Friedman | Greedy function approximation: a gradient boosting machine[END_REF]. We briefly recall it below.

Gradient Boosting in a Nutshell

Gradient boosting is an ensemble method that aims at learning a weighted majority vote over an ensemble of T weak predictors in a greedy way by learning one classifier per iteration. The final majority vote is of the form

∀x ∈ R d , sign H 0 (x) + T t=1 α t h a t (x) ,
where H 0 is an initial classifier fixed before the iterative process (usually set such that it returns the same value for every sample), and α t is the weight associated to the predictor h a t and is learned at the same time as the parameters a t of that classifier. Given a differentiable loss , the objective of the gradient boosting algorithm is to perform a gradient descent where the variable to be optimized is the ensemble and the function to be minimized is the empirical loss. The pseudo-code of gradient boosting is reported in Algorithm 1. First, the ensemble is constituted by only one predictor: the one that outputs a constant value minimizing the loss over the whole training set (line 1). Then at each

(x) + T t=1 α t K j=1 q t j cos ω t j • (x t -x) 1: H 0 ← H 0 (xi) = 1 2 ln 1+ 1 n n j=1 y j 1- 1 n n j=1 y j 2: for t = 1, . . . , T do 3: ∀i = 1, . . . , n, wi = exp(-yiH t-1 (xi)) 4: ∀i = 1, . . . , n, ỹi = yiwi 5: ∀j = 1, . . . , K, draw ω t j ∼ N (0, 2γ) d 6:
x t = argmin

x∈R d 1 n n i=1 exp -ỹi 1 K K j=1 cos(ω t j • (x -xi)) 7: ∀j = 1, . . . , K, q t j = 1 Z t exp -β √ n n n i=1 exp -ỹi cos ω t j • (x t -xi) 8: α t = 1 2 ln n i=1 1+y i K j=1 q t j cos ω t j •(x t -x i ) w i n i=1 1-y i K j=1 q t j cos ω t j •(x t -x i ) w i 9: ∀i = 1, . . . , n, H t (xi) = H t-1 (xi) + α t K j=1 q t j cos ω t j • (x t -xi) 10: end for
iteration, the algorithm computes for each training example the negative gradient of the loss (line 3), also called the residual and denoted by ỹi . The next step consists in optimizing the parameters of the predictor h a t that fits the best the residuals (line 4), before learning the optimal step size α t that minimizes the loss by adding h a t , weighted by α t , to the current vote (line 5). Finally, the model is updated by adding α t h a t (•) (line 6) to the vote.

Gradient Boosting with Random Fourier Features

Our main contribution takes the form of a learning algorithm which jointly optimizes a compact representation of the data and the model. Our method, called GBRFF1, leverages both Gradient Boosting and RFF. We describe its pseudo-code in Algorithm 2 which follows the steps of Algorithm 1. The loss function at the core of our algorithm is the exponential loss:

H T = 1 n n i=1 exp -y i H T (x i ) . (5)
Given H T , line 1 of Algorithm 1 amounts to setting the initial learner as

∀i ∈ {1, . . . , n}, H 0 (x i ) = 1 2 ln 1 + 1 n n j=1 y j 1 -1 n n j=1 y j . ( 6
)
The residuals of line 3 are defined as ỹi = -

∂ y i , H t-1 (x i ) ∂H t-1 (x i ) = y i e -yiH t-1 (xi) .
Line 4 of Algorithm 1 tends to learn a weak learner that outputs exactly the residuals' values by minimizing the squared loss; but, this is not well suited in our setting with the exponential loss (Equation ( 5)). To benefit from the exponential decrease of the loss, we are rather interested in weak learners that output predictions having a large absolute value and being of the same sign as the residuals. Thus, we aim at favoring parameter values minimizing the exponential loss between the residuals and the predictions of the weak learner as follows:

a t = argmin a 1 n n i=1 exp -ỹi h a (x i ) . (7)
Following the RFF principle, we can now define our weak learner as

h a t (x i ) = K j=1 q t j cos ω t j • (x t -x i ) , (8)
where its parameters are given by a t =({ω t j } K j=1 , x t , q t ). Instead of using a predefined set of landmarks [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF], we build this set iteratively, i.e., we learn one landmark per iteration. To benefit from the closed form of Equation ( 4), we propose the following greedy approach to learn the parameters a t . At each iteration t, we draw K vectors {ω t j } K j=1 ∼p K with p the Fourier transform of a given kernel (as defined in Equation ( 3)); then we look for the optimal landmark x t . Plugging Equation ( 8) into Equation ( 7) and assuming a uniform prior distribution over the random features, x t is learned to minimize

x t = argmin x∈R d f (x) = 1 n n i=1 exp -ỹi 1 K K j=1 cos(ω t j • (x -x i )) . (9)
Even if this problem is non-convex due to the cosine function, we can still compute its derivative and perform a gradient descent to find a possible solution. The partial derivative of Equation ( 9) with respect to x is given by

∂f ∂x (x)= 1 Kn n i=1   ỹi K K j=1 sin(ω t j • (x-x i ))   exp   - ỹi K K j=1 cos(ω t j • (x-x i ))   K j=1 ω t j .
According to Letarte et al. [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF], given the landmark x t found by gradient descent, we can now compute the weights of the random features q t as ∀j ∈ {1, . . . , K}, q t j =

1 Z t exp -β √ n n n i=1 exp -ỹi cos ω t j • (x t -x i ) , (10)
with β ≥ 0 a parameter to tune and Z t the normalization constant.

The last step concerns the step size α t . It is computed so as to minimize the combination of the current model H t-1 with the weak learner h t , i.e.,

α t = argmin α n i=1 exp -y i (H t-1 (x i )+αh t (x i )) = argmin α n i=1 w i exp -y i αh t (x i ) ,
where w i = exp(-y i H t-1 (x i )). In order to have a closed-form solution of α, we use the convexity of the above quantity and the fact that h t (x i ) ∈ [-1, 1] to bound the loss function to optimize. Indeed, we get

n i=1 w i e -yiαh t (xi) ≤ n i=1 1-y i h t (x i ) 2 w i e α + n i=1 1+y i h t (x i ) 2 w i e -α .
This upper bound is strictly convex. Its minimum α t can be found by setting to 0 the derivative w.r.t. α of the right-hand side of the previous equation. We get

n i=1 1 -y i h t (x i ) 2 w i e α = n i=1 1 + y i h t (x i ) 2 w i e -α ,
for which the solution is given by

α t = 1 2 ln n i=1 (1 -y i h t (x i ))w i n i=1 (1 + y i h t (x i ))w i .
The same derivation can be used to find the initial predictor H 0 . As usually done in the RFF literature [START_REF] Agrawal | Data-dependent compression of random features for large-scale kernel approximation[END_REF][START_REF] Rahimi | Random features for large-scale kernel machines[END_REF][START_REF] Sinha | Learning kernels with random features[END_REF] we use the RBF kernel k γ (x, x )=e -γ x-x 2 with as Fourier transform vectors of d numbers each drawn from the normal law with zero mean and variance 2γ that we denote N (0, 2γ) d .

Refining GBRFF1

In GBRFF1, the number of random features K used at each iteration has a direct impact on the computation time of the algorithm. Moreover ω t is drawn according to the Fourier transform of the RBF kernel and thus is not learned. The second part of our contribution is to propose two refinements. First, we bring to light the fact that one can drastically reduce the complexity of GBRFF1 by learning a rough approximation of the kernel, yet much simpler and still very effective, using K=1. In this scenario, we show that learning the landmarks boils down to finding a single real number in [-π, π]. Then, to speed up the convergence of the algorithm, we suggest to optimize ω t after a random initialization from the Fourier transform. We show that a simple gradient descent with respect to this parameter allows a faster convergence with better performance. These two improvements lead to a variant of our original algorithm, called GBRFF2 and presented in Algorithm 3.

Cheaper landmark learning using the periodicity of the cosine. As we set K=1, the weak learner h a t (x) is now simply defined as Draw ω ∼ N (0, 2γ) d 6:

h a t (x) = cos ω t • (x t -x i ) ,
b t = argmin b∈[-π,π] 1 n n i=1 exp -ỹi cos ω • xi -b) 7: ω t = argmin ω∈R d λ ω 2 2 + 1 n n i=1 exp -ỹi cos ω • xi -b t ) .
8:

α t = 1 2 ln n i=1 1+y i cos ω t •x i -b t w i n i=1 1-y i cos ω t •x i -b t w i 9: ∀i = 1, . . . , n, H t (xi) = H t-1 (xi) + α t cos ω t • xi -b t 10: end for
where its parameters are given by a t = (ω t , x t ). This formulation allows us to eliminate the dependence on the hyper-parameter K. Moreover, one can also get rid of β, because learning the weights q t j (line 7 of Algorithm 2) is no more necessary. Instead, since K=1, we can see α t learned at each iteration as a surrogate of these weights. As our weak learner is based on a single random feature, the objective function (line 6) to learn the landmark at iteration t becomes

x t = argmin x∈R d f ω t (x) = 1 n n i=1 exp -ỹi cos(ω t • (x -x i )) .
Let c ∈ 1, d be the index of the c-th coordinate of the landmark x t . We can rewrite the objective function as

f ω t (x t ) = 1 n n i=1 e -ỹi cos(ω t •x t -ω t •xi) = 1 n n i=1 e -ỹi cos(ω t c x t c + j =c ω t j x t j -ω t •xi) .
We leverage the periodicity of the cosine function along each direction to find the optimal c-th coordinate of the landmark

x t c ∈ [ -π ω t c , π ω t c
] that minimizes f ω t (x t ) by fixing all the other coordinates. Figure 1 illustrates this phenomenon on the twomoons dataset when applying GBRFF1 with K=1. The plots in the first row show the periodicity of the loss represented as repeating diagonal green/yellow stripes (light yellow is associated to the smallest loss). There is an infinite number of landmarks giving such a minimal loss at the middle of the yellow stripes. Thus, by setting one coordinate of the landmark to an arbitrary value, the algorithm is still able at any iteration to find along the second coordinate a value that minimizes the loss (the resulting landmark at the current iteration is depicted by a white cross). The second row shows that such a strategy allows us to get an accuracy of 100% on this toy dataset after 10 iterations. By generalizing this, instead of learning a landmark vector x t ∈ R d , we fix all but one coordinate of the landmark to 0, and then learn a single scalar b t ∈ [-π, π] that minimizes

f ω t (b t ) = 1 n n i=1 exp -ỹ i cos ω t • x i -b t .
Learning ω t for faster convergence. The second refinement concerns the randomness of the RFF due to vector ω t . So far, the latter was drawn according p and then used to learn b t . We suggest instead to fine-tune ω t by doing a gradient descent with as initialization the vector drawn from p. Supported by the experiments performed in the following, we claim that such a strategy allows us to both speed up the convergence of the algorithm and boost the accuracy. This update requires to add a line of code, just after line 6 of Algorithm 2, expressed as a regularized optimization problem:

ω t = argmin ω∈R d λ ω 2 2 + 1 n n i=1 exp -ỹi cos ω • x i -b t ) , its derivative being ∂f ω ∂ω (ω) = 2λω + 1 n n i=1 x i ỹi sin(ω•x i -b t ) e -ỹi cos(ω•xi-b t ) .

Experimental Evaluation

The objective of this section is three-fold: first, we aim to bring to light the interest of learning the landmarks rather than fixing them as done in Letarte et al. [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF]; second we study the impact of the number K of random features; lastly, we perform an extensive experimental comparison of our algorithms. The Python code of all experiments and the data used are publicly available4 .

Setting

For GBRFF1 and GBRFF2, we select by cross-validation (CV) the hyperparameter γ ∈ 2 {-2,...,2}

d

. For GBRFF2, we also tune λ ∈ {0, 2 {-5,...,-2} }. We compare our two methods with the following algorithms.

• LGBM [START_REF] Ke | Lightgbm: A highly efficient gradient boosting decision tree[END_REF] is a state-of-the-art gradient boosting method using trees as base predictors. We select by CV the maximum tree depth in {1, . . . , 10} and the L2 regularization parameter λ ∈ {0, 2 {-5,...,-2} }.

• BMKR [START_REF] Wu | Boosting based multiple kernel learning and transfer regression for electricity load forecasting[END_REF] is a Multiple Kernel Learning method based on gradient boosting with least square loss. It selects at each iteration the best kernel plugged inside an SVR to fit the residuals among 10 RBF kernels with γ ∈ 2 {-4,...,5} and the linear kernel k(x, x ) = x x . We select by CV the SVR parameter C ∈ 10 {-2,...,2} .

• GFC [START_REF] Oglic | Greedy feature construction[END_REF] is a greedy feature construction method based on functional gradient descent. It iteratively refines the representation learned by adding a feature that matches the residual function defined for the least squared loss. We use the final representation to learn a linear SVM where C ∈ 10 {-2,...,2} is selected by CV.

• PBRFF Fig. 2. Mean test accuracy over 20 train/test splits over the 16 datasets. We train the three methods using from 1 to 50 landmarks.

We consider 16 datasets coming mainly from the UCI repository that we binarized as described in Table 1. We generate for each dataset 20 random 70%/30% train/test splits. Datasets are pre-processed such that each feature in the training set has 0 mean and unit variance; the factors computed on the training set are then used to scale each feature in the test set. All parameters are tuned by 5-fold CV on the training set by performing a grid search.

Influence of Learning the Landmarks

We present in Figure 2 the behavior of the three methods that make use of landmarks and RFF, that is PBRFF, GBRFF1 and GBRFF2. With more than 25 landmarks, PBRFF and GBRFF1 show similar mean accuracy and reach about 87.5% after 50 iterations. However, for a small set of landmarks (in particular smaller than 25) GBRFF1 is consistently superior by about 1 point higher than PBRFF, showing the interest of learning the landmarks. But the certainly most striking result comes from the performance of our variant GBRFF2 which outperforms the two competing methods. This is particularly true for a small amount of landmarks. Notice that GBRFF2 is able to reach its maximum with about 20 landmarks, while GBRFF1 and PBRFF require more iterations without reaching the same performance. This definitely shows the benefit of learning the random features compared to drawing them randomly.

Influence of the Number of Random Features

A key parameter of GBRFF1 is K, the number of random features used at each iteration. To highlight its impact, we report in Figure 3 the mean test accuracy of GBRFF1 with K ∈ {1, 5, 10, 20} across all datasets and over the 20 train/test splits. To have a fair study, the comparison is performed according to the same total number of random features after the whole boosting process, that is T ×K with T the number of iterations. First of all, we observe that with a total of 1, 000 random features, K does not have a big impact on the performance. However, when decreasing the value of T × K, it becomes much more interesting in terms of accuracy to set K to a small value. This shows that the more we want a compact final representation, the more we need to refine the random features: it is better to weight each of the features greedily with α t (line 8 of Algorithm 3) rather than using the closed-form solution of Equation ( 10) (line 7 of Algorithm 2) to weight them all at once. Even if in the usual context of RFF it is desirable to have a large K value to approximate a kernel, this series of experiments shows that a simple rough approximation with K=1 along with a sufficient number of iterations allows the final ensemble to mimic the approximation of a new kernel suited for the task at hand.

Influence of the Number of Samples on the Computation Time

The specificities of GBRFF2 come from the number of random features K set to 1 at each iteration and the learning of ω t . We have already shown in Figure 2 that this allows us to get better results. We study in this section how GBRFF2 scales compared to the other methods. To do so, we consider artificial datasets with an increasing number of samples (generated with scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in python[END_REF] library's make classification function). The initial size is set to 150 samples, and we successively generate datasets with a size equal to the previous dataset size multiplied by 1.5. Here, we do not split the datasets in train and test as we are not interested in the accuracy. We report the time in seconds necessary to train the models and to predict the labels on the whole datasets. The parameters are fixed as follows: C = 1 for the methods using SVM or SVR; the tree depth is set to 5 for LGBM; K = 10, γ = 1 d , and β = 1 for PBRFF and GBRFF1; γ = 1 d and λ = 0 for GBRFF2. All the methods are run with 100 iterations (or landmarks) and are not run on datasets requiring more than 1000 seconds of execution time (because larger datasets requiring more than 1000 seconds by the fastest method do not fit in the RAM memory of the computer used for the experiments). We report the results in Figure 4.

We first recall that GBRFF2 learns at each iteration a random feature and a landmark while GBRFF1 only learns the landmark and PBRFF draws them randomly. Thus, GBRFF1 should present higher computation times compared to PBRFF. However, for datasets with a number of samples larger than 20, 000, GBRFF1 becomes cheaper than PBRFF. This is due to the fact that the SVM classifier learned by PBRFF does not scale as well as gradient boosting-based methods. The two-step method GFC is in addition also slower than GBRFF1. This shows the computational advantage of having a one-step procedure to learn both the representation and the final classifier. When looking at the time limit of 1000 seconds, both GBRFF1 and GBRFF2 are the fastest kernel-based methods compared to BMKR, GFC and PBRFF. This shows the efficiency of learning kernels in a greedy fashion. We also see that GBRFF2 performs faster than GBRFF1 for any number of samples. At the limit of 1000 seconds, it is able to deal with datasets that are 10 times larger than GBRFF1, due to the lower complexity of the learned weak learner used in GBRFF2. Finally, GBRFF2 is globally the second-fastest method behind the gradient boosting method LGBM that uses trees as base classifiers.

Performance Comparison Between All Methods

Table 2 presents for each dataset the mean results over the 20 splits using 100 iterations/landmarks for each method. Due to the size of the dataset "bankmarketing", we do not report the results of the algorithms that do not converge in time for this dataset, and we compute the average ranks and mean results over the other 15 datasets. In terms of accuracy, GBRFF2 shows very good results compared with the state-of-the-art as it obtains the best average rank among the six methods and on average the best mean accuracy leaving apart "bankmarketing". Interestingly, our method is the only kernel-based one that scales well enough to be applied to this latter dataset.

Comparison of LGBM and GBRFF2 on Toy Datasets

In this last experiment, we focus on LGBM and GBRFF2 which have been shown to be the two best performing methods in terms of accuracy and execution time. Even if BMKR is among the three best methods in terms of accuracy, we do not consider it for this experiment due to its poor execution time. Learning a classifier based on non-linear kernels through GBRFF2 has the advantage of being able to capture non-linear decision surfaces, whereas LGBM is not well suited for this because it uses trees as base learner. To illustrate this advantage, we consider three synthetics 2D datasets with non-linearly separable classes. The first one, called "swiss", represents two spirals of two classes side by side. The second one, namely "circles", consists of four circles with the same center and an increasing radius by alternating the class of each circle. The third dataset, called "board", consists of a four by four checkerboard with alternating classes in each cell. Here, both LGBM and GBRFF2 are run for 1000 iterations to ensure their convergence and parameters are tuned by CV as previously. Figure 5 gives evidence that GBRFF2 is able to achieve better results than LGBM using only a small amount of training examples, i.e., 500 or less. The performances are asymptotically similar for both methods on the board and circle datasets with a faster rate of convergence for GBRFF2. Furthermore, if we look at the decision boundaries and their associated performances at train and test time, we can see that LGBM is prone to overfit the training data compared to our approach, showing a drastic drop in performance between learning and testing. The learned decision boundaries are also smoother with GBRRF2 than with LGBM. These experiments show the advantage of having a non-linear weak learner in a gradient boosting approach.

Conclusion and Perspectives

In this paper, we take advantages of two machine learning approaches, gradient boosting and random Fourier features, to derive a novel algorithm that jointly learns a compact representation and a model based on random features. Building on a recent work [START_REF] Letarte | Pseudo-Bayesian learning with kernel Fourier transform as prior[END_REF], we learn a kernel by approximating it as a weighted sum of RFF [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF]. The originality is that we learn such kernels so that the representation and the classifier are jointly optimized. We show that we can benefit from a performance boost in terms of accuracy and computation time by considering each weak learner as a single trigonometric feature and learning the random part of the RFF. The experimental study shows the competitiveness of our method with state-of-the-art boosting and kernel learning methods.

The optimization of the random feature and of the landmark at each iteration can be computationally expensive when the number of iterations is large. A promising future line of research to speed-up the learning is to derive other kernel approximations where these two parameters can be computed with a closed-form solution. Other perspectives regarding the scalability include the use of standard gradient boosting tricks [START_REF] Ke | Lightgbm: A highly efficient gradient boosting decision tree[END_REF] such as sampling or learning the kernels in parallel.

Algorithm 2 :

 2 GBRFF1 Inputs : Training set S = (xi, yi) n i=1 ; Number of iterations T ; K number of random features; Parameters γ and β Output : sign H 0

Algorithm 3 :j 2 :

 32 GBRFF2 Inputs : Training set S = (xi, yi) n i=1 ; Number of iterations T ; Parameters γ and λ Output : sign H 0 (x) + T t=1 α t cos ω t • xi -b t 1: H 0 ← H 0 (xi) for t = 1, . . . , T do 3: ∀i = 1, . . . , n, wi = exp(-yiH t-1 (xi)) 4: ∀i = 1, . . . , n, ỹi = yiwi 5:

Fig. 1 .

 1 Fig. 1. GBRFF1 with K=1 on the two-moons dataset at different iterations. Top row shows the periodicity of the loss (light yellow indicates the minimal loss). Bottom row shows the resulting decision boundaries between the classes (blue & red) by fixing arbitrarily one coordinate of the landmark and minimizing the loss along the other one.

Table 1 .

 1 Description of the datasets (n: number of examples, d: number of features, c: number of classes) and the classes chosen as negative (-1) and positive (+1).

[ 6 ]

 6 that (1) draws with replacement n L landmarks from the training set; (2) learns a representation of n L features where each feature is computed using Equation (2) based on K=10 vectors drawn like our methods from N (0, 2γ) d ; (3) learns a linear SVM on the new representation. We select by CV its parameters γ ∈ 2 {-2,...,2} d , β ∈ 10 {-2,...,2} and the SVM parameter C ∈ 10 {-2,...,2} .

Fig. 3 .

 3 Fig. 3. Mean results the 16 datasets w.r.t. the same total number of random features T ×K for K ∈ {1, 5, 10, 20}, with T the number of boosting iterations.

Fig. 5 .

 5 Fig. 5. Comparison of LGBM and GBRFF2 on three synthetic datasets in terms of classification accuracy and decision boundaries (upper part of the figure) and in terms of performance w.r.t. the number of examples (last row of plots).

  Computation time in seconds required to train and test the six methods with fixed parameters on an artificial dataset having an increasing number of samples. The whole dataset is used for training and testing, and a method requiring more than 1000 seconds at a given step is not trained on the larger datasets.

		BMKR	GFC	PBRFF	GBRFF1	GBRFF2	LGBM
		1000
	Time in seconds	200 400 600 800
		1, 70 2 3, 82 9 8, 61 4 19 ,3 81 43 ,6 06 98 ,1 13 22 0, 75 3 49 6, 69 3 0	1, 11 7, 55 8 2, 51 4, 50 5 5, 65 7, 63 5 12 ,7 29 ,6 78 28 ,6 41 ,7 75 64 ,4 43 ,9 93
				Number of samples
	Fig. 4.

Table 2 .

 2 Mean test accuracy ± standard deviation over 20 random train/test splits. A '-' in the last row indicates that the algorithm did not converge in time on this dataset. Average ranks and mean results are computed over the first 15 datasets.

	Dataset	BMKR	GFC	PBRFF	GBRFF1	LGBM	GBRFF2
	wine	99.5 ± 1.0	99.3 ± 1.1	98.1 ± 2.1	98.3 ± 1.5	96.6 ± 3.2	98.5 ± 1.6
	sonar	78.8 ± 7.2	76.6 ± 3.2	76.7 ± 5.2	81.8 ± 3.5	82.4 ± 4.3	83.0 ± 5.0
	newthyroid	96.5 ± 1.7	96.5 ± 2.1	96.5 ± 1.5	95.3 ± 2.2	94.8 ± 2.9	96.9 ± 2.1
	heart	85.6 ± 4.0	79.4 ± 4.5	85.4 ± 3.5	83.6 ± 4.0	83.0 ± 3.5	83.1 ± 4.0
	bupa	68.1 ± 4.9	64.7 ± 3.2	69.0 ± 4.2	70.3 ± 4.9	72.0 ± 3.3	71.2 ± 4.5
	iono	94.2 ± 1.4	91.5 ± 2.3	94.2 ± 1.8	88.2 ± 2.3	93.3 ± 2.5	89.2 ± 2.1
	wdbc	96.1 ± 1.2	95.8 ± 1.3	96.5 ± 1.1	96.8 ± 1.1	95.8 ± 1.5	97.3 ± 1.2
	balance	96.0 ± 1.2	95.1 ± 2.0	98.9 ± 1.1	97.7 ± 0.7	93.5 ± 2.6	97.7 ± 0.6
	australian	85.9 ± 2.0	80.9 ± 2.4	84.6 ± 2.3	86.7 ± 1.7	85.5 ± 1.9	86.9 ± 1.9
	pima	76.4 ± 2.0	68.7 ± 2.6	76.1 ± 2.5	76.5 ± 2.7	75.5 ± 2.7	77.1 ± 2.5
	vehicle	96.6 ± 1.3	95.9 ± 0.8	96.5 ± 1.4	96.3 ± 1.2	96.7 ± 1.0	97.1 ± 1.0
	german	72.3 ± 1.8	64.3 ± 2.8	72.4 ± 1.4	73.7 ± 1.6	73.5 ± 1.7	74.0 ± 1.3
	splice	87.5 ± 1.0	87.0 ± 1.0	83.5 ± 0.7	83.9 ± 1.1	97.0 ± 0.5	92.4 ± 0.8
	spambase	93.5 ± 0.4	91.3 ± 0.6	91.6 ± 0.7	90.7 ± 0.7	95.6 ± 0.4	92.8 ± 0.6
	occupancy	99.3 ± 0.1	98.9 ± 0.7	98.9 ± 0.1	98.8 ± 0.1	99.3 ± 0.1	98.9 ± 0.1
	Mean	88.4 ± 2.1	85.7 ± 2.0	87.9 ± 2.0	87.9 ± 2.0	89.0 ± 2.1	89.1 ± 2.0
	Average Rank	2.88	4.94	3.75	3.81	3.44	2.19
	bankmarketing	-	-	-	89.7 ± 0.2	90.8 ± 0.2	90.0 ± 0.2

The code is available here: https://leogautheron.github.io

Acknowledgements. Work supported in part by French projects APRIORI ANR-18-CE23-0015, LIVES ANR-15-CE23-0026 and IDEXLYON ACADEMICS ANR-16-IDEX-0005, and in part by the Canada CIFAR AI Chair Program.


A Rosuel

P Vallet

P Loubaton

X Mestre

ON THE FREQUENCY DOMAIN DETECTION OF HIGH DIMENSIONAL TIME SERIES

Keywords: Spectral analysis, Detection Tests, High Dimensional Statistics, Random Matrix Theory

In this paper, we address the problem of detection, in the frequency domain, of a M -dimensional time series modeled as the output of a M × K MIMO filter driven by a K-dimensional Gaussian white noise, and disturbed by an additive M -dimensional Gaussian colored noise. We consider the study of test statistics based of the Spectral Coherence Matrix (SCM) obtained as renormalization of the smoothed periodogram matrix of the observed time series over N samples, and with smoothing span B. To that purpose, we consider the asymptotic regime in which M, B, N all converge to infinity at certain specific rates, while K remains fixed. We prove that the SCM may be approximated in operator norm by a correlated Wishart matrix, for which Random Matrix Theory (RMT) provides a precise description of the asymptotic behaviour of the eigenvalues. These results are then exploited to study the consistency of a test based on the largest eigenvalue of the SCM, and provide some numerical illustrations to evaluate the statistical performance of such a test.

INTRODUCTION

The detection of a low rank multivariate signal corrupted by a spatially uncorrelated noise with unknown statistics is an important signal processing problem that is met in the context of array processing of uncalibrated sensor networks, see e.g. [START_REF] Boonstra | Gain calibration methods for radio telescope arrays[END_REF] for applications to radio-astronomy, [START_REF] Ramírez | Detection of rank-p signals in cognitive radio networks with uncalibrated multiple antennas[END_REF] and [START_REF] Sala-Alvarez | Multiantenna glr detection of rank-one signals with known power spectral shape under spatially uncorrelated noise[END_REF] motivated by detection of primary signals in the context of cognitive radio. We also notice that the underlying signal model, called the errors in variables model, plays an important role in other fields such as econometrics, and was intensively studied in the past (see e.g. [START_REF] Söderström | Errors-in-variables methods in system identification[END_REF] for a review).

References [START_REF] Ramírez | Detection of rank-p signals in cognitive radio networks with uncalibrated multiple antennas[END_REF] and [START_REF] Sala-Alvarez | Multiantenna glr detection of rank-one signals with known power spectral shape under spatially uncorrelated noise[END_REF] developed GLRT tests which require solving rather difficult optimization problems to be implemented, in particular if the underlying errors in variables model is dynamic in the sense that the useful signal coincides with the output of a MIMO filter driven by a low-dimensional white noise and that the noise on each sensor is correlated in time (see [START_REF] Ramírez | Detection of rank-p signals in cognitive radio networks with uncalibrated multiple antennas[END_REF]Sec. V]).

The goal of this paper is to address the detection of the useful signal in the context of a dynamic errors in variables model when the observations dimension M is large, the number of available samples N is limited, and the rank K of the spectral density of the useful signal is much smaller than M and N . This timely context is in particular motivated by the considerable development of large sensor networks which tend to produce high-dimensional multivariate signals. While the high dimensionality of the observations poses a number of new statistical problems, it sometimes allows to simplify the performance analysis of traditional statistical inference schemes as shown in the present paper.

An important class of high-dimensional models, called the generalized dynamic linear factor models, were introduced at the end of the nineties in econometrics, see e.g. [START_REF] Forni | The generalized dynamic-factor model: Identification and estimation[END_REF][START_REF] Forni | The generalized dynamic factor model consistency and rates[END_REF]. It is assumed that M and N both converge towards +∞ in such a way that M N remains bounded ( M N may also converges towards 0) and that K remains fixed. The fundamental assumption formulated in these works is that the K non zero eigenvalues of the spectral density matrix of the useful signal converge towards +∞ while the spectral density of the noise remains bounded. Under this regime, the K largest eigenvalues of an estimate of the spectral density of the observation converge towards +∞ almost surely, thus leading to a consistent detection scheme. Moreover, it is possible to retrieve the frequency domain principal components of the useful signal. While this regime is certainly justified in the econometrics field, it is not the most relevant in the context of high dimensional array processing where relevant algorithms have the potential to produce a Signal to Noise Ratio (SNR) ratio gain of the order of the number of sensors. Therefore, regimes in which the SNR ρ, before applying such algorithms, is O(1 M ) term, are of special interest. Unfortunately, the generalized dynamic linear factor models assumptions lead to a larger order of magnitude SNR.

In contrast, large random matrix methods allow to obtain interesting results in situations where ρ = O( 1 M ). Previous related works addressed static models (also known as narrowband models in array processing) where a number of papers used the so-called spiked large random matrices defined as the sum of a low rank matrix due to the signal with a full rank random matrix representing the additive noise (see e.g. [START_REF] Benaych-Georges | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF][START_REF] Bianchi | Performance of Statistical Tests for Single-Source Detection Using Random Matrix Theory[END_REF][START_REF] Nadakuditi | Sample eigenvalue based detection of high-dimensional signals in white noise using relatively few samples[END_REF][START_REF] Vallet | Performance Analysis of an Improved MUSIC DoA Estimator[END_REF][START_REF] Vinogradova | Statistical inference in large antenna arrays under unknown noise pattern[END_REF]). The most complete results were obtained when the additive noise is temporally and spatially white. It appears that the standard detection test comparing the largest eigenvalue of the sample covariance matrix to a threshold is consistent if and only if M ρ 1 is larger than the ratio M N . In other words, the high-dimensionality of the observations produces a threshold effect on standard detectors. While the existing large random matrix results allow to address detection in static models, a considerable work is still needed to consider dynamic models (also known as wideband models in array processing) in regimes where ρ = O( 1 M ). In this paper, we propose to develop frequency domain methods where the sample covariance matrix is replaced by a frequency smoothed estimate Ĉ(ν) of the spectral coherence matrix of the observation at a Fourier frequency ν. Under some reasonable assumptions, we establish that at each Fourier frequency ν, Ĉ(ν) behaves as the sample covariance matrix of a static spiked large random matrix model whose noise part is temporally and spatially white. Using this useful result, we characterize the conditions under which the test comparing maxν λmax( Ĉ(ν)) to a certain threshold is consistent.

General notations. Vectors and matrices are denoted respectively as bold lower case and bold upper case. If A is a matrix, A represents its spectral norm and A F its Frobenius norm. If a1, . . . , aM ∈ C, diag(a1, . . . , aM ) represents the diagonal matrix with diagonal elements a1, . . . , aM , whereas if A is a matrix, diag(A) is the diagonal matrix representing the diagonal part of A. For any M × M Hermitian matrix A, λ1(A) ≥ . . . ≥ λM (A) denote the eigenvalues of A sorted in decreasing order. The notation N C M (0, I) is used for the standard M -dimensional complex Gaussian distribution.

MODEL AND ASSUMPTIONS

In this section, we introduce the signal model and some definitions and assumptions that will be used throughout the next sections.

We consider a M -dimensional observed signal (yn)n∈Z modeled as yn = un + vn

where (un)n∈Z is a useful signal defined as the output of an unknown causal and stable M ×K MIMO filter driven by a N C (0, I) white noise ( n)n∈Z, that is almost surely (a.s.),

un = +∞ k=0 H k n-k
and where the additive noise (vn) n∈Z is modeled as a M -dimensional stationary complex Gaussian time series whose components time series (v1,n)n∈Z, . . . , (vM,n)n∈Z are mutually independent. For each ν ∈ [0, 1], we define

H(ν) = +∞ k=0 H k e -i2πνk
, and for each m ∈ {1, . . . , M }, we denote by (rm(k)) k∈Z the covariance sequence of (vm,n)n∈Z and by sm(ν) the corresponding spectral density.

Denote by

ξ y (ν) = 1 √ N N n=1 yne -i2πν(n-1)
the finite Fourier transform of (yn)n∈Z over the sample window n = 1, . . . , N , and for an even integer B < N referred to as smoothing span,

Ŝy(ν) = 1 B + 1 B/2 b=-B/2 ξ y (ν + b/N )ξ y (ν + b/N ) *
the classical frequency smoothed estimate of the spectral density Sy(ν) of (yn)n∈Z, where

Sy(ν) = H(ν)H(ν) * + Sv(ν)
with Sv(ν) = diag(s1(ν), . . . , sM (ν)) the spectral density of (vn)n∈Z. We recall (see e.g. the classical reference [START_REF] Brillinger | Time series: data analysis and theory[END_REF]) that in the classical large sample size regime where B, N → ∞ while

M, K are fixed, E[ Ŝy(ν)] = Sy(ν) + O( B 2 N 2 ) and E[ Ŝy(ν) - E[ Ŝy(ν)] 2 ] = O( 1 B
) and the bias-variance compromise is achieved by choosing B → ∞ such that B N → 0. In the following, we study the statistical behaviour of the Spectral Coherence Matrix (SCM)

Ĉy(ν) = diag Ŝy(ν) -1 2 Ŝy(ν) diag Ŝy(ν) -1 2 (2)
under the following high-dimensional asymptotic regime. We assume that M = M (N ) and B = B(N ) are both functions of N such that

M ∼ N α and M B ----→ N →∞ c ∈ (0, 1) (3)
for α ∈ (0, 1), while K remains fixed with respect to N . Choosing B in such a way that M/B → 0 (B = O(N β ), β > α) would make Ŝy(ν) a consistent estimator of Sy(ν). However, in practice, for finite values of M and N , it is not necessarily possible to choose B in such a way that M B and B N , thus making the regime M/B → 0 inaccurate. Assumption 3 appears thus relevant.

Under this asymptotic regime, it is necessary to precise how certain quantities previously defined evolve with respect to M (and thus N ). We first consider the following assumption regarding the noise part (vn)n∈Z.

Assumption 1

The following hold lim sup

M →∞ max m=1,...,M k∈Z 1 + |k| 2 |rm(k)| < ∞ (4)
and

lim inf M →∞ min m=1,...,M inf ν∈[0,1] sm(ν) > 0 (5)
Note that condition (4) trivially implies that the noise spectral densities s1, . . . , sM are twice continuously differentiable and that lim sup

M →∞ max m=1,...,M sup ν∈[0,1] |s (i) m (ν)| < +∞ for i ∈ {0, 1, 2}, with s (i)
m the i-th order derivative of sm. The condition ( 5) is intuitively expected as our objective is to study the behaviour of the spectral coherence matrix [START_REF] Ramírez | Detection of rank-p signals in cognitive radio networks with uncalibrated multiple antennas[END_REF], which involves a renormalization by the inverse of the spectral density estimates of each time series (noise whitening). Note also that condition (5) together with (4) implies that the total power of the noise satisfies

0 < lim inf M →∞ 1 M E vn 2 ≤ lim sup M →∞ 1 M E vn 2 < ∞
The next assumption is related to the signal part (un)n∈Z.

Assumption 2

The following hold lim sup

M →∞ +∞ k=0 (1 + k) H k < ∞ (6)
and if h1(ν), . . . , hM (ν) denote the rows of H(ν), then lim

M →∞ max m=1,...,M sup ν∈[0,1] hm(ν) = 0 (7)
Since K is assumed fixed with respect to N , condition [START_REF] Forni | The generalized dynamic factor model consistency and rates[END_REF] implies that E un 2 = k≥0 H k 2 F = O(1) so that the "observed SNR" vanishes as M → ∞ at rate 1/M , that is ρ = E un 2 E vn 2 = O 1 M . We notice that this regime is in accordance with a number of works studying e.g. the behaviour of large sensor array processing techniques [START_REF] Vallet | Performance Analysis of an Improved MUSIC DoA Estimator[END_REF]: taking advantage of large number M of sensors to increase the observed SNR (or SNR before matched filtering in this context) ρ by a factor O(M ), it is reasonable to expect that despite a low SNR, reliable information can be extracted on the useful signal. In order to explain the significance of condition [START_REF] Benaych-Georges | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF], we notice that the power E|um,n| 2 of the contribution of un on sensor m can be written as E|um,n| 2 = hm(ν) 2 dν. Condition (7) thus implies that the SNR

E|um,n| 2
E|vm,n| 2 on each sensor m converges towards 0 when M → +∞. This condition in particular holds if the useful signal powers received on the various sensors are of the same order of magnitude, in which case

E|um,n| 2 E|vm,n| 2 = O( 1 M ).

ASYMPTOTIC BEHAVIOUR OF THE SCM Ĉy(ν)

In this section, we study the SCM Ĉy(ν) under the asymptotic regime (3), and show that this matrix can be approximated by a standard model from RMT. To that purpose, we study separately the pure noise case (yn = vn in (1)) and the noise free case (yn = un in ( 1)), and in the following, we denote by

ξ v (ν) = 1 √ N N n=1 vne -i2πν(n-1) ξ u (ν) = 1 √ N N n=1 une -i2πν(n-1)
the finite Fourier transforms of (vn)n∈Z and (un)n∈Z respectively, and also use the notation VN = {0, 1 N , . . . , N -1 N } for the set of Fourier frequencies. Due to space constraints, the proofs of the results are omitted.

We have the following approximation result regarding the finite Fourier transforms of (vn)n∈Z.

Theorem 1 Let Σv(ν) = 1 √ B+1 [ξ v ν -B 2N , . . . , ξ v ν + B 2N
]. Then under Assumption 1, for all ν ∈ VN , there exists a M ×(B+1) random matrix Z(ν) with i.i.d. NC(0, 1) entries such that

max ν∈V N Σv(ν) - 1 √ B + 1 Sv(ν) 1 2 Z(ν) a.s. ----→ N →∞ 0
Using standard concentration bound on the norm of Wishart matrices (see e.g. [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF]), we deduce that maxν∈V N Z(ν) = O( √ B) a.s. and

max ν∈V N Ŝv(ν) - 1 B + 1 Sv(ν) 1 2 Z(ν)Z(ν) * Sv(ν) 1 2
a.s.

----→

N →∞ 0 (8)
where

Ŝv(ν) = 1 B + 1 B/2 b=-B/2 ξ v (ν + b/N )ξ v (ν + b/N ) *
Thus, Theorem 1 shows that the smoothed periodogram of (vn)n∈Z behaves (asymptotically in operator norm) as a Wishart matrix with scale matrix coinciding with the spectral density Sv(ν).

We now turn to the study of the finite Fourier transforms of the signal part (un)n∈Z.

Theorem 2 Let Σu(ν) = 1 √ B+1 [ξ u ν -B 2N , . . . , ξ u ν + B 2N ] and Σ (ν) = 1 √ B+1 [ξ ν -B 2N , . . . , ξ ν + B 2N ]. Then under Assumption 2, it holds that max ν∈V N Σu(ν) -H(ν)Σ (ν) a.s.
----→ N →∞ 0 Note that the type of approximation stated in Theorem 2 is wellknown in the classical large sample size regime in which M, K, B are fixed while N → ∞ since in that case [START_REF] Brillinger | Time series: data analysis and theory[END_REF]Th. 4.5.2] sup

ν∈[0,1] Σu(ν) -H(ν)Σ (ν) = O log(N ) N a.s.
Of course, in the high dimensional regime in which M and B also converge to infinity, the result of Theorem 2 cannot be deduced from [START_REF] Brillinger | Time series: data analysis and theory[END_REF] and requires a specific study.

Using Theorems 1 and 2, we directly obtain as for [START_REF] Bianchi | Performance of Statistical Tests for Single-Source Detection Using Random Matrix Theory[END_REF], that for all ν ∈ VN there exists a M × (B + 1) random matrix X(ν) with i.i.d. NC(0, 1) entries such that

max ν∈V N Ŝy(ν) - 1 B + 1 Sy(ν) 1 2 X(ν)X(ν) * Sy(ν) 1 2
a.s.

----→

N →∞ 0 (9)
Moreover, using Assumption 2 (condition ( 7)), we can show that

max ν∈V N diag Ŝy(ν) -Sv(ν) a.s. ----→ N →∞ 0 (10)
Equipped with ( 9) and ( 10), we are now in position to study the behaviour of the SCM Ĉy(ν).

Corollary 1 Under Assumptions 1 and 2, there exists a M ×(B+1) random matrix X(ν) with i.i.d. NC(0, 1) entries such that

max ν∈V N Ĉy(ν) -Ξ(ν) 1 2 X(ν)X(ν) * B + 1 Ξ(ν) 1 2
a.s.

----→

N →∞ 0 (11)
where Ξ(ν) = Sv(ν) -1 2 H(ν)H(ν) * Sv(ν) -1 2 + I. Let us make two observations regarding Corollary 1.

First, notice that the operator norm approximation [START_REF] Vinogradova | Statistical inference in large antenna arrays under unknown noise pattern[END_REF] implies, thanks to Weyl's inequality, that

max m=1,...,M max ν∈V N λm Ĉy(ν) - λm Ξ(ν) 1 2 X(ν)X(ν) * B + 1 Ξ(ν) 1 2 a.s.
----→

N →∞ 0 (12)
which shows that any linear spectral statistic [START_REF] Bai | CLT for linear spectral statistics of large-dimensional sample covariance matrices[END_REF] based on the spectrum of Ĉy(ν) behaves (at first order) as the same linear spectral statistic applied to the spectrum of the corresponding Wishart matrix Ξ(ν)

1 2 X(ν)X(ν) * B+1 Ξ(ν) 1 2
. Second, an interpretation of the results of Theorems 1 and 2, and Corollary 1 can be given in the light of array processing application. Indeed, the time domain model ( 1) is usually referred to as wideband model, and in particular the useful signal contribution (un)n∈Z is not necessarily confined to a low-dimensional subspace of C M due to the filtering induced by (H k ) k≥0 . In the frequency domain, Theorems 1 and 2 show that we retrieve in some sense a narrowband model, and that the useful signal contribution is now confined to K-dimensional subspace of C M , which opens the possibility to use standard narrowband techniques for e.g. detecting the presence of the useful signal (un)n∈Z, using test statistics based on the eigenvalues of Ŝy(ν) or Ĉy(ν) at multiple Fourier frequencies ν ∈ VN .

We take profit of these two observations to study the consistency of a certain spectral detection test in the next section.

APPLICATION TO SPECTRAL DETECTION

As we have seen in the previous section, the SCM Ĉy(ν) behaves asymptotically in operator norm as a Wishart matrix, with associated scale matrix Ξ(ν) = Sv(ν) -1 2 H(ν)H(ν) * Sv(ν) -1 2 + I being a fixed rank K perturbation of the identity matrix. We can therefore exploit the well-known results on spiked models in the RMT literature [START_REF] Baik | Eigenvalues of large sample covariance matrices of spiked population models[END_REF]. We define for the remainder the function

φ(γ) = (γ+1)(γ+c) γ if γ > √ c (1 + √ c) 2 if γ ≤ √ c
and fix a frequency ν * N ∈ VN such that

ν * N ∈ argmax ν∈V N λ1 Sv(ν) -1 2 H(ν)H(ν) * Sv(ν) -1 2
We consider the following purely technical additional assumption.

Assumption 3 For all k ∈ {1, . . . , K}, there exists γ k > 0 such that

λ k Sv(ν * N ) -1 2 H(ν * N )H(ν * N ) * Sv(ν * N ) -1 2 ----→ N →∞ γ k
From ( 12) and the results of [START_REF] Baik | Eigenvalues of large sample covariance matrices of spiked population models[END_REF], we deduce the following indivual behaviour of the eigenvalues of the SCM.

(ν * N ) -1 2 H(ν * N )H(ν * N ) * Sv(ν * N ) -1 2
are sufficiently large:

γ k = lim N →∞ λ k Sv(ν * N ) -1 2 H(ν * N )H(ν * N ) * Sv(ν * N ) -1 2 > √ c
The result of Corollary 2 can be exploited to obtain consistent test statistics based on the eigenvalues of Ĉy, for detecting the presence of the useful signal in model (1), i.e. considering the hypothesis test H0 : yn = vn vs H1 : yn = un + vn Since the intrinsic dimensionality K of the useful signal is not necessarily known in practice, we consider the detection test solely based on the largest eigenvalue of the SCM, i.e.

T = 1 ((1+ √ c) 2 + ,+∞) max ν∈V N Ĉy(ν)
where is some threshold. We thus have the following consistency result. where Pi is the underlying probability measure under hypothesis Hi.

In practical situations, the condition γ1 > √ c means that the quantity Sv(ν * N ) -1 H(ν * N )H(ν * N ) * should be larger than M B . We now illustrate numerically the above asymptotic approximations when K = 1. The noise is chosen as a MA(1) process with standard Gaussian innovation (zn)n∈Z, i.e. vn = zn + θ1zn-1, whereas for the useful signal, we choose H k = CSNR 1 √ M β k (1, . . . , 1) T with β = 10 11 , and CSNR a factor controlling the SNR defined below.

Figures 1 and2 represent the probability of detection as a function of the probability of false alarm (ROC curve) for the test statistic T , for different values of the threshold and a total of 10 5 draws.

In Figure 1, the ROC curve is plotted for different values of SNR defined in the frequency domain as follows

SNR freq = sup ν∈[0,1] M m=1 hm(ν) 2 sm(ν)
As the SNR increases, the performance of the test improves as well. In Figure 2, we compute the ROC curve for increasing values of N

Corollary 2 2

 22 Under Assumptions 1, 2 and 3, for all k = 1, . . . , K, λ k Ĉy(ν * Consider the set of indexes K = {k ∈ {1, . . . , K} : γ k > √ c}. Corollary 2 implies that each "signal" eigenvalue λ k ( Ĉy(ν * N )) of the SCM at frequency ν * N for which k ∈ K asymptotically splits from the "noise" eigenvalues λK+1( Ĉy(ν * N )), . . . λM ( Ĉy(ν * N )), which concentrate in a neighborhood of the interval [(1-√ c) 2 , (1+ √ c) 2 ]. Likewise, the signal eigenvalues λ k ( Ĉy(ν * N )) for which k ∈ K are asymptotically absorbed in a neighborhood of [(1 -√ c) 2 , (1 + √ c) 2 ]. This phase transition phenomenon for the kth eigenvalue of Ĉy(ν * N ) thus occurs when the eigenvalues of Sv

Theorem 3

 3 Under Assumptions 1, 2 and 3, and if γ1 > √ c then for all ∈ 0, φ(γ1) -(1 + √ c) 2 , and i ∈ {0, 1} Pi lim N →∞ T = i = 1

Fig. 1 .

 1 Fig. 1. ROC curve of T with varying SNR with B = N 0.7 and B = 2M . As N grows, the performance of the test improves as expected since we get closer to the high dimensional asymptotic regime (3).

Fig. 2 .

 2 Fig. 2. ROC curve of T with varying N

If K = 1, M ρ represents the signal to noise after spatial matched filtering.

The authors were supported by the grant ANR-17-CE40-0003 of the French National Research Agency ANR (project HIDITSA).


M W Argyle

Sala

M Scardia

J.-L Prieur

L Pansecchi

R W Argyle

M Sala

Speckle observations with PISCO in Merate: VIII. Astrometric measurements of visual binaries in 2007 and new orbits of the multiple system Zeta Aqr

Keywords: close -binaries, visual -astrometry -techniques, interferometric -individual (Zeta Aqr, ADS 15971

Speckle observations with PISCO in Merate. VIII. Astrometric measurements of visual binaries in 2007 and

new orbits of the multiple system Zeta Aqr

Introduction

This paper deals with the results of speckle observations of visual binary stars made in Merate (Italy) during the second semester of 2007 with the Pupil Interferometry Speckle camera and COronagraph (PISCO) on the 102 cm Zeiss telescope of INAF -Osservatorio Astronomico di Brera (OAB, Brera Astronomical Observatory). It is the eighth of a series (Scardia et al. 2005, 2006, 2007, 2008b, Prieur et al. 2008, Scardia et al. 2009, Prieur et al., 2009, herein: Papers I to VII), whose purpose is to contribute to the determination of binary orbits. PISCO was developed at Observatoire Midi-Pyrénées (France) and first used at Pic du Midi from 1993 to 1998. It was moved to Merate in 2003 and used there since.

We briefly describe our observations in Sect. 2. Then we present and discuss the astrometric measurements in Sect. 3. Finally in Sect. 4 we propose new revised orbits for the multiple system Zeta Aqr and derive estimate values for the component masses.

Observations and description of the sample

The observations were carried out with the PISCO speckle camera with the ICCD (Intensified Charge Coupled Device) detector belonging to Nice University (France). This instrumentation is presented in Prieur et al. (1998) and our observing procedure is described in detail in Paper VI. For the present observations, thanks to an improvement of our software and a faster computer, we were also able to compute in real time the restricted triple correlation (Aristidi et al., 1997) that we use for resolving the 180 • ambiguity (see Sect. 3.2).

The description of our sample can be found in our previous papers (e.g., Paper VI). It basically includes all the visual binaries for which new measurements are needed to improve their orbits, that are accessible with our instrumentation.

The distribution of the angular separations measured in this paper is displayed in Fig. 1a and shows a maximum for ρ ≈ 0 ′′ .7. The largest separation of ρ = 4 ′′ .4 was obtained for ADS 11178. The smallest separation was measured for KUI 93, with ρ = 0 ′′ .17. This is close to the diffraction limit ρ d = λ/D ≈ 0 ′′ .13 with the R filter (i.e. λ = 650 nm) and the Zeiss telescope whose aperture is D = 1.02 m.

The distribution of the apparent magnitudes m V of the binaries measured in this paper is presented in Fig. 1b and the difference of magnitudes ∆m V between the two components in Fig. 1c. The telescope aperture and detector sensitivity lead to a limiting magnitude of about m V = 9. A careful handling of the image processing enabled us to measure a few couples with ∆m V > 3 (Fig. 1c).

Using the Hipparcos parallaxes, we were able to construct the HR diagram of those binaries that is displayed in Fig. 2. We only plotted the objects for which the relative uncertainty on the parallax was smaller than 50%. It can be seen that a large part of the HR diagram is covered by our sample. In the future we would like to acquire a more sensitive detector in order to observe fainter (and cooler) main sequence stars.

Astrometric measurements

The astrometric measurements of the observations made during the second semester of 2007 are displayed in Table 1. The format of the first ten columns of this table is the same as for the previous papers of this series. It is described for instance in Paper VI. In Col. 11, we report some information about the secondary peaks of the auto-correlation files: diffuse (Diff.) or elongated (Elon.). The last three columns are dealing with the residuals and will be presented in Sect. 3.3.

The characteristics of the PISCO R and V filters can be found in Table 1 of Paper III. Some objects were observed with no filter because they were too faint. This is indicated with W (for "white" light) in the filter column (Col. 5). The corresponding bandpass is 420 nm corresponding to that of the ICCD detector, with a central wavelength of about 650 nm, close to that of the R filter.

The average values of the errors of the 283 measurements reported in this table are 0 ′′ .014 ± 0 ′′ .009 and 0 • .56 ± 0 • .35 for ρ and θ, respectively. Our error determination procedure is described in detail in Paper III.

We have noticed some peculiarities in the autocorrelations of ADS 504, 11579, and 16011, with the presence of small symmetric patches (that may be due to detector noise ?). We will re-observe those objects in order to verify that those peculiarities are artefacts.

The triple system ADS 11344

ADS 11344 is a triple star whose components A, B and C give rise, in the literature, to three different systems: HU 66AB, HU 66BC and STT 351AC (also reported as STT 351AB-C). In fact the C component is the brightest star and B the faintest. As the separation of the AB system is very small, many observers are unable to resolve it and confine themselves to measure AB, as it was a single star (A), with respect to C. This is a source of confusion in the literature1 .In many cases the reported measurements of STT 351AC are in fact those of STT 351AB-C, when the AB system was not resolved. This also happened to us for our first observation of ADS 11344 made in 2004.690 with (ρ = 0 ′′ .802, θ = 25. • 7), reported in Paper II, and wrongly attributed to the STT 351AC system.

The duplicity of the auto-correlation secondary peaks is clear in the observation of ADS 11344 made in 2007.671 (see Fig. 3). We measured them separately and found that the innermost peaks correspond to STT 351AC and the outermost ones to HU 66BC. Finally, each couple composed by the big and small secondary peaks correspond to HU 66AB. This identification is validated with the small residuals that we obtain with Novakovic (2008b)'s orbit of HU 66BC, and with Seymour (2002)'s orbit of HU 66AB (see Table 1).

Quadrant determination

As our measurements were obtained from the symmetric auto-correlation files, the θ values first presented a 180 • ambiguity. To try resolving this ambiguity and determine the quadrant containing the companion, we have used Aristidi et al. (1997)'s method by computing and analysing restricted triple correlation files. For the couples with the largest separations, a straightforward determination could be done when the companions were visible in the long integration files.

As a result, in Table 1, we are able to give the non-ambiguous (or "absolute") position angles of 213 out of 283 measurements, i.e. 76% of the total. They are marked with an asterisk in Col 9. Otherwise, our angular measurements were reduced to the quadrant reported in the "Washington Double Star Catalogue" (Mason et al, 2009, hereafter WDS Catalogue).

Our "absolute" θ values are consistent with the values tabulated in the WDS Catalogue for all objects except for ADS 10, 10905, 11558, 15954, 16204, 16307, 16958 and 17122. In all those cases, however, a clear and high contrast between the two secondary peaks of the triple correlation is observed, which is a good indication of the validity of our quadrant determination. In the following, we discuss each case using the usual convention of numbering the quadrants (Q) from 1 to 4 to indicate the North-East, South-East, South-West and North-West quadrants, respectively.

ADS 10: Our value of Q=2 is consistent both with the space-based Tycho measurement of 1991.75 [START_REF] Høg | Tycho-2 Cat[END_REF] and with Douglass et al. (1997)'s measurement of 1992.878. It is likely that the quadrant Q=4 reported in the WDS Catalogue comes from the other 12 observations listed in the "Fourth Catalogue of Interferometric Measurements of Binary Stars" (Hartkopf et al. 2009, hereafter IC4). The two components of ADS 10 have a similar brightness, with ∆m V = 0.16, which well accounts for the difficulty of determining the quadrant. Furthermore, most of those observations of IC4 1, for which Hipparcos parallaxes were obtained with a relative error smaller than 50% (184 objects).

(O-C) ∆θ(O-C) (mm) (") (") ( • ) ( • ) (") ( • ) 16153+0416 STF2027 9980
(mm) (") (") ( • ) ( • ) (") ( • )

Fig. 3

Auto-correlation of the triple star ADS 11344.

were made in V , whereas we observed with no filter, which corresponds to a bandpass close to that of the R filter (see Sect. 2). There can exist a quadrant inversion between V and R, since the spectral type of the primary component is A0 and is probably associated to a cooler companion.

ADS 10905: Our quadrant determination Q=2 in R is consistent with the values we found in 2004.690 in Merate (Paper II) and in 1998.654 at the Pic du Midi with a 2m-telescope (Scardia et al., 2000a). It is however in disagreement with the quadrant Q=4 reported for the other 58 observations listed in IC4, made in V . Like ADS 10, this discrepancy can be due to a quadrant inversion that is likely since ∆m V = 0.12, with a spectral class of A0 III for the primary.

ADS 11558: Our value of Q=2 in R is inconsistent with the fourth quadrant given for the 27 measures listed in IC4, all made in V . The Wilson spectral classes are A1n, for the primary and G0 for the secondary (WDS Catalogue notes). The small difference in magnitude, ∆m V = 0.14, might account for the quadrant inversion since the G0 component could appear as the brightest one in R. We already observed this object in 2004.657 (Paper II) but the seeing conditions were not good enough for enabling the quadrant determination.

ADS 15954: Our quadrant Q=4 obtained with no filter is in agreement with Tycho measurement of 1991.69 [START_REF] Høg | Tycho-2 Cat[END_REF], and in disagreement with the other 9 observations reported in IC4, made in V . For this object, the components have a larger brightness difference, with ∆m V = 0.4, and the spectral type of the primary is O9V.

ADS 16204: Here also, our value of Q=3 obtained with no filter is in agreement with Tycho measurement of 1991.76 [START_REF] Høg | Tycho-2 Cat[END_REF] and in disagreement with the other 7 observations reported in IC4, all made in V . The very small value of ∆m V = 0.06 well explains the difficulty of determining the good quadrant. Furthermore, we noticed an inversion in the magnitudes of the components tabulated in the WDS Catalogue (the faintest magnitude being placed in the column of the primary component and vice-versa), which might be another source of error. We observed this binary also on 2006.951 (Paper VI) but the low signal-to-noise ratio did not allow us to determine the quadrant. ADS 16307: Our quadrant Q=1 determined in R is in agreement with Tycho determination of 1991.72 [START_REF] Høg | Tycho-2 Cat[END_REF], but in disagreement with the other 5 observations reported in IC4, and made in V . This discrepancy is difficult to explain since ∆m V = 0.81. The primary spectral class is A0.

ADS 16958: Here conversely, our quadrant Q=1 determined without filter is consistent with all the observations reported in IC4 (all made in V ) but not with the Tycho measurement (Q=3) of 1992.07 [START_REF] Høg | Tycho-2 Cat[END_REF] whose quadrant was adopted by the WDS Catalogue. The spectral class is F2 and the difference in magnitude of the components is significant: ∆m V = 0.33.

ADS 17122: Our determination of Q=3 made with no filter is in agreement with Tycho measurements of 1991.67 [START_REF] Høg | Tycho-2 Cat[END_REF], but not with the first quadrant (adopted by the WDS Catalogue) of the other 7 observations reported in IC4, all made in V. The difference in magnitude of the components is ∆m V = 0.17 and the spectral class A5.

Comparison with published ephemerides

The (O -C) (Observed minus Computed) residuals of the measurements for the 65 systems with a known orbit in Table 1 are displayed in Cols. 13 and 14 for the separation ρ and position angle θ, respectively. The orbital elements used for computing the ephemerides were retrieved from the "Sixth Catalogue of Orbits of Visual Binary Stars" (Hartkopf & Mason, 2009, hereafter OC6). The corresponding authors are given in Col. 12, using the style of the OC6 references.

The residuals reported in this table were computed with the most recent orbits found in OC6, but for some objects (i.e., for ADS 293 and 11635), we also give the O -C values relative to old orbits found in the previous issues of OC6, when they are still valid. For STF 2909 AB, we also give the residuals obtained with our new orbits presented in Sect. 4, for comparison.

For most of the orbits used for computing the residuals reported in Table 1, the equinox of the node Ω was 2000.0, and the correction to be applied to the computed position angle was negligible. Such a correction ADS 14573 ADS 15971 (Olevic & Cvetkovic, 2004) ADS 10429

Fig. 4

Residuals of the measurements of Table 1 from the published orbits.

was only necessary for ADS 9378 (equinox: 1950.0) and for ADS 9969, 16314 and 17149 (equinox: 1900.0).

Fig. 4 shows that the residuals are well centered around the origin, with a rather large scatter that can be explained by the (old) age of many orbits. The mean values computed with the residuals of Table 1

are < ∆ρ O-C >= 0 ′′ .00±0 ′′ .05 and < ∆θ O-C >= 0 • .2±1 • .7.
The small values obtained for those offsets provide an additional validation of our calibration made with a grating mask (see Paper III), which is in good agreement with the measurements made by the other observers. In the following, we examine the cases of ADS 10429, 14573 and 15971 that appear with the largest residuals in Fig. 4.

ADS 10429: this couple is difficult to measure because of a large brightness difference between the two components: ∆m V = 2.5. This may explain why there are only 3 speckle measurements in IC4: two made in 1995 by the US Naval Obs. in Washington, and one obtained by our group in 2004.659 with ρ = 0 ′′ .867 ± 0 ′′ .006 and θ = 10 • .3 ± 0 • .9 (Paper II). All the other measurements were obtained visually. The revision of Olevic et al (1993)'s orbit would be too premature, because it is still satisfactory in terms of position angle, and also because the orbital arc covered by the companion is very small.

ADS 14573: here also the orbital arc covered by the companion is too small to justify a new computation: the resulting orbit would be completely undetermined (grade 5).

ADS 15971: the recent orbit of STF 2909 AB computed by Olevic & Cvetkovic (2004) leads to rather large residuals with both our measurement and the last ones reported in IC4. The companion is moving closer to the primary than expected by Olevic's orbit. We propose a new orbit in the next section.

4 Revised orbits of the multiple system Zeta Aqr -ADS 15971 4.1 Visual orbit AB of Zeta Aqr ADS 15971 -STF 2909 AB -HIP 110960 -Zeta Aqr, is a famous binary star, with more than 1200 measurements made since its discovery in 1777, by Christian Mayer in Mannheim with a Bird mural quadrant. It was first reported in Mayer's catalogue (actually the first catalogue of binary stars) that contained 89 couples and was published in the "Astronomische Jahrbuch für 1784". The first real measures of ADS 15971 were made by W. Herschel in 1779 with a 6-inch telescope.

This object was observed with PISCO at the Pic du Midi in 1998and in Merate in 2004, 2005and 2007, and two of us have also contributed to this effort with micrometer observations. MS used a 23-cm refractor in Merate in 1982and 1984, and the 38-cm GPO in La Silla in 1992. RWA obtained an annual series of measures between 1990 and 2008 using the 20-cm refractor at Cambridge University Observatories. The whole list of measurements used for computing the new orbits was kindly provided by B. Mason from the WDS data base maintained by the United States Naval Observatory. We also added a recent measurements made with a micrometer by RWA with the 20-cm Cambridge Observatory refractor (Argyle, 2009).

We first computed the preliminary orbital elements with the analytical method of [START_REF] Kowalsky | Procès-verbaux de l'Université Imperiale de Kasan Ling[END_REF]. We then used them as initial values for the least-squares method of Hellerich (1925), which leads to an improvement of the orbital elements with an estimation of their errors.

The total list of 1211 measurements was obtained with a wide range of telescopes (with apertures from 7-cm to 4-m). To obtain a more robust solution, we rejected the aberrant measurements for which the residuals were larger than 1.7-σ, i.e., with ∆ρ O-C > 0 ′′ .31 and ∆θ O-C > 4 • 6. This resulted in a selection of 1189 observations. The final orbital elements are presented in Table 2 with the errors obtained by Hellerich's method. For all the tables of this section, we use the same format as for the previous papers. A full description can be found, for instance in Paper VI. The mean standard deviations of the residuals are 0 ′′ .22 and 2 • .9 for ρ and θ, respectively. This scatter is explained by the large range of instruments and epochs, and also by the perturbation of a third invisible body (see below).

The apparent orbit is shown in Fig. 5 as a solid line and the observational data plotted as small crosses or, in the case of PISCO observations, as small circles.

Perturbation orbit Bb-P of Zeta Aqr

The zoomed part of the orbit in Fig. 5b on the most recent observations exhibits an oscillatory trend which betrays the presence of a third star perturbing the orbital motion of the AB couple. This property was first noticed by Strand (1942), who published the first orbits for the close and wide pairs simultaneously. A series of subsequent studies (Strand 1942, Franz 1958, Harrington 1968, Heintz 1984c) precised the parameters of the perturbation orbit with a period of 25.5 years and a semi-major axis of 0 ′′ .06. The mass of the third star was estimated at 0.22 M ⊙ by Harrington (1968) or at 0.4 M ⊙ by Heintz (1984c).

It is difficult to locate this companion, since the masses of the A and B components are similar and the presence of a companion "a" close to A or another companion "b" close to B would produce similar effects on the AB orbit. On the basis of parallax observations, Franz (1958) linked it to a possible Bb pair, whereas Heintz (1984) favoured a hypothetical Aa pair.

At optical wavelengths, [START_REF] Weigelt | [END_REF] found a close companion to Zeta Aqr A with a separation of 0 ′′ .064±0.005 in 1978.964, which is much too small to correspond to the perturbation orbit. In near infrared, a companion was detected close to Zeta Aqr B at 0 ′′ .17 by [START_REF] Mccarthy | [END_REF] with the Kitt Peak 3.8-m telescope. The angular separation and the estimated mass were both in good agreement with Harrington's orbit. We will assume in this paper that this observation is valid and that the perturbation is due to the Bb couple detected by McCarthy et al. According to the convention used by the OC6, the perturbation orbit will be called Bb-P, since its corresponds to the orbit of B relative to the center of mass of the outer Bb system. Using all the photographic and speckle measurements published since 1923, we computed the residuals of our orbit of ADS 15971 AB in Cartesian coordinates, from which we derived the elements of the perturbation orbit Bb-P and its uncertainties with Hellerich (1925)'s least-squares method. To our knowledge, the uncertainties on those elements had never been computed before. The results are presented in Tables 2, 3 and in Fig. 6.

The (O -C) residuals obtained with the combination of orbits AB and Bb-P, restricted for reasons of space to the observations made after 2002, are given in Table 4. The corresponding mean standard deviations are 0 ′′ .11 and 1 • .25 for ρ and θ, respectively, which is about half of the values obtained with the orbit AB alone (see Sect. 4.1). The combined ephemerides for 2009-2020 are presented in Table 5.

Our two orbits are still temporary since only half of the AB orbit has been monitored yet. The orbital elements of the main orbit are still poorly known (for instance the period has an uncertainty of 10%).

Discussion: nature of the perturbing component

In the SIMBAD data base, the two components A and B of Zeta Aqr are identified separately as HD 213052 (F3V, m V = 4.42) and HD 213051 (F6IV, m V = 4.5). Some physical parameters derived for this system from our orbital elements are presented in Table 6. Both the linear size of a (Col. 5) and M total (Col. 6) were computed from our orbital elements and using the Hipparcos parallax revised by van Leeuwen (2007) and reported in Col. 3. Using the relation we established for IV-type stars (Scardia et al, 2008c, p 380) we derived a total mass of 3.59 ± 0.96 M ⊙ . The theoretical value for two stars F3V+F6IV is M total = 2.64 M ⊙ (Straizys & Kuriliene, 1981). The mass of the invisible companion could thus be as large as 0.95±0.96 M ⊙ . Unfortunately the corresponding uncertainty prevents us from drawing any conclusion about the nature of the invisible companion.

The mass of the b component can also be estimated with the measurement of angular separation of Bb made by [START_REF] Mccarthy | [END_REF], which was ρ = 0. ′′ 17 in 1981.94. The ephemerides for this epoch with our orbit Bb-P lead to ρ ′ = 0 ′′ .056. Since ρ/ρ ′ = (M B +M b )/M b , it comes: M B /M b = 2.0 and M b = 0.65 M ⊙ . Hence the b component could be a late-type main sequence star of type K3-K5. [START_REF] Mccarthy | [END_REF] measured its flux in near infrared and found K = 6.0 mag with a strong infrared excess K -L = 1.3 ± 0.2 mag, redder than that expected of even the coolest main-sequence stars and corresponding to a colour temperature of 1440 K. It is thus expected to be much fainter than the A, B components, (∆m V > 4 mag.) which well explains the absence of detection with visual observations.

Conclusion

In the second semester of 2007, we obtained 283 new measurements of 279 visual binaries with PISCO in Merate, with an average accuracy of 0 ′′ .014 for the angular separation and 0 • .6 for the position angles. The total number of measurements made in Merate since 2004 now approaches 1600. This work is thus a good contribution to the continuing monitoring of long period visual binary systems, which is important for refining systemic stellar masses.

We obtained new orbital elements of the perturbation orbit Bb-P of Zeta Aqr and were able, for the first time, to compute their errors. Infrared observations with a 4-m class telescope would be very useful to constrain this orbit, or more basically to provide a clear confirmation of the presence of the invisible companion at the right position.

  Fig. 1Distribution of the angular separations of the 283 measurements of Table1(a), of the visual magnitudes of the corresponding objects (b) and of the differences of magnitude between the two components of those binaries (c).

  Fig.2HR diagram of the binaries measured in Table1, for which Hipparcos parallaxes were obtained with a relative error smaller than 50% (184 objects).

  Notes: In column 9, the exponent * indicates that the position angle θ could be determined without the 180 • ambiguity. Orbit references in column 12: Bdl2007a =Brendley & Hartkopf (2007), Cou1973b = Couteau (1973), Cve2007c = Cvetkovic (2007), Doc1984b = Docobo & Costa (1984), Doc1990a = Docobo & Costa (1990a), Doc1990b = Docobo & Costa (1990b), Doc2007d = Docobo & Ling (2007), Doc2008c = Docobo & Ling (2008), Hei1963a = Heintz (1963), Hei1965c = Heintz (1965), Hei1975b = Heintz (1975), Hei1976 = Heintz (1976), Hei1981b = Heintz (1981), Hei1984b = Heintz (1984), Hei1991 = Heintz (1991), Hei1993b = Heintz & Strom (1993), Hei1996a = Heintz (1996), Hei1997 = Heintz (1997), Hei2001 = Heintz (2001), Hop1964b = Hopmann (1964), Hrt1989 = Hartkopf et al. (1989), Hrt2001b = Hartkopf & Mason (2001), Hrt2008 = Hartkopf et al. (2008), Lin2008a = Ling (2008), Nov2006e = Novakovic & Todorovic (2006), Nov2007d = Novakovic (2007), Nov2008a = Novakovic (2008a), Nov2008b = Novakovic (2008b), Ole1993 = Olevic et al. (1993), Ole2001 = Olevic & Jovanovic (2001), Ole2004a = Olevic & Cvetkovic (2004), Pop1969b = Popovic (1969), Pop2000a = Popovic et al. (2000), Pbx2000b = Pourbaix (2000), Sca1983a = Scardia (1983), Sca2000b = Scardia et al. (2000b), Sca2001d = Scardia (2001), Sca2002d = Scardia et al. (2002), Sca2008a = Scardia et al. (2008a), Sca2008c = Scardia et al. (2008c), Sey2002 = Seymour (2002), Sod1999 = Söderhjelm (1999), Sta1977b = Starikova (1977), Sta1982b = Starikova (1982), Tok1984 = Tokovinin (1984), Val1981d = Valbousquet (1981), WSI2004a = Mason et al. (2004a), WSI2004b = Mason et al. (2004b), WSI2006b = Mason et al. (2006), Zul1977b = Zulevic (1977).

  Fig. 6 New orbit of ADS 15971: combination of the AB and Bb-P orbits (solid line) compared to all the photographic and speckle measurements (crosses): (a) position angle vs epoch, (b) angular separation vs epoch, and (c) orbit in the plane of the sky.

Table 1

 1 Measurements of binaries with PISCO and residuals with published orbits (begin).

	WDS	Name	ADS Epoch Fil. Eyep.	ρ	σρ	θ	σ θ Notes Orbit	∆ρ(O-C) ∆θ(O-C)
						(mm) (")	(")	( • ) ( • )		(")	( • )
	00029+4715 A 800	10 2007.967 W	20 1.584 0.010 113.3 * 0.9
	00093+7943 STF 2	102 2007.969 R	20 0.825 0.008 18.0 * 0.3	Hei1997	-0.02	1.1
	00116+5558 STF 7	143 2007.967 R	20 1.311 0.024 211.4 0.6
	00214+6700 STT 6AB	293 2007.969 R	20 0.653 0.008 154.4 * 0.6	Sca2000b	-0.02	-0.1
	"	"	"	"	"	"	"	"	" "	Ole2001	0.04	1.1
	00287+3718 A 1504AB	382 2007.972 R	20 0.593 0.008 42.9 1.0
	00366+5609 A 914	504 2007.969 R	20 0.440 0.008 22.6 1.0	Nov2008a	-0.02	-2.8
	00373+5801 BU 1097	515 2007.969 R	20 0.510 0.009 72.5 * 0.5
	00444+3337 STF 55	618 2007.969 W	20 2.212 0.011 329.7 * 0.3
	00480+5127 STF 59AB	659 2007.969 W	20 2.251 0.013 148.0 * 0.3
	00550+2338 STF 73AB	755 2007.972 R	20 1.011 0.008 319.8 * 0.6	Doc1990b	-0.02	-0.9
	01041+2635 COU 351	-2007.972 W	20 0.797 0.014 243.7 0.4
	01089+4512 AC 13AB	936 2007.964 W	20 0.595 0.009 263.3 * 0.9
	01213+1132 BU 4AB	1097 2007.964 R	20 0.594 0.016 107.7 * 1.0	Sca2001d	0.06	-1.0
	01360+0739 STF 138AB	1254 2007.972 R	20 1.717 0.012 58.2 * 0.3
	01501+2217 STF 174	1457 2007.964 R	20 2.841 0.041 164.6 * 0.3
	02211+2956 A 962	1792 2007.964 W	20 0.863 0.008 66.6 * 0.3
	14455+4223 STT 285AB	9378 2007.542 R	20 0.472 0.016 92.5 0.5	Cou1973b	-0.01	-0.2
	14484+2422 STF1884	9389 2007.542 R	20 2.161 0.011 54.9 * 0.3
	14497+0759 A 1110AB	9400 2007.543 R	20 0.683 0.008 243.5 * 0.8
	14497+4843 STF1890	9406 2007.537 R	20 2.692 0.013 45.1 * 0.3
	15018-0008 BU 348AB	9480 2007.545 R	10 0.487 0.003 108.0 * 0.3
	15038+4739 STF1909	9494 2007.534 R	20 1.842 0.009 57.9 * 0.3	Sod1999	0.01	0.1
	15056+1138 STF1907	9498 2007.537 W	20 0.913 0.008 350.6 0.3
	15075+0914 STF1910	9507 2007.567 R	20 4.038 0.031 211.2 * 0.3
	15116+1007 A 1116	9530 2007.545 W	20 0.808 0.008 50.9 * 0.5
	15183+2650 STF1932AB	9578 2007.534 R	20 1.621 0.010 262.1 * 0.5	Hei1965c	-0.01	-0.4
	15210+2104 HU 146	9600 2007.534 R	20 0.677 0.039 125.0 * 1.7
	15232+3017 STF1937AB	9617 2007.534 R	10 0.510 0.003 139.6 * 0.3	WSI2006b	-0.01	-1.3
	15245+3723 STF1938Ba,Bb 9626 2007.515 R	20 2.278 0.011	6.5 * 0.4	Sod1999	0.03	0.4
	15246+5413 HU 149	9628 2007.543 R	20 0.663 0.008 270.5 * 0.3
	15261+1810 STF1940	9634 2007.567 R	20 0.415 0.008 331.3 1.0
	15264+4400 STT 296AB	9639 2007.545 R	20 2.091 0.011 275.4 * 0.5
	15277+0606 STF1944	9647 2007.545 R	20 0.679 0.013 299.1 0.6
	15278+2906 JEF 1	-2007.515 R	10 0.288 0.003 144.2 * 0.3	Tok1984	0.01	3.1
	15300+2530 STF1950	9675 2007.570 R	20 3.297 0.016 90.7 * 0.5
	15329+3122 COU 610	-2007.515 R	10 0.815 0.004 199.1 * 0.3
	15348+1032 STF1954AB	9701 2007.551 R	20 4.025 0.020 172.8 * 0.3	WSI2004a	0.03	-0.1
	15360+3948 STT 298AB	9716 2007.551 R	20 0.979 0.008 175.4 * 0.3	Sod1999	-0.02	0.5
	15361+4849 HU 652	9718 2007.515 R	20 1.106 0.014 184.3 0.6
	15396+7959 STF1989	9769 2007.551 R	20 0.670 0.008 25.3 * 0.3	Hrt2008	0.04	1.3
	15405+1840 A 2076	9742 2007.537 R	20 0.716 0.008 184.1 0.7
	15413+5959 STF1969	9756 2007.543 W	20 0.955 0.008 27.4 * 0.4	Hei1975b	0.09	-0.7
	15427+2618 STF1967	9757 2007.537 R	10 0.697 0.006 112.9 * 0.3	Hrt1989	-0.03	0.2
	15432+1340 BU 619	9758 2007.570 R	20 0.671 0.009	1.3 * 0.9 Diff.
	15487+8337 STF2034	9853 2007.551 R	20 1.106 0.009 109.4 * 0.4
	15498+4431 BU 621	9802 2007.583 R	20 0.672 0.010 28.1 * 1.0
	15509+1911 A 2078	9809 2007.581 W	20 1.060 0.008 162.6 * 0.4
	15568+1229 STF1988	9850 2007.581 R	20 1.885 0.010 250.6 * 0.6
	16009+1316 STT 303AB	9880 2007.515 R	20 1.533 0.011 172.1 * 0.3
	16030+1359 STF2000	9904 2007.534 R	20 2.564 0.013 227.1 * 0.3
	16133+1332 STF2021AB	9969 2007.551 R	20 4.107 0.021 355.6 * 0.3	Hop1964b	0.02	0.9
	16135+7147 A 1136	9985 2007.570 W	20 0.727 0.011	8.6 * 1.0
	16137+4638 A 1642	9975 2007.534 R	20 0.753 0.009 182.1 * 0.5	Hrt2001b	0.03	-0.5
	16145+0531 STF2023	9974 2007.534 R	20 1.925 0.010 223.8

* 0.3

Table 4

 4 ADS 15971: O-C residuals obtained with the combination of the orbits AB and Bb-P. The symbol P indicates PISCO measurements. Measurements with "Ary" in Col. 4 were made by RWA with a 20 cm refractor.

	Epoch	∆ρ (O-C) ∆θ (O-C) Observer
		( ′′ )	( • )
	2002.531	0.039	0.511 WSI
	2002.800	0.006	-0.031 Ary
	2003.740	-0.049	0.621 Alz
	2003.804	0.020

Table 5

 5 New ephemeris of ADS 15971 -STF 2909 AB obtained with the combination of the AB and Bb-P orbits. .113 171.2 2015.0 2.313 164.7 2010.0 2.151 169.9 2016.0 2.337 163.9 2011.0 2.188 168.7 2017.0 2.358 163.0 2012.0 2.223 167.6 2018.0 2.375 162.1 2013.0 2.256 166.6 2019.0 2.389 161.2 2014.0 2.286 165.6 2020.0 2.401 160.3

	Epoch	ρ	θ	Epoch	ρ	θ
		( ′′ )	( • )		( ′′ )	( • )
	2009.0 2

Table 6

 6 Physical parameters of Zeta Aqr AB (a and M total ) derived from the new orbital elements.

	mV	Spec. type πHIP	a	M total
		(mas)	( ′′ ) (AU)	M⊙
	3.65 F3V+F6IV 35.68	3.38 94.7	3.6
		±1.24 ±0.02 ±3.4	±1.0

The last reported observation of STT 351AC in IC4 made in

2005.594 by WSI2006b, with ρ = 0 ′′ .58 and θ = 28. • 5 seems suspicious. It is inconsistent with the trend of the previous observations, especially for ρ.

www.an-journal.org

Acknowledgements. We thank the members of the United States Naval Observatory, Washington DC, for kindly sending on request some lists of measurements of visual binaries.

Table 3

New ephemeris of the perturbation orbit Bb-P of ADS 15971. Here ρ ′ is the angular separation between B and the center of mass of the Bb system. The angular separation of Bb is estimated at 3.0 ρ ′ (see text).

Epoch