# A Hybrid Approach for Learning to Shift and Grasp with Elaborate Motion
Primitives

Zohar Feldman1, Hanna Ziesche1, Ngo Anh Vien1, Dotan Di Castro1 1 Bosch Center
for Artificial Intelligence (BCAI)

###### Abstract

Many possible fields of application of robots in real world settings hinge on
the ability of robots to grasp objects. As a result, robot grasping has been
an active field of research for many years. With our publication we contribute
to the endeavor of enabling robots to grasp, with a particular focus on bin
picking applications. Bin picking is especially challenging due to the often
cluttered and unstructured arrangement of objects and the often limited
graspability of objects by simple top down grasps. To tackle these challenges,
we propose a fully self-supervised reinforcement learning approach based on a
hybrid discrete-continuous adaptation of soft actor-critic (SAC). We employ
parametrized motion primitives for pushing and grasping movements in order to
enable a flexibly adaptable behavior to the difficult setups we consider.
Furthermore, we use data augmentation to increase sample efficiency. We
demonstrate our proposed method on challenging picking scenarios in which
planar grasp learning or action discretization methods would face a lot of
difficulties.

##  I Introduction

The most fundamental atomic task of robotic manipulation is grasping. When
successful, it enables the robot to do more complex manipulation tasks such as
pick and place or bin picking [1]. In the bin picking task, several objects
are placed inside a bin, and the goal is to remove all the objects from the
bin and place them at a target position. Bin picking has enormous implications
in many areas of robotic manipulation, e.g., industrial assembly, logistics
automation, domestic robotics, pick and place tasks, and more. In addition,
this task is among one of the most complex robotic manipulation problems,
since it exhibits many challenges, yet to be addressed, such as noise and
occlusions in perception, object obstructions and collision in motion
planning. Therefore, a resilient and robust approach to pick an object out of
a bin is needed.

There are two approaches to robotic manipulation. The classical approach is
based on planning [2, 3]. In this approach, a domain expert carefully plans
what the robot should do in order to achieve the task at hand. The second,
more recent approach consists in learning such skills using machine learning
techniques. The planning in this approach is reduced to the necessary
requirements and the system learns what to do by adopting trial and error
methods such as reinforcement learning [4]. At one extreme side of the
automation-learning spectrum, many recent approaches attempt to solve these
tasks using an end-to-end learning approach based on visual input (i.e., ”from
pixels”; e.g., [5, 6]). There are however also methods, that lie somewhere
between purely planning and purely learning based approaches. One of these
methods is learning based on _motion primitives_. These methods define a set
of required primitives to solve the problem at hand [7] where for each
primitive there are parameters defining how this primitive is executed.

Figure 1:  The robot is tasked with the goal of emptying one bin to another,
shifting the objects (top-right) when necessary in order to facilitate the
grasping (bottom-right). The inferred motion primitives have high flexibility
to be able to cope in challenging conditions

In this work we propose a method based on continuously parametrized motion
primitives. The parameters defining these motion primitives are provided as
output of a deep neural network architecture and can be learned in an end-to-
end fashion. The advantages for this approach are twofold. _Firstly_ , we do
not need to apply discretization to the continuous action of the robot. This
enables computation and memory efficiency during training, and generalization
between similar situations. _Secondly_ , our approach allows us to add more
parameters to each primitive while avoiding the curse of dimensionality of
discretization. Thus we are able to work with actions of a higher degree of
freedom. In other words, in our approach the output of the network scales
linearly in the dimensionality of the actions, while in discretized settings,
the size of the output space increases exponentially.

Specifically, in this work we focus on two primitives: shifting and grasping
as depicted in Fig. 1. For grasping we propose a parameterization with two
discrete parameters (2D position in x-y plane in a heightmap image) and three
continuous parameters (yaw, pitch, and gripper width) while for shifting we
propose a parameterization with two discrete parameters (2D position in x-y
plane in a heightmap image) five parameters (yaw, pitch, roll, shifting
direction and shifting distance). While our algorithm now involves both
discrete and continuous actions, we do not need to employ a joint hybrid
formulation as proposed in [8]. Due to the fact that our continuous actions
depend on the choice of a discrete action, we can instead resort to a
hierarchical reinforcement learning (RL) formulation and policy optimization.
Our formulation shares some similarity with recent work [9] in RL. However
these methods have not been extended to robotic applications, which possess
more modeling challenges, e.g. high-dimensional input spaces.

Our main contributions in this work are as follows:

  1. 1.

We present an algorithm based on Soft Actor Critic (SAC [10]) for learning
complex grasping and shifting policies directly from pixels.

  2. 2.

We devise a bin picking methodology featuring low sample complexity to
overcome the physical limitations for successful grasping. This is
specifically relevant, when using parallel-jaw grippers on objects in hard to
grasp locations.

  3. 3.

We propose a novel network architecture, which is designed to model hybrid
action spaces, i.e. discrete-continuous actions.

The paper is organized as follows. In Section II we review related work for
this paper. In Section III we describe the mathematical framework that we use
throughout the paper. In Section IV we describe the methodology used whereas
in Section V we show the applicability of our work in simulation. We conclude
in Section VI and point out future work directions.

##  II Related work

As a prerequisite for many robot manipulation tasks, grasping objects has been
a long-standing subject of interest in robotic research and a large body of
literature exists on the topic. While early approaches were mainly based on
analytic methods, performing constrained optimization on certain criteria
defining stable grasp (see [11] for a comprehensive review), data-driven
approaches dominate the field since at least around the year 2000. The survey
by Bogh et al. [12] provides an excellent overview of these ’classical’ data-
driven approaches. With the success of deep neural networks, in particular in
the domain of computer vision, also the robotic grasping community has more
and more turned towards deep learning methods. A substantial amount of
publications from the recent years follow the paradigm of supervised learning
and cast the problem of robot grasping as one of predicting oriented
rectangles from visual inputs. As noted first by [13], given a depth image
five degrees of freedom (DoF) are sufficient to specify a rectangle that fully
defines a 7 DoF robot grasping pose (3D position, a 3D rotation, gripper
opening width).

Grasping via supervised learning: Distinct approaches differ in how they
tackle the huge search space of possible rectangles. Lenz et al. [13] suggest
a two step approach relying on one fast but inaccurate step, performing an
exhaustive search over the space of possible rectangles, followed by a slower
accurate step to re-rank the top 100 rectangles. In contrast, Redmon et. al.
[14] and Kumra et al. [15] propose a faster one-step approach, superseding the
exhaustive search step. More recent approaches [16, 17] regress rectangles in
a generative fashion by predicting pixel-based grasp angles, gripper width and
grasp success scores and using these images as input to a final network
predicting the grasp rectangles. Instead of using manually labeled images like
in those work, Mahler et al. [18] combined a grasp sampling heuristics with an
analytic method to determine force closure based on wrench space analysis in
order to generate a large data base of object meshes annotated with grasp
poses and corresponding grasp quality metrics [19]. A similar approach was
proposed by Schmidt et al. [20]. More recent approaches [21, 22, 23] propose a
grasp sampling network that is based on a variational autoencoder, and trained
with simulated or annotated data.

Grasping via self-supervised learning: Another solution to the problem of
training data is provided by direct self-supervision, i.e. by running a
prescribed grasping policy either in simulation or on a real robot to collect
training data annotated automatically by the observed grasp success. This
method was pioneered on a large scale by Pinto et al. [24]. Levine et al. [25]
combined this strategy with a vision-based servoing-mechanism to guide the
robot to successful grasps and scaled the approach further to a multi-robot
distributed setup. In a subsequent publication, the approach was instead
formulated as a RL problem [26]. Khansari et al. [27] built upon this work and
extended it with an image-based action representation, which they found to
increase sample efficiency. Furthermore, they trained purely in simulation and
used data augmentation to successfully transfer to the real world. In a recent
work, Song et al. [28] also propose an extension of [29] to full 6D pose
grasping. In contrast they however do no longer consider non-prehensile
manipulations. Furthermore they rely on a manually collected grasping dataset,
while our approach is trained in a fully self-supervised manner.

Grasping via non-prehensile manipulation: Although some of the methods
mentioned above have been demonstrated or even explicitly trained in cluttered
scenes, grasping actions may not always be enough in the presence of multiple
objects, occlusions or workspace limitations. Therefore, another stream of
research has focused on combining grasping with non-prehensile actions. Dogar
et al. [30] propose a fully planning-based approach using non-prehensile
motion primitives to extract an object from clutter. Boularis et al. [31]
formulate a finite dimensional Markov decision process (MDP) and use tabular
Q-learning to train a policy combining pushing and grasping. Their approach
relies on hand-crafted features and kernel methods to model the value
function, reward and transition function. In [32] the authors investigate
multi-task supervised learning to perform grasping, pushing and poking actions
and find that the multi-task model shows better performance than single-task
models trained on the same amount of data. A hierarchical approach to combine
pushing and grasping primitives has been considered by Danielczuk et al. [33],
who propose a high-level planner that queries the low-level motion primitive-
based policies and chooses actions based on the returned quality estimates. In
a subsequent publication [33] the authors train a visuomotor policy to uncover
objects in clutter in order to increase graspability. To facilitate the
difficult exploration problem, the authors define several teacher heuristics.
Most related to our work are the approaches proposed by Zeng et al. [29] and
by Berscheid et al. [34]. Both define motion primitives for pushing and
grasping actions and learn them using Q-learning on a discretized action space
consisting of a 2D Cartesian grasping position in the (x,y)-image plane and a
yaw rotation around the z-axis. Both methods, however, are restricted to top-
down grasps with discrete yaw rotations. A recent work from Berscheid et al.
[35] proposes a hybrid approach of learning 6 DoF grasping. However this
method is still simplified to planar grasping in which lateral DoF ($z$-depth,
roll, and pitch angles) are inferred from a model-based controller. In our
work, we built on these approaches but extend them to allow for full end-to-
end 6DoF grasp learning and more flexible motion primitives. We achieve this
by giving up the strict discretization in favour of a hybrid discrete-
continuous action space.

##  III Background

Recent years have shown great improvements in model-free deep reinforcement
learning in several domains. In particular SAC (Soft Actor-Critic) [10] has
gained popularity for continuous control settings due to its sample efficiency
and stability and its superior performance both in simulation benchmarks and
in real world settings [36]. For these reasons and due to the fact that SAC
also has shown good performance when trained on pixel inputs [37], we use SAC
as backbone algorithm for our approach. SAC is an off-policy actor-critic
algorithm and as such jointly trains a pair of state-action value functions
$Q^{\pi}_{\phi_{i}}$, $i=1,2$ and an stochastic policy $\pi_{\theta}$. As SAC
is based on the paradigm of maximum-entropy RL, the actor is trained to
maximize the cumulative expected return while at the same time maximizing its
entropy, i.e. acting as stochastic as possible. In standard SAC, the actor is
parametrized as a Gaussian policy $\pi_{\theta}$ and is trained on the
following loss function:

$\displaystyle\mathcal{L}(\theta)=$
$\displaystyle\mathbb{E}_{a\sim\pi_{\theta}}\left[Q^{\pi}(s,a)-\alpha\log\pi_{\theta}(a|s)\right],$

where $Q^{\pi}(s,a)=\underset{i=1,2}{\mathrm{min}}Q_{\phi_{i}}^{\pi}(s,a)$,
and $\alpha$ is a coefficient defining the weight between entropy and reward.
The critics $Q_{\phi_{i}}$ are trained via Double Deep Q-learning [38] with
targets provided by corresponding temporally delayed target networks
$Q_{\bar{\phi}_{i}}$, i.e the critic loss is given by

$\displaystyle\mathcal{L}(\phi_{i})=\mathbb{E}_{\overset{s,a,s^{\prime},r\sim\mathcal{D}}{a^{\prime}\sim\pi_{\theta}(\cdot|s^{\prime})}}\left[\left(Q_{\phi_{i}}(s,a)-\left(r+\gamma\,y_{t}(s^{\prime},a^{\prime})\right)\right)^{2}\right]$

