Li Liu

Da Chen

Laurent D Cohen

Jiasong Wu

Michel Paques

Huazhong Shu

Anisotropic Tubular Minimal Path Model with Fast Marching Front Freezing Scheme

Keywords: Minimal path model, Anisotropy enhancement, Riemannian metric, Path feature, Tubular structures

In this work, we introduce an anisotropic minimal path model based on a new

Riemannian tensor integrating the crossing-adaptive anisotropic radius-lifted tensor field and the front freezing indicator by appearance and path features.

The non-local path feature only can be obtained during the geodesic distance computation process by the fast marching method. The predefined criterion derived from path feature is able to steer the front evolution by freezing the point causing high bending of the geodesic to solve the shortcut problem. We performed qualitative and quantitative experiments on synthetic and real images (including retinal vessels, rivers and roads) and compare with the minimal path models with classical anisotropic Riemannian metric and dynamic isotropic metric, which demonstrated the proposed method can detect desired targets from complex tubular tree structures.

Introduction

Tubular structure extraction is a crucial task in many computer vision applications. Various tubular structure extraction algorithms and techniques have been studied, which can be roughly divided into two categories: supervised and unsupervised methods [START_REF] Kirbas | A review of vessel extraction techniques and algorithms[END_REF][START_REF] Zhang | Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores[END_REF]. Supervised approaches usually achieve better performance due to the prior information from the training components to decide whether a pixel should belong to a vessel or not. These models construct feature vectors in order to train the respective classifiers [START_REF] Ricci | Retinal blood vessel segmentation using line operators and support vector classification[END_REF][START_REF] Soares | Retinal vessel segmentation using the 2-d gabor wavelet and supervised classification[END_REF] or apply the neural network [START_REF] Sironi | Multiscale centerline detection[END_REF][START_REF] Zhu | Retinal vessel segmentation in colour fundus images using extreme learning machine[END_REF] with the help of manual label. However, it is very time-consuming to find discrimination features for constructing powerful models by training. Besides, the obtained models always meet for particular requirement. Those unsupervised methods, without manual labeling information, usually have lower computation complexity and more simple procedure [START_REF] Zhang | Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores[END_REF], such as the pattern recognition techniques [START_REF] Armande | Thin nets extraction using a multi-scale approach[END_REF][START_REF] Hoover | Locating blood vessels in retinal images by piecewise threshold probing of a matched filter response[END_REF], model-based approaches [START_REF] Lam | General retinal vessel segmentation using regularization-based multiconcavity modeling[END_REF][START_REF] Vermeer | A model based method for retinal blood vessel detection[END_REF], tracking-based algorithms [START_REF] Bekkers | A multiorientation analysis approach to retinal vessel tracking[END_REF][START_REF] Tolias | A fuzzy vessel tracking algorithm for retinal images based on fuzzy clustering[END_REF] and so on.

The minimal path technique [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF] has been adopted to detect tubular structures thanks to its efficiency and global optimality [START_REF] Deschamps | Fast extraction of minimal paths in 3d images and applications to virtual endoscopy[END_REF][START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF][START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF][START_REF] Péchaud | Extraction of tubular structures over an orientation domain[END_REF]. It is often taken as an interactive tool such that users are allowed to provide manual intervention such as source points, which constrains the minimal path to delineate the tubular structures. In essence, a geodesic path is a curve linking two points derived by globally minimizing a weighted curve length [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]. It is defined by integrating a geodesic metric P : Ω × R d → R + , where Ω ⊂ R d stands for the open and bounded domain and d is the dimension of Ω, along a regular path γ:

L(γ) = 1 0 P(γ(t), γ (t))dt, (1)
where γ : [0, 1] → Ω ⊂ R d is a regular curve of Lipschitz continuity. The minimal path models have been extensively exploited for tubular structure segmentation, where both the centerline positions and thickness of the tubular structures are required [START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF][START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF][START_REF] Péchaud | Extraction of tubular structures over an orientation domain[END_REF]. The isotropic tubular minimal path model [START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF] only makes use of the curve position information, which may increase the risk of shortcuts.

More general anisotropic Riemannian metric [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF][START_REF] Chen | Vessel extraction using anisotropic minimal paths and path score[END_REF] is proposed to segment tubular structures in conjunction of the pre-detected vessel anisotropy features.

However, this extended tubular minimal path model still has difficulties when dealing with complicated situation such as extract a weak vessel from a tree structure. Moreover, some researchers try to improve the classical tubular minimal path models by reducing the user intervention [START_REF] Chen | Vessel tree extraction using radiuslifted keypoints searching scheme and anisotropic fast marching method[END_REF][START_REF] Rouchdy | Geodesic voting for the automatic extraction of tree structures. methods and applications[END_REF][START_REF] Kaul | Detecting curves with unknown endpoints and arbitrary topology using minimal paths[END_REF]. For these models, only one (or several) initial root point(s) is (are) necessary to extract a tubular tree structure. Recently, the higher order properties of the geodesic curves are exploited for shortest path estimation as well as tubular structure segmentation. Specifically, Ulen et al. [START_REF] Ulen | Shortest paths with higher-order regularization[END_REF] proposed a curvature and torsion regularized shortest path model by computing the curvature and the torsion.

Based on the Eikonal equation framework, the curvature-penalized minimal path models [START_REF] Chen | Global minimum for a finsler elastica minimal path approach[END_REF][START_REF] Bekkers | Nilpotent approximations of sub-Riemannian distances for fast perceptual grouping of blood vessels in 2D and 3D[END_REF][START_REF] Duits | Optimal paths for variants of the 2D and 3D Reeds-Shepp car with applications in image analysis[END_REF] are established over an orientation-lifted space such that the orientation dimension can be used to represent the curve tangents.

Controlling the course of the fast marching front propagation is a solution for the above problem. Criterions are defined to make the front propagation stop in undesired directions or promote in specific directions [START_REF] Lo | Vessel tree extraction using locally optimal paths[END_REF][START_REF] Liao | Automatic 3d segmentation and quantification of lenticulostriate arteries from high-resolution 7 tesla mra images[END_REF]. In contrast to the minimal path models using static geodesic metrics which are fixed during fast marching front propagation, the dynamic minimal path models allow to update the respective geodesic metrics in the course of geodesic distance estimation.

In [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF], the authors proposed a new fast marching front propagation scheme in order to address the shortcut and short branches combination problems. The basic idea is to freeze the front points of which the respective local path features violate the pre-defined criteria. In this case the fast marching fronts will not pass through the points corresponding to undesired geodesic paths. Chen et al. [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF] proposed an adaptively Riemannian metric using a dynamic metric construction way in order to penalize the appearance feature coherence property. In this case the obtained geodesic paths favour to pass by the region of slow-varying tubular features. Geoffrey et al. [START_REF] Parker | Estimating distributed anatomical connectivity using fast marching methods and diffusion tensor imaging[END_REF] created a connectivity metric to determine the reasonable pathways of connection. Windheuser et al. [START_REF] Windheuser | Beyond connecting the dots: A polynomial-time algorithm for segmentation and boundary estimation with imprecise user input[END_REF] estimated the curvature information by computing the angles between each pair of adjacent edges. Krueger et al. [START_REF] Krueger | Robust and efficient object segmentation using pseudo-elastica[END_REF] incorporated the curvature estimation into a greedy core algorithm to improve the efficiency. However, these models track geodesic paths only in the spatial space where the vessel thickness information cannot be extracted.

In this paper, we focus on the dynamic minimal path-based tubular structure segmentation, for which the geodesic metrics are established over a radius-lifted space. The existing tubular minimal path models [START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF][START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF][START_REF] Péchaud | Extraction of tubular structures over an orientation domain[END_REF] often suffer from the shortcut and short branches combination problems, especially for the case that the target is weakly-defined and crosses or parallels with strong ones, see Fig. 1. In this figure, the target vessel has low gray level contrast, while its neighbour vessels are strongly-defined (i.e., high gray level contrast). In Fig. 1(c), we show the minimal path extraction result obtained by the anisotropic tubular minimal path model [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF]. One can observe short branches combination occur.

