# Fixed Point Method to Analyze Differences between Hipparcos and ICRF2

Authors: María José Martínez Usó; Francisco J. Marco Castillo; José Antonio López Ortí
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

From 1998, the International Astronomical Union (IAU) adopted a new Celestial Reference System: the International Celestial Reference System (ICRS). The first optical materialization was the Hipparcos catalogue, defining the Hipparcos Celestial Reference Frame (HCRF). The compilation of subsequent radio sources catalogues culminated in the current representation of the ICRF, the ICRF2 catalogue that is not sufficiently dense to cover all astrometrical purposes. Linking Hipparcos and ICRF2 is essential to uniformize the reference regardless of whether it is visible (HCRF) or not (ICRF). Many working groups provide their own complementary catalogs, some of whose sources are also in the ICRF2, with different reduction processes for positions. The point is that they provide information in more than one reference for a small number of objects. Some of these projects have been used by us to study the Hipparcos-ICRF2 differences: a certain number of couples of catalogs can be interrelated using a set of parameters. With these couples, we build a closed cycle with the same ending and departure couple. The parameters obtained from each couple affect the next; thus we have an iterative process whose fixed point is the solution that stabilizes it, providing a preliminary link for Hipparcos-ICRF2.


---

## Body

## 1. Introduction

A Reference System, particularly a Celestial Reference System, is a set of prescriptions and conventions together with the modelling required to define at any time a triad of axes. Until the end of the 20th century, Celestial Reference Systems were very much linked to the different earth parameters; thus it was necessary to update them every few years, but this situation changed when radio sources began to be used to define a really inertial, nearly invariant in time, system.Along the 20th century successive FK catalogs were compiled, motivated for two main causes: first, the technical improvements that provide more accurate positions; second, the need to remove different systematic errors without known physical origin. The final result was the FK5 [1] (there was also a FK6 [2], but it was never used as a materialization for any Celestial Reference System). As astronomers observed, FK5 presented also some characteristics that set it aside from the searched inertiality (see, e.g., [3–5]).It was for these reasons that the International Astronomical Union (IAU) decided in 1997, at its 23rd general assembly, that the IAU Celestial Reference System was the International Celestial Reference System (ICRS) from January 1, 1998, in replacement of the FK5. The first materialization arose from Hipparcos mission that provided the catalogue of the same name defining the Hipparcos Celestial Reference Frame (HCRF) [6]. To this aim, several radio sources were considered to align the HCRF (see [7] and its references for more details), approximately, with the theoretical ICRS. It is necessary to highlight that this catalogue contains approximately 118000 visible wavelengths stars. In order to assure the continuity of the Celestial Reference Systems, different authors have carried out the study of the relations between HCRF and the FK5 (see [3–5, 8] among others).As previously stated, the new system is expected to be inertial (but it is not) and it is for this reason that several observational programs have been developed on the basis of Very Large Baseline Array (VLBA) in order to obtain the positions of celestial objects without significant proper motions. These programs have provided the ICRF series of nonvisible wavelengths catalogues [9], of which the ICRF2 [10] is currently in use. The ICRF series catalogs are not dense enough to cover all astrometrical purposes and, for this reason, many working groups provide their own complementary catalogs, some of whose positions are also in the ICRF2 but do not coincide numerically due to several reasons such as the measure methods, the different instruments employed for the compilation, and the method used to reduce the positions.Linking Hipparcos and ICRF2 is essential in order to uniformize the reference independently of whether the source is visible (HCRF) or not (ICRF). In the future, the Gaia project (seehttp://sci.esa.int/gaia/) is expected to provide a practical solution to this problem but, for the present time, different efforts have been carried out to obtain a link between both references. Every observational program has its own reduction processes for positions, but the question is that each of them provides information in more than one reference for a small number of objects. Currently the method to relate two catalogs with common sources is to consider the differences in the data (residuals) and the search of an adjustment for such residuals. With the word “couple,” we refer to this process. Next, we list the couples that we use along the present work.(a)
Assafin et al. [11] provide 300 radio sources, with positions observed in their Rio Survey (Rio henceforth) program, together with the data that such radio sources have in the ICRF2 (Couple Rio-ICRF2).
(b)
We take 82 of the radio sources from Rio that are also contained in Hipparcos (Couple Rio-Hipparcos).
(c)
The comparison of ICRF2-Zach14, where Zach14 refers to N. Zacharias and M. I. Zacharias [12], provides 682 positions obtained after reducing the observations with respect to the Tycho-2 [6] (i.e., linked to, but much denser than the HCRF) and using the UCAC4 ([13]).
(d)
The case UCAC4-Hipparcos is different from the three previous ones, because it contains high quality positions in visible wavelength. We have used 104077 common sources. That is why this couple of catalogs has the greatest influence in the level of reliability of our results.The aim of this paper is the study of the current situation of the link between HCRF and ICRF2, using the four couples of catalogs (a), (b), (c), and (d) previously defined to obtain the links among them. Every relation affects the next pair of catalogues, of which we obtain the corresponding relations. This process will continue until the closed cycle made up by the couples has been completed:(Rio-ICRF2)-(Rio-Hipparcos)-(Hipparcos-UCAC4)-(ICRF2-Zach14)where we identify UCAC4 and Zach14 (see [12]). The departure point is the couple Rio-ICRF2 and the closed cycle is completed when we reach Rio-ICRF2 again using the procedures that we explain later.In each of the substeps, the relations are given by infinitesimal rotations and deformations between the references defined by each set of data. Since the process is a closed circle, as soon as the parameters (infinitesimal rotations and deformations) are fixed for one of the couples, the process finishes and the fixed point of the problem has been found. We give a much more detailed description of the method in Section2.All problems come from the fact that, although we search the parameters of a certain functional development (we work with vector fields whose components belong to theL
2
(
S
2
) Hilbert space), in practice for each couple of references we have only a discrete set of unnecessary homogeneously distributed points on the sphere. This implies that the method of the discrete least squares (DLS henceforth), which is currently used to determine the parameters, is highly unstable (due to an ill conditioned normal matrix. Equivalently, orthogonality is not preserved when we go from the theoretic-continuous case to the practical-discrete) and inefficient (if a higher order of development is needed, all the coefficients of lower order must be recomputed). See, for example, [8, 14]. We consider these topics in Section 3.In Section4, we implement the method described in Section 2, using the techniques of Section 3 and we present the results obtained, that is to say, the fixed point of the problem.We finish the paper with a brief discussion about the results and perspectives of the present work.


## 2. Searching for the Fixed Point for the Hipparcos-ICRF2 Corrections

The fixed point implementation applied to the different couples of catalogs that have been listed in Section1 is described in Table 1 and runs as shown in Figure 1 in a closed circle until the relationships stabilize.Table 1
Fixed point iterative process (see explanation in text).
Rio-ICRF2
Rio-Hip
Hip-UCAC4
ICRF2-Zach14
(Starting point)Θ
R
0
Rio +Θ
R
0 → Θ
H
0
Hip +Θ
H
0 → Θ
U
0
Zach14 +Θ
U
0 → Θ
I
0
ICRF2 +Θ
I
0 → Θ
R
1
Rio +Θ
R
1 → Θ
H
1
Hip +Θ
H
1 → Θ
U
1
Zach14 +Θ
U
1 → Θ
I
1
⋯
⋯
⋯
⋯
Θ
R
f
Rio +Θ
R
f → Θ
H
f
Hip +Θ
H
f → Θ
U
f
Zach14 +Θ
U
f → Θ
I
fFigure 1
Fixed point iterative process (see explanation in text).The complete process could be subdivided in the following steps.Step
1
(
a
). Initialization of the process from Rio-ICRF2: we use u to represent the nonparametrical adjustment for the differences in Rio-ICRF2; V
F
u is the first order vector field development, with the unknowns Θ
R
0. The integral on the sphere is a g
u
(
Θ
) function and we obtain its minimum:
(1)
min
⁡
Θ
g
u
(
Θ
)
=
min
⁡
Υ
,
Θ
∫
S
2
[
u
0
(
α
,
δ
)
-
V
F
u
0
(
Θ
,
α
,
δ
)
]
2
d
σ
⟶
Θ
R
0
.
Let G
u be the application of the previously obtained adjustment on one of the catalogs of the couple. Depending on which catalog we apply the correction, this would be used with a plus or a minus sign.Step
1
(
b
). Initialization for Río-Hipparcos: analogously, we denote the adjustment and the field by v and V
F
v, respectively. In v we have applied G
u to the Rio positions. Then g
v
(
Θ
) is the corresponding integral on the sphere whose minimum is searched in this step
(2)
min
⁡
Θ
g
v
(
Θ
)
=
min
⁡
Θ
∫
S
2
[
v
0
(
α
,
δ
)
-
V
F
v
0
(
Θ
,
α
,
δ
)
]
2
d
σ
⟶
Θ
H
0
.
We denote by G
v the application of this adjustment on one of the catalogs, in a similar way to Step 1
(
a
). Analogous procedure is carried out for the rest of couples.Step
1
(
c
). Let us return to Rio-ICRF2 taking into account all the previous G
u
,
G
v
,
… already considered and used. For simplicity, we denote by G the composition of all the transformations; then u
1
=
G
(
u
) and we proceed with the adjustment by
(3)
min
⁡
Θ
f
(
Θ
)
=
min
⁡
Θ
∫
S
2
[
u
1
(
α
,
δ
)
-
V
F
u
1
(
Θ
,
α
,
δ
)
]
2
d
σ
⟶
Θ
H
1
.
Step2. Iteration continues until the coefficients stabilize.


## 3. Methods and Mathematical Models

The nonhomogeneity in the spatial distribution causes functional orthogonality not to turn into algebraical orthogonality in the discrete case if DLS is used. This problem is especially serious when we estimate high order harmonics.In previous papers ([8, 15]), we have highlighted these potential problems using DLS and we have proposed an alternative method, called mixed method (MM henceforth) that overcomes these errors. In brief, the MM involves the following steps:(a)
a selection of the estimation off
~ at any point by means of a kernel nonparametrical method (KNP) (see Section 3.1). Other methods of estimation could also be employed;
(b)
the computation off
~ values on an equispaced selected grid;
(c)
the application ofa
i
=
〈
f
,
P
i
〉
/
〈
P
i
,
P
i
〉 discretizing the integrals to obtain the searched coefficients of a parametrical model; this parametrical model is given by a series development on a complete and orthogonal basis in a certain Hilbert space; see Section 3.2.In addition, our MM has the following advantages (for a more complete exposition see [15]).(i)
We can calculate the estimation at equally spaced points on the sphere. These computations are carried out using the initial data (obtained from the catalogs and unnecessary homogeneously distributed). The grid of equal spaced points employed in the KNP is used in the functional adjustment.
(ii)
The (i) property allows an increase in the order of harmonics without having to recalculate the coefficients of lower orders again because they are fixed. This enables a sequential order increase.
(iii)
Any additional assignation of weights is not required because they are included in the method itself, via the properties of the kernel.
(iv)
The calculation of each coefficient requires the computation of integrals on the sphere which can be numerically carried out with sufficient accuracy and low computational cost (because orthogonality is preserved).


### 3.1. Kernel Regression Methods

In this subsection we expose in detail the kernel regression method in the one-dimensional case. Then, we extend the formulas (already without providing demonstrations) for the spherical surface and the solid sphere. This later case is used for the differences in Hip-UCAC4 described in Section3.2.


#### 3.1.1. Unidimensional Kernel Regression

The DLS regression on discrete data(
x
i
,
y
i
)
,
1
≤
i
≤
n is carried out as a finite and linear up to order r combination of orthogonal functions ϕ
j in a certain domain with respect to an inner product <
,
>. It implies obtaining the α
j
^
,
1
≤
j
≤
r values minimizing
(4)
∑
i
=
1
n
[
y
i
-
∑
j
=
1
r
α
j
ϕ
j
(
x
i
)
]
2
.
Alternatively, we propose building upon the general definition of a one-dimensional regression curve from data (see [16]) and then suppose that a generic, although unknown, relationship is fulfilled:
(5)
y
i
=
m
(
x
i
)
+
ε
i
,
where the errors ε
i are N
(
0
,
σ
). By definition, the regression curve is given by
(6)
m
(
x
)
=
E
[
Y
∣
X
=
x
]
=
∫
y
f
Y
∣
X
=
x
(
y
∣
X
=
x
)
d
y
=
∫
y
f
(
X
,
Y
)
(
x
,
y
)
f
X
(
x
)
d
y
,
where f
X
(
x
), f
(
x
,
y
), and f
(
y
∣
x
) are the marginal density of X, the joint density of X and Y, and the conditional density of Y given X, respectively. A kernel estimate of f
X
(
x
) is
(7)
f
^
X
(
x
)
=
1
n
h
x
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
,
and a kernel estimate of f
(
x
,
y
) is
(8)
f
^
(
x
,
y
)
=
1
n
h
x
h
y
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
K
y
(
y
-
y
i
h
y
)
,
where K
(
x
) is called a kernel function if K
≥
0, ∫
K
(
u
)
d
u
=
1, ∫
u
K
(
u
)
d
u
=
0, and ∫
u
2
K
(
u
)
d
u
<
∞. Applying these kernel properties together with formulas (7) and (8), we obtain a linear function of y called the Nadaraya-Watson kernel estimator:
(9)
m
^
N
W
(
x
)
=
∑
i
=
1
n
K
x
(
(
x
-
x
i
)
/
h
x
)
y
i
∑
j
=
1
n
K
x
(
(
x
-
x
j
)
/
h
x
)
≡
∑
i
=
1
n
ϖ
i
y
i
,
where w
i are the weights given by
(10)
ϖ
i
=
1
n
h
x
K
x
(
(
x
-
x
i
)
/
h
x
)
f
^
X
(
x
)
.
This expression has been obtained by direct application of the properties of a general kernel K and the general properties of the continuous density functions. Further, this kind of method is called kernel nonparametric (KNP henceforth) method.At this point, we have an estimatorm
^
(
x
) of the unknown function that may be computed on any point of the domain. On this basis, it is possible to obtain an estimate for each coefficient α
j provided that, on the one hand, the function to adjust is approximated by means of the m
^
(
x
) continuous estimation and, on the other hand, this function can be developed in series of the orthogonal functions {
ϕ
j
}
1
≤
j
≤
r:
(11)
α
^
j
=
〈
m
^
,
ϕ
j
〉
〈
ϕ
j
,
ϕ
j
〉
.
The integrals, included in the previous formula as inner product, are computed by means of a numerical method, using the grid of equally spaced points. On these same points of the grid the values of the function to adjust are calculated using the regression m
^
(
x
). If the real data correspond accidentally to equispaced abscises, the results obtained from (11) using MM are similar to the results obtained applying DLS from (4).To implement the method on a certain set of discrete points, all we need to do is to choose a kernel, for instance, the Epanechnikov kernel, defined by(12)
K
(
x
)
=
{
3
4
(
1
-
x
2
)
|
x
|
≤
1
0
|
x
|
>
1
.
The selection of this kernel is justified on the basis that this is the most efficient among a large range of kernels (see [16] again). Finally, the selection of the h value, called bandwidth, is done using expressions that minimize the asymptotic mean integrate square error (AMISE) over the whole domain (see [16–18] or [19], for a more detailed exposition of the regression nonparametrical methods).


#### 3.1.2. KNP on the Sphere and Solid KNP

Next we describe the two-dimensional KNP method on the sphere. Along this paragraph, we use(
Δ
α
cos
⁡
δ
)
i and Δ
δ
i for the numerical differences in positions for the common object i of each couple of catalogs, α being the right ascension and δ the declination. The use of Δ
α
cos
⁡
δ and Δ
δ for the calculation of the coefficients of the model is commonly carried out employing the n individual residuals to estimate the parameters of the selected adjustment models m
α, m
δ:
(13)
∑
i
=
1
n
{
[
(
Δ
α
cos
⁡
δ
)
i
-
m
α
(
α
i
,
δ
i
)
]
2
+
[
(
Δ
δ
)
i
-
m
δ
(
α
i
,
δ
i
)
]
2
}
.
If necessary, it is also possible to introduce weights depending on the statistical characteristics of the data.Depending on the characteristics of the searched function it may be required to meet certain mathematical hypothesis of regularity. In this case, the vector field should have integrable square on the sphereL
2
(
S
2
). Thus, we want to find
(14)
min
⁡
∫
{
[
Δ
α
cos
⁡
δ
(
α
,
δ
)
-
m
α
(
α
,
δ
)
]
2
n
n
n
n
n
n
n
+
[
Δ
δ
(
α
,
δ
)
-
m
δ
(
α
,
δ
)
]
2
}
d
S
.
As in the one-dimensional case, we deal with random variables and their corresponding mathematical expectations (replacing means) and variances, calculated as integrals. To this aim, the use of the probability density function of each random variable is required. As these densities are unknown, the employment of a kernel estimator is especially useful.In this sense in formula (14) Δ
α
cos
⁡
δ, Δ
δ are the nonparametrical kernel adjustment computations of the random variables given by their names.Let us remember that nonparametric kernel adjustments compute the conditional expectation of a certain random variable that depends on others. For example, ifX is the random variable (
Δ
α
cos
⁡
δ, Δ
δ
) the KNP method consists of finding
(15)
E
[
X
∣
(
α
,
δ
)
]
=
∫
D
x
f
(
x
∣
(
α
,
δ
)
)
d
x
=
∫
D
x
f
(
x
,
α
,
δ
)
f
(
α
,
δ
)
(
α
,
δ
)
d
x
,
where D is the domain of X, f
(
x
,
α
,
δ
) is the joint density function of the three random variables, and f
(
α
,
δ
)
(
α
,
δ
) is the marginal density. All these elements may be unknown so they may be approximated using
(16)
f
^
(
x
,
α
,
δ
)
=
1
n
h
x
h
α
h
sin
δ
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
×
K
α
(
α
-
α
i
h
α
)
K
δ
(
sin
δ
-
sin
δ
i
h
sin
δ
)
,
and the condition
(17)
1
4
π
μ
(
D
)
∫
D
∫
S
2
f
^
(
x
,
α
,
δ
)
cos
⁡
δ
d
x
d
α
d
δ
=
1
must be fulfilled. We proceed analogously for the marginal density. Taking the same kernel for all the random variables and considering their properties we reach a similar expression to the one-dimensional Nadaraya-Watson one, but on the sphere:
(18)
m
X
(
α
,
δ
)
=
∑
i
=
1
n
w
i
x
i
,
w
i
=
K
α
(
(
α
-
α
i
)
/
h
α
)
K
δ
(
(
sin
δ
-
sin
δ
i
)
/
h
sin
δ
)
∑
j
=
1
n
K
α
(
(
α
-
α
j
)
/
h
α
)
K
δ
(
(
sin
δ
-
sin
δ
j
)
/
h
sin
δ
)
,
where K is given in (12). It can be demonstrated that the lesser h is the more steeply the adjustment is. A large h provides smoother results. The theoretical optimum values have been studied in [16].Analogously, if we introduce a radial component, we obtain an estimation of the residuals on the solid sphere that generalizes formula (18).


### 3.2. Vector Spherical Developments

Vector spherical harmonics have been widely used in several astronomical lines of work. See in particular [5] within the context of the Hipparcos-FK5 comparison. This issue was also studied by us in [8], where we truncate the VSH development up to first order, which is enough to determine the rotational and deformation parameters. VSH was also studied by us in the context of the Hipparcos spin in [3].Let us consider the vector field over the celestial sphere (in the case of a surface vector field, we take a constant radiusr
=
1 and the last summand disappears):
(19)
Δ
X
≡
V
(
r
,
α
,
δ
)
=
V
α
(
r
,
α
,
δ
)
e
α
+
V
δ
(
r
,
α
,
δ
)
e
δ
+
V
r
(
r
,
α
,
δ
)
e
r
=
Δ
α
cos
⁡
δ
e
a
+
Δ
δ
e
d
+
Δ
r
e
r
,
V
α, V
δ, and V
r being the scalar fields of the vector field and e
α, e
δ the unitary vectors in the tangent plane and in the directions of the right ascension and declination, respectively, and e
r is the unitary vector position.On the other hand, provided that we are on the surface of the unitary sphere, the only vector spherical harmonics involved are theS
l
,
m spheroidal spherical harmonics and the T
l
,
k toroidal spherical harmonics. Anyway, in general, there is an additional radial component denoted by R
l
,
k. Thus, we suppose that the vector fieldV has a mathematical development:
(20)
V
(
r
,
α
,
δ
)
=
∑
k
≥
1
∑
l
=
-
k
k
[
r
l
,
m
R
l
,
m
+
t
l
,
m
T
l
,
m
+
s
l
,
m
S
l
,
m
]
.
The definitions of the vectors of the orthogonal and complete basis for the L
2(S
2) Hilbert space of the vector fields on the sphere are
(21)
R
n
,
m
=
r
r
Y
n
,
m
;
S
n
,
m
=
r
∇
Y
n
,
m
;
T
n
,
m
=
-
r
×
∇
Y
n
,
m
,
where Y
n
,
m
(
n
≥
0
,
-
n
≤
m
≤
n
) are the usual surface spherical harmonics. As we have previously stated, if we are interested only in the spheroidal and toroidal components, then, the estimations may be computed, due to the functional orthogonality, as
(22)
t
^
l
,
m
=
∫
S
2
V
·
T
k
,
l
d
σ
∥
T
k
,
l
∥
2
;
s
^
l
,
m
=
∫
S
2
V
·
S
k
,
l
d
σ
∥
S
k
,
l
∥
2
.
The denominators are exactly calculated, whereas for the numerators an estimation is obtained using the proposed MM: for the calculation of the components of the vector fieldV at regularly spaced points on the sphere we can use the kernel regression method, which is computationally efficient and, in addition, it is rather accurate for the problem that we are discussing. It is important to emphasize that, once the adjustment has been established forV on the set of points of the sphere, this same set can be used for the numerical integration of the numerators up to any order of the development. Thus, we can easily calculate independently the estimations for higher orders of harmonics.Truncation of (20) up to the first order, together with formula (19), requires (see, e.g., [8]) fulfilling the following conditions (23) in order to verify the compatibility of the system:
(23)
Δ
α
cos
⁡
δ
=
[
t
1,0
cos
⁡
δ
-
t
1,1
sin
α
sin
δ
-
t
1
,
-
1
cos
⁡
α
sin
δ
]
+
[
s
1,1
cos
⁡
α
-
s
1
,
-
1
sin
α
]
,
Δ
δ
=
[
-
t
1,1
cos
⁡
α
+
t
1
,
-
1
sin
α
]
+
[
-
s
1,1
sin
α
sin
δ
-
s
1
,
-
1
cos
⁡
α
sin
δ
+
s
1,0
cos
⁡
δ
]
,
where it is evident that t
1
,
-
1
=
ε
x, t
1
,
1
=
ε
y, and t
1
,
0
=
ε
z, ε
x, ε
y, and ε
z being the infinitesimal rotations around the corresponding axis.


### 3.3. Average for the Solid Hipparcos-UCAC4 Residuals

In order to use Hipparcos and UCAC4, we need the data for the common stars that we denote by (r
i, α
i, δ
i), 1
≤
i
≤
N. Notice that we also use the module of the radius vector because good parallax data are available. To reduce this solid case to the surface case, we proceed in the following way: we classify these differences in concentric zones with radii:
(24)
R
°
<
R
¹
<
R
²
<
⋯
<
R
k
.
In this way, we have for each I
(
j
)
=
[
R
(
j
)
,
R
(
j
+
1
)
) interval the (
Δ
α
r
cos
⁡
δ
r
,
Δ
δ
r
)
(
j
) differences, where superindex (j) denotes the slice where the differences are. For each I
(
j
) we have a vector field of differences projected on the R
-
j
,
j
+
1 point which is the arithmetical mean between R
(
j
) and R
(
j
+
1
):
(25)
V
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
=
V
α
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
e
α
+
V
δ
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
e
δ
,
e
α, e
δ being the unitary vectors in the directions of RA and DEC, respectively, that determine the plane tangent to the sphere on each point. Let us call averaged vector field to the field obtained by numerically integrating
(26)
V
(
α
,
δ
)
=
∫
R
-
0,1
R
-
k
-
1
,
k
[
V
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
]
d
r
.
Once this field is obtained, we proceed as in Section 3.2.


## 3.1. Kernel Regression Methods

In this subsection we expose in detail the kernel regression method in the one-dimensional case. Then, we extend the formulas (already without providing demonstrations) for the spherical surface and the solid sphere. This later case is used for the differences in Hip-UCAC4 described in Section3.2.


### 3.1.1. Unidimensional Kernel Regression

The DLS regression on discrete data(
x
i
,
y
i
)
,
1
≤
i
≤
n is carried out as a finite and linear up to order r combination of orthogonal functions ϕ
j in a certain domain with respect to an inner product <
,
>. It implies obtaining the α
j
^
,
1
≤
j
≤
r values minimizing
(4)
∑
i
=
1
n
[
y
i
-
∑
j
=
1
r
α
j
ϕ
j
(
x
i
)
]
2
.
Alternatively, we propose building upon the general definition of a one-dimensional regression curve from data (see [16]) and then suppose that a generic, although unknown, relationship is fulfilled:
(5)
y
i
=
m
(
x
i
)
+
ε
i
,
where the errors ε
i are N
(
0
,
σ
). By definition, the regression curve is given by
(6)
m
(
x
)
=
E
[
Y
∣
X
=
x
]
=
∫
y
f
Y
∣
X
=
x
(
y
∣
X
=
x
)
d
y
=
∫
y
f
(
X
,
Y
)
(
x
,
y
)
f
X
(
x
)
d
y
,
where f
X
(
x
), f
(
x
,
y
), and f
(
y
∣
x
) are the marginal density of X, the joint density of X and Y, and the conditional density of Y given X, respectively. A kernel estimate of f
X
(
x
) is
(7)
f
^
X
(
x
)
=
1
n
h
x
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
,
and a kernel estimate of f
(
x
,
y
) is
(8)
f
^
(
x
,
y
)
=
1
n
h
x
h
y
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
K
y
(
y
-
y
i
h
y
)
,
where K
(
x
) is called a kernel function if K
≥
0, ∫
K
(
u
)
d
u
=
1, ∫
u
K
(
u
)
d
u
=
0, and ∫
u
2
K
(
u
)
d
u
<
∞. Applying these kernel properties together with formulas (7) and (8), we obtain a linear function of y called the Nadaraya-Watson kernel estimator:
(9)
m
^
N
W
(
x
)
=
∑
i
=
1
n
K
x
(
(
x
-
x
i
)
/
h
x
)
y
i
∑
j
=
1
n
K
x
(
(
x
-
x
j
)
/
h
x
)
≡
∑
i
=
1
n
ϖ
i
y
i
,
where w
i are the weights given by
(10)
ϖ
i
=
1
n
h
x
K
x
(
(
x
-
x
i
)
/
h
x
)
f
^
X
(
x
)
.
This expression has been obtained by direct application of the properties of a general kernel K and the general properties of the continuous density functions. Further, this kind of method is called kernel nonparametric (KNP henceforth) method.At this point, we have an estimatorm
^
(
x
) of the unknown function that may be computed on any point of the domain. On this basis, it is possible to obtain an estimate for each coefficient α
j provided that, on the one hand, the function to adjust is approximated by means of the m
^
(
x
) continuous estimation and, on the other hand, this function can be developed in series of the orthogonal functions {
ϕ
j
}
1
≤
j
≤
r:
(11)
α
^
j
=
〈
m
^
,
ϕ
j
〉
〈
ϕ
j
,
ϕ
j
〉
.
The integrals, included in the previous formula as inner product, are computed by means of a numerical method, using the grid of equally spaced points. On these same points of the grid the values of the function to adjust are calculated using the regression m
^
(
x
). If the real data correspond accidentally to equispaced abscises, the results obtained from (11) using MM are similar to the results obtained applying DLS from (4).To implement the method on a certain set of discrete points, all we need to do is to choose a kernel, for instance, the Epanechnikov kernel, defined by(12)
K
(
x
)
=
{
3
4
(
1
-
x
2
)
|
x
|
≤
1
0
|
x
|
>
1
.
The selection of this kernel is justified on the basis that this is the most efficient among a large range of kernels (see [16] again). Finally, the selection of the h value, called bandwidth, is done using expressions that minimize the asymptotic mean integrate square error (AMISE) over the whole domain (see [16–18] or [19], for a more detailed exposition of the regression nonparametrical methods).


### 3.1.2. KNP on the Sphere and Solid KNP

Next we describe the two-dimensional KNP method on the sphere. Along this paragraph, we use(
Δ
α
cos
⁡
δ
)
i and Δ
δ
i for the numerical differences in positions for the common object i of each couple of catalogs, α being the right ascension and δ the declination. The use of Δ
α
cos
⁡
δ and Δ
δ for the calculation of the coefficients of the model is commonly carried out employing the n individual residuals to estimate the parameters of the selected adjustment models m
α, m
δ:
(13)
∑
i
=
1
n
{
[
(
Δ
α
cos
⁡
δ
)
i
-
m
α
(
α
i
,
δ
i
)
]
2
+
[
(
Δ
δ
)
i
-
m
δ
(
α
i
,
δ
i
)
]
2
}
.
If necessary, it is also possible to introduce weights depending on the statistical characteristics of the data.Depending on the characteristics of the searched function it may be required to meet certain mathematical hypothesis of regularity. In this case, the vector field should have integrable square on the sphereL
2
(
S
2
). Thus, we want to find
(14)
min
⁡
∫
{
[
Δ
α
cos
⁡
δ
(
α
,
δ
)
-
m
α
(
α
,
δ
)
]
2
n
n
n
n
n
n
n
+
[
Δ
δ
(
α
,
δ
)
-
m
δ
(
α
,
δ
)
]
2
}
d
S
.
As in the one-dimensional case, we deal with random variables and their corresponding mathematical expectations (replacing means) and variances, calculated as integrals. To this aim, the use of the probability density function of each random variable is required. As these densities are unknown, the employment of a kernel estimator is especially useful.In this sense in formula (14) Δ
α
cos
⁡
δ, Δ
δ are the nonparametrical kernel adjustment computations of the random variables given by their names.Let us remember that nonparametric kernel adjustments compute the conditional expectation of a certain random variable that depends on others. For example, ifX is the random variable (
Δ
α
cos
⁡
δ, Δ
δ
) the KNP method consists of finding
(15)
E
[
X
∣
(
α
,
δ
)
]
=
∫
D
x
f
(
x
∣
(
α
,
δ
)
)
d
x
=
∫
D
x
f
(
x
,
α
,
δ
)
f
(
α
,
δ
)
(
α
,
δ
)
d
x
,
where D is the domain of X, f
(
x
,
α
,
δ
) is the joint density function of the three random variables, and f
(
α
,
δ
)
(
α
,
δ
) is the marginal density. All these elements may be unknown so they may be approximated using
(16)
f
^
(
x
,
α
,
δ
)
=
1
n
h
x
h
α
h
sin
δ
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
×
K
α
(
α
-
α
i
h
α
)
K
δ
(
sin
δ
-
sin
δ
i
h
sin
δ
)
,
and the condition
(17)
1
4
π
μ
(
D
)
∫
D
∫
S
2
f
^
(
x
,
α
,
δ
)
cos
⁡
δ
d
x
d
α
d
δ
=
1
must be fulfilled. We proceed analogously for the marginal density. Taking the same kernel for all the random variables and considering their properties we reach a similar expression to the one-dimensional Nadaraya-Watson one, but on the sphere:
(18)
m
X
(
α
,
δ
)
=
∑
i
=
1
n
w
i
x
i
,
w
i
=
K
α
(
(
α
-
α
i
)
/
h
α
)
K
δ
(
(
sin
δ
-
sin
δ
i
)
/
h
sin
δ
)
∑
j
=
1
n
K
α
(
(
α
-
α
j
)
/
h
α
)
K
δ
(
(
sin
δ
-
sin
δ
j
)
/
h
sin
δ
)
,
where K is given in (12). It can be demonstrated that the lesser h is the more steeply the adjustment is. A large h provides smoother results. The theoretical optimum values have been studied in [16].Analogously, if we introduce a radial component, we obtain an estimation of the residuals on the solid sphere that generalizes formula (18).


## 3.1.1. Unidimensional Kernel Regression

The DLS regression on discrete data(
x
i
,
y
i
)
,
1
≤
i
≤
n is carried out as a finite and linear up to order r combination of orthogonal functions ϕ
j in a certain domain with respect to an inner product <
,
>. It implies obtaining the α
j
^
,
1
≤
j
≤
r values minimizing
(4)
∑
i
=
1
n
[
y
i
-
∑
j
=
1
r
α
j
ϕ
j
(
x
i
)
]
2
.
Alternatively, we propose building upon the general definition of a one-dimensional regression curve from data (see [16]) and then suppose that a generic, although unknown, relationship is fulfilled:
(5)
y
i
=
m
(
x
i
)
+
ε
i
,
where the errors ε
i are N
(
0
,
σ
). By definition, the regression curve is given by
(6)
m
(
x
)
=
E
[
Y
∣
X
=
x
]
=
∫
y
f
Y
∣
X
=
x
(
y
∣
X
=
x
)
d
y
=
∫
y
f
(
X
,
Y
)
(
x
,
y
)
f
X
(
x
)
d
y
,
where f
X
(
x
), f
(
x
,
y
), and f
(
y
∣
x
) are the marginal density of X, the joint density of X and Y, and the conditional density of Y given X, respectively. A kernel estimate of f
X
(
x
) is
(7)
f
^
X
(
x
)
=
1
n
h
x
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
,
and a kernel estimate of f
(
x
,
y
) is
(8)
f
^
(
x
,
y
)
=
1
n
h
x
h
y
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
K
y
(
y
-
y
i
h
y
)
,
where K
(
x
) is called a kernel function if K
≥
0, ∫
K
(
u
)
d
u
=
1, ∫
u
K
(
u
)
d
u
=
0, and ∫
u
2
K
(
u
)
d
u
<
∞. Applying these kernel properties together with formulas (7) and (8), we obtain a linear function of y called the Nadaraya-Watson kernel estimator:
(9)
m
^
N
W
(
x
)
=
∑
i
=
1
n
K
x
(
(
x
-
x
i
)
/
h
x
)
y
i
∑
j
=
1
n
K
x
(
(
x
-
x
j
)
/
h
x
)
≡
∑
i
=
1
n
ϖ
i
y
i
,
where w
i are the weights given by
(10)
ϖ
i
=
1
n
h
x
K
x
(
(
x
-
x
i
)
/
h
x
)
f
^
X
(
x
)
.
This expression has been obtained by direct application of the properties of a general kernel K and the general properties of the continuous density functions. Further, this kind of method is called kernel nonparametric (KNP henceforth) method.At this point, we have an estimatorm
^
(
x
) of the unknown function that may be computed on any point of the domain. On this basis, it is possible to obtain an estimate for each coefficient α
j provided that, on the one hand, the function to adjust is approximated by means of the m
^
(
x
) continuous estimation and, on the other hand, this function can be developed in series of the orthogonal functions {
ϕ
j
}
1
≤
j
≤
r:
(11)
α
^
j
=
〈
m
^
,
ϕ
j
〉
〈
ϕ
j
,
ϕ
j
〉
.
The integrals, included in the previous formula as inner product, are computed by means of a numerical method, using the grid of equally spaced points. On these same points of the grid the values of the function to adjust are calculated using the regression m
^
(
x
). If the real data correspond accidentally to equispaced abscises, the results obtained from (11) using MM are similar to the results obtained applying DLS from (4).To implement the method on a certain set of discrete points, all we need to do is to choose a kernel, for instance, the Epanechnikov kernel, defined by(12)
K
(
x
)
=
{
3
4
(
1
-
x
2
)
|
x
|
≤
1
0
|
x
|
>
1
.
The selection of this kernel is justified on the basis that this is the most efficient among a large range of kernels (see [16] again). Finally, the selection of the h value, called bandwidth, is done using expressions that minimize the asymptotic mean integrate square error (AMISE) over the whole domain (see [16–18] or [19], for a more detailed exposition of the regression nonparametrical methods).


## 3.1.2. KNP on the Sphere and Solid KNP

Next we describe the two-dimensional KNP method on the sphere. Along this paragraph, we use(
Δ
α
cos
⁡
δ
)
i and Δ
δ
i for the numerical differences in positions for the common object i of each couple of catalogs, α being the right ascension and δ the declination. The use of Δ
α
cos
⁡
δ and Δ
δ for the calculation of the coefficients of the model is commonly carried out employing the n individual residuals to estimate the parameters of the selected adjustment models m
α, m
δ:
(13)
∑
i
=
1
n
{
[
(
Δ
α
cos
⁡
δ
)
i
-
m
α
(
α
i
,
δ
i
)
]
2
+
[
(
Δ
δ
)
i
-
m
δ
(
α
i
,
δ
i
)
]
2
}
.
If necessary, it is also possible to introduce weights depending on the statistical characteristics of the data.Depending on the characteristics of the searched function it may be required to meet certain mathematical hypothesis of regularity. In this case, the vector field should have integrable square on the sphereL
2
(
S
2
). Thus, we want to find
(14)
min
⁡
∫
{
[
Δ
α
cos
⁡
δ
(
α
,
δ
)
-
m
α
(
α
,
δ
)
]
2
n
n
n
n
n
n
n
+
[
Δ
δ
(
α
,
δ
)
-
m
δ
(
α
,
δ
)
]
2
}
d
S
.
As in the one-dimensional case, we deal with random variables and their corresponding mathematical expectations (replacing means) and variances, calculated as integrals. To this aim, the use of the probability density function of each random variable is required. As these densities are unknown, the employment of a kernel estimator is especially useful.In this sense in formula (14) Δ
α
cos
⁡
δ, Δ
δ are the nonparametrical kernel adjustment computations of the random variables given by their names.Let us remember that nonparametric kernel adjustments compute the conditional expectation of a certain random variable that depends on others. For example, ifX is the random variable (
Δ
α
cos
⁡
δ, Δ
δ
) the KNP method consists of finding
(15)
E
[
X
∣
(
α
,
δ
)
]
=
∫
D
x
f
(
x
∣
(
α
,
δ
)
)
d
x
=
∫
D
x
f
(
x
,
α
,
δ
)
f
(
α
,
δ
)
(
α
,
δ
)
d
x
,
where D is the domain of X, f
(
x
,
α
,
δ
) is the joint density function of the three random variables, and f
(
α
,
δ
)
(
α
,
δ
) is the marginal density. All these elements may be unknown so they may be approximated using
(16)
f
^
(
x
,
α
,
δ
)
=
1
n
h
x
h
α
h
sin
δ
∑
i
=
1
n
K
x
(
x
-
x
i
h
x
)
n
n
n
n
n
n
n
n
n
n
n
n
n
n
n
×
K
α
(
α
-
α
i
h
α
)
K
δ
(
sin
δ
-
sin
δ
i
h
sin
δ
)
,
and the condition
(17)
1
4
π
μ
(
D
)
∫
D
∫
S
2
f
^
(
x
,
α
,
δ
)
cos
⁡
δ
d
x
d
α
d
δ
=
1
must be fulfilled. We proceed analogously for the marginal density. Taking the same kernel for all the random variables and considering their properties we reach a similar expression to the one-dimensional Nadaraya-Watson one, but on the sphere:
(18)
m
X
(
α
,
δ
)
=
∑
i
=
1
n
w
i
x
i
,
w
i
=
K
α
(
(
α
-
α
i
)
/
h
α
)
K
δ
(
(
sin
δ
-
sin
δ
i
)
/
h
sin
δ
)
∑
j
=
1
n
K
α
(
(
α
-
α
j
)
/
h
α
)
K
δ
(
(
sin
δ
-
sin
δ
j
)
/
h
sin
δ
)
,
where K is given in (12). It can be demonstrated that the lesser h is the more steeply the adjustment is. A large h provides smoother results. The theoretical optimum values have been studied in [16].Analogously, if we introduce a radial component, we obtain an estimation of the residuals on the solid sphere that generalizes formula (18).


## 3.2. Vector Spherical Developments

Vector spherical harmonics have been widely used in several astronomical lines of work. See in particular [5] within the context of the Hipparcos-FK5 comparison. This issue was also studied by us in [8], where we truncate the VSH development up to first order, which is enough to determine the rotational and deformation parameters. VSH was also studied by us in the context of the Hipparcos spin in [3].Let us consider the vector field over the celestial sphere (in the case of a surface vector field, we take a constant radiusr
=
1 and the last summand disappears):
(19)
Δ
X
≡
V
(
r
,
α
,
δ
)
=
V
α
(
r
,
α
,
δ
)
e
α
+
V
δ
(
r
,
α
,
δ
)
e
δ
+
V
r
(
r
,
α
,
δ
)
e
r
=
Δ
α
cos
⁡
δ
e
a
+
Δ
δ
e
d
+
Δ
r
e
r
,
V
α, V
δ, and V
r being the scalar fields of the vector field and e
α, e
δ the unitary vectors in the tangent plane and in the directions of the right ascension and declination, respectively, and e
r is the unitary vector position.On the other hand, provided that we are on the surface of the unitary sphere, the only vector spherical harmonics involved are theS
l
,
m spheroidal spherical harmonics and the T
l
,
k toroidal spherical harmonics. Anyway, in general, there is an additional radial component denoted by R
l
,
k. Thus, we suppose that the vector fieldV has a mathematical development:
(20)
V
(
r
,
α
,
δ
)
=
∑
k
≥
1
∑
l
=
-
k
k
[
r
l
,
m
R
l
,
m
+
t
l
,
m
T
l
,
m
+
s
l
,
m
S
l
,
m
]
.
The definitions of the vectors of the orthogonal and complete basis for the L
2(S
2) Hilbert space of the vector fields on the sphere are
(21)
R
n
,
m
=
r
r
Y
n
,
m
;
S
n
,
m
=
r
∇
Y
n
,
m
;
T
n
,
m
=
-
r
×
∇
Y
n
,
m
,
where Y
n
,
m
(
n
≥
0
,
-
n
≤
m
≤
n
) are the usual surface spherical harmonics. As we have previously stated, if we are interested only in the spheroidal and toroidal components, then, the estimations may be computed, due to the functional orthogonality, as
(22)
t
^
l
,
m
=
∫
S
2
V
·
T
k
,
l
d
σ
∥
T
k
,
l
∥
2
;
s
^
l
,
m
=
∫
S
2
V
·
S
k
,
l
d
σ
∥
S
k
,
l
∥
2
.
The denominators are exactly calculated, whereas for the numerators an estimation is obtained using the proposed MM: for the calculation of the components of the vector fieldV at regularly spaced points on the sphere we can use the kernel regression method, which is computationally efficient and, in addition, it is rather accurate for the problem that we are discussing. It is important to emphasize that, once the adjustment has been established forV on the set of points of the sphere, this same set can be used for the numerical integration of the numerators up to any order of the development. Thus, we can easily calculate independently the estimations for higher orders of harmonics.Truncation of (20) up to the first order, together with formula (19), requires (see, e.g., [8]) fulfilling the following conditions (23) in order to verify the compatibility of the system:
(23)
Δ
α
cos
⁡
δ
=
[
t
1,0
cos
⁡
δ
-
t
1,1
sin
α
sin
δ
-
t
1
,
-
1
cos
⁡
α
sin
δ
]
+
[
s
1,1
cos
⁡
α
-
s
1
,
-
1
sin
α
]
,
Δ
δ
=
[
-
t
1,1
cos
⁡
α
+
t
1
,
-
1
sin
α
]
+
[
-
s
1,1
sin
α
sin
δ
-
s
1
,
-
1
cos
⁡
α
sin
δ
+
s
1,0
cos
⁡
δ
]
,
where it is evident that t
1
,
-
1
=
ε
x, t
1
,
1
=
ε
y, and t
1
,
0
=
ε
z, ε
x, ε
y, and ε
z being the infinitesimal rotations around the corresponding axis.


## 3.3. Average for the Solid Hipparcos-UCAC4 Residuals

In order to use Hipparcos and UCAC4, we need the data for the common stars that we denote by (r
i, α
i, δ
i), 1
≤
i
≤
N. Notice that we also use the module of the radius vector because good parallax data are available. To reduce this solid case to the surface case, we proceed in the following way: we classify these differences in concentric zones with radii:
(24)
R
°
<
R
¹
<
R
²
<
⋯
<
R
k
.
In this way, we have for each I
(
j
)
=
[
R
(
j
)
,
R
(
j
+
1
)
) interval the (
Δ
α
r
cos
⁡
δ
r
,
Δ
δ
r
)
(
j
) differences, where superindex (j) denotes the slice where the differences are. For each I
(
j
) we have a vector field of differences projected on the R
-
j
,
j
+
1 point which is the arithmetical mean between R
(
j
) and R
(
j
+
1
):
(25)
V
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
=
V
α
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
e
α
+
V
δ
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
e
δ
,
e
α, e
δ being the unitary vectors in the directions of RA and DEC, respectively, that determine the plane tangent to the sphere on each point. Let us call averaged vector field to the field obtained by numerically integrating
(26)
V
(
α
,
δ
)
=
∫
R
-
0,1
R
-
k
-
1
,
k
[
V
(
j
)
(
R
-
j
,
j
+
1
,
α
,
δ
)
]
d
r
.
Once this field is obtained, we proceed as in Section 3.2.


## 4. Numerical Results

The fixed point method, together with the use of VSH and the mixed method, needs only four steps to converge. The results are listed in Table2.Table 2
Results after four steps (in mas.).
ε
x
ε
y
ε
z
s
1,0
s
1,1
s
1
,
-
1
Δ
R
0.25
−0.26
−0.27
0.03
0.22
0.19
Δ
H
1.32
−3.52
−28.72
−0.61
15.69
35.38
Δ
Z
−1.37
−12.92
−25.60
−9.83
14.03
6.28
Δ
I
−0.26
−0.80
−1.18
−1.29
0.14
0.44The ideal Hipparcos-ICRF2 differences are given by the −Δ
R values from Table 2. On the contrary, we could say that currently the relations obtained with the available data are given by Δ
H
+
Δ
Z
-
Δ
I (see Table 3). The ideal relations have been obtained from a set of radio sources common to HCRF and ICRF2. In this case, the number of radio sources employed in the adjustment is low, but the positions are very accurate. On the contrary, the current relations have been obtained using the part of the cyclical process where the most massive catalogs are used. It is interesting to note that the massive catalogs provide more positions but with lower accuracy. These current relations should be obtained again when more accurate and massive catalogs are available.Table 3
Current relations obtained for the Hipparcos-ICRF2 (in mas.).
ε
x
ε
y
ε
z
s
1,0
s
1,1
s
1
,
-
1
Hip-ICRF2
0.21
−15.64
−53.14
−9.15
29.58
41.22The qualitative conclusion that we reach from these results is that new reductions of the Hipparcos positions must be carried out on the basis of modern observations such as the Sloan Digital Sky Survey (SDSS) Project (seehttp://www.sdss.org/) or the use of techniques similar to the employed in [19] to compile Zach14.


## 5. Conclusions

The aim of this paper is the study of the current situation of the link between HCRF and ICRF2, using the four couples of catalogs (Rio-ICRF2)-(Rio-Hipparcos)-(Hipparcos-UCAC4) -(ICRF2-Zach14) in a closed fixed point process with departure point in the couple Rio-ICRF2.In each of the steps, the relations are given by infinitesimal rotations and deformations between the references that each set of data defines. We have also pointed out the problems that appear in each step and in the proposed method. All of the problems come from the fact that, although we search the parameters of a certain functional development, we work with vector fields whose components belong to theL
2
(
S
2
) Hilbert space, but in practice for each couple of references we have only a discrete set of unnecessary homogeneously distributed points on the sphere. That has led us to employ the MM instead of the current DLS and a kernel nonparametrical regression.The fixed point method, together with the use of VSH and the MM, has a very quick convergence and it has provided us with the ideal Hipparcos-ICRF2 differences, but with the data available nowadays the relations obtained are very much different. We conclude from this that new reductions of the Hipparcos positions must be carried out on the basis of modern observations or the employ of techniques similar to the employed in [12] to compile Zach14.



---

*Source: 917583-2014-10-14.xml*