where
$y_{t}(s^{\prime},a^{\prime})=\underset{i=1,2}{\mathrm{min}}Q_{\bar{\phi}_{i}}(s^{\prime},a^{\prime})-\alpha\log\pi_{\theta}(a^{\prime}|s^{\prime})$.
Here, states $s$, actions $a$, next states $s^{\prime}$ and rewards are
sampled from a replay buffer which is continually populated as the training
progresses. The actions $a^{\prime}$ in state $s^{\prime}$ are sampled from
the current policy. In their subsequent publication [36] Haarnoja et al. also
propose a method to automatically tune the hyperparameter $\alpha$, which
implicityly controls the the policy’s exploration.

##  IV Method

We model the Bin Picking task as a finite-horizon Markov Decision Process
(MDP; [39]) $\left(\mathcal{S},\mathcal{A},\mathcal{T},r,\gamma,H\right)$ with
state space $\mathcal{S}$, action space $\mathcal{A}$, transition probability
function $\mathcal{T}$, reward function $r$, discount factor $\gamma$, and
horizon of $H$ steps. In each step $t=1,\ldots,H$, a state
$s_{t}\in\mathcal{S}$ is observed, and an action $a_{t}\in\mathcal{A}$ is
chosen according to a policy $\pi(a_{t}|s_{t})$. Upon the application of
$a_{t}$ in $s_{t}$, a reward $r\left(s_{t},a_{t}\right)$ is received, and the
system transitions into a new state $s_{t+1}$ according to $\mathcal{T}$. We
represent the state $s_{t}$ by a heightmap image with four channels, the color
(RGB) and height (Z) from the table on which the robot and bin are installed.
The heightmap is computed from an image captured by RGB-D camera overlooking
the bin area. Given the camera intrinsic and extrinsic parameters, the image
is transformed to a color point-cloud in the robot coordinate system, whose
origin is conveniently located at the base link attached to the table and
z-axis points in the inverse direction of gravity. The point-cloud is then
projected orthographically onto a 2-dimensional grid with granularity of
$5\times 5$mm in the xy-plane containing the bin. The action $a_{t}$ is a
motion primitive corresponding to either a grasp or a shift, each defined by a
distinct set of parameters, as detailed in Section IV-A. A reward $r_{t}$ of
$1$ is received in step $t$ if $a_{t}$ results in a successful grasp,
otherwise it is $0$. We use RL to train the policy
$\pi\left(a_{t}|s_{t}\right)$ to maximize the $Q$-function defined as

$Q(s_{t},a_{t})\triangleq\mathbb{E}_{s^{\prime}\sim\mathcal{T}(s^{\prime}|s,a)}\left[\sum_{i=t}^{H}\gamma^{i}r(s_{i},a_{i})\right].$

The Bellman equation,

$Q_{t}(s_{t},a_{t})=\mathbb{E}_{\pi,\mathcal{T}}\left[r(s_{t},a_{t})+Q_{t+1}(s_{t+1},a_{t+1})\right],$
(1)

provides a means to evaluate the policy $\pi$ by calculating the $Q$-function
recursively, and lies in the core logic of most (off-policy) RL algorithms,
including the algorithm presented here.

###  IV-A Motion Primitives

In each step the policy outputs a primitive type
$\phi\in\left\\{\textbf{g}\text{(rasp)},\textbf{s}\text{(hift)}\right\\}$ as
well as the primitive parameters sets $a^{\phi}$ that together define the
maneuver to be executed. Each primitive is executed essentially in open-loop
as follows.
Grasp: The gripper is first oriented according to the inferred Euler angles
$\left(i^{g},j^{g},k^{g}\right)$. Subsequently, the gripper tool-center-point
(TCP) is moved to a point above the target position defined by the Cartesian
coordinates $\left(x^{g},y^{g},z^{g}\right)$ and then lowered to the target
height $z^{g}$ having the width between the fingers set to $w^{g}$ cm. Upon
reaching the target pose or detecting a collision, the gripper is closed and
lifted up $20$cm at which point the gripper is again signaled to close. The
grasp is considered successful if the distance read between the fingers
exceeds a threshold set to a value slightly below the smallest size of all
objects in consideration. The parameter set
$a^{\text{g}}=\left(x^{g},y^{g},j^{g},k^{g},w^{g}\right)$ includes all the
above parameters, except $z^{g}$ which is extracted from the heightmap in
position $\left(x^{g},y^{g}\right)$, and the roll angle $i^{g}$ which is set
to $0$ to ensure equal finger heights for stable top-down grasping.
Shift: The gripper TCP is moved to a target pose
$\left(x^{s},y^{s},z^{s},i^{s},j^{s},k^{s}\right)$ with clamped fingers, after
which it is moved $\overrightarrow{d^{s}}$ cm in a horizontal direction
defined by a rotation angle $\overrightarrow{k^{s}}$ around the z-axis. Here,
$a^{\text{s}}=\left(x^{s},y^{s},i^{s},j^{s},k^{s},\overrightarrow{d^{s}},\overrightarrow{k^{s}}\right)$
whereas, similarly to the grasp primitive, $z^{s}$ is extracted from the
heightmap.

###  IV-B Reinforcement Learning for Bin Picking

####  IV-B1 Architecture

We use a Fully Convolutional Network (FCN) to infer the parameter set
$a^{\phi}$ and approximate the value $Q^{\phi}(s,a^{\phi})$ for each primitive
type $\phi$ and heightmap image $s$. The underlying algorithmic framework and
network architecture follows what can be schematically viewed as a combination
of convolutional SAC for continuous actions and $Q$-learning for discrete
actions. Namely, the network infers an embedding for each of the heightmap
pixels using a Pixel Encoder module. By viewing each pixel embedding as a
distinct state, an Actor module convolutes over these ”states” and infers a
Gaussian action for each of them, resulting in an action map $A^{\phi}$.
Actions are concatenated with their corresponding pixel embedding and
evaluated by a Critic module, resulting in a $Q$-value map $Q^{\phi}$. Figure
2 depicts the network architecture which is shared between both primitives.
Notably, each pixel $(h^{\prime},w^{\prime})$ in the action map $A^{\phi}$
induces a fully specified parameters set $a^{\phi}$. The spatial parameters
$(x^{\phi},y^{\phi})$ are derived from the location association of the
heightmap pixels, whereas the rest of the parameters are given by
$A^{\phi}[h^{\prime},w^{\prime},:]$, the values in the action map
corresponding to the pixel coordinates along all channels. Similarly,
$Q^{\phi}[h^{\prime},w^{\prime},1]$ represents the $Q$-value of the state-
action pair $(s,a^{\phi})$. The map $Q^{\phi}$ thus represents a $Q$-function
$Q^{\phi}(s,a^{\phi})$ for a discrete set of actions corresponding to the
heightmap pixels, and as such, it is trained using a $Q$-learning schema for
discrete actions. The actor, on the other hand, is trained using a soft actor-
critic schema. Figure 2 depicts this action selection process. Note that for
training the Q-network, the Actor module needs not be executed. Instead, we
can concatenate action parameters $a^{\phi}$ stored in a replay buffer with
the pixel embedding of the corresponding state $s$ and then evaluate it using
the critic to predict a $Q$-value.

####  IV-B2 Inference

Inferring an action for a given state $s$ entails computing all the network
modules end-to-end, producing the action map $A^{\phi}$ and $Q^{\phi}$-map for
both primitives. The primitive type is selected as
$\phi^{\ast}=\arg\max_{\phi}\max_{h^{\prime},w^{\prime}}Q^{\phi}[h^{\prime},w^{\prime},1]$.
The primitive parameters set is calculated by first selecting the best pixel
for $\phi^{\ast}$,
$(h^{\ast},w^{\ast})=\arg\max_{h^{\prime},w^{\prime}}Q^{\phi^{\ast}}[h^{\prime},w^{\prime},1]$
and extracting $a^{\phi^{\ast}}$ from $(h^{\ast},w^{\ast})$ as described
above.

####  IV-B3 Training

We collect data throughout the experiments and store it in a replay memory.
When sampling mini-batches for training, we use data augmentation in order to
increase the sample efficiency [37]. Particularly, we create task-invariant
versions of the sampled experiences $\left(s_{t},a_{t},r_{t},s_{t+1}\right)$
by rotating the heightmap image $s_{t}$ by a random angle, as well as rotating
the relevant angles in $a_{t}$ by the same angle. In both primitives, we
rotate the yaw angle $k^{\phi}$, while in the shift primitive we also rotate
the shift direction $\overrightarrow{k^{s}}$. By that, we create new samples
that induce similar results, that is, similar $r_{t}$ and $s_{t+1}$. We train
each of our models using the following loss functions:
Critic Loss:

$\mathcal{L}_{\text{critic}}=\begin{cases}\text{BCE}\left(Q_{t}^{\phi}(s_{t},a_{t}^{\phi}),y_{t}\right)&t=H\\\
\text{MSE}\left(Q_{t}^{\phi}(s_{t},a_{t}^{\phi}),y_{t}\right)&\text{otherwise}\end{cases},$

where $y_{t}=r_{t}+\gamma\max_{\phi,a}Q_{t+1}^{\phi}\left(s_{t},a\right)$, BCE
is the Binary Cross Entropy loss, and MSE is the Mean Squared Error loss. $H$
denotes a fixed predefined horizon. Note that we use the BCE loss for the
final step only , since the corresponding ground label
$y_{H}=r_{H}\in\\{0,1\\}$. For all other steps $t<H$, $y_{t}$ is a real number
and we thus use the MSE loss. The critic loss is minimized with respect to the
network parameters of the Pixel Encoder, Pixel Action Encoder, and Critic
modules.
Actor Loss:

$\mathcal{L}_{\text{actor}}=Q_{t}^{\phi}(s_{t},a_{t}^{\phi})-\alpha\log\pi_{t}^{\phi}(a_{t}^{\phi}|s_{t})$

The actor loss is minimized with respect to the parameters of the Pixel
Encoder and Actor modules.

###  IV-C Implementation Details

##### Finite Horizon

In finite-horizon MDPs, the Q-function is time-dependent and, respectively,
the Q-functions in the different steps should be approximated by distinct
networks to be theoretically sound. This, however, entails training $H$ neural
networks, which may impose a large computational burden. Most RL
implementations circumvent this issue by treating the MDP as infinite-horizon,
regardless of the actual case, and using the discount factor to moderate the
effect of future steps. Here, we choose to use distinct networks for the
different steps and instead take different relaxing measures. Regardless of
the number of steps allowed to empty the bin, we use a fixed and small horizon
of $H=2$ and train two networks, $\Phi_{0}$ and $\Phi_{1}$. While this choice
helps alleviating the hurdles mentioned above which are even amplified by
having sparse rewards and large action space, it can also be motivated by the
observation that bin picking typically does not benefit much from looking more
than a few steps ahead. Indeed, looking ahead beyond the current state is
mostly beneficial when shifting is required to allow a consecutive grasp, and
most likely a single shift manoeuvre should suffice in this case. Respectively
of this relaxation, we use $\Phi_{0}$ to infer an action in step $t<H$ and
$\Phi_{1}$ for $t=H$. During training, we use all the recorded experiences for
updating the networks of all steps independently of the actual step they were
performed within the episode.

##### Exploration Heuristics

In order to increase the chances of seeing an effective result for either a
grasp attempt or a shift attempt when taking exploration steps, we use a
change detection procedure as described in [40] to locate pixels that
correspond to objects. Very roughly, we subtract the point-cloud of the
current image state from a reference point-cloud of an image with empty bin
and mask pixels in which there is a notable difference. We then sample a
portion of the exploration actions from a uniform distribution over all the
object pixels. We also assume to have at our disposal a bounding box of the
bin, which is either known or can be obtained by some external perception
utility. We define sites of interest on the gripper and wrist mounted camera,
and transform these points according to any candidate target pose to verify
its feasibility by checking if the transformed points fall inside the bounding
box of the bin. If there is at least one point that falls outside the bin, we
know that a reach attempt would result in a collision and therefore discard
the attempt in advance. We also use this computation as complementary
exploration heuristic for searching feasible orientations for a given
translation.

##### Mini-Batch Sampling