In other words, the obtained minimal path prefers to pass through strong segments belonging to different vessels. In order to simultaneously obtain the centerlines and thickness of tubular structures as well as to overcome the shortcut and short branches combination problems, we proposed an anisotropic minimal path model for tubular structure segmentation via a dynamic radius-lifted anisotropic metric and state-of-the-art anisotropic fast marching method [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using lattice basis reduction[END_REF].

The motivation of this work is that the anisotropic tubular minimal path model [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF] cannot take into account path features for distance estimation. Using the path features detected during the front propagation for minimal path computation is able to obtain desired results even in complicated situation. The main contribution in this paper lies at the introduction of a dynamic anisotropic Riemannian metric in the radius-lifted space using back-tracked geodesic paths. In Riemannian metric in a radius-lifted space and its numerical implementation in conjunction with the fast marching method. Section 4 shows the experimental results and the final section 5 gives the conclusion.

Minimal Paths for Tubular Structure Segmentation

In this section, we give the background on anisotropic radius-lifted minimal path model [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF] for tubular structure segmentation.

Tubular Geometry Features Computation

In this section, we first consider the extraction of the local tubular geometry features from the image I : Ω → R through the optimally oriented flux (OOF)

filter [START_REF] Law | Three dimensional curvilinear structure detection using optimally oriented flux[END_REF]. Note that the local geometry features can also be extracted using other steerable filters such as [START_REF] Law | Three dimensional curvilinear structure detection using optimally oriented flux[END_REF][START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF][START_REF] Hannink | Crossing-preserving multi-scale vesselness[END_REF]. The output of the OOF filter is established in a multi-scale space defined as Ω := Ω × [r min , r max ], where [r min , r max ] is the radius space1 . The output of the OOF filter at the position x and scale r is a symmetric matrix of size d × d as follows:

F(x, r) = 1 r (∂ i,j G σ ) * I * χ Sr (x, r), (2)
where (∂ i,j G σ ) is the Hessian matrix of the Gaussian kernel G σ with variance σ. χ Sr is a characteristic function of a disk with radius r.

The eigenvalues λ i (i = 1, 2, • • • , d) extracted using the OOF filter are the values of the oriented flux along the corresponding eigenvectors v i as following:

λ i (x, r) = v T i (x, r) F(x, r)v i (x, r). (3)
Without loss of generality, we assume that

λ 1 (•) ≤ λ 2 (•) ≤ • • • λ d (•).
In this paper, we only consider the 2D tubular segmentation problem, which means that d = 2. For the case that the tubular structures have lower gray levels than background, the vector v 1 (x, r * ) indicates the orientation of the tubular structure at the point x inside the vessel structures, where r * is the optimal vessel scale of point x.

Radius-lifted Minimal Paths for Tubular Segmentation

We denote by S + d the set of symmetric positive definite matrices of size d × d (d = 2, 3) and let Lip([0, 1], Ω) be the set of Lipchitz continuous curves γ : [0, 1] → Ω. The radius-lifted minimal path model is proposed in [START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF] for tubularity segmentation. Benmannsour and Cohen [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF] generalize this isotropic model [START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF] to the anisotropic case. The basic idea for both models is to lift a spatial geodesic to the radius-lifted space Ω by adding one abstract dimension to represent the vessel radii.

The path length L for such a radius-lifted curve γ ∈ Lip([0, 1], Ω) can be measured through a Riemannian tensor field M scale : Ω → S + 3 ,

L(γ) = 1 0 γ (u), M scale (γ(u))γ (u) du, (4)
where •, • represents the Euclidean scalar product. The tensor field M scale is constructed [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF] by the eigenvalues λ i and eigenvectors v i ,

M scale (x, r) =   M aniso (x, r) 0 0 P scale (x, r)   , (5)
where M aniso : Ω → S + 2 is a tensor field associated to the spatial anisotropy with size 2×2 and P scale : Ω → R + is a scalar function. Both can be constructed by combining the eigenvalues and eigenvectors of the OOF response F in Eq. [START_REF] Zhang | Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores[END_REF] as follows:

M aniso (x) = exp(αλ 1 (x))v 2 (x)v T 2 (x) + exp(αλ 2 (x))v 1 (x)v T 1 (x), (6)
P scale (x) = β exp(αλ 2 (x)), (7)
where α ∈ R and β ∈ R + are two parameters which control the regularization 105 of the spatial and radius dimensions, respectively. Note that the value of α is negative for dark-vessel and bright-background case and α > 0, otherwise.

In the Riemannian minimal path model, once the metric tensor is given, one can define the relevant geodesic distance U ŝ(x) as the minimal curve length between x and the fixed source point ŝ as follows

U ŝ(x) = inf γ∈Lip([0,1], Ω) {L(γ); γ(0) = ŝ, γ(1) = x}, (8)
where the values of U ŝ can be regarded as the arrival times of a front propagating with oriented velocity related to the metric tensor M -1 scale . The geodesic distance map U ŝ from the source point ŝ satisfies the anisotropic Riemannian Eikonal equation

∇U ŝ(x) M -1 scale (x) = 1, ∀x ∈ Ω\{ŝ}, (9)
with boundary condition U ŝ(ŝ) = 0, where u M = √ u T M u for any matrix M ∈ S + 3 . A geodesic path Ĉx,ŝ linking x to the source point ŝ can be tracked from the point x by solving the following ordinary differential equation (ODE) till the source point ŝ is reached

C x,ŝ (s) = - M -1 scale ( Cx,ŝ (s))∇U ŝ( Cx,ŝ (s)) M -1 scale ( Cx,ŝ (s))∇U ŝ( Cx,ŝ (s)) (10)
with Cx,ŝ (0) = x. The geodesic Cx,ŝ is parameterized by its arc-length. The 110 final geodesic C ŝ,x ∈ Lip([0, 1], Ω) joining from the source point ŝ to x can be obtained by reversing and reparameterizing Cx,ŝ .

The New Anisotropic Geodesic Metric with Nonlocal Information

in Radius-lifted Space

Overview. The main goal in this section is to establish a new anisotropic geodesic metric in a radius-lifted space in order to avoid the shortcut and short branches combination problems as shown in Fig. 1. In contrast to [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF] for which only the local geometry information are considered to establish the geodesic metrics, the front freezing-based minimal path model [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF] constructs the geodesic metrics by taking into account the nonlocal information such as the local path.

During the geodesic distance front propagation, the front points for which the features do not satisfy the given criterion will be frozen, where the features for each front point are extracted through two extra points. This scheme proposed in [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF] is able to avoid the shortcut problem in some extent. However, it is restricted to the case of centerline detection and cannot take advantages of the path orientation due to the isotropic nature of the used metrics in [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF]. In this section, we extend this front frozen scheme to a radius-lifted space through an anisotropic Riemannian metric. This is done by constructing a new Riemannian tensor field M dyn : Ω → S + 3 during the geodesic distance computation which is actually carried out in a front advancing procedure. The invoked tensor field M dyn consists of two ingredients: the crossing-adaptive anisotropic radius-lifted tensor field M adap : Ω → S + 3 and the front freezing indicator δ : Ω → {1, ∞}. It can be expressed by

M dyn (x, r) = M adap (x, r) × δ(x, r), ∀(x, r) ∈ Ω. ( 11
)
We will explain the construction details for the tensor field M adap in Section 3.1

115
and the scalar indicator δ in Section 3.2, respectively.

Computation of the Crossing-Adaptive Tensor Field

By applying the OOF filter to an image involving vessel tree structures, one can obtain an optimal direction v 1 (x, r) at each point (x, r) belonging to the tree structures through eigen-decomposition to the response of the filter OOF

F(x, r) = 2 i=1 λ i (x, r)v i (x, r), s. t. λ 1 (•) ≤ λ 2 (•) (12)
We suppose that inside the tubular structures the intensities are lower than background, thus the corresponding eigenvalue λ 2 can be used to characterize the appearance feature of the vessels, by which we define a vesselness map

ζ : Ω → R + 0 ζ(x) = max max r λ 2 (x, r), 0 . (13)
The value of ζ(x) is derived from λ 2 at the optimal scale η(x) ∈ [r min , r max ],

where the η : Ω → [r min , r max ] is a map defined by

η(x) = arg max r∈[rmin,rmax] {λ 2 (x, r)} . (14)
Based on the optimal scale map η and the eigenvectors v 1 , we can define the tubular feature vector p : Ω → R 2 that indicates the orientation for a tubular structure

p(x) = v 1 (x, η(x)). (15)
The tubular direction p(x) well describes the vessel orientation at non-crossing point x. However, for crossing points, denoted by y, there will be at least two vessels across one another. This means that the highly anisotropic tensors 120 M aniso (y, •) defined in Eq. ( 5) may lead to incorrect geodesic tangents at these which can be done by invoking the orientation scores [START_REF] Zhang | Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores[END_REF] or by structure tensors [START_REF] Law | Gradient competition anisotropy for centerline extraction and segmentation of spinal cords[END_REF]. In this section, we construct our crossing-adaptive radius-lifted tensor field through the tool of structure tensor.

The shortcut and short branches combination problems usually occur where the vessel crosses its neighbours, since the eigenvector of a crossing point usually indicates the orientation of the vessel with strong appearance feature. So the speed computed from the anisotropic geodesic metric is slower along the weak vessel direction than along the strong one. To solve the shortcut and short branches combination problems, we make use of a crossing-adaptive tensor field M adap : Ω → S + 3 established in the radius-lifted space Ω, the anisotropy of which is kept at single-vessel region and removed or reduced at crossing points.

It can be constructed by

M adap (x, r) = diag exp(-αλ 2 (x, r))T smooth (x, r), P scale (x, r) , ( 16
)
where α is a positive constant which has been used in Eqs. ( 6) and [START_REF] Armande | Thin nets extraction using a multi-scale approach[END_REF]. The tensor field T smooth : Ω → S + 2 is computed via a Gaussian kernel G p with standard derivative p and the identity matrix

I d of size 2 × 2 T smooth (x, •) = (G p * T )(x) (G p * )(x) + I d -1 , T (x) = (x)p(x)p T (x), (17)
where ∈ R + is a sufficiently small constant to generate a regular matrix M adap (•). Note that for any radius values r 1 , r 2 ∈ [r min , r max ] where r 1 = r 2 , one always has T smooth (x, r 1 ) = T smooth (x, r 2 ). The scalar-valued function

: Ω → R + 0 is a weighted function over the image domain Ω. In this paper, we propose three ways to compute the function , where the first way is to invoke the vesselness map ζ (see Eq. ( 13)).

(x) = ζ(x), ∀x ∈ Ω. (18)
Alternatively, the weighted function can be computed as a binary-valued function, which depends on a threshold value Th ∈ R + and also on the vesselness

map ζ (x) =      1, if ζ(x) > Th 0, otherwise. (19)
The last method utilizes the skeleton of the vessel mask. The function is used in order to reduce the influence from the regions outside the vessel structures.

In other words, the information inside the vessel region dominates the tensor field T smooth .

In Eq. ( 18), the function is set to be the vesselness map which has small values outside the vessel regions but large values inside, where an example can 135 be found in the first row of column (a) of Fig. 3. This is an adaptive way for the computation of the tensor field T smooth comparing to the other two ways which set as a binary-valued function. Considering the vessel segmentation has been studied for a long time, the vessel mask is easy to be computed by many existing approaches and the skeleton can be obtained from the vessel mask. For the second method expressed in Eq. ( 19), a point x for which (x) = 1 implies that x is located inside a vessel. An example is shown in the second row of column (a) in Fig. 3. In the third method, (x) = 1 means that x is a centreline point, see the third row of column (a) of Fig. 3 for an example. In this case, the information on the centreline positions dominates the tensor field T smooth . For the crossing point, the tensor field T smooth is impacted greater by the vessel with bigger radius. Note that the main difference between T smooth and the structure tensor field used in [START_REF] Law | Gradient competition anisotropy for centerline extraction and segmentation of spinal cords[END_REF] lies at the existence of the function . As discussed above, it can reduce or avoid the effects derived from the non-vessel regions.

This means that we are able to make use of a Gaussian filter with more flexible variance values p.

We express the tensor filed T -1 smooth by

T -1 smooth (x, •) = Ω G p (y) (x -y)p(x -y)p T (x -y)dy G p (y) (x -y)dy + I d . (20)
One can see that the Gaussian filtering is operated on the feature vectors p(y)

(indicating the vessel orientations) at the vicinity of a point x. If x is a vessel segment point and the vectors p(y) change slowly, then the eigenvector of the tensor (G p * T )(x), which corresponds to the largest eigenvalue, will approximate to the feature vector p(x). For the vessel points x nearby a crossing structure, the feature vectors potentially vary fast, leading to the tensors (G p * T )(x) which are nearly isotropic. As shown in Fig. 3, the anisotropy of the metric at point (e.g. the red point) located on the crossing section is removed apparently by the Gaussian filtering, however, the anisotropy of the metric of the single vessel point (e.g. cyan point) is contained adaptively.

Computation of the Front-Freezing Indicator

In this section, we construct the proposed front-freezing indicator δ with non-local path features derived from the local truncated geodesic paths to solve the shortcut and short branches combination problems. The non-local feature is described by the bending measure of the local path computed through two extra points located on it. Moreover, we first present the principle for the construction of the indicator δ and the numerical implementation method based on the fast marching method will be described in Section 3.3.

Here we give the method to identify the extra two points on a local geodesic path linking a target point x = (x, r x ) to the source point ŝ = (s, r s ). First of all, the geodesic path can be computed by solving the gradient descent ODE (see in Eq. ( 10)) on the obtained geodesic distance map U ŝ, which has already been computed within the region where the front visited [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF]. During the front propagation, one can easily track a radius-lifted geodesic path Cx,ŝ parameterized by its arc-length by solving the corresponding gradient descent ODE [START_REF] Vermeer | A model based method for retinal blood vessel detection[END_REF]. Since Cx,ŝ lies in the radius-lifted space, we have Cx,ŝ = ( Gx,s , η) where Gx,s (•) ∈ Ω is the projected path and η(•) ∈ [r min , r max ]. We first detect a point q = (q, r q ) from Cx,ŝ such that its projected path |G x,q | = Γ, where |G x,q | stands for the Euclidean length of G x,q and Γ is a positive constant. Note that q is the first extra point and the second one, referred to as m = (m, r m ), is defined as the middle point of the part of the path Cx,q . The process to search the two extra points are shown in Fig. 4(a). The bending measure K is estimated by the angle between two vectors related to the two extra points q, m ∈ Ω as follows [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF]

K(x, r) = m -q, x -m m -q x -m , (21)
where •, • denotes the scalar product and ||•|| represents the norm of the vector.

170

The values of K are ranged at [-1, 1]. The estimation of K is done during the fast marching front propagation, which will be described in Section 3.3. We mention that a small value of K implies that the corresponding geodesic path likely appears with a sharp turning.

Algorithm 1 Fast Marching Method with Dynamic Riemannian metric Output: Minimal action map U ŝ.

Initialization: ∀x = (x, r) ∈ Ω \ {ŝ}, set U ŝ(x) ← ∞, L(x) ← F ar; U ŝ(ŝ) ← 0, L(ŝ) ← T rial ; F(x) ← 0.
1: while stopping criterion is not reached do 2:

Find xmin = (x min , r min ), the T rial point which minimizes U ŝ;

3:

Set L(x min ) ← Accepted;

4:

if F(x min ) = 0 then 5:

Backtrack the local geodesic path from xmin ;

6:

Compute the bending measure K(x min , r min ) via Eq. ( 21);

7:

if K(x min , r min ) ≤ K 0 then 8:

Set U ŝ(x min ) ← ∞; Front point freezing procedure.

9:

Set F(x min ) ← 1;

10:

else Regular fast marching scheme.

11:

for all ŷ ∈ N * (x min ) and L(ŷ) = Accepted do 12:

Set L(ŷ) ← T rial and update U ŝ(ŷ);

13:

end for Set U ŝ(x min ) ← ∞; Front point freezing procedure.

δ(x) =      1 K(x) > K 0 , +∞ otherwise (22)
where K 0 is a given threshold. The bending measure K(x) > K 0 which yields δ = 1 means the fast marching front is propagated as usual. By using the indicator δ, the shortcut and short branches combination problems can be avoided as shown in Fig. 1. Our method with an appropriated threshold is able to detect the desired vessel although the neighbouring vessel has stronger appearance feature.

Fast Marching Implementation

The fast marching method estimates the geodesic distance map on a discretization orthogonal grid of the image domain Ω. The single-pass scheme of the fast marching method makes it possible to estimate the indicator δ(x) in the course of the geodesic distance computation for each point x ∈ Ω.

In the course of the fast marching front propagation, all the grid points are labeled as three classes: Far, Trial and Accepted. The geodesic distance will be estimated by solving the Hopf-Lax operator defined in [START_REF] Mirebeau | Anisotropic fast-marching on cartesian grids using lattice basis reduction[END_REF]. The stencil N (x) determines the neighbourhood for each point x, and its inverse neighbourhood using the projected path Gxmin,q . In Figs. If F(x min ) = 0, the local geodesic path Cx min, ŝ is tracked and the indicator δ(x min ) is computed. Then the tensor M dyn (x min ) can be constructed as in Eq. ( 11) and the geodesic distance value can be estimated. According to the front-freezing indicator δ(x min ), if K ≤ K 0 , the value of F(x min ) is set to 1 which means the front propagation needs to be frozen at xmin . The front propagation scheme can be terminated till the end point has been tagged as Accepted or all the grid points have been tagged as Accepted. The pseudo-code for our fast marching-based implementation procedure can be found in Algorithm 1.

N * (x)
The fast marching front propagation process during the computation of the geodesic distance map is displayed in Fig. 5. The geodesic distance map obtained in the radius-lifted space is projected to the domain Ω by minimizing it along radius dimension as follows:

H s (x) = min r∈[rmin,rmax] U ŝ(x, r) (23)
The obtained geodesic map H s (x) is shown on the original vessel images as in

Discussion

The proposed method in essence constructs the geodesic metric in a dynamic way and we discuss the differences between our method and two state-of-the-art dynamic minimal path approaches [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF][START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF].

• In work [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF], authors utilize the dynamic isotropic metric established only in the image domain. They introduced a new front freezing fast marching scheme based on the path features which can be detected through the appearance feature of the path or the geometric feature such as the bending measure of the paths. However, This dynamic minimal path approach only invokes isotropic geodesic metrics established in the spatial domain (the image domain), so it cannot benefit from the vessel anisotropy information and also can only detect the vessel centerlines. In order to address the issues mentioned above, we propose a dynamic anisotropic and crossing-adaptive Riemannian metric a in radius-lifted domain based on the geometric features of the local path. As a consequence, the proposed model is able to take advantages of the vessel anisotropy information and also to detect the vessel centerlines and boundaries simultaneously.

• In the work [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF], the authors proposed a minimal path model associated with an adaptively Riemannian metric embedded with an appearance feature coherence penalty. The key differences between the model [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF] and the proposed model are mainly described as three aspects. Firstly, the method [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF] penalizes the geodesic metric by a weighted scalar-valued function derived from the coherence property of the appearance features, instead of freezing the front points. Sometimes one needs to carefully control the importance between the tubular confidence information and the coherence penalization. Secondly, the dynamic metric in [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF] is still built in the spatial domain3 while our method can detect both centerlines and thickness in one step. Finally, we analyze the implementation difference on the adaptive establishment. In [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF], the authors build the dynamic metric using a reference point form a truncated geodesic path vector to distinguish the correct vessel direction in crossing points. In contrast, the basic idea in the the proposed work for the considered crossing-adaptive anisotropy is to reduce or remove the anisotropy property at the crossing points and to keep the anisotropy property at the single vessel points, done by the tool of structure tensor.

Experimental Results and Discussion

In this section, we conduct the numerical experiments on qualitative and quantitive comparison between state-of-the-art minimal path models and the Basic Setting. The front-freezing indicator δ in Eq. ( 22) is computed based on the bending measure operator in Eq. ( 21) of the local truncated geodesic. In numerical implementation, the local geodesic length Γ and the bending threshold K 0 should be given. In this work, we conduct an experiment to study the sensitivity of the segmented results to the two crucial parameters Γ and K 0 on synthetic images, the dataset of which will be described detailed in Section 4.1. The segmentation accuracy and running time are shown in Fig. 6. We set Γ = 8 and K 0 = 0.9 which are reasonable choice for the experiment, since the segmentation accuracy is high and running time is relative small. Besides, the parameter = 0.05 in Eq. ( 17) is used to determine the anisotropic property of the structure tensor. The parameter α = 5 in Eq. ( 16) is related to the influence of the appearance features. The parameter β = 0.5 controls the radius speed. The experiments are performed on a standard Intel Core i5 of 2.7GHz architecture with 8Gb RAM.

For the proposed dynamic minimal path model, we consider three types of the scalar functions to construct the crossing-adaptive anisotropic radius-lifted Riemannian metrics, respectively relying on the vesselness map (dArR), vessel mask (dArR-M) and vessel skeleton (dArR-S), see Section 3.1. In order to show the advantages of the proposed model in tubular structure extraction, we compare with the ArR metric [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF], the variant of the isotropic dynamic minimal path model (dIsoM) [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF] and the appearance feature coherence-penalized Riemannian metric with adaptively anisotropy (AFC) [START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF]. Note that for the sake of fair comparison, the dIsoM metric-driven minimal paths are extracted by submitting the scalar-valued function P scale defined in Eq. ( 7) to the geodesic distance propagation scheme with front freezing procedure as depicted in Algorithm 1.

Evaluation Score. In order to evaluate the proposed dynamic minimal path models quantitatively, we consider an accuracy measurement to validate the segmentation results. Let S represent the segmentation region from the considered models and G denotes the region from ground truth data. In addition, # • be the number of grid points within the set. Thus the measurement is defined as

R = # S ∩ G # S . (24)
One can point out that the measurement R is ranged within the interval [0, 1],

where R = 1 means that the vessel segmentation is exactly identical to the ground truth. When the measurement R is used in the vessel centerlines comparison, the ground truth region G is the dilated centerline from the ground truth by a disk of radius 2 (in grid points) and S represents the set of grid points located in the detected centerlines.

Experiments on Synthetic Images

In this section, we compare the proposed dArR metric to the ArR metric [START_REF] Benmansour | Tubular structure segmentation based on minimal path method and anisotropic enhancement[END_REF]. The crossing-adaptive tensor field for this dArR metric is computed by the associated vesselness map. The experiments are conducted on synthetic images in order to study the robustness of the dArR metric. These synthetic images include high curvature lines and crossing lines. All images are blurred by Gaussian noise with 5 different standard deviation levels ranging from 0.01 to 0.09 to evaluate the robustness against noise. Each level of the noise is added independently 3 times. The ground truth is a binary segmentation of all pixels corresponding to the structures.

The segmentation results in terms of the accuracy measurement R are shown in Table . 1. The extracted centerlines and boundaries are demonstrated in Fig. 7. From both the qualitative results in Fig. 7 and the quantitative results from the image data. Thus in case the target vessels are smooth (without sharp turning), the dArR metric can obtain desired minimal paths.

Experiments on Retinal Images

In this section, we validate the proposed minimal path model on 40 retinal patches obtained from the DRIVE dataset [START_REF] Staal | Ridge-based vessel segmentation in color images of the retina[END_REF] and 40 retinal patches from the IOSTAR dataset [START_REF] Zhang | Robust retinal vessel segmentation via locally adaptive derivative frames in orientation scores[END_REF]. Each patch includes a retinal artery vessel which is near a stronger retinal vein vessel or crosses with it. The average size of the patches from the DRIVE and the IOSTAR is 160 × 277 and 223 × 326, respectively. Our goal is to extract an artery vessel between two given points. In order to show the advantages of using the crossing-adaptive tensor field, we directly apply the radius-lifted tensor field M adap defined in Eq. ( 16) for minimal path extraction under a regular anisotropic fast marching algorithm. For simplicity, we refer to M adap as caArR in the remaining of this paper.

In Fig. 9, we illustrate the effect of the caArR metric for which the anisotropy is reduced at the crossing points. In column (a), the yellow lines indicate the anisotropy feature vectors of the tubular structures detected by the OOF filter.

We can observe that the anisotropy feature vectors located at the cross section are approximately collinear to the directions of the vessel of strong appearance feature. This may lead to incorrect geodesic tangents on these crossing points results from the long vessel crossing or near with strong ones. The third row shows the detect results from the image with long and tortuosity vessels. The last two rows denote detected geodesics of vessel regions with strong tortuosity.

We observe that the classical ArR metric suffered from short combination problem. In the last two columns, the proposed dArR metric can detect the desired vessel regions by the front frozen scheme.

We quantitatively compare the proposed dynamic models including the dArR, dArR-M and dArR-S metrics, to the existing ArR, dIsoM and AFC metrics on those patches of retina images. The quantitative evaluation is computed with racy results for the vessel regions with respect to different metrics are shown in Table . 2. The detection accuracy of the vessel centerlines is shown in Table . 3.

(a) (b) (c) (d)
From these quantitive results, we can see that the proposed metrics indeed overcomes the shortcuts and short branches combination problems and outperform 380 all the other compared metrics.

The dIsoM 4 and dArR metrics (also including the dArR-M and dArR-S metrics) are able to take into account the curve bending measure to freeze front points so as to avoid shortcuts and short branches combination problems. From Tables. 2 and 3, we can see that these metrics achieve better vessel segmenta-385 4 As mentioned above, dIsoM metric is a variant of the one proposed in [START_REF] Liao | Progressive minimal path method for segmentation of 2d and 3d line structures[END_REF]. tion results than the static ones on those retinal image patches. However, the proposed metrics differ to dIsoM by encompassing the vessel anisotropy features to build an anisotropic (radius-lifted) Riemannian metric. During the fast marching front propagation and within the vessel regions, the vessel anisotropy features can greatly help the fast marching fronts to propagate faster along the direction of the vessel. As a consequence, in conjunction with the crossingadaptive tensor filed, the performance of the proposed method in the vessel tracking tasks can be improved by the anisotropy enhancement.

In Fig. 12, we show the minimal paths obtained respectively using the AFC metric and the dArR metric. The experimental results are conducted on retina image patches, showing that our proposed dArR metric can detect the correct vessels even if its neighboring or crossing vessels have similar appearance features.

Finally, in Fig. 13, a failure example of vessel extraction for the proposed metric is given, in which parts of the vessel region are mistaken and the vessel centerline is not correct. The vessel directions detected by the OOF do not

Experiments on Satellite Images

In this section, we conduct the experiments on satellite images with rivers, bridges and streets from Google Earth to evaluate the proposed metric.

Conclusion

In this paper, there are mainly three contributions of the proposed method to overcome the shortcut and short branches combination problems in tubular structure detection. Firstly, we propose a crossing-adaptive anisotropic tensor field in the radius-lifted space, which is used to reduce the in favourable effects from the incorrect vessel anisotropy features at the crossing points. This is done by reducing the anisotropy in the crossing points adaptively and simultaneously keeping the high anisotropy in the single vessel points. Secondly, we propose a new way to reduce the influence from the region outside the vessel structure.

Finally, we show a way for exploiting front freezing scheme for the fast marching algorithm in a radius-lifted space. We show the practical application of the proposed dynamic minimal path model on synthetic images, retinal images and satellite images. The results demonstrate that our method indeed outperforms classical minimal path methods with static geodesic metrics and the dynamic isotropic metric. In the future work, the proposed method will be extended to 3D tubular structure detection application and try to reduce the manually intervention.

Figure 1 :

 1 Figure 1: (a)The yellow and cyan dots incidate the source and end points. (b) The red dash line indicates the target tubular structure centerline. (c) The geodesic path (blue line) obtained by the anisotropic tubular minimal path model [16], where the short branches combination problem occurs. (d) The geodesic path (blue line) derived from the proposed model. (Yellow lines denote the vessel boundaries)

Fig. 1 (

 1 Fig. 1(d), we present the result obtained by the proposed model. One can point out that the blue line indeed correctly indicates the target. The paper outline is shown as follows. In Section 2, we briefly introduce the background for the anisotropic radius-lifted Riemannian metric. Section 3 describes the main contribution of this work including the proposed dynamic

Figure 2 :

 2 Figure 2: (a) A retina image with a crossing point denoted as cyan dot. (b) The red curves illustrate the vessel direction passing through the crossing point. Yellow curve denotes the vessel direction detected by the OOF filter.

Figure 3 :

 3 Figure 3: (a) A crossing structure with two sampled points. Columns (b) and (c): The unit balls for the proposed M adap with p = 1 and p = 7 at the single vessel point denoted as cyan dot. Columns (d) and (e): The unit balls for the proposed M adap with p = 1 and p = 7 at the crossing point displayed as a red dot.

Figure 4 :

 4 Figure 4: (a) The extra points q and m identified from the projected geodesic path Gx,s. (b) and (c) Projected path Gx min ,q (red lines) with small and large bending measures, respectively. The red dot in figure (b) (or figure (c)) indicates the spatial component x min of the latest T rial point xmin = (x min , r min ).

  end while Now we give the computation method for the front-freezing indicator δ : Ω → {1, ∞}, which takes into account the non-local feature to freeze the front points causing a geodesic path of high values of K. The criterion is based on the fact that the orientations of the tubular structure such as blood vessels usually varies smoothly without sharp turning. The indicator δ is constructed to determine the bending range of the local tubular structure:

  can be described as N * (x) := {ŷ ∈ Z 3 ; x ∈ N (ŷ)}. The fast marching front is comprised of all the grid points tagged as Trial. In each geodesic distance updating step, we apply the method presented in Section 3.2 by setting that x := xmin to compute the value of the indicator K for the latest Trial point 2

Figure 5 :

 5 Figure 5: The first row shows the detected vessel centerlines (blue lines). The source and end points are respectively denoted by red and cyan dots. The following rows display the course of the fast marching fronts propagation. Columns (a) and (c): Geodesic distance maps derived from the ArR metric [16]. Columns (b) and (d): Geodesic distance maps computed by the proposed dArR metric.

  4(b) and (c), we illustrate the spatial component x min of the latest Trial point xmin = (x min , q) and the corresponding 195 projected path Gxmin,q . Moreover, the path in Fig. 4(b) corresponds to a small bending measure is given . If K(x min ) ≤ K 0 which yields δ(x) = ∞ such that the criterion is violated and the propagation at xmin should be frozen, see Fig. 4(c) for an example of the large bending measure. The proposed Riemannian tensor field M dyn is updated during the geodesic 200 distance computation which is actually carried out by computing the front-freezing indicator δ(x, r) in a front advancing procedure as described in Section 3.2. Moreover, we can observe if the propagation should be frozen or not is based on the indicator δ(x, r). In implementation, once a point (x, r) is frozen, all the points at position x with different radius scale [r min , r max ] are frozen and labeled as Accepted. This scheme is helpful to save calculation consuming. It is achieved by defining a new flag map, i.e. a scalar-valued function, F : Ω → {0, 1}, where F(x) = 1 denotes that the point at position x is frozen. During fast marching fronts propagation, once the latest Trial point xmin = (x min , r min ) on the front is chosen, the F(x min ) value can be checked. If F(x min ) = 1, the point should be frozen and labeled as Accepted ;

Fig. 5 .

 5 Fig. 5. Columns (b) and (d) illustrate the geodesic distance maps computed by the our proposed dynamic metric with vesselness map as the scalar function (dArR). The black points in columns (b) and (d) indicate the points which have been frozen due to the violation of the criterion (i.e. δ(•) = ∞). The geodesic distance maps obtained by the anisotropic tubular minimal path method with anisotropic radius-lifted Riemannian metric (ArR) [16] are shown in columns (a)

Figure 6 :

 6 Figure 6: (a) Segmentation accuracy R and (b) Running time corresponding to the local geodesic length Γ and the bending threshold K 0 .

  (a) Origin (b) δ = 0.05 (c) δ = 0.09 (d) δ = 0.05 (e) δ = 0.09

Figure 7 :Figure 8 :

 78 Figure 7: Boundaries (red lines) and centerlines (blue lines) detected from synthetic images. Source and end points are indicated by red and cyan dots, respectively. (a) Original images. (b) and (c): Segmentation results by the proposed dArR metric. (d) and (e): The segmentation results form the ArR metric.

Figure 9 :

 9 Figure 9: Qualitative comparison between the ArR and caArR metrics. The geodesic paths are indicated by yellow lines. The source and end points are indicated by red and cyan dots. Column (a): Vessel anisotropy feature vector filed indicated by lines. Columns (b) and (c): Geodesic paths extracted by the ArR and caArR metrics, respectively. The pictures in row 2 respectively illustrate the details in the regions within the red rectangles of Row 1.

Figure 10 :

 10 Figure 10: Qualitative comparison results between the ArR and caArR metrics on retinal image patches. Blue lines represent the centerlines. Yellow and cyan dots represent the source and end points. (a) Patches from retinal images. (b) and (c): Geodesic paths extracted by the ArR and caArR metrics. (d) The segmented vessel regions.

Figure 11 :

 11 Figure 11: Qualitative comparison results between the ArR and dArR metrics on retinal images. Blue line represent the centerline. Blue lines represent the centerline. Yellow and cyan dots represent the source and end points. (a) Patches from retinal images. (b) and (c): Geodesic paths extracted by the ArR and dArR metrics. (d) The segmented vessel regions.

Figure 12 :

 12 Figure 12: Qualitative comparison results for the AFC and dArR metrics on retinal images. The blue lines represent the extracted geodesic paths. The blue and cyan dots represent the source and end points. (a): Patches from retina images. (b) and (c): Geodesic paths from the AFC [29] and dArR metrics, respectively.

Figure 13 :

 13 Figure 13: (a) Ground truth. (b) Image patch with two given points (yellow and cyan dots). (c) The extracted vessel centerline and boundary by the proposed metric. (d) The segmented vessel region.

  4.

  410

. 14 ,

 14 we illustrate the qualitative comparison for the ArR metric, the caArR metric and the dArR metric, which are shown from the second to the fourth columns. The first row shows the segmentation results from the image with river crossing several bridges whose intensities are similar to these in the background region. The ArR and caArR metrics-driven minimal paths miss the river direction from the middle located bridge, where the bending of the geodesic path varies strongly. The bending constraint is able to effectively handle this problem. The second row illustrates the segmentation results of a satellite image with two neighbouring rivers and complex background. We aim to detect the smaller one by given two points. We observe that the correct river region cannot be recognized by the ArR metric. The dArR-based minimal paths are capable of tracking the favourable river correctly without being effected by the stronger river.

Figure 14 :

 14 Figure 14: Qualitative comparison results between the ArR and dArR metrics on satellite images. The Red region describes the detected road or river and the blue line represents the target centerline. The yellow and cyan dots represent the source and end points. (a) Image patches. (b) to (d): The segmentation results derived from the ArR, caArR and dArR metrics, respectively.

Table 1 :

 1 Quantitative comparison results on synthetic images.

				Noise level δ
	Method
		0.01	0.03	0.05	0.07	0.09
	ArR Metric	0.828 0.836 0.843 0.834 0.842
	dArR Metric 0.875 0.861 0.863 0.860 0.867
	proposed dynamic metrics.

Table 3 :

 3 Quantitative comparison results of vessel centerlines detection on retinal images from the DRIVE and IOSTAR datasets.

	Dataset ArR caArR dIsoM dArR dArR-M dArR-S
	DRIVE	0.425 0.557	0.803	0.849	0.842	0.850
	IOSTAR	0.683 0.730	0.832	0.862	0.859	0.865

Table 4 :

 4 Computation time (in seconds) of each part of the proposed method on retinal images from the DRIVE and IOSTAR datasets.

	Dataset Vesselness Metric	FM
	DRIVE	1.977s	0.118s 2.080s

In the space Ω, a point x ∈ Ω is an ordered pair (x, r) where x ∈ Ω ⊂ R

and r ∈ [r min , rmax].

the point of smallest geodesic distance value among all the Trial points

Indeed, the authors[START_REF] Chen | Minimal paths for tubular structure segmentation with coherence penalty and adaptive anisotropy[END_REF] proposed to use one more step to obtain the vessel thickness with a region-constrained minimal path technique.

the A-V ground truth by the accuracy measurement R. The segmentation accu-

Acknowledgments

We would like to thank the editor and the anonymous referees for their valuable comments, which greatly improved this manuscript. This work was supported in part by the National Key Research and Development Program of China (2017YFC0107903, 2017YFC0109202, 2018ZX10201002-003), by the National Natural Science Foundation of China (Nos. 61902224, 61876037, 61271312, 61401085, 31571001, 31640028, 31400842, 61572258 and 11301074), by the Shortterm Recruitment Program of Foreign Experts (WQ20163200398), by Roche pharma (project AMD short) and by a grant from the French Agence Nationale de la Recherche ANR-16-RHUS-0004 (RHU TRT cSVD). The first author thanks China ScholarShip Council (CSC) for the financial support to finish this work.


Sohini Upadhyay

Mikhail Yurochkin

Mayank Agarwal

Yasaman Khazaeni

Djallel Bouneffouf

Graph Convolutional Network Upper Confident Bound: Online Semi-Supervised Learning with Bandit Feedback

We formulate a new problem at the intersection of semi-supervised learning and contextual bandits, motivated by several applications including clinical trials and ad recommendations. We demonstrate how Graph Convolutional Network (GCN), a semisupervised learning approach, can be adjusted to the new problem formulation. We also propose a variant of the linear contextual bandit with semisupervised missing rewards imputation. We then take the best of both approaches to develop multi-GCN embedded contextual bandit. Our algorithms are verified on several real world datasets.

Introduction

We formulate the problem of Online Partially Rewarded (OPR) learning. Our problem is a synthesis of the challenges often considered in the semi-supervised and contextual bandit literature. Despite a broad range of practical cases, we are not aware of any prior work addressing each of the corresponding components. Next we justify each of the keywords and give motivating examples.

• Online: data is often naturally collected over time and systems are required to make predictions (take an action) before they are allowed to observe any response from the environment.

• Partially: oftentimes there is no response available, e.g. a missing label or environment not responding to system's action.

• Rewarded: in the context of online (multiclass) supervised learning we assume that the environment will provide the true label -however in many practical systems we can only hope to observe feedback indicating whether our prediction is good or bad (1 or 0 reward), the latter case obscuring the true label for learning.

Practical scenarios that fall under the umbrella of OPR range from clinical trials to dialog orchestration. In clinical trials, reward is partial, as patients may not return for followup evaluation. When patients do return, if feedback on their treatment is negative, the best treatment, or true label, remains unknown and the only available information is a reward of 0 for the treatment administered. In dialog systems, a user's query is often directed to a number of domain specific agents and the best response is returned. If the user provides negative feedback to the returned response, the best available response is uncertain and moreover, users can choose to not provide feedback at all. In many applications, obtaining labeled data requires a human expert or expensive experimentation, while unlabeled data may be cheaply collected in abundance. Learning from unlabeled observations to improve prediction performance is the key challenge of semi-supervised learning [START_REF] Chapelle | Semi-supervised learning[END_REF]. One of the possible approaches is the continuity assumption, i.e. points closer to each other in the feature space are more likely to share a label [Seeger, 2000]. When the data has a graph structure, another approach is to perform node classification using graph Laplacian regularization, i.e. penalizing difference in outputs of the connected nodes [START_REF] Zhu | [END_REF]. The latter approach can also be applied without the graph under the continuity assumption by building similarity based graph. We note that the problem of online semi-supervised leaning is rarely considered, with few exceptions [Yver, 2009;Valko et al., 2012;Bouneffouf et al., 2020b]. In our setting, the problem is further complicated by the bandit-like feedback in place of labels, rendering the existing semi-supervised approaches inapplicable. We will however demonstrate how one of the recent approaches, Graph Convolutional Networks (GCN) [START_REF] Kipf | [END_REF], can be extended to our setting.

The multi-armed bandit problem provides a solution to the exploration versus exploitation trade-off, informing a player how to pick within a finite set of decisions while maximizing cumulative reward in an online learning setting. Optimal solutions have been developed for a variety of problem formulations [Auer et al., 2002a;Auer et al., 2002b;[START_REF] Allesiardo | [END_REF][START_REF] Bouneffouf | Djallel Bouneffouf and Raphaël Féraud. Multi-armed bandit problem with known trend[END_REF]Balakrishnan et al., 2018;[START_REF] Bouneffouf | Optimal exploitation of clustering and history information in multi-armed bandit[END_REF][START_REF][END_REF]Bouneffouf et al., 2017a;[START_REF] Yu | [END_REF]Bouneffouf et al., 2017b]. These approaches do not take into account the relationship between context and reward, potentially inhibiting overall performance. In Linear Upper Confidence Bound (LINUCB) [START_REF] Li | A contextual-bandit approach to personalized news article recommendation[END_REF]Chu et al., 2011;[START_REF] Bouneffouf | Djallel Bouneffouf and Irina Rish. A survey on practical applications of multi-armed and contextual bandits[END_REF]Bouneffouf et al., 2020a;Bouneffouf, 2020] and in Contextual Thompson Sampling (CTS) [Agrawal and Goyal, 2013], the authors assume a lin-ear dependency between the expected reward of an action and its context; the representation space is modeled using a set of linear predictors. These algorithms assume that the bandit can observe the reward at each iteration, which is not the case in our setting. Several authors have considered variations of partial/corrupted rewards [START_REF] Bartók | [END_REF][START_REF] Gajane | [END_REF], however the case of entirely missing rewards has not been studied to the best of our knowledge.

In this paper we focus on handling the problem of online semi-supervised learning with bandit feedback. We first review some existing methods in the respective domains and propose extensions to each of them to accommodate our problem setup. Then we proceed to combine the strengths of both approaches to arrive at an algorithm well suited for the Online Partially Rewarded learning as demonstrated with experiments on several real datasets.

Preliminaries

In this section we review two approaches coming from the respective domains of semi-supervised learning and contextual bandits, emphasizing their relevance and shortcomings in solving the OPR problem.

Graph Convolutional Networks

Neural networks have proven to be powerful feature learners when classical linear approaches fail. Classical neural network, Multi Layer Perceptron (MLP), is dramatically overparametrized and requires copious amounts of labeled data to learn. On the other hand, Convolutional Neural Networks are more effective in the image domain [START_REF] Krizhevsky | [END_REF], partially due to parameter sharing exploiting relationships between pixels. Image structure can be viewed as a grid graph where neighboring pixels are connected nodes. This perspective and the success of CNNs inspired the development of convolution on graphs neural networks [START_REF] Henaff | [END_REF][START_REF] Defferrard | [END_REF][START_REF] Bronstein | [END_REF] based on the concept of graph convolutions known in the signal processing communities [Shuman et al., 2013]. Though all these works are in the realm of classical supervised learning, the idea of convolving signal over graph nodes is also widely applied in semi-supervised (node classification) learning [Belkin et al., 2006], where the graph describes relationships among observations (cf. grid graph of features (pixels) in CNNs). [START_REF] Kipf | [END_REF] proposed Graph Convolutional Network (GCN), an elegant synthesis of convolution on graphs ideas and neural network feature learning capability, which significantly outperformed prior semi-supervised learning approaches on several citation networks and knowledge graph datasets.

To understand the GCN method, let X ∈ R T ×D denote a data matrix with T observations and D features and let A denote a positive, sparse, and symmetric adjacency matrix A of size T × T . The GCN embedding of the data with one hidden layer of size L is g(X) =  ReLU( ÂXW ) ∈ R T ×L , where  is degree normalized adjacency with self connections:  = (D + I T ) -1/2 (A + I T )(D + I T ) -1/2 and D ii = T j=1 A ij is the diagonal matrix of node degrees. W ∈ R D×L is a trainable weight vector. Resulting embedded data goes into the softmax layer and the loss for backpropagation is computed only on the labeled observations. The product ÂX gives the one-hop convolution -signal from a node is summed with signals from all of its neighbours achieving smooth transitions of the embeddings g(X) over the data graph. Although a powerful semi-supervised approach, GCN is not suitable for the Online and Rewarded components of OPR. It additionally requires a graph as an input, which may not be available in some cases.

Contextual Bandit

Following [START_REF] Langford | [END_REF], the contextual bandit problem is defined as follows. At each time t ∈ {1, ..., T }, a player is presented with a context vector x t ∈ R D , x t ≤ 1 and must choose an arm k ∈ {1, ..., K}. r t,k ∈ [0, 1] is the reward of the action k at time t, and r t ∈ [0, 1] K denotes a vector of rewards for all arms at time t. We operate under the linear realizability assumption, i.e., there exist unknown weight vectors

θ * k ∈ R D with θ * k 2 ≤ 1 for k = 1, . . . , K so that ∀k, t : E[r t,k |x t ] = θ * k x t
Hence, the r t,k are independent random variables with expectation x t θ * k . One solution to the contextual bandit problem is the LIN-UCB algorithm [START_REF] Li | A contextual-bandit approach to personalized news article recommendation[END_REF] where the key idea is to apply online ridge regression to incoming data to obtain an estimate of the coefficients θ k for k = 1, . . . , K. At each step t, the LINUCB policy selects the arm with the highest upper confidence bound of the reward k

(t) = argmax k (µ k + σ k ), where µ k = θ k x t is the expected reward for arm k, σ k = α x t A -1
k x t is the standard deviation of the corresponding reward scaled by exploration-exploitation trade-of parameter α (chosen a priori) and A k is the covariance of the k-th arm context. LINUCB requires a reward for the chosen arm, r t,k(t) , to be observed to perform its updates. In our setting reward may not be available at every step t, hence we need to adjust the LINUCB algorithm to learn from data with missing rewards.

Proposed algorithms

In this section we formally define Online Partially Rewarded (OPR) problem and present a series of algorithms, starting with natural modifications of GCN and LINUCB to suit the OPR problem setting and conclude with an algorithm building on strengths of both GCN and LINUCB.

Problem setting

We now formally define each of the OPR keywords:

• Online: at each step t = 1, . . . , T we observe observation x t and seek to predict its label ŷt using x t and possibly any information we had obtained prior to step t.

• Partially: after we made the prediction ŷt , environment may not provide any feedback (we will use -1 to encode absence of feedback) and we have to proceed to step t+1 without knowledge of the true y t .

• Rewarded: suppose there are K possible labels y t ∈ {1, . . . , K}. The environment at step t, if it responds to our prediction ŷt , will not provide true y t , but instead a response h t ∈ {-1, 0, 1}, where h t = 0 indicates ŷt = y t and h t = 1 indicates ŷt = y t (-1 indicates missing response).

Note on absence of environment response. We assume that there is no dependence on x t in whether environment will respond or not. This is a common setting in semi-supervised learning [START_REF] Chapelle | Semi-supervised learning[END_REF] -we have access to limited samples from the joint distribution of data and labels P(x, y) and larger amount of samples from the data marginal P(x) with the goal to infer P(y|x) using both. This assumptions is justified in some applications of interest, e.g. whether user will provide feedback to the dialog agent is independent of what the user asked.

Rewarded Online GCN

There are three challenges to be addressed to formulate Rewarded Online GCN (ROGCN): (i) online learning; (ii) the environment only responds with 0 or 1 to our predictions and (iii) the potential absence of graph information. As we shall see, there is a natural path from GCN to ROGCN. Suppose there is a small portion of data and labels available from the start, X 0 ∈ R T0×D and y 0 ∈ {-1, 1, . . . , K} T0 , where D is the number of features, K is the number of classes and T 0 is the size of initially available data. When there is no graph available we can construct a k-NN graph (k is a parameter chosen a priori) based on similarities between observations -this approach is common in convolutional neural networks on feature graphs [START_REF] Henaff | [END_REF][START_REF] Defferrard | [END_REF] and we adopt it here for graph construction between observations X 0 to obtain graph adjacency A 0 . We provide details in Section 4. Now that we have X 0 , y 0 , A 0 , we can train GCN with L hidden units (a parameter chosen a priori) to obtain initial estimates of hidden layer weights W 1 ∈ R D×L and softmax weights W 2 ∈ R L×K .

Next we start to observe the stream of data -as new observation x t arrives, we add it to the graph and data matrix, and append -1 (missing label) to y. Then we run additional training steps of GCN and output a prediction to obtain environment response h t ∈ {-1, 0, 1}. Here 1 indicates correct prediction, hence we include it to the set of available labels for future predictions; 0 indicates wrong prediction and -1 an absence of a response, in the later two cases we continue to treat the label of x t as missing. This procedure is summarized in Algorithm 1.

Bounded Imputation LINUCB

Contextual multi-armed bandits offer a powerful approach to online learning when true labels are not available and the environment's response to a prediction is observed at every observation instead. However, in our OPR problem setting, the environment may not respond to the agent for every observation. Classic bandit approach such as Linear Upper Confidence Bound (LINUCB) [START_REF] Li | A contextual-bandit approach to personalized news article recommendation[END_REF]] may be directly applied to OPR, however it would not be able to learn from observations without environment response. We propose to

Algorithm 1 ROGCN 1: Input: W1, W2, X0, y0, Â0 2: Set X = X0, y = y0, Â = Â0 3: for t = T0 + 1 to T do 4: Append xt to X, -1 to y 5:

Update  with new edges if graph information is available or build k-NN similarity graph from X to obtain  6:

Update W1 and W2 through GCN backpropagation with inputs X, Â, y 7:

Retrieve GCN prediction ŷt and observe environment response ht ∈ {-1, 0, 1} for ŷt 8:

if ht = 1 then 9:

Replace last entry of y with ŷt 10:

end if 11: end for combine LINUCB with a user defined imputation mechanism for the reward when environment response is missing. In order to be robust to variations in the imputation quality, we only allow imputed reward to vary within agent's beliefs. To make use of the context in the absence of the reward, we consider a user defined imputation mechanism I(•) : R D → ∆ K-1 , which is expected to produce class probabilities for an input context x t to impute the missing reward. Typically any imputation mechanism will have an error of its own, hence we constrain the imputed reward to be within one standard deviation of the expected reward for the chosen arm:

r t (k, x t ) = max(µ k -σ k , min(I(x t ) k , µ k + σ k )). (1)
As in ROGCN, we can take advantage of small portion of data to initialize bandit parameters b k = for all k ∈ K do 5:

θ k ← A -1 k * b k θ k ← θ k / θ k 2 6: σ k ← α x t A -1 k xt µ k ← θ k xt 7:
end for 8:

Predict ŷt = argmax k (µ k + σ k ), and observe environment response ht ∈ {-1, 0, 1} 9:

if ht = 1 then 10:

A k ← A k + xtx t for k = 1, . . . ,

Multi-GCN embedded UCB

We have presented two algorithms for OPR learning, however both approaches pose some limitations: ROGCN is unable to learn from missclassified observations and has to treat them as missing labels, while BILINUCB assumes linear relationship between data features and labels and even with perfect imputation is limited by the performance of the best linear classifier. Note that the bandit perspective allows one to learn from missclassfied observations, i.e. when the environment response h t = 0, and the neural network perspective facilitates learning better features such that linear classifier is sufficient. This observation motivates us to develop a more sophisticated synthesis of GCN and LINUCB approaches, where we can combine advantages of both perspectives.

To begin, we note that if K = 2, a h t = 0 environment response identifies the correct class, hence the OPR reduces to online semi-supervised learning for which GCN can be trivially adjusted using ideas from ROGCN. To take advantage of this for K > 2 we can consider a suite of GCNs for each of the classes, which then necessitates a procedure to decide which of the GCNs to use for prediction at each step. We propose to use a suite of class specific GCNs, where prediction is made using contextual bandit with context of k-th arm coming from the hidden layer representation of k-th class GCN and, when missing, reward is imputed from the corresponding GCN.

We now describe the multi-GCN embedded Upper Confidence Bound (GCNUCB) bandit in more details. Let g(X) (k) = Â ReLU( ÂXW (k) 1 ) denote the k-th GCN data embedding and let g(X)

(k) t denote the embedding of observation x t . We will use this embedding (additionally normalized to unit l 2 norm) as context for the corresponding arm of the contextual bandit. The advantage of this embedding is its graph convolutional nature coupled with expressive power of neural networks. We note that as we add new observation x t+1 to the graph and update weights of the GCNs, the embedding of the previously observed x 1 , . . . , x t evolves. Therefore instead of dynamically updating bandit parameters b k and A k as it was done in BILINUCB, we maintain set of indices for each of the arms C k = {t : ŷt = k or h t = 1}. At any step we can compute corresponding bandit context covariance and weight estimate using current embedding:

A k = t∈C k g(X) (k) t g(X) (k) t (2) θ k = A -1 k t∈C k r t,k g(X) (k) t , θ k = θ k / θ k 2 (3)
where r t,k is the reward that was observed or imputed at step t for arm k (recall that we are imputing using prediction of the binary GCN corresponding to the arm chosen by the bandit). Now we can compute expected value and standard deviation for the reward on each arm. The prediction is made based on the upper confidence bounds for the rewards of the arms:

µ k = θ k g(X) (k) t σ k = α g(X) (k) t A -1 k g(X) (k) t ŷt = argmax k (µ k + σ k ).
(4)

Then we observe the environment response h t ∈ {-1, 0, 1}.

Unlike ROGCN, GCNUCB is able to learn from mistakes, i.e. when h t = 0 -although as before we don't know the true class, we can be sure that it was not ŷt , hence we can use this information to improve GCN corresponding to the class ŷt . We summarize GCNUCB in Algorithm 3. Similarly to ROGCN and BILINUCB we can use a small amount of data X 0 and labels y 0 converted to binary labels y

(k) 0

∈ {-1, 0, 1} T0 (as before -1 encodes missing label) for each class k to initialize GCNs weights

W (k) 1 , W (k) 2
for k = 1, . . . , K and index sets C k for each of the arms k = 1, . . . , K. Adjacency matrix if not given is obtained as in ROGCN.

Algorithm 3 GCNUCB 1: Input: W (k) 1 , W (k) 2 , C k , r •,k , y (k) 0 ∀k, X0, Â0, α 2: Set y (k) = y (k) 0 k = 1, . . . , K, X = X0, Â = Â0 3: for t = T0 + 1 to T do 4:
Append xt to X, -1 to each of y (1) , . . . , y through GCN backpropagation with inputs X, Â, y (k) for k = 1, . . . , K 7:

Retrieve embeddings g(X)

(k) t ∀k 8:

Compute A k (Eq. ( 2)) and θ k (Eq. ( 3)) ∀k 9:

Make prediction ŷt using Eq. ( 4) and observe environment response ht 10:

if ht = 1 then 11:

For each k replace last entry of y (k) with 1 if ŷt = k and 0 otherwise 12:

Append t to each C k and 1 to r •,k if ŷt = k and 0 otherwise 13:

else if ht = 0 (learning from mistakes) then 14:

Replace last entry of y (ŷ t ) with 0 15:

Append t to C ŷt and 0 to r

Experiments

In this section we compare baseline method LINUCB which ignores the data with missing rewards to ROGCN, BILIN-UCB and GCNUCB -algorithm proposed in this paper. We consider four different datasets: CNAE-9 and Internet Advertisements from the the UCI Machine Learning Repository 1 , Cora2 , and Warfarin [START_REF] Sharabiani | [END_REF]. Cora is naturally a graph structured data which can be utilized by ROGCN, BILINUCB with ROGCN based imputation and GCNUCB. For other datasets we use a 5-NN graph built online from the available data as follows.

Suppose at step t we have observed data points x i ∈ R D for i = 1, . . . , t. Weights of the similarity graph computed as follows:

A ij = exp x i -x j 2 2 σ 2 . ( 5
)
As it was done by [START_REF] Defferrard | [END_REF] we set σ =

1 t t i=1 d(i, i k )
, where d(i, i k ) denotes L 2 distance between observation i and its k-th nearest neighbour indexed i k . The k-NN adjacency A is obtained by setting all but k (excluding itself) corresponding closest entries of A ij , i, j = 1, . . . , t to 0 and symmetrizing. Then, as in [START_REF] Kipf | [END_REF], we add self connections and row normalize  = (D+I T ) -1/2 (A+I T )(D+I T ) -1/2 , where D ii = T j=1 A ij is the diagonal matrix of node degrees.

For pre-processing we discarded features with large magnitudes (3 features in Internet Advertisements and 2 features in Warfarin) and row normalized all observations to have unit l 1 norm.

For all of our algorithms that use GCN we use default parameters of the GCN and Adam optimizer [START_REF] Kingma | [END_REF]. Default parameters are as follows: 16 hidden units, learning rate of 0.01, 0.0005 weight decay, and dropout of 0.5.

To emulate the OPR setting we randomly permute the order of the observations in a dataset and remove labels for some portion (we experiment with three settings: 25%, 50% and 75% missing labels) of the observations chosen at random. For all methods we consider initial data X 0 and y 0 to represent a single observation per class chosen randomly (T 0 = K). At a step t = T 0 + 1, . . . , T each algorithm is given a feature vector x t and is ought to make a prediction ŷt . The environment response h t ∈ {-1, 0, 1} is then observed and algorithms moves onto step t + 1. To compare performance of different algorithms at each step t we compare ŷt to true label y t available from the dataset (but concealed from the algorithms themselves) to evaluate running accuracy. Defined as such, accuracy is inversely proportional to regret. Imputation Methods. We test two different imputation functions I(•) for BILINUCB -a ROGCN and simple kmeans clustering with 10 clusters. Henceforth, we denote these two approaches as BILINUCB-GCN and BILINUCB-KMeans. In BILINUCB-GCN we update ROGCN with incoming observations and use the softmax class prediction to impute missing reward when needed. In BILINUCB-KMeans, we use the mini-batch k-means algorithm to cluster incoming observations online and impute missing reward with the average non-missing reward of all observations in the corresponding cluster.

Running accuracy results. We noticed that BILINUCB with both imputation approaches and GCNUCB are more robust to data ordering when we use baseline LINUCB for first 300 steps and then proceed with the corresponding algorithm (see Figure 1a where aforementioned algorithms and LINUCB perform the same until step 300 and then have individual running accuracies). For all LINUCB based approaches we used exploration-exploitation trade-off parameter α = 0.25. Results are summarized in Table 1. Since ordering of the data can affect the problem difficulty, we performed 10 data resampling for each setting to obtain error bars.

GCNUCB outperforms the LINUCB baseline and our other proposed methods in all of the experiments, validating our intuition that a method synthesizing the exploration capabilities of bandits coupled with the effective feature representation power of neural networks is the best solution to the OPR problem. We see the greatest increase in accuracy between GCNUCB and the alternative approaches on the Cora dataset which has a natural adjacency matrix. This suggests that GCNUCB has a particular edge in OPR applications with graph structure. Such problems are ubiquitous. Consider our motivating example of dialog systems -for dialog systems deployed in social network or workplace environments, there exists graph structure between users, and user information can be considered alongside queries for personalization of responses.

Role of bounding the imputed reward. Notice that on average, a BILINUCB method outperforms LINUCB and ROGCN. To understand the role of these imputation bounds, we analyze the effects of random imputation. We denote this use of random imputation as BILINUCB-Random and ILINUCB-Random as the same without bounding the imputed reward.

We define BILINUCB-KMeans and ILINUCB-KMeans similarly. We summarize the overall accuracy of each method on CNAE-9 in Table 2.

As the purpose of these bounds is to correct for errors in the imputation method, we expect to see its impact the most when imputation is inaccurate. This is exactly what we see in Table 2 and Figure 1a. When we use a reasonable imputation method, k-means, the imputation bounds do not improve, or only make slight improvements to the results. The improvement gain is much more evident with random imputation, and across both imputation methods, the bounds have a larger impact when there is more reward missing.

Visualizing GCNUCB context space. Recall that the context for each arm of GCNUCB is provided by the corresponding binary GCN hidden layer. The motivation for using binary GCNs to provide the context to LINUCB is the ability of GCN to construct more powerful features using graph convolution and neural networks expressiveness. To see how this procedure improves upon the baseline LINUCB utilizing input features as context, we project the context and the corresponding bandit weight vectors, θ 1 , . . . , θ K , for both LINUCB and GCNUCB to a 2-dimensional space using t-SNE [START_REF] Maaten | [END_REF]. In this experiment we analyzed CNAE-9 dataset with 25% missing labels. Recall that the bandit makes prediction based on the upper confidence bound of the regret: argmax k (θ k x k,t + σ k ) and that x k,t = x t ∀k = 1, . . . , K for LINUCB and x k,t = g(X)

(k) t for GCNUCB. To better visualize the quality of the learned weight vectors, for this experiment we set α = 0 and hence σ k = 0 resulting in a greedy bandit, always selecting an arm maximizing expected reward θ k x t,k . In this case, a good combination of contexts and weight vectors is the one where observations belonging to the same class are well clustered and corresponding bandit weight vector is directed at this cluster. For LINUCB (Figure 1b, 68% accuracy) the bandit weight vectors mostly point in the direction of their re- spective context clusters, however the clusters themselves are scattered, thereby inhibiting the capability of LINUCB to effectively distinguish between different arms given the context. In the case of GCNUCB (Figure 1c, 77% accuracy) the context learned by each GCN is tightly clustered into two distinguished regions -one with context for corresponding label and binary GCN when it is the correct label (points with bolded colors), and the other region with context for the label and GCN when a different label is correct (points with faded colors). The tighter clustered contexts allow GCNUCB to effectively distinguish between different arms by assigning higher expected reward to contexts from the correct binary GCN than others, thereby resulting in better performance of GCNUCB than other methods.

Conclusion and Discussion

We have defined and studied the problem of Online Partially Rewarded (OPR) learning, which combines challenges from semi-supervised learning and multi-armed contextual bandits. We have developed ROGCN and BILINUCB -exten-sions of popular algorithms in the corresponding domains to solve the OPR problem. Our main contribution, GCNUCB algorithm, is the efficient synthesis of the strengths of the two approaches. Our experiments show that GCNUCB, which combines feature extraction capability of the graph convolution neural networks and natural ability of contextual bandits to handle online learning with reward (instead of labels), is the best approach for OPR across a LINUCB baseline and other algorithms that we proposed.

  =-1 x t x t for k = 1, . . . , K. Bounded Imputation LINUCB (BILINUCB) is summarized in Algorithm 2.Algorithm 2 BILINUCB 1: Input: α, b, A, I(•) 2: for t = T0 + 1 to T do 3:Update I(•) with xt and retrieve I(xt) ∈ ∆ K-1 4:

Table 1 :

 1 Total average accuracy

			25% Missing labels
		CNAE-9	Internet Ads	Warfarin	Cora
	LINUCB	67.57 ± 2.90	90.08 ± 0.64	53.70 ± 0.77	38.06 ± 3.45
	ROGCN	64.73 ± 2.67	88.22 ± 1.73	47.72 ± 9.40	48.57 ± 7.75
	BILINUCB-GCN	67.27 ± 2.79	89.91 ± 0.73	53.70 ± 0.77	37.66 ± 3.92
	BILINUCB-KMeans	67.69 ± 4.30	90.37 ± 0.63	52.53 ± 4.83	39.11 ± 2.68
	GCNUCB	77.10 ± 1.89	93.14 ± 0.39	55.19 ± 3.40	66.01 ± 1.35
			50% Missing labels
		CNAE-9	Internet Ads	Warfarin	Cora
	LINUCB	64.25 ± 3.55	88.62 ± 0.67	51.87 ± 5.12	38.85 ± 2.74
	ROGCN	65.96 ± 3.69	88.38 ± 1.93	49.37 ± 8.29	47.71 ± 9.25
	BILINUCB-GCN	63.52 ± 3.31	88.40 ± 0.73	51.75 ± 5.32	38.08 ± 2.97
	BILINUCB-KMeans	67.37 ± 5.18	89.95 ± 0.66	54.20 ± 0.30	39.20 ± 1.76
	GCNUCB	74.55 ± 1.82	92.62 ± 0.37	56.51 ± 3.43	63.47 ± 2.26
			75% Missing labels
		CNAE-9	Internet Ads	Warfarin	Cora
	LINUCB	61.67 ± 3.16	86.66 ± 0.99	52.99 ± 2.61	33.92 ± 0.04
	ROGCN	65.67 ± 5.28	88.31 ± 1.81	47.48 ± 5.41	49.63 ± 5.06
	BILINUCB-GCN	61.36 ± 3.79	86.68 ± 1.04	50.04 ± 11.44	32.21 ± 5.99
	BILINUCB-KMeans	57.16 ± 3.57	88.21 ± 0.99	51.21 ± 7.12	32.51 ± 4.98
	GCNUCB	70.82 ± 2.33	91.45 ± 0.89	53.31 ± 2.98	58.29 ± 2.80

Table 2 :

 2 CNAE-9 total average accuracy

	% Reward Missing ILINUCB-Random BILINUCB-Random
	25	67.29 ± 4.15	67.65 ± 4.30
	50	65.67 ± 5.20	67.19 ± 5.37
	75	49.77 ± 4.68	56.36 ± 3.71
	% Reward Missing ILINUCB-KMeans BILINUCB-KMeans
	25	67.92 ± 3.98	67.69 ± 4.30
	50	67.14 ± 4.84	67.37 ± 5.18
	75	56.62 ± 4.40	57.16 ± 3.57

https://archive.ics.uci.edu/ml/datasets.html

https://people.cs.umass.edu/ mccallum/data.html