To overcome the impact of sparse rewards on training, particularly as early
stages, we sample mini-batches that are comprised of half success samples and
half failures.

Figure 2:  Action selection network: Both grasping and shifting primitives
$\phi$ share this architecture. The Pixel Encoder FCN module encodes the
spatial features of each of the image pixels, whereas all other FCN modules
act as Fully Connected Layers in a context of a single pixel using kernels of
size $1$.

##  V Experiments

###  V-A Simulation environment

In order to assess the performance of our algorithm in a binpicking scenario,
we set up a suitable MuJoCo [41] environment, as depicted in Fig. 1. This
environment contains a Franka Emika Robot with a Realsense d415 camera mounted
on its wrist. Two bins of particularly small size $30\,\text{cm}\times
40\,\text{cm}$ and with a height of 20 cm are placed in front of the robot. We
purposely create this challenging setting to evaluate the limit of different
approaches, and demonstrate the benefit of high DoF grasp learning methods
like ours. First, the small size of the bin renders direct grasping infeasible
for a large fraction of the area inside the bin. Secondly, the high bin walls
preclude top-down grasps for many of the remaining graspable object positions,
as these would lead to collisions with the walls. Lastly, the wrist mounted
camera adds to the difficulty of avoiding collisions. We use this setup both
for training the algorithm and for evaluation. During all experiments we spawn
4 objects from a predefined set at a random position above one bin at the
beginning of each episode. In each step we execute the action returned from
the policy $\pi$ as described in Section IV. We use ROS for control and
execute all motion commands using a Cartesian impedance controller.

###  V-B Experimental setup

We conduct several experiments to evaluate our algorithm and compare it to one
following baseline, called _Learning shifting for grasping_ proposed by
Berscheid et al. [42]. Similar to Zeng et. al. [7], this baseline resorts to
action discretization and is restricted to planar grasps, i.e. allows only for
rotations around the axis perpendicular to the table top. With our experiments
we aim to answer primarily the following questions:

  * •

Do the _non-prehensile shifting actions_ improve the performance of our bin
picking algorithm, in particular in challenging settings?

  * •

Do our flexible and continuously parametrized motion primitives show an
advantage over _fully discrete parametrizations_?

  * •

How well can our learned solution _generalize to unseen objects_?

In order to answer these questions, we set up two different scenarios. All
scenarios share the same general setup described above and differ only in the
way the objects are placed into the bin. In particular, we consider the
following two settings: (1) Randomized bin: We spawn objects at random
positions above the bin and let them drop into it. This scenarios is also used
for training, and (2) Close-to-wall bin: We spawn objects only at random
positions close to the bin walls. This is a very challenging setting as the
spawn objects are likely very close to the bin wall, so grasping such objects
often requires shifting. In other cases if an object can be grasped without a
support of shifting, its grasp pose might be harder to find as it involves
more constraints, e.g. stable grasp pose, collision avoidance to the bin walls
and nearby objects.

The first setting allows us a general evaluation of our algorithm in a
challenging setup and a comparison to the chosen baselines. The second setting
is designed to test our algorithm in situations where direct grasps and top-
down grasps are impossible and therefore flexible shifting and grasping
motions are necessary. For each of the test scenarios and algorithms we
execute 100 runs and compare average performance according to the following
metrics, proposed by Zeng et al. [7]:
Clearance: percentage % of objects cleared from the bin
Completion rate: percentage % of runs in which the robot succeeds in emptying
the bin without executing 10 successive unsuccessful actions.
Success rate: average percentage % of successful grasps per episode.
Action efficiency: defined as $\frac{\text{\\# objects in scenario}}{\text{\\#
actions before completion}}$ this metric measures the capability of the policy
to avoid unnecessary actions.
Another important metric is the compute time required for inferring the best
action from the camera input. Our method takes $\approx 0.01s$ when run on a
GeForce RTX 2080. It is thus real time capable and would even allow for a
closed-loop setup.

TABLE I: Experimental Results for the Random Setting (Mean $\%$) Method | Clearance | Completion | Grasp Succ. | Action Eff.
---|---|---|---|---
[42] | $72$ | $32$ | $30$ | $70$
Ours | 96 | 85 | 80 | 76
[42] (G) | $68.3$ | $21.7$ | $27.4$ | $64.5$
Ours (G) | 85 | 51 | 55 | 85
[42] (U) | $31.9$ | $1.5$ | $43.5$ | $41.9$
Ours (U) | 82.4 | 53 | 46 | 73
TABLE II: Experimental Results for the Challenging Setting (Mean $\%$) Method | Clearance | Completion | Grasp Succ. | Action Eff.
---|---|---|---|---
[42] | $48$ | $1$ | $15$ | 66
Ours | 88 | 66 | 40 | $34$
[42] (G) | $23.8$ | $1.6$ | $7.7$ | $62.3$
Ours (G) | 72 | 30 | 33 | 72
[42] (U) | $33.7$ | $0.9$ | $11.6$ | 80
Ours (U) | 77 | 49 | 39 | $62$

###  V-C Evaluation results

Tables I and II summarize the results of the random and challenging settings,
respectively. In the basic experiments, the robot is presented with the same
object types that were used for training, cubes and rods. In addition, we test
each algorithm in modifications. In the first (G), we allow only grasp
actions, while in the second (U), the robot is given unseen objects,
specifically cans, elongated boxes, balls, and tubes. As expected, in the
challenging setting the performance was worse than the in the random setting
emphasizing the fact that most failure cases result from objects being placed
at location from which they cannot be grasped with direct top-down grasping.
It is also apparent that the method we propose gave better results than [42].
Here, the relatively tall bin renders top-down grasping and shifting
practically impossible in a large area of the bin, as opposed to the setting
described in [42]. Therefore grasping and shifting attempts of objects near
the bin walls were mostly ineffective. We also observed that most of these
hard cases were solved by applying 6DoF grasps which allow the gripper to
reach objects without colliding with the bin wall, and much less by shifting.
The typical failure cases were concerned with either not being able to learn
the right shift approach, having repeated attempts in misleading cases (e.g.
trying to grasp to close rods that appear as one), and in some less common
occasions having spawned the objects outside the bin. It is also apparent that
our solution is preferable in respect to generalization. While both algorithms
show performance degradation in presence of unseen objects, the impact is
significantly lower on our algorithm that the baseline. Notably baseline
sometimes shows better action efficiency in the challenging setting. This is
due to fact that only completed scenes contribute to this metric and that the
baseline is not able to solve most of the scenes that require many shifts and
therefore necessarily come with a low action efficiency.

##  VI Conclusion

This paper proposes a unified actor-critic framework for learning continuous
and high-DoF pushing and grasping policies in robotic bin-picking. Bin-picking
on industrial applications requires non-trivial grasping skills, e.g. avoiding
collisions with the bin wall and nearby objects. It additionally requires non-
prehensible actions, e.g. shifting. Existing methods often resort to a
discretization of the action space. As a result, they scale poorly to high-
dimensional action spaces and cannot easily be extended to 6DoF grasp learning
or even higher DoF shifting skill learning. The main contributions of our
proposed framework can be summarized as follows. _First_ , this framework
easily allows more DoFs for the action space based on a well-known actor-
critic RL framework. Hence it enables continuous actions efficiently. _Second_
, we made an important contribution on the choice of the network architecture
with which we can easily model hybrid action spaces, i.e. discrete-continuous
actions. The discrete actions are primitives (grasping, shifting), and
grasp/shift centers ($xy$-planar), while the continuous actions are
grasp/shift configurations.

We demonstrate the proposed framework on a challenging setting which is common
in many industrial applications, e.g. featuring a small bin, a high wall bin,
ubiquitous collisions etc. This setting requires both high DoF grasp/shift
skill configurations, and non-prehensible shifting actions. We showed that our
approach is able to complete the tasks with higher rates in comparison to
existing baselines. While our method can show great benefits and potentials,
there are open research questions yet to be addressed. First, we plan to
evaluate how it can handle sim2real transfer in which a policy trained
extensively in simulation can be evaluated in real bin-picking scenarios.
Second, beside the two grasp and shift primitives we would like to see how
more primitives, e.g. hand-over or throwing, can also be integrated and learnt
in an end-to-end fashion using this framework.

## References

  * [1] I. Lenz, H. Lee, and A. Saxena, “Deep learning for detecting robotic grasps,” _The International Journal of Robotics Research_ , vol. 34, no. 4-5, pp. 705–724, 2015.
  * [2] H. B. Maynard and K. B. Zandin, _Maynard’s industrial engineering handbook_ , 2001.
  * [3] G. Chryssolouris, _Manufacturing systems: theory and practice_. Springer Science & Business Media, 2013.
  * [4] R. S. Sutton and A. G. Barto, _Reinforcement learning: An introduction_. MIT press, 2018.
  * [5] S. Levine, C. Finn, T. Darrell, and P. Abbeel, “End-to-end training of deep visuomotor policies,” _The Journal of Machine Learning Research_ , vol. 17, no. 1, pp. 1334–1373, 2016.
  * [6] T. Zhang, Z. McCarthy, O. Jow, D. Lee, X. Chen, K. Goldberg, and P. Abbeel, “Deep imitation learning for complex manipulation tasks from virtual reality teleoperation,” in _ICRA_. IEEE, 2018, pp. 1–8.
  * [7] A. Zeng, S. Song, S. Welker, J. Lee, A. Rodriguez, and T. Funkhouser, “Learning synergies between pushing and grasping with self-supervised deep reinforcement learning,” in _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2018, pp. 4238–4245.
  * [8] M. Neunert, A. Abdolmaleki, M. Wulfmeier, T. Lampe, T. Springenberg, R. Hafner, F. Romano, J. Buchli, N. Heess, and M. Riedmiller, “Continuous-discrete reinforcement learning for hybrid control in robotics,” in _Conference on Robot Learning_. PMLR, 2020, pp. 735–751.
  * [9] O. Delalleau, M. Peter, E. Alonso, and A. Logut, “Discrete and continuous action representation for practical rl in video games,” _arXiv preprint arXiv:1912.11077_ , 2019.
  * [10] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, “Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 1861–1870.
  * [11] A. Bicchi and V. Kumar, “Robotic grasping and contact: a review,” _Proceedings 2000 ICRA. Millennium Conference. IEEE International Conference on Robotics and Automation. Symposia Proceedings (Cat. No.00CH37065)_ , vol. 1, pp. 348–353 vol.1, 2000.
  * [12] J. Bohg, A. Morales, T. Asfour, and D. Kragic, “Data-driven grasp synthesis—a survey,” _IEEE Transactions on Robotics_ , vol. 30, no. 2, pp. 289–309, 2014.
  * [13] I. Lenz, H. Lee, and A. Saxena, “Deep learning for detecting robotic grasps,” in _1st International Conference on Learning Representations, ICLR 2013, Scottsdale, Arizona, USA, May 2-4, 2013, Workshop Track Proceedings_ , Y. Bengio and Y. LeCun, Eds., 2013. [Online]. Available: http://arxiv.org/abs/1301.3592
  * [14] J. Redmon and A. Angelova, “Real-time grasp detection using convolutional neural networks,” 2015.
  * [15] S. Kumra and C. Kanan, “Robotic grasp detection using deep convolutional neural networks,” 2017.
  * [16] D. Morrison, P. Corke, and J. Leitner, “Learning robust, real-time, reactive robotic grasping,” _The International Journal of Robotics Research_ , vol. 39, no. 2-3, pp. 183–201, 2020. [Online]. Available: https://doi.org/10.1177/0278364919859066
  * [17] S. Kumra, S. Joshi, and F. Sahin, “Antipodal robotic grasping using generative residual convolutional neural network,” _CoRR_ , vol. abs/1909.04810, 2019\. [Online]. Available: http://arxiv.org/abs/1909.04810
  * [18] J. Mahler, J. Liang, S. Niyaz, M. Laskey, R. Doan, X. Liu, J. A. Ojea, and K. Goldberg, “Dex-net 2.0: Deep learning to plan robust grasps with synthetic point clouds and analytic grasp metrics,” 2017.
  * [19] J. Mahler, F. T. Pokorny, B. Hou, M. Roderick, M. Laskey, M. Aubry, K. Kohlhoff, T. Kröger, J. Kuffner, and K. Goldberg, “Dex-net 1.0: A cloud-based network of 3d objects for robust grasp planning using a multi-armed bandit model with correlated rewards,” in _IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2016, pp. 1957–1964.
  * [20] P. Schmidt, N. Vahrenkamp, M. Wächter, and T. Asfour, “Grasping of unknown objects using deep convolutional neural networks based on depth images,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_ , 2018, pp. 6831–6838.
  * [21] A. Mousavian, C. Eppner, and D. Fox, “6-dof graspnet: Variational grasp generation for object manipulation,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 2901–2910.
  * [22] H. Fang, C. Wang, M. Gou, and C. Lu, “Graspnet-1billion: A large-scale benchmark for general object grasping,” in _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020_. IEEE, 2020, pp. 11 441–11 450.
  * [23] M. Sundermeyer, A. Mousavian, R. Triebel, and D. Fox, “Contact-graspnet: Efficient 6-dof grasp generation in cluttered scenes,” _arXiv preprint arXiv:2103.14127_ , 2021.
  * [24] L. Pinto and A. Gupta, “Supersizing self-supervision: Learning to grasp from 50k tries and 700 robot hours,” _CoRR_ , vol. abs/1509.06825, 2015. [Online]. Available: http://arxiv.org/abs/1509.06825
  * [25] S. Levine, P. Pastor, A. Krizhevsky, and D. Quillen, “Learning hand-eye coordination for robotic grasping with deep learning and large-scale data collection,” _CoRR_ , vol. abs/1603.02199, 2016. [Online]. Available: http://arxiv.org/abs/1603.02199
  * [26] D. Kalashnikov, A. Irpan, P. Pastor, J. Ibarz, A. Herzog, E. Jang, D. Quillen, E. Holly, M. Kalakrishnan, V. Vanhoucke, and S. Levine, “Qt-opt: Scalable deep reinforcement learning for vision-based robotic manipulation,” _CoRR_ , vol. abs/1806.10293, 2018. [Online]. Available: http://arxiv.org/abs/1806.10293
  * [27] M. Khansari, D. Kappler, J. Luo, J. Bingham, and M. Kalakrishnan, “Action image representation: Learning scalable deep grasping policies with zero real world data,” 2020.
  * [28] S. Song, A. Zeng, J. Lee, and T. Funkhouser, “Grasping in the wild: Learning 6dof closed-loop grasping from low-cost demonstrations,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4978–4985, 2020.
  * [29] A. Zeng, S. Song, S. Welker, J. Lee, A. Rodriguez, and T. A. Funkhouser, “Learning synergies between pushing and grasping with self-supervised deep reinforcement learning,” _CoRR_ , vol. abs/1803.09956, 2018. [Online]. Available: http://arxiv.org/abs/1803.09956
  * [30] M. Dogar and S. Srinivasa, “A planning framework for non-prehensile manipulation under clutter and uncertainty,” _Autonomous Robots_ , vol. 33, no. 3, pp. 217–236, June 2012.
  * [31] A. Boularias, J. A. Bagnell, and A. Stentz, “Learning to manipulate unknown objects in clutter by reinforcement,” in _AAAI_ , ser. AAAI’15. AAAI Press, 2015, p. 1336–1342.
  * [32] L. Pinto and A. Gupta, “Learning to push by grasping: Using multiple tasks for effective learning,” _CoRR_ , vol. abs/1609.09025, 2016. [Online]. Available: http://arxiv.org/abs/1609.09025
  * [33] M. Danielczuk, A. Kurenkov, A. Balakrishna, M. Matl, D. Wang, R. Martín-Martín, A. Garg, S. Savarese, and K. Goldberg, “Mechanical search: Multi-step retrieval of a target object occluded by clutter,” _CoRR_ , vol. abs/1903.01588, 2019. [Online]. Available: http://arxiv.org/abs/1903.01588
  * [34] L. Berscheid, P. Meißner, and T. Kröger, “Robot learning of shifting objects for grasping in cluttered environments,” _CoRR_ , vol. abs/1907.11035, 2019. [Online]. Available: http://arxiv.org/abs/1907.11035
  * [35] L. Berscheid, C. Friedrich, and T. Kröger, “Robot learning of 6 dof grasping using model-based adaptive primitives,” _CoRR_ , vol. abs/2103.12810, 2021. [Online]. Available: https://arxiv.org/abs/2103.12810
  * [36] T. Haarnoja, A. Zhou, K. Hartikainen, G. Tucker, S. Ha, J. Tan, V. Kumar, H. Zhu, A. Gupta, P. Abbeel, and S. Levine, “Soft actor-critic algorithms and applications,” _CoRR_ , vol. abs/1812.05905, 2018. [Online]. Available: http://arxiv.org/abs/1812.05905
  * [37] I. Kostrikov, D. Yarats, and R. Fergus, “Image augmentation is all you need: Regularizing deep reinforcement learning from pixels,” _arXiv preprint arXiv:2004.13649_ , 2020.
  * [38] S. Fujimoto, H. van Hoof, and D. Meger, “Addressing function approximation error in actor-critic methods,” in _ICML_ , ser. Proceedings of Machine Learning Research, J. G. Dy and A. Krause, Eds., vol. 80. PMLR, 2018, pp. 1582–1591.
  * [39] M. L. Puterman, _Markov Decision Processes_. Wiley and Sons, 1994.
  * [40] R. Finman, T. Whelan, M. Kaess, and J. J. Leonard, “Toward lifelong object segmentation from change detection in dense rgb-d maps,” in _Proceedings of European Conference on Mobile Robots (ECMR ’13)_ , September 2013, pp. 178 – 185.
  * [41] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ International Conference on Intelligent Robots and Systems_. IEEE, 2012, pp. 5026–5033.
  * [42] L. Berscheid, P. Meißner, and T. Kröger, “Robot learning of shifting objects for grasping in cluttered environments,” _arXiv preprint arXiv:1907.11035_ , 2019.


# Constraints on Cosmographic Functions of Cosmic Chronometers Data Using
Gaussian Processes

A. M. Velasquez-Toribio<EMAIL_ADDRESS>Nucleo Cosmo-ufes & Departamento
de Física, Universidade Federal do Espirito Santo, 29075-910 Vitǿria - ES,
Brasil  Júlio C. Fabris<EMAIL_ADDRESS>Nucleo Cosmo-ufes &
Departamento de Física, Universidade Federal do Espirito Santo, 29075-910
Vitǿria - ES, Brasil  National Research Nuclear University MEPhI, Kashirskoe
sh. 31, Moscow 115409, Russia

###### Abstract

We study observational constraints on the cosmographic functions up to the
fourth derivative of the scale factor with respect to cosmic time, i.e., the
so-called snap function, using the non-parametric method of Gaussian
Processes. As observational data we use the Hubble parameter data. Also we use
mock data sets to estimate the future forecast and study the performance of
this type of data to constrain cosmographic functions. The combination between
a non-parametric method and the Hubble parameter data is investigated as a
strategy to reconstruct cosmographic functions. In addition, our results are
quite general because they are not restricted to a specific type of functional
dependency of the Hubble parameter. We investigate some advantages of using
cosmographic functions instead of cosmographic series, since the former are
general definitions free of approximations. In general, our results do not
deviate significantly from $\Lambda CDM$. We determine a transition redshift
$z_{tr}=0.637^{+0.165}_{-0.175}$ and $H_{0}=69.45\pm 4.34$. Also assuming
priors for the Hubble constant we obtain $z_{tr}=0.670^{+0.210}_{-0.120}$ with
$H_{0}=67.44$ (Planck) and $z_{tr}=0.710^{+0.159}_{-0.111}$ with
$H_{0}=74.03$(SH0ES). Our main results are summarized in table 2.

###### pacs:

98.80.-k, 95.36.+x, 98.80.Es

##  I Introduction

The accelerated expansion of the Universe Riess et al. 1998 ; Perlmutter et
al. 1999  is currently one of the main problems of cosmology. The simplest
explantion is given by the $\Lambda CDM$ model making use of cold dark matter
and a cosmological constant which matches well with different observational
data types: supernovae Ia, cosmic microwave background (CMB), baryons acoustic
ocillations (BAO), etc, see e.g. Scolnic et al. 2018 ; Planck Collaboration et
al. (2018) ; Alam et al. (2016) .

The $\Lambda CDM$ model includes a component of dark energy which is
responsible for the accelerating expansion of the universe. However, from a
theoretical point of view, one of the big problems that affect this model is
the difference between the measured value of the cosmological constant and
that required by quantum theory (depending on the energy scale), which can
reach a part in $10^{120}$. This gives rise to the cosmological constant fine
tuning problem. Other theoretical problem is the cosmological coincidence
problem, i.e., why does the universe start to accelerate when the structures
reach the non-linear regime? Or why the present values of the densities of
dark energy and dark matter are of the same order of magnitude? For a review
of these and another theoretical problems on the $\Lambda CDM$ model, see
Peebles and Ratra (2003) , Perivolaropoulos (2008) and Bullock and Boylan-
Kolchin, (2017).

On the other hand, there are many theoretical models that try to explain the
accelerating expansion of the Universe, for example, among the most popular we
have: the scalar field models (quintessence) Horndeski 1974 , k-essence models
Horndeski 1974 , models including dissipationHorndeski 1974 , models inspired
by the renormalization group Horndeski 1974 , modified gravity models and the
Horndeski’s model that have the most general scalar-tensor action yielding up
to second-order equations of motion (Horndeski 1974 ), and so forth. For a
review see Weinberg, et al. (2013) .

In addition, another way to understand the accelerated expansion is the
phenomenological approach, which consists in assuming a specific functional
form for the equation of state and studying its implications. The prototype of
this approach is the parameterization $CPL$ Chevallier and Polarski 2001 ;
Linder 2003 . Other examples are: the $q(z)$ parametrization Ishida et al.
2007 , the Wetterich parameterizations Wetterich, 2004 , among others e.g.
Copeland et al. 2006 .

However, all the mentioned models to explain the accelerated expansion are
dependent on the theory of gravitation. In this regard there is an approach
that avoids it by assuming only the metric structure, without considering the
theory of gravitation and is called the model-independent approach. We assumed
this approach in our investigation. It is also common to call this approach as
the cosmographic approach and derive it from a Taylor series of the scale
factor or, equivalently, of the Hubble parameter e.g. see references Weinberg
(1972) ; Chiba and Nakamura 1998 ; Visser (2003) ; Visser 2004 ; Caldwell 2004
.

It is worth noting that the Taylor series expansion implies fundamental
difficulties with the convergence and the truncation of the series. Thus we
can not use observational data for redshift greater than one, $z>1$. For a
detailed discussion of this question, see the reference Cattone and Visser
(2007) . Consequently, in the present paper, to avoid this problem we
calculate the cosmographic functions instead of the cosmographic parameters
derived from the Taylor series. The cosmographic functions are assumed as
kinematic definitions on the scale factor or, equivalently, the Hubble
parameter. In this case the cosmographic parameters correspond to the
cosmographic functions with $z=0$, i.e, currently evaluated.

In addition, a statistical method complementary to the model-independent
approach is the method of Gaussian Processes($GP$), which is a non-parametric
statistical method. This means that the reconstruction of the Hubble function
and its derivatives are performed without assuming a specific theoretical
model. Also as observational data we used the Hubble parameter data, since
from the theoretical point of view the cosmographic functions are directly
dependent on the Hubble parameter and its derivatives. On the other hand, a
reconstruction that uses, for example, measurements of the luminosity distance
of supernovae Ia involves the computation of the derivative of the supernovae
Ia distance data to determine the Hubble parameter and its derivatives. Thus
this process adds a source of error propagation. This is the main reason
because we only use Hubble parameter data in our research.

Our main objective is determine observational constraints in model-independent
form using a non-parametric statistical technique, such as Gaussian processes,
and observational data from the Hubble parameter. In this way we avoid using
the Taylor series expansion. We also discuss how to use these results to
derive observational constraints on specific models by constructing a simple
correspondence between the parameters of a given model and the cosmographic
parameters. We discuss the implications of this methodology.

The present paper is organized as follows. We present the basic definitions in
the section 2. The observational data and the simulated data are discussed in
the section 3. We summarize the method of Gaussian Processes in section 4. We
present the results of the reconstruction of our cosmographic functions in
section 5 and in section 6 we give our conclusions.

##  II Cosmographic Functions

We assume that our universe is homogenous and isotropic and therefore can be
described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric,

$\displaystyle
ds^{2}=-c^{2}dt^{2}+a(t)^{2}\left[\frac{dr^{2}}{1-kr^{2}}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right)\right],$ (1)

where $a(t)$ is the scale factor and the $k$ is a constant that can have the
value $=+1,0,-1$ and represent the three-dimensional curvature. Using this
metric it is possible to introduce the cosmographic functions defined directly
from the scale factor,

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{1}{a}\frac{da(t)}{dt},$ $\displaystyle q(t)$
$\displaystyle=$
$\displaystyle-\frac{1}{a}\frac{d^{2}a(t)}{dt^{2}}\frac{1}{H^{2}},$
$\displaystyle j(t)$ $\displaystyle=$
$\displaystyle\frac{1}{a}\frac{d^{3}a(t)}{dt^{3}}\frac{1}{H^{3}},$ (2)
$\displaystyle s(t)$ $\displaystyle=$
$\displaystyle\frac{1}{a}\frac{d^{4}a(t)}{dt^{4}}\frac{1}{H^{4}},$
$\displaystyle...$ $\displaystyle=$ $\displaystyle...$ $\displaystyle m_{n}$
$\displaystyle=$
$\displaystyle\frac{1}{a}\frac{d^{n}a(t)}{dt^{n}}\frac{1}{H^{n}}.$

If we evaluate these functions at the current time ($z=0$), then they are
called the cosmographic parameters and are usually referred as the Hubble
constant ($H_{0}$), the deceleration parameter ($q_{0}$), the jerk parameter
($j_{0}$), the snap parameter ($s_{0}$), etc.111Curiously in the literature
there are proposals of names for the fifth and sixth derivatives as crackle
and pop (Dunajski and Gibbons 2008 ).

In the literature the traditional way to introduce cosmographic parameters is
by using the Taylor expansion for the Hubble parameter see, for example,
Weinberg (1972)  and Visser (2003) . However, the use of observational data
and the Taylor series leads to two main inconsistencies. The first arises when
we use observational data with $z>1$, since from the mathematical point of
view the series has a convergence radius at the most $|z|=1$, thus it is only
justified to use data with redshift $z<<1$. Otherwise the errors would
increase and the series would not be useful to represent the original
function. The other inconsistency arises when a new redshift variable is
introduced to transform the domain of redshift values and to be able to use
data with $z>1$, the problem is that there is no one-way to do this.
Furthermore, a given redshift parametrization contributes to error
propagation. In what follows we briefly review these two inconsistencies.

  * •

Taylor series and their error in cosmological models: the Hubble’s function
can be expressed in Taylor series as a function of the redshift using the
relation $1+z=\frac{a_{0}}{a(t)}$ Weinberg (1972) ,

$\displaystyle
H(z)=H_{0}+\frac{dH}{dz}\Bigg{|}_{0}z+\frac{1}{2!}\frac{d^{2}H}{dz^{2}}\Bigg{|}_{0}z^{2}+\frac{1}{3!}\frac{d^{3}H}{dz^{3}}\Bigg{|}_{0}z^{3}+O(z^{4}).$
(3)

The series is written around the current redshift, $z=0$, and mathematically
represents the approximate Hubble’s law for an arbitrary order. 222On the
generality of the equation above it is important to mention that for
gravitation theories of high derivative one can not reconstruct his equation
of state without making use of an external hypothesis, because these theories
include new degrees of freedom. For details see references de la Cruz-Dombriz,
et al. (2016) and Busti, et al., (2015). In this context, we consider the
truncation error. To briefly illustrate this question we have used three
cosmological models with Hubble’s law exact, so that we can easily determine
its Taylor series. We consider the flat $\Lambda CDM$ model, the $CPL$ model
and a $q(z)$ parametrization. In all case we consider the flat universe.
Thus, for the case of the flat cosmological constant model the Hubble
parameter can be written,

$\displaystyle H(z)$ $\displaystyle=$ $\displaystyle
H_{0}\left[\Omega_{m0}(1+z)^{3}+(1-\Omega_{m0})\right]^{1/2},$ (4)

where $\Omega_{m0}$ is the the matter density parameter.
On the other hand, for the CPL (Chevallier-Polarski-Linder) parametrization
Chevallier and Polarski 2001  and Linder 2003  we have,

$\displaystyle
H(z)=H_{0}[\Omega_{m0}(1+z)^{3}+(1-\Omega_{m0})(1+z)^{3(1+w_{0}+w_{1})}e^{3w_{1}(\frac{1}{1+z}-1)}]^{1/2},$
(5)

where $w_{0}$ and $w_{1}$ are the parameters of the model and although it does
not behave very well for high redshift it allows reproduce the Hubble law and
the distances of different models of dark energy with very good precision and
in the literature it has been intensively used, see for example Padmanabhan
and Choudhury (2003)  and Scherrer. (2015)  and references therein. This model
represents dynamic dark energy.

We are also consider another model of dynamic dark energy with four free
parameters Ishida et al. 2007  for which have

$\displaystyle
H(z)=H_{0}\left[(1+z)^{(1+q_{i})}\left(\frac{q_{i}\left(\frac{1+z_{t}}{1+z}\right)^{1/\tau}-q_{i}}{q_{i}(1+z_{t})^{2(1+q_{i})}-q_{i}}\right)^{\tau(q_{i}-q_{f})}\right],$
(6)

where $q_{i}$ and $q_{f}$ are the initial and final deceleration parameter
respectively and $z_{t}$ is the transition redshift from a decelerated
universe to an accelerated universe and $\tau$ is the width of the transition.
This model is quite generic and allows to study in greater detail the
transition for an accelerated universe. It is also interesting to know that
this parametrization includes several other models of dark energy, for details
and observational constraints of the parameters see Dos santos et al. (2016) .

Figure 1: In the top left the Hubble function for $\Lambda CDM$, in the middle
left for the CPL model and in the botton left for the $q(z)$ model and in the
right side the quantification of the error respectively for each model. In all
case we use $\Omega_{m0}=0.32$ and we consider two values for the today Hubble
parameter: the highest value corresponds to the SH0ES Collaboration and the
lowest value corresponds to the Planck Collaboration.

All these models have an analytical expression for the Hubble parameter and
have been constrained observationally. We calculate the Taylor series for each
of these models up to the third power in the redshift variable.

In Figure 1, we show the exact Hubble function versus the third order
approximation. We can observe that up around $z=0.5$ the values between the
exact function and the third order approximation are more or less equivalent.
To quantify explicitly this difference we plot the magnitude of the error
defined as:

$\displaystyle Error$ $\displaystyle=$
$\displaystyle\frac{H(z)_{exact}-H(z)_{approximate}}{H(z)_{exact}}.$ (7)

Also in figure 1, right column, we can see that the error between the exact
function and the third-order function for the range of redshift, $0<z<1$, is
relatively small, but for higher values of the redshift the error increases
considerably. Therefore, the use of the Taylor series to compute observational
constraints for data with $z>1$ leads to considerable truncation errors. One
way to circumvent this problem could be to select observational data with
$z<1$, this would give consistent results within the truncation error of the
Taylor series.

  * •

New redshift variable: a widely used proposal in the literature on the
truncation problem is to use a new variable for the redshift. The initial and
best motivated proposal is given by the reference Cattone and Visser (2007)
and establishes the following relationship

$\displaystyle y$ $\displaystyle=$ $\displaystyle\frac{z}{1+z},$ (8)

where $y$ is the new redshift variable. This new redshift variable allows
mapping the redshift values with $z>1$ to a region with $y<1$ and can be
defined as the change of the emission wavelength divided by the observed
wavelength,

$\displaystyle 1-y$ $\displaystyle=$
$\displaystyle\frac{\lambda_{e}}{\lambda_{0}}=\frac{1}{z}.$ (9)

The correspondence between the values of redshift $z$ and the values of the
variable $y$ are: in the past $z\in(0,\infty)$ and corresponding the values
$y\in(0,1)$; in the future $z\in(-1,0)$ and corresponding to the values
$y(\infty,0)$.

However, this redshift variable transform is not unique in the literature
there are other proposals such as Aviles, et.al; 2012 ,

$\displaystyle y_{1}$ $\displaystyle=$ $\displaystyle\arctan(z),$ (10)
$\displaystyle y_{2}$ $\displaystyle=$ $\displaystyle\arctan(\frac{z}{z+1}),$
(11) $\displaystyle y_{3}$ $\displaystyle=$ $\displaystyle\frac{z}{1+z^{2}}.$
(12)

In this sense, the cosmological constraints depend on the parametrization
assumed for the redshift. When comparing the results using the variable $y$
with those obtained using the original redshift variable inconsistencies
arise. It is quite instructive to show this explicitly. For example, we can
consider the expressions for the luminosity distances in the variables $y$ and
$z$ respectively Weinberg (1972)  and Cattone and Visser (2007) ,

$\displaystyle d_{L}(y)$ $\displaystyle=$
$\displaystyle\frac{c}{H_{0}}y[1-\frac{1}{2}(-3+q_{0})y+\frac{1}{6}(12-5q_{0}+3q_{0}^{2}+6(j_{0}+\Omega_{0}))y^{2}+O(y^{3})],$
(13) $\displaystyle d_{L}(z)$ $\displaystyle=$
$\displaystyle\frac{c}{H_{0}}z[1+\frac{1}{2}(1-q_{0})z+\frac{1}{6}(1-q_{0}-3q_{0}^{2}+j_{0})z^{2}+O(z^{3})].$
(14)

Now we can use these two expressions to determine observational constraints on
the parameter space $(q_{0}$, $j_{0})$. To determine these constraints we can
use a robust data set from the point of view of error control, particularly
systematic errors. We use data from type Ia Supernovae, specifically the
Pantheon sample with 1048 supernovae Scolnic et al. 2018 . In appendix A we
briefly describe the method used to estimate parameters using these data.

In figure 2 we show the results of the confidence contours for the parameters
space $(q_{0},j_{0})$ using the type Ia Supernovae data. On the left side we
present the results for the original redshift, $z$, with the best fit
$q_{0}=-0.579\pm 0.103$ and $j_{0}=0.977\pm 0.232$ with $1\sigma$. On the
right side we show the results for the redshift, $y$, with the best fit
$q_{0}=-1.107\pm 0.103$ and $j_{0}=0.650\pm 0.232$ with $1\sigma$. As we can
see, the results are incompatible since that the confidence contours do not
overlap.333Similar results are obtained in a recent paper by Capozziello et
al. 2020  also using type Ia Supernovae.

Therefore, to avoid the problems mentioned about the Taylor series we propose
to use directly the definitions of cosmographic functions as given above.
Thus, for our analysis it is convenient to rewrite the cosmographic functions
explicitly as a function of the Hubble parameter and its derivatives with
respect to the redshift,

$\displaystyle q(z)$ $\displaystyle=(1+z)\frac{H^{\prime}}{H}-1$ (15)
$\displaystyle j(z)$
$\displaystyle=\frac{H^{\prime\prime}}{H}(1+z)^{2}+\left(\frac{H^{\prime}}{H}\right)^{2}(1+z)^{2}-2\frac{H^{\prime}}{H}(1+z)+1$
(16) $\displaystyle s(z)$
$\displaystyle=-\frac{1}{2}(1+z)^{2}\frac{[H^{2}(z)]^{\prime\prime\prime}}{H^{2}(z)}-\frac{1}{2}(1+z)^{2}\frac{[H^{2}(z)]^{\prime\prime}H^{\prime}(z)}{H^{3}(z)}$
(17)
$\displaystyle+\frac{(1+z)^{2}[H^{2}(z)]^{\prime\prime}}{H^{2}}+\frac{(1+z)^{2}[H^{2}(z)]^{\prime}H^{\prime}(z)}{H^{3}(z)}+\frac{3(1+z)H^{\prime}(z)}{H(z)}+1,$

where the prime represents derivative with respect to redshift. Furthermore,
using the snap parameter definition we have $\ddddot{a}=saH^{4}$, from the
jerk parameter we have $\dddot{a}=jaH^{3}$ and of the deceleration parameter
we have $\ddot{q}=-qaH^{2}$, thus we can write a useful relationship between
these cosmographic functions given by the expression jerk ,

$\displaystyle s(z)$ $\displaystyle=$
$\displaystyle-(1+z)\frac{dj(z)}{dz}-j(z)\left(2+3q(z)\right).$ (18)

It is also useful to write the cosmographic parameters for the flat $\Lambda
CDM$ model,

$\displaystyle q_{0}$ $\displaystyle=$
$\displaystyle-1+\frac{3}{2}\Omega_{m0}$ (19) $\displaystyle j_{0}$
$\displaystyle=$ $\displaystyle 1$ (20) $\displaystyle s_{0}$ $\displaystyle=$
$\displaystyle 1-\frac{9}{2}\Omega_{m0}.$ (21)

As we can see for the flat $\Lambda CDM$ model these cosmographic parameters
are defined by the value of the matter parameter today. These results are
useful to make a comparison with our results. We use the set of equations (16)
- (18) to reconstruct the cosmographic functions, i.e, $q(z)$, $j(z)$ and
$s(z)$ using observational data from the Hubble parameter and simulated data.

Figure 2: In the left we shown the observational constraints on the parameters
$q_{0}$ and $j_{o}$ using the redshift (z) and in the right shows the
observational constraints using the $y$-redshift. We can see that both results
are incompatible. The data used to determine these observational constraints
have been data from Type Ia supernovae, for details see appendix A.

##  III Data

As we have observed in the previous section the cosmographic functions can be
written directly as a function of the Hubble parameter and its derivatives.
This allows us to discern that the most appropriate data to minimize error
propagation are the Hubble data. Therefore, we focus on this class of data.

###  III.1 Observational Data

There are two efficient and widely used forms to obtain Hubble parameter
measurements.

####  III.1.1 Cosmic Chronometers(CC)

This method is based on the expression of the differential age of the universe
as a function of redshift,

$\displaystyle H(z)$ $\displaystyle=$
$\displaystyle-\frac{1}{1+z}\frac{dz}{dt}.$ (22)

This method was proposed by Jimenez and Loeb (2002) and consists in directly
measuring the amount $dz/dt$ and, consequently, the Hubble parameter. The most
used data to measure this amount have been passively evolving galaxies with
high-resolution spectroscopic data along with synthetic catalogs to limit the
age of the oldest stars in the galaxy. A complete description of this
methodology can be reviewed for the $SDSS$ and $GDS$ data in the reference
Jimenez et al. (2002).

####  III.1.2 The Radial BAO Size Method

The method is based on measurements of the acoustic scale of $BAO$, which is
more accurate with respect at the $CC$ method Gaztanaga et al. 2009 . This
accuracy is understandable because $BAO$ mainly depends on a spatial
measurement compared to the first method where a time measurement is required
which increases systematic errors. However, this method of $BAO$ requires
assuming a prior in the radius of the sound horizon. Hence, we have

$H(z)=-\frac{r_{bao}}{r_{cmb}}H_{fiducial}.$ (23)

This method depends on the fiducial model which is assumed to be the flat
$\Lambda CDM$ model.

On the other hand, in the literature there are different compilations of
samples of the Hubble parameters data. We use the sample presented by Zheng et
al. (2016)  what includes data of both $CC$ and $BAO$, but we not including
the points $z=0.4497$ and $z=2.34$. The first point is excluded because it
overlaps with another data and has a negligible effect on the results. The
other point has a strong influence on the estimation of parameters and because
it is very restrictive we have excluded it. Therefore our sample has 34
points. This sample has been used in Velásquez-Toribio and Magnago (2020)
showing results compatible with the literature.

In figure 3 we show the data used. We can see that the errors of the $CC$ data
can reach up to a value of $\Delta H\approx 62Km\>s^{-1}Mpc^{-1}$ and the
errors of the BAO data have a mean of about $7Km\>s^{-1}Mpc^{-1}$.

Before using these data it is important to verify their Gaussian nature
Podariu et al. 2001 . To verify this we determine the factor $N$ which must be
small to be consistent with the gaussianity hypothesis. Briefly the
methodology used is shown in Appendix B. The results are presented in Table 1
and we can see that the data are compatible with the gaussianity hypothesis.
Therefore, we have a suitable sample to use the Gaussian process method.

###  III.2 Simulation Data

Our simulation follows the prescriptions established in Ma and Zhang (2011) :
we choose as background cosmology the flat $\Lambda CDM$ model using the
values $\Omega_{m0}=0.315\pm 0.001$ of the Planck Collaboration et al. (2018)
and $H_{0}=69.45\pm 4.34$( best fit of $H(z)$ observational data). Then, with
these values we determine the $H_{fiducial}$ using the equation (5). We assume
that the deviation of the simulated value with respect to fiducial value is
calculated as: $\Delta H=H_{sim}-H_{fiducial}$, where the $\Delta H$ can be
derived from a Gaussian distribution of the form $N(0,\tilde{\sigma})$, where
the $\tilde{\sigma}$ is a random number which also can be drawn from the
following Gaussian distribution $N\left(\sigma_{0},\varepsilon(z)\right)$. In
this Gaussian distribution we introduce the information about the error as
follows: we define the $\varepsilon(z)$ as
$\varepsilon(z)=\frac{\sigma_{+}(z)-\sigma_{-}(z)}{4}$ and the parameter
$\varepsilon(z)$ is chosen so that the probability of $\tilde{\sigma}$ falling
within the strip of $95.4\%$. We also define the
$\sigma_{0}=(\sigma_{+}+\sigma_{-})/2$. We assume a linear function for the
$\sigma$ that allows us to model the simulated errors. This can be deduced
from figure 3 where we can see that, except for some special cases, the error
increases with the redshift. Thus, in general we define
$\sigma_{\pm}=a_{\pm}z+b_{\pm}$ where $a$ and $b$ are chosen according to the
errors of the Hubble data and the $+$ sign indicates the upper line and the
$-$ sign the lower line as shown in the lower part of figure 3. For example,
for the case of $CC$ we assume the straight lines (orange in the figura 3):
$\sigma_{+}=19.999z+11.105$ and $\sigma_{-}=3.650z+1.65$ and in the case of
$BAO$ the blue lines.

Table 1: Gaussian nature of the data $N$ | $\bar{z}$ | $\bar{H}$ | $\bar{\sigma}$ | $N_{\sigma}$
---|---|---|---|---
2 | $0.084$ | $69$ | $10.2342$ | $1.41$
5 | $0.184$ | $75.939$ | $2.878$ | $1.443$
4 | $0.332$ | $81.785$ | $6.309$ | $1.861$
5 | $0.417$ | $83.232$ | $4.824$ | $1.429$
4 | $0.516$ | $88.687$ | $5.341$ | $0.926$
3 | $0.624$ | $91.191$ | $4.544$ | $0.624$
3 | $0.759$ | $102.146$ | $5.696$ | $0.164$
3 | $0.965$ | $132.079$ | $14.120$ | $0.34$
5 | $1.445$ | $160.331$ | $8.717$ | $0.088$
2 | $2.333$ | $221.328$ | $6.933$ | $0.427$

Figure 3: We shown the Hubble data set (CC+BAO) used in this investigation and
also the lines used to calculate the error of simulated data. The orange lines
represent the lines for the CC data and the blue lines represent the lines for
the BAO data.

##  IV Gaussian Processes

The gaussian processes ($GP$) are a nonparametric statistical method allowing
the reconstruction of a given function directly from the manipulation of data.
They can be thought as a extension to a multivariate normal distribution, that
is, a $GP$ is a infinity collection of random variables, but where every
finite set have a normal distribution. A $GP$ is a distribution over
functions, but in the practice a $GP$ allows us to derive a posterior
distributions by simply considering a finite set of points that are associated
with the observational data. For details see Rasmussen and Williams (2006) . A
$GP$ can be written as:

$\displaystyle f(x)\sim GP\left(\mu(x),k(x,\tilde{x})\right),$ (24)

where the value of $f$ when evaluated at a point $x$ is a gaussian random
variable with mean $\mu(x)$. Additionally, the value of the function $f$ at
the point $x$ is not independent of the value of the function $f$ at some
other point nearby $\tilde{x}$, but is related by the covariance function
$k(x,\tilde{x})$. Considering observational data $(x_{i},y_{i})$, assuming
that the errors are gaussian and that $y_{i}=f(x_{i})\pm\sigma_{i}$, where
$\sigma_{i}$ are the $1\sigma$ error bars and $i=1,...,N$, we can to
reconstruct the function $f$ at chosen points. This function is denoted by
$f^{*}$. This reconstruction can be done through the joint distribution
between $f^{*}$ and $y_{i}$. For details of the exact expressions see the
reference Seikel et al. (2012a). In general, this reconstructed function has a
reconstructed mean given by,

$\displaystyle\bar{f}^{*}(x)=\sum_{i,j=1}^{N}{k(x,x_{i})(M^{-1}})_{i,j}(f(x_{j})-\mu(z))$
(25)

in our reconstruction we chosen as prior mean function $\mu(z)$ a constant
value and for the standard deviation:

$\displaystyle\sigma(x)=k(x,x)-\sum_{i,j=1}^{N}{k(x,x_{i})(M^{-1}})_{i,j}k(x_{j},x)$
(26)

where $M_{ij}=k(x_{i},x_{j})+c_{ij}$ and $c_{ij}$ represents the covariance of
the input data.

For our calculations, we assume the exponential function as a covariance
function which is given by

$\displaystyle
k(x,\tilde{x})=\sigma_{f}^{2}\exp\left(-\frac{(x-\tilde{x})^{2}}{2\ell^{2}}\right),$
(27)

where $\sigma_{f}$ and $\ell$ are called hyperparameters which are determined
by maximizing the log marginal likelihood, see Rasmussen and Williams (2006) ,

$\displaystyle\ln\mathcal{L}=-\frac{1}{2}\sum{\left[f(x_{i})-\mu(x_{i})\right](M^{-1})_{i,j}\left[f(x_{j})-\mu(x_{j})\right]}-\frac{1}{2}\ln|M|-\frac{1}{2}N\ln
2\pi,$ (28)

where $M$ represents the determinant of $M_{ij}$. It is important to note that
instead of optimizing these hyperparameters, they can be marginalized by
using, for example, the MCMC method. However, a simple method of testing the
validity of the optimization process is to vary the initial value of the
hyperparameters and see if the values obtained change significantly. We have
done this process and we determine that our hyperparameters do not change
significantly. Therefore, we are confident in using the optimization process.
On the other hand, the $MCMC$ method should be used to obtain high precision
restrictions associated with high quality data.

The question of how to choose a suitable covariance function is an important
problem for GP. In the present investigation we show some evidence to consider
the Gaussian nature of the data, see table I. However, we consider that this
evidence is not conclusive. Some papers as Rasmussen and Williams (2006) ;
Rasmussen and Williams (2006)  have used other covariance functions to
estimate $H_{0}$ and have shown that within $1\sigma$ different choices are
consistent. Therefore, in the present paper we only use the squared
exponential covariance function and we leave the specific study of the effects
of choosing different covariance functions for further investigation.

Additionally, the GP allows to reconstruct the derivative of the data. To
implement this method, we use the public package $GaPP$ Seikel et al. (2012a).
For applications of GP in cosmology consider references Seikel et al. 2012b
and Velásquez-Toribio et al. 2019  and for another GP methodology to see
Holsclaw et al. (2010b) ; 2010a ; Holsclaw et al. (2011) .

##  V Results

In figure 4 we present the reconstruction of the Hubble parameter using the
$GaPP$ code. In the top we present the model-independent reconstruction of
$H(z)$, where we have used the $CC+BAO$ data to determine the Hubble constant
($H_{0}=69.45\pm 4.34$). In the other two figures we reconstructed $H(z)$
using as prior for the Hubble constant the values of the Planck and SH0ES
Collaborations. In all cases, the reconstruction is done with $1\sigma$ of
uncertainty. In figure 5 we can observe the reconstruction of the cosmographic
functions $q(z)$, $j(z)$ and $s(z)$ without using a prior value for the Hubble
constant. In particular the deceleration parameter shows that the transition
redshift, $z_{tr}$, from a decelerated universe to an accelerated universe,
defined as $q(z_{tr})=0$, is into the region $z<1$ with $2\sigma$.

We obtain the values: $z_{tr}=0.637^{+0.165}_{-0.175}$ (model-independent) and
$z_{tr}=0.670^{+0.210}_{-0.120}$ with $H_{0}=67.44$ (Planck) and
$z_{tr}=0.710^{+0.159}_{-0.111}$ with $H_{0}=74.03$(SH0ES). In accordance with
these results we can conclude that $z_{tr}\in\left[0.550,0.870\right]$. These
values for, $z_{tr}$, are in accordance with other results determined in the
literature, for example, see Farooq and Ratra (2013a) ; farooq et al. (2013b)
; Farooq et al. (2017) ; Yu et al. (2018) ; Mamon (2018) ; Velásquez-Toribio
et al. (2011) .

On the other hand, for values $z>1$ the function $q(z)$ shows a trend to
remain in the region $q>0$ associated with a phase dominated by matter.
However, for high redshift values we see a abrupt drop which can be associated
with the few observational data for this region (see figure 3).

In figure 6 we can see that the reconstructed function $j(z)$ does not exclude
the flat $\Lambda CDM$ model with $2\sigma$, however the average value is
different from the flat $\Lambda CDM$ model for all the reconstruction. The
figure shows a trend for $j<0$ when $z>1$, this is due to the propagation of
errors for high redshift of the derivatives of the Hubble parameter. For the
case of the Hubble constant of the Planck collaboration the value of $j_{0}$
is close to the value of the flat $\Lambda CDM$ model.

With respect to the snap parameter if we consider the value of the matter
density parameter measured by the Planck collaboration, $\Omega_{m0}=0.315$,
then for the flat $\Lambda CDM$ model implies a value of $s_{0}=-0.4175$. Our
results with $1\sigma$ are marginally in agreement with this value (see values
in table 2). Secondly when $z>1$ the error increases remarkably because
reconstructing the parameter $s(z)$ implies reconstructing third-order
derivatives for the Hubble parameter.

In figure 7 we show the same results as figure 6 but for simulated data. On
the left side, we use simulated data with the error prescription from the $CC$
technique and on the right side we use the simulated data compatible with the
error from the $BAO$ technique. In both cases we observe that the
reconstructed functions behave similarly to the observational data, but with
less propagation of errors.

##  VI Discussion and Conclusions

In this paper we presented observational constraints on some of the
cosmographic functions, specially the deceleration parameter, the jerk
parameter and the snap parameter using a model-independent approximation. The
statistical method used has been the Gaussian processes through the $GaPP$
public package. The data used for this study are observational and simulated
data from the Hubble parameter. We directly estimate the value of the Hubble
constant from the data used, $H_{0}=69.45\pm 4.34$, and use this value to
perform our simulations. However, in order to observe the effect of priors
values of the Hubble constant on the functions $q(z)$, $j(z)$ and $s(z)$ we
use the values from the Planck collaboration and from the $SH0ES$
collaboration SHOES Collaboration et al. (2018) .

On the other hand, even though it is not our objective to discuss the Hubble
tension, it is interesting to note that our result marginally includes both
the Planck result and the SH0ES result, see figure 6. However, this is
primarily due to the large dispersion of the $CC$ measures. But what is
interesting is to observe that following the same methodology presented here
and using data from future observational projects such as the $LSST$ the $CC$
measurements can become excellent discriminators for the Hubble tension
problem.

We propose to use the cosmographic functions given by the equations (15-17)
instead of the cosmographic parameters derived from Taylor’s expansion
commonly used in the literature. By using cosmographic functions we avoid the
Taylor’s expansion truncation problem and consequently we can use data with
$z>1$ without restrictions. Also if we use cosmographic functions it is not
necessary a new redshift variable which can lead to inconsistent results.

In general we can observe that the current observational data of the Hubble
parameter do not allow to exclude the model $\Lambda CDM$ which still fits the
reconstruction of the cosmographic functions. On the other hand, the simulated
data, in particular, the Hubble parameter data obtained using the BAO
technique seems promising to distinguish between alternative models to the
$\Lambda CDM$ model.

An interesting application of our reconstruction is the study of specific dark
energy models by writing the parameters of a given dark energy model as a
function of the cosmographic parameters.

To implement this methodology we can proceed as follows: given a specific
Hubble function we can use our system of equations (15-17) to determine a
system of algebraic equations between the cosmographic parameters
$(q_{0},j_{0},s_{0})$ and the parameters of a given model $(\theta_{i})$.
Then, we can solve the model parameters as a function of the cosmographic
parameters: $\theta_{i}=\theta_{i}(q_{0},j_{0},s_{0})$.

As an instructive example in the present work we are going to apply this
methodology to a two-parameter cosmological model, specifically the $CPL$
model. In this case the Hubble function is given by equation (5) so that the
parameters of the model are $\theta_{i}=(w_{0},w_{1},\Omega_{m0})$, where
initially we assume the parameter $\Omega_{m0}$ as known. This means that to
determine the parameters $(w_{0},w_{1})$ we only need equations (15) and (16).
Thus, we get an algebraic system of equations that provide us with the
following solutions:

$\displaystyle w_{0}$ $\displaystyle=$
$\displaystyle-\frac{1}{(1-\Omega_{m0})}+\frac{2(q_{0}+1)}{3(1-\Omega_{m0})},$
(29)

$\displaystyle w_{1}$ $\displaystyle=$
$\displaystyle-\frac{2}{3}\frac{(j_{0}-1)}{(\Omega_{m0}-1)}-3w_{0}-3w_{0}^{2}.$
(30)

If we consider the best fits and the value $\Omega_{m0}=0.315$ (Planck) we
obtain for the parameters: $w_{0}=-1.132,w_{1}=-0.461$.

Additionally, if we consider the equation for the snap parameter, equation
(17), we can also determine the parameter $\Omega_{m0}$ as a function
explicitly of the cosmographic parameters so after some manipulations we
obtain,

$\displaystyle\Omega_{m0}$
$\displaystyle=\frac{-1+9j_{0}+6q_{0}+10j_{0}q_{0}-10q_{0}^{2}+2q_{0}^{3}}{2(-2q_{0}+(3+q_{0})(2j_{0}+q_{0}^{2})+s_{0})}$
(31)
$\displaystyle+\frac{(1-2q_{0})}{2(-2q_{0}+(3+q_{0})(2j_{0}+q_{0}^{2})+s_{0})}$
$\displaystyle\times\left(1+9j_{0}^{2}+2j_{0}(-9-4q_{0}+8q_{0}^{2})+4q_{0}(q_{0}(9+q_{0}(5+2q_{0}))+2(-1+s_{0}))-4s_{0}\right)^{1/2}$
$\displaystyle+\frac{s_{0}}{2(-2q_{0}+(3+q_{0})(2j_{0}+q_{0}^{2})+s_{0})}.$

Now we can use this equation to derive the value of $\Omega_{m0}$ directly
from the cosmographic parameters and then we can use this value to determine
the parameters $(w_{0},w_{1})$ using the equations (29-30). If we do this we
get: $w_{0}=-1.036,w_{1}=-0.119,\Omega_{m0}=0.247$.

In principle we can use this methodology for other cosmological models. If we
add more cosmographic functions, then we can include models with more free
parameters. In a subsequent investigation we intend to extend this methodology
to other models of dynamic dark energy. Therefore, it is important to note
that in this way we can infer the values of the free parameters of a given
model using the cosmographic approximation without using the Taylor expansion.

It is important to mention that recently the reference Munoz and Escamilla-
Rivera (2020)  have proposed a similar method but in another context. They
determine the cosmographic parameters as function of the parameters of
different dark energy models and use observational data to estimate the values
of the parameters of dark energy and, with these results, the cosmographic
parameters.

As we mentioned, our methodology is very useful to avoid using the Taylor
expansion, but in compensation it depends on the data of the Hubble parameter,
which are still a small sample when compared to supernovae Ia, for example.
Furthermore we need to reconstruct derivatives of observational data. The fact
that a small noise in the measurement data can cause a large error in the
derivatives is a difficult problem to analyze. For this reason, we hope to
investigate in a later work the use of different statistical methods to
reconstruct the derivatives as well as to discuss the importance of data
quality.

Figure 4: In the top the reconstruction of the Hubble Parameter using
observational data of Hubble parameter in a model-independent form. In the
middle the reconstruction of the Hubble Parameter with a prior $H_{0}=74.03$
($SH0ES$ collaboration). In the bottom the reconstruction using a prior
$H_{0}=67.44$ (Planck’s collaboration).

Figure 5: Observational constraints on the cosmography functions using the
observational data of the Hubble parameter in a model-independent form. Figure
6: Observational constraints on the Hubble constant, $H_{0}$, where we show
our result in comparison with the main results published in the literature,
for example, compiled in divalentino . The orange vertical band corresponds to
the SH0ES team and the blue vertical band corresponds to the Planck
Collaboration. The figure was made using the public program available online
at github.com/lucavisinelli/H0TensionRealm. divalentino .

The figure was made using the public program available at

Figure 7: Observational constraints on the cosmography functions using the
observational data of Hubble parameter for the two values of the Hubble
constant. In left side for $H_{0}=74.03$ and in right side for $H_{0}=67.44$.

Figure 8: In the column on the left we used the simulation that uses the $CC$ prescription and in the right column we have used the $BAO$ prescription. In both cases we have used the value of the Hubble parameter estimated directly from the observational data. The value used is in Table 2.  Table 2: The best-fitting parameter values and uncertainties with $1\sigma$ error bars for the cosmographic parameters from the observational and simulated data. The value $69.45$ $\pm$ $4.34$ is the best-fit for all observational data (CC + BAO) and is assumed as input for our simulations. Data |  $H_{0}$($1\sigma$) |  $q_{0}$($1\sigma$) |  $j_{0}$($1\sigma$) |  $s_{0}$($1\sigma$)
---|---|---|---|---
CC+BAO |  $69.45$ $\pm$ $4.34$ | $-0.672^{+0.265}_{-0.223}$ | $1.200^{+0.300}_{-0.390}$ | $-0.670^{+0.250}_{-0.390}$
CC+BAO |  $67.44$ $\pm$ $0.58$ (Value of Planck) | $-0.670^{+0.210}_{-0.120}$ | $0.985^{+0.414}_{-0.380}$ | $-0.600^{+0.379}_{-0.300}$
CC+BAO |  $74.03$ $\pm$ $1.42$ (Value of SH0ES) | $-0.710^{+0.159}_{-0.111}$ | $1.714^{+0.593}_{-0.592}$ | $-0.700^{+0.670}_{-0.567}$
Simulation CC |  $69.45$ $\pm$ $4.34$ | $-0.657^{+0.123}_{-0.187}$ | $1.155^{+0.501}_{-0.509}$ | $-0.899^{+0.485}_{-0.458}$
Simulation BAO |  $69.45$ $\pm$ $4.34$ | $-0.690^{+0.092}_{-0.085}$ | $1.291^{+0.285}_{-0.263}$ | $-0.716^{+0.294}_{-0.285}$

## Acknowledgements

A.M.V.T appreciate the computational facilities of the UFES to develop the
present investigation and also the authors thank Célia Escamilla-Rivera for
helpful discussions on the interpretation of cosmographic parameters. J. C.
Fabris thanks Fundação de Amparo Pesquisa e Inovação do Espirito Santo (FAPES,
project number 80598935/17) and Conselho Nacional de Desenvolvimento
Científico e Tecnológico (CNPq, grant number 304521/2015-9) for partial
support.

## References

  * (1) Riess. A. G. et al. (Supernova Search Team), Astron.J.116, 1009 (1998), astro-ph/9805201.
  * (2) Perlmutter.S. et al. (Supernova Cosmology Project), Nature 391, 51 (1998), astro-ph/9712212.
  * (3) Scolnic, D. M., et al., Astrophys.J. 859, 2, 101, (2018).
  * (4) Planck Collaboration et al., arXiv e-prints , arXiv:1807.06209 (2018), 1807.06209.
  * (5) Alam S., et al., 2017, Mon. Not. Roy. Astron. Soc., 470, 2617
  * (6) Peebles, P. J. E. and Ratra, B., Rev. Mod. Phys.75(2003) 559 [arXiv:astro-ph/0207347]
  * Perivolaropoulos (2008) Perivolaropoulos, L., arXiv:0811.4684
  * Bullock and Boylan-Kolchin, (2017) Bullock, J. S. and Michael Boylan-Kolchin, Annual Review of Astronomy and Astrophysics, vol. 55, pp. 343-387 (2017)
  * (9) Caldwell, R.R.; Dave, R.; Steinhardt, P.J., Phys. Rev. Lett. 80 (8), 1998 arXiv:astro-ph/9708069
  * (10) Almendariz.P, Mukanov.V and Steinhardt, P., Phys.Rev.D63, 103510,(2001)
  * (11) Winfried Zimdahl, Josep Triginer, Diego Pavon, Phys.Rev. D54 (1996) 6101-6110
  * (12) Shapiro, I.L and Sola, T., J.Phys.A40:6583-6593,2007
  * (13) Horndeski, G. W., Int. J. Theor.Phys. 10, 363 (1974)
  * (14) Weinberg, D., et al., Physics Reports, Volume 530, Issue 2, p. 87-255, (2013)
  * (15) Chevallier, M. and D. Polarski, Int. J. Mod. Phys. D 10,213 (2001).
  * (16) Linder, E. V., Phys. Rev. Lett.90, 091301 (2003).
  * (17) Ishida, E. E. O, Reis R. R. R, Toribio A. V. and Waga I, Astroparticle Physics 28, 7 (2007)
  * (18) Wetterich, C., Phys. Lett. B, 594, 17, 2004
  * (19) Copeland, E. J., Sami, M. and Tsujikawa, S., Int. J. Mod. Phys. D 15, 1753 (2006) [arXiv:hep-th/0603057];
  * (20) Weinberg, S., Gravitation and cosmology: Principles and applications of the general theory of relativity , Wiley, New York, (1972).
  * (21) Chiba, T. and Nakamura,T., Prog. Theor. Phys. 100 (1998) 1077 [arXiv:astro-ph/9808022].
  * (22) Visser, M., Class. Quant. Grav.21 (2004) 2603 [arXiv:gr-qc/0309109].
  * (23) Visser, M., Gen. Rel. Grav.37(2005) 1541 [arXiv:gr-qc/0411131].
  * (24) Caldwell, R and Kamionkowski, M., JCAP 09 (2004) 009 • e-Print: astro-ph/0403003
  * (25) Cattone, C and Visser, M. Class.Quant.Grav., 24, 5985-5998, (2007).
  * (26) Dunajski, M and Gibbons, G., Classical and Quantum Gravity, 25 235012, (2008)
  * de la Cruz-Dombriz, et al. (2016) de la Cruz-Dombriz, et al., JCAP 1612 12, 042, (2016).
  * Busti, et al., (2015) Busti, C. et al., (Phys.Rev. D92 (2015) no.12, 123512
  * (29) Padmanabhan T. and Choudhury,T.R., Mon. Not. Roy. Astr. Soc. 344, 823 (2003).
  * (30) Scherrer R. J., Phys. Rev. D, 92, 043001, (2015).
  * (31) Dos santos, V. M. Reis, R.R.R and Waga, I., JCAP 02, 066, (2016).
  * (32) Aviles, A., Gruber C., Luongo,O. and Quevedo, H., Phys. Rev.D,86, 123516, (2012).
  * (33) Capozziello, S. R. D’Agostino and O. Luongo, Mon. Not. Roy. Astron. Soc. 494 (2020) no.2, 2576
  * (34) Poplawski, N.J., Class.Quant.Grav.24:3013-3020,2007
  * Jimenez and Loeb (2002) Jimenez, R. and Loeb, A., ApJ 573, 37 (2002).
  * Jimenez et al. (2002) Jimenez, R., Verde, L., Treu, T and Stern, D., ApJ 593, 622 (2003).
  * (37) Gaztanaga E., Miquel R., Sanchez E., 2009, Phys. Rev. Lett., 103, 091302
  * (38) Zheng X., Ding X., Biesiada M., Cao S., Zhu Z., Astrophys.J. 825, 17, (2016). [arXiv:1604.07910]
  * (39) Velásquez-Toribio, A.M.; Magnago, A., EPJC, 80, 562, (2020).
  * (40) Podariu S., Souradeep T., Gott III J. R., Ratra B., Vogeley M. S., ApJ, 559, 9, (2001).
  * (41) Ma, C. and Zhang, T-J., Astrophys.J. 730, 74, (2011).
  * (42) Rasmussen, C. E. and Williams, C. K. I. , Gaussian Process for Machine Learning, 2006
  * (43) Busti, V.C., Clarkson, C and Seikel, M. MNRAS, 441,1, 11 (2014)
  * (44) Hai Yu, Ratra, B and Wang, F.Y, ApJ, 856, 3, (2018).
  * Seikel et al. (2012a) Seikel M., Clarkson C., Smith M.,JCAP, 6, 36 (2012).
  * (46) Seikel, M. and Clarkson, C., Phys. Rev. D 86, 083001 (2012).
  * (47) Velásquez-Toribio, A.M.; Machado, M.M. and Fabris, J., EPJC, 79, 1010, (2019).
  * (48) Holsclaw T et al., Phys. Rev. Lett.105 (2010) 241302 [arXiv:1011.3079];
  * (49) Holsclaw T et al., Phys. Rev. D 82 (2010) 103502 [arXiv:1009.5443]
  * (50) Holsclaw T et al., Phys. Rev. D 84 (2011) 083501 [arXiv:1104.2041].
  * (51) Farooq, O and Ratra,B., ApJ 766, L7 (2013). arXiv:1301.5243.
  * (52) Farooq, O., Crandall, S., and Ratra, B., Phys. Lett. B 726, 72 (2013).
  * (53) Farooq, O.,et al. Astrophys.J. 835, 1, 26 (2017).
  * (54) Yu, H., Ratra,B., and Wang,F.-Y., ApJ 856(1), 3 (2018). arXiv:1711.03437
  * (55) Mamon, A. A., Mod.Phys.Lett.A 33, 1850056 (2018).
  * (56) Velasquez-Toribio, A.M., and Bedran, M.L., Braz. J. Phys, 41, 59 (2011)
  * (57) SHOES Collaboration, ApJ, 855, 136, (2018).
  * (58) Munõz, C. Z. and Escamilla-Rivera, C. [arXiv:2005.02807]
  * (59) Tripp, R., A&A, 331, 815, (1998)
  * (60) Conley, A et al., ApJS, 192, 1, (2011).
  * (61) Arjona, R., Cardona, W. and Nesseris, S., Phys. Rev. D99, 043516 (2019), [arXiv:1811.02469].
  * (62) Di Valentino, E., et al. (2021), ArXiv:2103.01183

##  Appendix A Supernovae Ia

In this study we use the data from Supernovas Ia called "Pantheon" sample
(Scolnic et al. 2018 ) which is the largest combined sample of SNIa and
consists of $1048$ data with redshifts in the range $0.01<z<2.3$. It is a
collection of SNe Ia discovered by the Pan-STARRS1 (PS1) Medium Deep Survey
and SNe Ia from Low-z, SDSS, SNLS and HST surveys. This supernova Ia
compilation uses The SALT 2 program to transform light curves into distances
using a modified version of the Tripp formula Tripp 1998 ,

$\mu=m_{B}-M+\alpha x_{1}-\beta c+\Delta_{M}+\Delta_{B},$ (32)

where $\mu$ is the distance modulus, $\Delta_{M}$ is a distance correction
based on the host-galaxy mass of the SNIa and $\Delta_{B}$ is the distance
correction based on predicted bias from simulations. Also $\alpha$ is the
coefficient of the relation between luminosity and stretch, $\beta$ is the
coefficient of the relation between luminosity and color and $M$ is the
absolute $B$-band magnitude of a fiducial SNIa with $x_{1}=0$ and $c=0$. Also
$c$ is the color and $x_{1}$ is the light-curve shape parameter and $m_{B}$ is
the log of the overall flux normalization. An uncertainty matrix $\bf{C}$ is
defined such that,

$\chi_{SNIa}^{2}=\Delta\vec{\mu}^{T}.\mathbf{C}^{-1}.\Delta\vec{\mu},$ (33)

where $\Delta\vec{\mu}=\vec{\mu}_{obs}-\vec{\mu}_{model}$ and
$\vec{\mu}_{model}$ is a vector of distance modulus from a given cosmological
model and $\vec{\mu}_{obs}$ is a vector of observational distance modulus. The
distance module is defined as $\vec{\mu}=\vec{m}-M$, where $M$ is the absolute
magnitude and $\vec{m}$ is the apparent magnitude which is given by the
expression

$\displaystyle\vec{m}_{model}$
$\displaystyle=M+5Log_{10}(D_{L})+5Log_{10}(\frac{c/H_{0}}{1Mpc})+25$ (34)
$\displaystyle=\bar{M}+25+5Log(D_{L}).$ (35)

where $D_{L}=\frac{H_{0}}{c}d_{L}$ and $\bar{M}=M+5Log(\frac{c/H_{0}}{1Mpc})$
is a nuisance parameter, which depends on the Hubble constant $H_{0}$ and the
absolute magnitude $M$. To minimize with respect to the nuisance parameter we
follow a process similar at the references Conley et al. (2011)  and Arjona et
al. (2019) . Therefore the $\chi^{2}_{\bar{M}marg}$ is,

$\chi^{2}_{\bar{M}marg}=a+\log\frac{e}{2\pi}-\frac{b^{2}}{e},$ (36)

where,

$\displaystyle a$ $\displaystyle=$ $\displaystyle\Delta\vec{m}^{T}\cdot
C^{-1}\cdot\Delta\vec{m},$ (37) $\displaystyle b$ $\displaystyle=$
$\displaystyle\Delta\vec{m}^{T}\cdot C^{-1}\cdot\mathrm{I},$ (38)
$\displaystyle e$ $\displaystyle=$
$\displaystyle\mathrm{I}^{T}.C^{-1}.\mathrm{I}$ (39)

where $\Delta\vec{m}=\vec{m}_{obs}-\vec{m}_{model}$ and $\mathrm{I}$ is the
identity matrix.

##  Appendix B The Gaussianity of the data

An important test of Gasussianity is to determine the parameter $N_{\sigma}$
as shown in table 1. Then the weighted mean for the Hubble parameter is given
by Podariu et al. 2001 ,

$\bar{H}=\frac{\sum_{i=1}^{N}H_{i}(z_{i})/\sigma_{i}^{2}}{\sum_{i=1}^{N}1/\sigma_{i}^{2}}.$
(40)

We also compute the weighted bin redshift, it is given by,

$\bar{z}=\frac{\sum_{i=1}^{N}z_{i}/\sigma_{i}^{2}}{\sum_{i=1}^{N}1/\sigma_{i}^{2}},$
(41)

The associated weighted error is given by

$\bar{\sigma}=(\sum_{i=1}^{N}1/\sigma_{i}^{2})^{-1/2}.$ (42)

Using the above definitions we determine a reduced goodness-of-fit $\chi^{2}$
for each bin as

$\chi_{\nu}^{2}=\frac{1}{N-1}\sum{\frac{H-\bar{H}}{\sigma_{i}^{2}}}.$ (43)

And thus the number of standard deviations with respect from unity is defined
as

$N_{\sigma}=\mid\chi_{\nu}-1\mid\sqrt{2(N-1)}.$ (44)