Xin Su
email: xin.su@inria.fr

Mira Rizkallah
email: mira.rizkallah@irisa.fr

Thomas Maugey
email: thomas.maugey@inria.fr

Christine Guillemot
email: christine.guillemot@inria.fr

Rate-Distortion Optimized Super-Ray Merging for Light Field Compression

Keywords: Super-Ray (SR) Merging, Rate-Distortion Minimization, Light Field (LF) Compression, Shape-Adaptive DCT (SA-DCT)

In this paper, we focus on the problem of compressing dense light fields which represent very large volumes of highly redundant data. In our scheme, view synthesis based on convolutional neural networks (CNN) is used as a first prediction step to exploit inter-view correlation. Super-rays are then constructed to capture the inter-view and spatial redundancy remaining in the prediction residues. To ensure that the super-ray segmentation is highly correlated with the residues to be encoded, the super-rays are computed on synthesized residues (the difference between the four transmitted corner views and their corresponding synthesized views), instead of the synthesized views. Neighboring super-rays are merged into a larger super-ray according to a rate-distortion cost. A 4D shape adaptive discrete cosine transform (SA-DCT) is applied per super-ray on the prediction residues in both the spatial and angular dimensions. A traditional coding scheme consisting of quantization and entropy coding is then used for encoding the transformed coefficients. Experimental results show that the proposed coding scheme outperforms HEVC-based schemes at low bitrate.

I. INTRODUCTION

Light fields (LF) are defined as the representation of radiance of light rays emitted along several directions by the different points in a 3D scene. Several devices have been developed for light fields capture, either based on camera arrays [START_REF] Wilburn | High performance imaging using large camera arrays[END_REF], on single moving cameras, or on arrays of microlenses [START_REF] Ng | Light field photography with a hand-held plenoptic camera[END_REF], etc. Light fields have recently gained in popularity due to the variety of potential applications in computational photography and computer vision, however they represent very large volumes of data with challenges in terms of storage, transmission and processing.

In this paper, we focus on the problem of compressing dense light fields captured by plenotic cameras. First methods for compressing synthetic light fields appeared late 90's essentially based on classical coding tools as vector quantization or using JPEG coding for each sub-aperture view, yielding however limited compression performances. It is only recently that compression solutions have been proposed for dense real light fields captured by plenoptic cameras. The proposed solutions can be classified into two categories: either coding the array of sub-aperture images extracted from the lenslet image as a pseudo video sequence in [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF], [START_REF] Liu | Pseudosequence-based light field image compression[END_REF], or directly encoding the lenslet images captured by plenoptic cameras In [START_REF] Liu | Pseudosequence-based light field image compression[END_REF]- [START_REF] Jia | Optimized Inter-View Prediction Based Light Field Image Compression With Adaptive Reconstruction[END_REF], with extensions of HEVC with dedicated prediction modes. Multiview video compression and disparity compensation techniques are considered in [START_REF] Li | Efficient intra prediction scheme for light field image compression[END_REF]. A homography-based low rank approximation [START_REF] Jiang | Light Field Compression With Homography-Based Low-Rank Approximation[END_REF] is used to exploit angular correlation of LF. Besides being represented and encoded as images or videos, the LF is represented by 4D Gaussian mixture models in [START_REF] Verhack | Steered mixture-of-experts for light field coding, depth estimation, and processing[END_REF] and by graphs containing minimum amount of color and disparity information in [START_REF] Su | Graph-based light fields representation and coding using geometry information[END_REF].

In this paper, we propose a compression scheme based on view synthesis. Four corner views of the LF are first encoded by HEVC-Inter and transmitted. The whole LF is then synthesized from the four corner views using the convolutional neural networks (CNN) based architecture proposed in [START_REF] Kalantari | Learning-based view synthesis for light field cameras[END_REF]. The prediction residues are then transformed using a 4D-shape adaptive Discrete Cosine Transform (4D SA-DCT) which exploits both spatial and angular correlation remaining in the residue signals. The support of the 4D SA-DCT is defined by a segmentation of the light field into super-rays. Super-rays can be seen as a set of super-pixels that are coherent across all light field views, taking into account disparity information. Note that local transforms have also been investigated for light fields compression in [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF], however, the support of the local transform was defined by co-located super-pixels in the different views, not taking into account disparity.

To ensure that the super-ray segmentation is highly correlated with the residual signals, the super-rays are computed on synthesized residues (the difference between the four transmitted corner views and their corresponding synthesized views). Neighboring super-rays with similar homogeneous residues are merged into a larger super-ray to have a better spatial energy compaction by optimizing a rate-distortion cost. Experimental results show that the proposed coding scheme yields rate-distortion gains at low bitrates (e.g. < 0.04 bpp corresponding to a PSNR quality up to 35 dB) compared with HEVC-based coding schemes, while being comparable or slightly worse at higher bitrates. CNN based synthesis method [START_REF] Kalantari | Learning-based view synthesis for light field cameras[END_REF], as shown in Fig. 1 (red arrows). To improve the quality of the synthesized light field, the residuals between the synthesized and original views are encoded by applying local super-ray based shape adaptive DCT (SA-DCT) (see Fig. 1, blue arrows).

II. LIGHT FIELD CODING SCHEME

A. Scheme Overview

1) The processing chain of the four corner views: At the decoder, the decompressed four corner views are used to synthesize the whole LF using the CNN-based architecture of [START_REF] Kalantari | Learning-based view synthesis for light field cameras[END_REF], as shown by the yellow region in Fig. 1. The first CNN is trained to model the disparity from the four input views, while the second CNN is used to estimate the color of the synthesized views. The synthesis quality depends on the QP value of the HEVC-inter coder.

As shown by the yellow region in Fig. 1, we compute the residual signals as the difference between the decompressed four corner views LF cor and their corresponding synthesized views LF cor . The four images of residues are then warped onto the other views using the disparity estimated by the CNN. The super-ray segmentation is then computed by applying the SLIC algorithm [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF] on the set of residue images, but also taking into account disparity when performing the clustering. Computing the segmentation on the residue images rather than on the synthesized views, similar residues signals are more likely to be grouped into one segment which can benefit the following energy compaction in transform domain.

2) The processing chain of the residues: The synthesized views LF and the super-ray segmentation are computed in the same manner at the encoder. We apply a local spatial SA-DCT on the residuals for each view using the super-ray as a support of the transform. Spatial SA-DCT coefficients of each superray corresponding to the same frequency form a U × V block in the angular domain, on which a second angular SA-DCT is applied to capture angular dependencies. A traditional coding scheme consisting of quantization and entropy coding is then used for encoding the transformed coefficients. At the decoder side, the decoded residuals are added to the synthesized views to obtain the final decompressed LF.

B. Super-Ray Segmentation on Residues

The concept of super-ray has been introduced in [START_REF] Hog | Superrays for Efficient Light Field Processing[END_REF] as an extension of super-pixels [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF] to group light rays coming from the same 3D object, i.e. to group pixels having similar color values and being close spatially in the 3D space. While the authors in [START_REF] Hog | Superrays for Efficient Light Field Processing[END_REF] estimate disparity only at the centroid of the super-rays, here we consider a scheme using dense disparity maps to synthesize the entire light field from a sparse set of views. The disparity maps used in the tests reported below, have been estimated by the first CNN of the view synthesis architecture of [START_REF] Kalantari | Learning-based view synthesis for light field cameras[END_REF] from the four corner views which are available at both the encoder and decoder. Having a dense disparity map for each view, the pixels in all the views are clustered using a method similar to SLIC [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF] with a weighted combination of a color distance, a spatial distance and in terms of depth.

Fig. 2 (a) shows that the super-rays computed on the set of synthesized color images are well aligned with the edges of the objects, but are not really suited for the residue images they are supposed to represent (see the edge of the building). Indeed the residues generally lie on both sides of object edges. Since the real residue signal is not available at the decoder, we generate an estimation of it (called synthesized residue), that is used to build a more accurate segmentation. For that purpose, the residue images at the four corner views are first computed as the difference between the coded/decoded corner views and their synthesized versions (see Fig. 1). The residue images at the other viewpoints are then obtained by warping the corner residue images to the other positions using the estimated disparity. We see in Fig. 2 (b) that the super-ray segmentation computed on these synthesized residue images has better correlation with the real residual signals, which is proved by the energy compaction comparison between Fig. 2 (a) and (b) in section IV-A. Since the synthesized residues are available at both the encoder and decoder, we can obtain the same super-ray segmentation on both sides, and do not need to transmit it.

C. 4D Shape Adaptive DCT (SA-DCT)

While 2D DCT applied on a square or rectangular support may fail to capture correlation at image edges, we consider here a separable 4D shape-adaptive DCT with a support defined by the super-ray segmentation. Fig. 3 illustrates how a 4D SA-DCT is applied on the i-th super-ray SR i in the LF. A spatial 2D SA-DCT is first applied per view on each superray (i.e. on the super-pixel of the view which belongs to the considered super-ray). In each view, the obtained coefficients form a rectangular block with non-zero coefficients only in the top-left area with the DC component located at the left-top corner, as shown in the middle of Fig. 3. These coefficients are sorted from low to high frequency by following a Zig-Zag order. Coefficients corresponding to the same frequency band but from different views form an U × V block in angular domain, as shown on the right of Fig. 3. Another SA-DCT, i.e. an angular SA-DCT, is then applied on this block. Note that some values may be missing in the U × V block, since the size of SR i varies in different views. Generally, the 4D DCT T(•) can be computed as

{X i,b } = T(SR i ) , T(•) = DCT u ⊗ DCT v ⊗ DCT y ⊗ DCT x Spatial SA-DCT Angular SA-DCT (1)
where

{X i,b } ∈ R U ×V ×N , N is the maximum size of SR i in different views, b = 1, 2, • • • , U × V × N .
For some values of b, X i,b may be missing due to the non-regular shape of SR i . The positions of missing elements in {X i,b } are available at both encoder and decoder, since the super-ray segmentation is known. DCT * ∈ R n×n corresponds to a n-point DCT (n is the number of elements of SR i in corresponding coordinate), and ⊗ denotes a Kronecker product operator.

D. Quantization and Entropy Coding

At the end of those two transform stages, coefficients are grouped into a 2-dimensional array X where X(i, b) is the b-th band in super-ray SR i . Using the observations on all the super-rays in a training dataset (Rose [START_REF] Jiang | Light Field Compression With Homography-Based Low-Rank Approximation[END_REF]), we can find the best ordering for quantization. We first sort the variances of coefficients with enough observations in decreasing order. We then split them into 64 groups g: {X i,b , ∀i, b ∈ g}. All the remaining coefficients with less observations will be considered in the last group. We use the zigzag ordering of the JPEG quantization matrix to assign the quantization step size Q for each. A simple rounding procedure then results in the quantized coefficients X q (i, b) = round(X(i, b)/Q(X(i, b)). that are further coded using an arithmetic coder.

III. RATE-DISTORTION OPTIMIZED SUPER-RAY MERGING

In order to increase compression performances, we improve the super-ray segmentation with a rate-distortion optimized super-ray merging. Four initial segmentations are performed using different initial numbers of clusters leading to different super-ray sizes as shown in Fig. 4.a. The resulting segmentations are referred to as layers. Note that we modify the super-ray segmentation in layer l respecting to the boundaries of super-rays in layer l -1, to make sure the boundaries are coherent at different layers. For instance, if super-ray SR l-1 i in layer l -1 is across two (or maybe more) superrays in layer l, one of these super-rays in layer l is enlarged to completely contain SR l-1 i , and the other super-rays are reduced correspondingly. We choose to enlarge the super-ray that initially contains most parts of SR l-1 i . The merging results {SR i } are initialized by the super-rays {SR l=0

i , ∀i} at layer 0 and the merging starts from layer 1 to layer 3. At each time, we only consider one super-ray SR l i at layer l, which consists of several super-rays {SR l-1 j } at layer l -1, i.e. SR l i = {SR l-1 j } = {SR l-1 j ∈ SR l i , ∀j}, as shown in Fig. 4 (b). {SR l-1 j } will be merged into SR l i , i.e. {SR l-1 j } ⇒ SR l i , if and only the rate-distortion cost J reduces after merging.

The rate cost of the 4D SA-DCT coefficients (after quantization) is modeled by their entropy computed per coefficient group as

R = g entropy ({Q(X i,b ), ∀i, b ∈ g}) . ( 2
)
The distortion is computed as the distortion of SA-DCT coefficients before and after quantization as

D = i b [X i,b -Q(X i,b )] 2 . ( 3
)
Thus, the rate-distortion cost is

J = D + λR , ( 4
)
where λ = 1 represents the relation between rate and distortion. The merging procedure based on the minimization of the rate-distortion cost J is detailed in Algorithm 1. Fig. 4 (c) gives an illustration of merging results using Flower1 dataset in [START_REF] Kalantari | Learning-based view synthesis for light field cameras[END_REF].

} = {SR l-1 j ∈ SR l i , ∀j} in layer l -1; Compute J with {{SR i } \ {SR l-1 j }} SR l i , i.e. {SR l-1 j } are replaced by SR l i Eq.4; if J > J then Merge, {SR l-1 j } ⇒ SR l i ; Update {SR i }; end end end

IV. EXPERIMENTS

We test our coding scheme on four real LF with 8 × 8 subaperture images of size (X = 536, Y = 376) from the dataset used in [START_REF] Kalantari | Learning-based view synthesis for light field cameras[END_REF], called Flower1, Flower2 and Cars.

A. Energy Compaction

We first evaluate the effectiveness of our 4D segmentation, by analyzing the compaction of the energy after transfor-mation. Therefore, 4D SA-DCT is applied on: 1) superpixels used in [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF] (computed on synthesized color without disparity compensation or merging) 2) super-rays computed on synthesized residues without disparity compensation or merging 3) super-rays computed on synthesized residues with disparity compensation but without merging 4) super-rays computed on synthesized residues with disparity compensation and merging. Fig. 5 shows the percentage of energy carried by a given percentage of coefficients obtained by SA-DCT on these segmentations. The blue curve is the baseline method [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF]. The red curves show the impact of using the synthesized residues to compute the segmentation. However, due the error in the synthesized residue, the improvement is limited. The Yellow curve shows the impact of using the disparity information, while the purple curve measures the effect of the merging. Thanks to the merging operation which compensates the errors in the super-ray segmentation, the proposed contributions bring a significant increase in terms of energy compaction (∼10%) with respect to a direct use of a super-pixel segmentation per view [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF].

Percentage of coefficients (%)

B. Rate-PSNR results

The final rate-distortion performance of the proposed scheme is evaluated in comparison with three baseline methods: 1) HEVC-lozenge [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF], the whole LF is considered as a video and compressed by HEVC with a lozenge sequence, 2) CNN+HEVC [START_REF] Rizkallah | Impact of light field compression on focus stack and extended focus images[END_REF], the same CNN based view synthesis is applied here, while the residues are compressed by HEVC, 3) CNN+SA-DCT (no merging, no disparity) [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF], our previous coding scheme presented in [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF] using the same CNN based view synthesis, however, there is no disparity compensation or super-merging strategy. Note that the coding methods using CNN based prediction are performed with best pairs of parameters (Q, QP ) where Q is the quality parameter used to compress the residues and QP is used in the HEVC intercoding of the four corners. The obtained rate-distortion curves are shown in Fig. 6.

The proposed CNN+SA-DCT coding scheme yields better or comparable rate-PSNR performance at low bitrate than HEVC based reference methods. The improvement of the proposed CNN+SA-DCT compared with the baseline method in [START_REF] Rizkallah | Graph-based Transforms for Predictive Light Field Compression based on Super-Pixels[END_REF] indicates the effectiveness of super-ray merging. It allows the proposed coding scheme to capture more information with fewer bits (at low bitrate), compared with HEVC based encoders. Tab. I shows the improvement in terms of bjontegaard metric at low bitrate (< 0.04 bpp corresponding to a PSNR quality up to 35 dB) obtained by our coding scheme. However, as shown in Tab. I, at bitrates higher than 0.04 bpp, the HEVC based encoders (HEVC lozenge and CNN+HEVC) generally outperform the proposed coding scheme at high bitrates. This is due to the fact that the proposed scheme does not have very complex and high quality prediction strategies in residue coding which is useful at high bitrate.

V. CONCLUSION

In this paper, we have presented a rate-distortion optimized super-ray merging to exploit the correlation in the spatial and angular dimensions of light fields. After the CNN-based view synthesis, the residue inside each super-ray is compacted into a few coefficients using 4D shape-adaptive DCT transform. The experimental results show that the proposed light field coding scheme can yield rate-distortion gains compared with HEVC based compression, especially at low bitrate.

Fig. 1

 1 Fig.1 depicts the proposed coding scheme. Let LF = {I u,v } denote a light field, where u = 1, . . . , U and v = 1, . . . , V are the view indices. Four views at the corners LF cor = {I 1,1 , I 1,V , I U,1 , I U,V } are encoded using HEVC-Inter and used to synthesize the whole light field with the

Fig. 1 .

 1 Fig. 1. Overview of proposed coding scheme.

  (a) Super-ray segmentation computed on synthesized color is shown with corresponding color image (left) and real residual image (right). (b) Super-ray segmentation computed on synthesized residues is shown with corresponding synthesized residual image (left) and real residual image (right).

Fig. 2 .

 2 Fig. 2. Super-ray segmentation computed on (a) synthesized color and (b) synthesized residues. Only the center view (u = 4, v = 4) has been shown here. Segments in (a) are well aligned with the edges of color image, however, not aligned with the discontinuities within the residue images, see the border of the building. Segments in (b) obtained using synthesized residues better follow the discontinuities of the real residue signals.

Fig. 3 .

 3 Fig. 3. An illustration of 4D SA-DCT performed on super-ray.

  (a) Super-rays in different layers. (b) Merging. (c) Merging result of Flower1 [15].

Fig. 4 .

 4 Fig. 4. An illustration of super-ray merging. (a) Super-rays at different layer l = 0, 1, 2, 3 from left to right. (b) Super-rays {SR l-1 j } = {SR l-1 j ∈ SR l i , ∀j} in layer l -1 are merged into SR l i , i.e. {SR l-1 j } ⇒ SR l i , if the rate-distortion cost J has been reduced. (c) A merging result of Flower1 [15].

Algorithm 1 :

 1 Merging based on rate-distortion minimization Data: Super-rays at different layer {SR l=0 i }, {SR l=1 i }, {SR l=2 i }, {SR l=3 i } Result: Merged super rays {SR i } Initialization: {SR i } = {SR l=0 i }; for Layer l = 1 to 3 do for Each super-ray SR l i in layer l do Compute J with {SR i } Eq.4; Find super-rays {SR l-1 j

Fig. 5 .

 5 Fig. 5. Energy compaction of the transformed coefficients of Flower1 and Cars using different super-ray computations. The super-rays computed on synthesized residues with merging and disparity compensation yield higher energy compaction.

  Fig. 6. Rate-distortion comparison. From left to right: Flower 1, Cars and Flower 2.

	40												38										40
		CNN+HEVC
		HEVC lozenge
		CNN+SA-DCT										36
		CNN+SA-DCT in [16]
	35												34										35
	Y PSNR (dB)		32 34								Y PSNR (dB)	30 32			32.5 33							Y PSNR (dB)				33 34
	30															32							30
			30									28			31.5											32
																31
			28												30.5
				0.005	0.01	0.015	0.02	0.025	0.03	0.035	0.04		26			0.01	0.015		0.02	0.025	0.03	0.035					31	0.01	0.015	0.02		0.025	0.03
	25																						25
	0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2	0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2	0	0.02	0.04	0.06	0.08	0.1	0.12	0.14	0.16	0.18	0.2
					Bitrate(bpp)										Bitrate(bpp)									Bitrate(bpp)

TABLE I BJONTEGAARD

 I COMPARISON (∆PSNR (DB)) AT low BITRATE (< 0.04 BPP) AND high BITRATE (> 0.04 BPP)

				Our CNN+SA-DCT vs
		CNN+HEVC	HEVC lozenge	CNN+SA-DCT in [16]
		Low	High	Low	High	Low	High
	Flower 1	0.22	0.05	1.92	-0.19	0.03	0.55
	Cars	0.21	-0.17	0.48	-0.50	0.17	1.08
	Flower 2	0.1	-0.26	1.77	-0.85	0.09	0.38

ACKNOWLEDGMENT

This work has been supported in part by the EU H2020 Research and Innovation Programme under grant agreement No 694122 (ERC advanced grant CLIM).


Emmanuel Paradis

Molecular dating of phylogenies by likelihood methods: a comparison of models and a new information criterion

Keywords: large-scale estimation, molecular dating, penalised likelihood, rate smoothing, relaxed molecular clock

Dating the divergence in a phylogenetic tree is a fundamental step in evolutionary analysis. Some extensions and improvements of the penalised likelihood method originally presented by Sanderson are introduced. The improvements are the introduction of alternative models, including one with non-correlated rates of molecular substitution ("relaxed" model), a completely reworked fitting algorithm that considers the high-dimensionality of the optimisation problem, and the development of a new information criterion for model selection in the presence of a penalised term. It is also shown that the strict clock model is a special case of the present approach. An extensive simulation study was conducted to assess the statistical performance of these improvements. Overall, the different estimators studied here appeared as unbiased though their variance varied depending on the fitted and the simulated models and on the number of calibration points. The strict clock model gave good estimates of branch lengths even in the presence of heterogeneous substitution rates. The correlated model gave the best estimates of substitution rates whatever the model used to simulate the data. These results, which are certainly the first from an extensive simulation study of a molecular dating method, call for more comparison with alternative methods, as well as further work on the developments introduced here.

Introduction

The molecular divergence among a set of molecular sequences is the product of the time they separated and the rate at which substitutions accumulated. The crux of the problem of phylogenetic molecular dating lies in the fact these two components are generally confounded, and untangling them can be done only with some assumptions. One of these assumptions is that molecular substitutions accumulate at a constant rate-the molecular clock. When this assumption is applied in its most radical version (i.e., a single rate for all taxonomic groups and constant in time), the exercise of molecular dating is considerably simplified since it can be accomplished even without a phylogenetic tree, for instance with pairwise genetic distances.

It is a fact that the molecular clock cannot be held generally, though this assumption may be of some value for recently diverged sequences [START_REF] Brown | Rate variation and estimation of divergence times using strict and relaxed clocks[END_REF]. A number of studies have characterized variation in substitution rates among taxonomic groups and through time (e.g., [START_REF] Pereira | A mitogenomic timescale for birds detects variable phylogenetic rates of molecular evolution and refutes the standard molecular clock[END_REF]. Thus, a rooted phylogenetic tree appears as appropriate to handle this problem because the internal nodes of such a tree might be interpreted as branching points through time. [START_REF] Sanderson | A nonparametric approach to estimating divergence times in the absence of rate constancy[END_REF] was one of the first authors to propose a general solution by assuming that substitution rates on two contiguous branches are likely to be similar (in other words, they are auto-correlated). By contrast, previous authors handled the same problem by applying strong restrictions on how rates vary in a tree (see citations in [START_REF] Sanderson | A nonparametric approach to estimating divergence times in the absence of rate constancy[END_REF].

In a subsequent paper, [START_REF] Sanderson | Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach[END_REF] refined his approach by introducing a likelihood component while keeping the assumption of auto-correlated rates. At the same time, Bayesian methods were developed with similar modeling assumptions though with different fitting approaches (e.g., [START_REF] Thorne | Estimating the rate of evolution of the rate of molecular evolution[END_REF][START_REF] Thorne | Divergence time and evolutionary rate estimation with multilocus data[END_REF]. Since then, the Bayesian approach to molecular dating has flourished in a diversity of refinements (e.g., [START_REF] Drummond | Relaxed phylogenetics and dating with confidence[END_REF][START_REF] Yang | Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds[END_REF][START_REF] Lepage | A general comparison of relaxed molecular clock models[END_REF][START_REF] Lartillot | PhyloBayes 3: a Bayesian software package for phylogenetic reconstruction and molecular dating[END_REF][START_REF] Guindon | Bayesian estimation of divergence times from large sequence alignments[END_REF][START_REF] Dos Reis | Approximate likelihood calculation on a phylogeny for Bayesian estimation of divergence times[END_REF][START_REF] Wilkinson | Dating primate divergences through an integrated analysis of palaeontological and molecular data[END_REF]. The Bayesian approach has not been the only one explored during the last decade. [START_REF] Britton | Phylogenetic dating with confidence intervals using mean path lengths[END_REF] proposed the mean path length to estimate divergence times from a non-ultrametric tree with branch lengths in number of substitutions. [START_REF] Xia | A distance-based least-square method for dating speciation events[END_REF] developed a method based on least squares to estimate a chronogram including the possibility to model auto-correlated rate variation. These two methods have the advantage of being easily and quickly run.

In the next section, I first review Sanderson's framework and show how it can be generalized to include other models of substitution rate variation. I then introduce a new information criterion that helps to decide which models of rate variation best describes the data. I also propose some improvements on model fitting by penalised likelihood. An extensive simulation was done to assess the statistical properties of the methods presented in this paper.

Methods

Sanderson's penalised likelihood framework

Sanderson's penalised likelihood function (denoted Ψ) is based on two components:

Ψ = ln L -Φ, (1)
where L is the likelihood function of a model of the branch lengths, and Φ is a function constraining rate variation. Essentially, L is a parametric component and Φ is a nonparametric one. [START_REF] Sanderson | Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach[END_REF] originally used a Poisson model for the first component and a constraint that rates on contiguous branches are more likely to be similar. Specifically, the likelihood function is:

L = ∏ i ζ x i i e -ζ i x i ! ,
with x i being the number of substitutions observed on branch i, and ζ i = r i t i where r i is the substitution rate and t i is the time length of the same branch. The x i 's are the input data usually from a non-ultrametric tree estimated by maximum likelihood, Bayesian estimation, or a distance-based method. Typically, they will be expressed as mean number of substitution per site, so it is not required to know the sequence length. The product is made over all branches of the tree. The log-likelihood function is thus:

ln L = ∑ i x i ln ζ i -ζ i -ln x i !. ( 2
)
The nonparametric part is:

Φ = ∑ k, j
(r k -r j ) 2 + var(r basal ),

where k and j denote two contiguous branches (i.e., consecutive in time), and the subscript 'basal' denotes the basal branches of the tree. Φ is further multiplied by a "smoothing" parameter λ which controls the trade-off between both components. In the present paper, I will leave aside the problem of estimating the smoothing parameter and will set, somehow arbitrarily, λ = 1 in the rest of the study.

It appears clearly that such a penalised likelihood approach can be a general framework which does not need to be restricted to the above assumptions on rate variation. [START_REF] Kim | Penalized likelihood phylogenetic inference: bridging the parsimony-likelihood gap[END_REF] discussed how this can be generalised in a broader context of phylogenetic inference in-cluding linking likelihood, parsimony, and Bayesian methods. In the next section I show how non-correlated relaxed clock models can be implemented using this framework in the context of molecular dating.

Non-correlated relaxed clock model

It is possible to relax the assumption of auto-correlation among rates; however, to permit estimation it is necessary to introduce another constraint on rate variation. Bayesian methods typically use a lognormal prior on substitution rates: such a distribution can be used to constrain the overall variation of these rates over the tree. Other distributions may be used, such as Γ, uniform, or Beta.

Let F be the theoretical cumulative density function of the distribution of the substitution rates.

The nonparametric component can now be defined as:

Φ = ∑ i F (r i ) -F(r i ) 2 ,
where F (r i ) is the empirical cumulative distribution function of the rates. This is then used in Equation (1) to estimate the parameters in the same way than for the auto-correlated model (see

Fitting Algorithm

). This penalty function can be seen as a special case of the general formulation from Kim and Sanderson (2008, their eq. 3).

How to choose a formulation for F? Typically, previous studies suggest that substitution rate variation follows a power-like distribution with most rates being small and few large. Several probabilistic functions may be used to model such distributions. The gamma law, Γ α,s , appears as a practical choice here since the shape parameter α can be estimated with the mean of the rates throughout the tree, and the scale parameter s can be fixed equal to one, resulting in an appropriate distribution. Thus the penalty term above can be calculated without the need to estimate additional parameters.

It must be noted that the expression 'relaxed clock' has been used with different meanings in the literature. To avoid confusion, in this paper I name "relaxed" model the model assuming that each branch of the tree has its own substitution rate with no assumption on correlation between contiguous branches.

Maximum likelihood approach with discrete rate variation

An issue in dating a phylogenetic tree with variable rates is how one decides which part of the tree evolves faster than others. [START_REF] Yang | A heuristic rate smoothing procedure for maximum likelihood estimation of species divergence times[END_REF] discussed the problem of assigning a rate value to a particular branch (see also a critique of earlier approaches by [START_REF] Sanderson | A nonparametric approach to estimating divergence times in the absence of rate constancy[END_REF]. On the other hand, the penalised likelihood approach does not assign rates to branches, though it requires to estimate a rate value to each branch which may result in the need to estimate a large number of parameters. [START_REF] Yang | A heuristic rate smoothing procedure for maximum likelihood estimation of species divergence times[END_REF] developed a method to assign rates to branches based on the initial (non-clock) tree.

The divergence times are then estimated by maximum likelihood.

Here I propose a method that is inspired from discrete rate category method used in phylogenetic inference [START_REF] Yang | Maximum likelihood phylogenetic estimation from DNA sequences with variable rates over sites: approximate methods[END_REF]. With this method, it is assumed that rates vary in a discrete way so that we can make categories of branches characterised by different rates. However, we do not know in which category each branch belongs to, but by estimating the frequency of each category, we can calculate the contribution of each branch to Equation ( 2) by summing the contribution of each rate weighted by its frequency. Let us denote as c the number of categories, so the number of parameters to be estimated, apart of the dates, will be 2c -1. Typically, c will be small (between 1 and 10) because if it is large then the present model tend to the "relaxed" one with the additional cost of estimating frequencies. Obviously, c = 1 implies a strict clock model.

By contrast to the two previous non-clock models, the present one is fitted by full maximum likelihood since there is no penalty term Φ.

A new information criterion

The recent literature has seen a number of contributions developing alternative models relaxing the molecular clock hypothesis. An important issue is how to select the correct model. Bayesian methods use criteria based on Bayes factors (e.g., [START_REF] Lepage | A general comparison of relaxed molecular clock models[END_REF][START_REF] Ho | Improved characterisation of among-lineage rate variation in cetacean mitogenomes using codon-partitioned relaxed clocks[END_REF]. In the present framework, a likelihood-based criterion must be used. A widely used model selection criterion is the Akaike information criterion (AIC, [START_REF] Akaike | Information theory and an extension of the maximum likelihood principle[END_REF]. We recall that AIC = -2 ln L + 2k with k being the number of parameters estimated from the data (sometimes called the free parameters). Here ln L is given by Equation 2, but we must also include the contribution of the penalty term Φ so that the smaller Φ the smaller the information criterion. In other words, we should favour models with predictions on rate variation that best conform to the data. A simple solution could be to compute -2 ln L + 2k + λΦ but this would ignore the fact that Φ is actually made of several elements. A solution is to make a multivariate decomposition of Φ which will summarize how its elements vary. Among the possibilities to do this, a singular value decomposition comes as a good choice since it does not require a special form of Φ (eigenvalue decomposition requires a square matrix). So the new criterion, which I denote as ΦIC, is defined as:

ΦIC = -2 ln L + 2k + λδ {Φ} ,
where δ x is the singular value of x. The notation {Φ} is to point out that the decomposition is done on a vector made with the elements of Φ and not on Φ itself. Like for other information-based criteria (see review in [START_REF] Konishi | Information Criteria and Statistical Modeling[END_REF], the model with the smallest value of ΦIC is selected as the best model describing the data.

Fitting algorithm

Molecular dating of phylogenies usually implies estimating a large number of parameters. For instance, fitting the penalised models to a fully dichotomous tree with n tips and a single calibration point will need to estimate 3n -4 parameters (2n -2 rates and n -2 dates) whereas fitting the discrete categories model will imply estimating n + 2c -3 parameters. For instance, if n = 100 (a medium-sized tree in today's standards) and c = 1, the number of free parameters will be 296 and 99, respectively. Giving the difficulties of optimising a function with so many parameters, Bayesian methods have clearly known a wide success. [START_REF] Sanderson | Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach[END_REF] reported several algorithmic difficulties when optimising his penalised likelihood function. He used a gradient-free method as well as quasi-Newton gradient-based methods.

Here, I use the PORT optimization routines [START_REF] Gay | Usage summary for selected optimization routines[END_REF] as implemented in R version 2.15.2 (R Development Core Team, 2012). The advantage of PORT is that they are relatively robust to irregularities in the objective function (by contrast to some other non-linear optimisation methods such as [START_REF] Schnabel | A modular system of algorithms for unconstrained minimization[END_REF]. Furthermore, they can handle infinite values, and can be used with or without gradients (i.e., first partial derivatives) of the objective function. If these gradients are not provided, the routine computes them numerically using the objective function. If this one has few parameters (less than ten), both methods (with or without gradients) have similar computing times; however, in the case of many parameters the gain in running times can be at least one order of magnitude. The gradients of the penalised likelihood functions used in this paper and their derivations are provided in Appendix A.

Another advantage of PORT is that bounds for parameter estimation can be set. This is used for dates known within intervals which are thus considered as additional free parameters to be estimated within these intervals. For the substitution rates, the lower and upper bounds are set to ε and 100 -ε where ε is by default set to 10 -8 (this can be changed by the user). A convergence diagnostic is output by PORT which makes possible to assess the quality of the solution. It is also possible to control the number of iterations used during optimisation which appears useful here because convergence may be very slow to achieve due to the high dimensionality of the problem; so, the user may assess whether increasing the number of iterations could improve convergence.

Another control on the optimization process is the size of the steps around the current solution to examine whether a better solution can be found. Setting the minimum and the maximum step size correctly improves substantially the optimisation. In the present problem, there are two kinds of parameters, the rates and the dates, which typically have distinct ranges of variation: 10 -8 -10 -2

for the former, and 1-100 for the latter. Therefore, I have implemented an algorithm of alternate optimizations which can be summarised as follows:

1. A first set of estimates is obtained by optimising the objective function over all parameters.

The step size is set to 1 (the default of PORT).

2. Optimise the objective function over the rates keeping the dates constant at the values of the current solution. The step size is allowed to vary between 10 -8 and 0.1.

3. Optimise the objective function over the dates keeping the rates constant at the values output at step 2. The step size is allowed to vary between 10 -3 and 0.5.

4. If the objective function value is improved, the solution found from the combination of steps 2 and 3 is taken as the new solution. These are repeated until convergence-or until a fixed number of alternate optimisations has been performed.

Finally, the unknown dates are initially set using a random algorithm which, when repeated several times, allows one to assess the importance of these initial values on the final results.

The methods presented in this paper have been implemented in ape [START_REF] Paradis | APE: analyses of phylogenetics and evolution in R language[END_REF]. The computing time of this new implementation was assessed with random trees of size 20, 50, 100, or 200 tips, simulated using the three models described in the next section, and fitting to each of them four different models: the "relaxed" model, the correlated model, the clock model (c = 1), and a discrete model with c = 10. The number of calibration points were 1, 3, 5, or 10 as described below. For each combination of these components, the simulation was replicated ten times resulting in 1920 recorded running times which were measured with the R function system.time on a modern laptop (hardware: CPU duo-core 2.93 GHz; software: Ubuntu 12.04 64-bit, R 2.15.2, ape 3.0-6). A regression analysis was performed to characterise the most important factors affecting computing times. The best model can then be used to predict the computing times for larger data sets. Note that this simulation protocol is almost similar to the one described below with the difference that one additional model was fitted and fewer replications were used.

Simulation study

Ultrametric trees were simulated with random topologies generated with the rtree function in ape [START_REF] Paradis | Analysis of Phylogenetics and Evolution with R (Second Edition)[END_REF]. In order to not bias branch lengths in a particular direction, the branching times of these trees were drawn from a uniform distribution between 1 and 50. The number of tips were equal to 20, 50, 100, or 200. This was replicated 100 times so that 400 initial ultrametric trees were generated. From each of them, three trees were simulated with different models of variation in the substitution rates: (i) the rates were generated from a Γ distribution with shape α = 0.5 and rate ρ = 1; (ii) the rates were simulated along the initial tree with a Brownian motion model with variance σ 2 = 0.01 and rescaled to have a similar range of variation than in the previous model; (iii) a clock model by adding to each branch length a small noise generated from a normal distribution with mean 0 and variance 0.5. In the first two models, the branch lengths of the initial trees were multiplied by the simulated rates. Thus, 1200 such trees were generated.

For each of the 400 initial trees, a set of ten calibration points was made by extracting the branching times of the root node (which was always included in the subsequent analyses) and nine randomly chosen nodes. The 1200 non-ultrametric trees were analysed by fitting three models: "relaxed", correlated, and clock (with c = 1) using 1, 3, 5, or 10 calibration points (which were considered as exactly known dates), so 4800 data sets were analysed (1200 trees × 4 sets of calibration points). In the end, 14,400 chronograms were estimated (4800 data sets × 3 models). All trees, chronograms, and sets of calibration points were saved for further analyses. For each estimated chronogram, the mean of the difference between the real and the estimated branch lengths and the mean of the difference between the real and the estimated substitution rates were computed.

Results

The 1920 computing times varied between 0.017 and 325.777 seconds (mean = 20.71, median = 2.42, SD = 45.28). Table 1 presents a summary of computing times for the considered tree sizes. The longest computing time (5 min 26 sec) was obtained fitting a clock model to a tree with n = 200 simulated with a correlated model using five calibration points. Remarkably, fitting the same model to the same tree with 10 calibration points was much faster (< 2 sec). However, the computing times had a very asymmetric distribution with a mean much higher than the median. A relatively complex linear model was found to fit well: Figure 1 shows one of the 400 sets of simulated trees. Among the 4800 simulated data sets, 174 could not be analysed completely because of lack of convergence of the model fitting procedure. Overall, the ΦIC selected the correct model in 61% of the cases; however, this percentage depended on the simulated model since when the simulated model was the "relaxed" one the correct model was selected in 82.7%, and the percentage reached 100% when the simulated model was the clock one (Table 2).

Overall, and somehow surprisingly, the results were not affected by the size of the tree, so they are presented below pooled over the different values of n. Figure 2 shows the distribution of the error on branch length estimation with respect to whether the correct model was used for fitting the chronogram and the number of calibration points. The eight boxes reveal no bias in these estimations even when the wrong model was fitted. On the other hand, using the correct model substantially decreased the dispersion of the estimates. Without surprise, increasing the number of calibration points led to better estimates even with the wrong model.

Figure 3 displays the same results broken down with respect to the nine combinations of simulated-fitted models. Overall, the clock model gave good estimates of branch lengths even if the trees were simulated with other models. When the data were simulated with a clock model, the fit with this model yielded almost unbiased branch length estimates with a very small variance (rightmost panel).

Like for the branch lengths, the mean error rate on the estimation of substitution rates was not affected by tree size. No bias was observed in these estimates; however, their dispersion was reduced when fitting the correct model and/or when the number of calibration points increased (Fig. 4). Figure 5 shows how this error varied with respect to the simulated and fitted models. In most cases, no bias was observed. The estimates from the "relaxed" model were slightly positively biased (i.e., overestimated) when the trees were simulated with another model (first box in second and third panels). Whatever the simulated model, a surprising result appears: the correlated model seems to estimate substitution rates better than the two other models, and it performed almost as well as the clock model when the data were clock-like. The dispersion of these errors were much larger when fitting a "relaxed" model compared to the two others though this was much reduced with the data simulated from this model in spite of a number of outliers with positive errors.

Discussion

Molecular dating is a fundamental step in evolutionary analysis. It is used to assess the tempo of phenotypic evolution [START_REF] Garland | Rate tests for phenotypic evolution using phylogenetically independent contrasts[END_REF], the rates of lineage diversification [START_REF] Nee | Tempo and mode of evolution revealed from molecular phylogenies[END_REF],

or to date the colonization of new areas or habitats [START_REF] Warren | Molecular phylogeography reveals island colonization history and diversification of western Indian Ocean sunbirds (Nectarinia: Nectariniidea)[END_REF]. Currently, two classes of methods are widely used: Bayesian methods, which are by far the most widespread with the computer programs Multidivtime [START_REF] Thorne | Divergence time and evolutionary rate estimation with multilocus data[END_REF], BEAST [START_REF] Drummond | BEAST: Bayesian evolutionary analysis by sampling trees[END_REF], and PAML [START_REF] Yang | PAML 4: Phylogenetic analysis by maximum likelihood[END_REF], and the penalised likelihood method as implemented in the program r8s [START_REF] Sanderson | r8s: inferring absolute rates of molecular evolution and divergence times in the absence of a molecular clock[END_REF]. Other methods have been proposed: mean path length [START_REF] Britton | Phylogenetic dating with confidence intervals using mean path lengths[END_REF], likelihood methods (e.g., [START_REF] Yang | A heuristic rate smoothing procedure for maximum likelihood estimation of species divergence times[END_REF], or least squares using distances as input data [START_REF] Xia | A distance-based least-square method for dating speciation events[END_REF].

In spite of this diversity of methods, there is apparently a confusion among practitioners about the relative merits of each method. Probably one reason for this is the scarcity of assessment of the statistical performance of these methods. [START_REF] Brown | Rate variation and estimation of divergence times using strict and relaxed clocks[END_REF] compared the performance of strict and relaxed clock methods by simulating phylogenies with n = 5, 10, or 20 species. To my knowledge, the present paper is the first extensive study of the statistical properties of a molecular dating method in realistic conditions of data analysis.

Computing time is clearly an issue in data analysis. In a recent paper, dos Reis and Yang (2011) stated "a typical Bayesian analysis for a phylogeny of <50 species might take several days." Such long calculations are at the expense of the time spent in the interpretation and understanding of the results. The computing times recorded in the present study were always less than six minutes, but this masks a great heterogeneity since 50% of these, regardless of the tree size, were less than or equal to 2.42 seconds. Furthermore, the linear regression analysis, characterising almost 80% of the variation among the recorded computing times, showed that they are affected by a number of factors including some interactions among them. The possibility to fit several models in reasonable times is clearly a strength for the present likelihood and penalised likelihood methods. Besides, short computing times make easier the assessment of statistical properties with simulations as done here.

Several results emerge from the present simulation study. First, using the correct model (i.e., the one used to simulate the data) gives better estimates of branch lengths. Second, using more calibration points results also in better estimation. Third, and this is more surprising, using the strict clock model gives good estimates of branch lengths whatever the model used to simulate the data. A possible explanation for this result may come from the way substitution rates were simulated. These rates were generated with stable distributions, so that in the estimation process replacing those rates by their mean value could result in easier estimation, especially considering that the clock model has less parameters to estimate than the two penalised models.

The fact that the number of tips n in the tree did not affect estimation may seem surprising at first because bigger trees have more observations (the x i 's in Eq. 2). However, the number of parameters to estimate for all the models considered here is proportional to n, so increasing the size of the tree does not actually increase the quantity of information available for estimation.

The simulation study reported here considered that the input phylogenetic trees are correctly estimated. In reality this may not be always the case, but the goal here was to focus on some issues that have not been studied in-depth until now. Moreover, the simulation approach used here is similar to the one used in other studies (e.g., [START_REF] Brown | Rate variation and estimation of divergence times using strict and relaxed clocks[END_REF]. Another assumption of the present simulations is that the calibration points are known exactly and without error. This also is not likely to be true in real situations. This point has been extensively discussed in the literature. For instance, [START_REF] Yang | Bayesian estimation of species divergence times under a molecular clock using multiple fossil calibrations with soft bounds[END_REF] have shown that, in the case of Bayesian estimation, incorrectly assumed "hard bounds" (exactly known ages) will lead to wrong results whatever the quantity of sequence data.

A point that will need future attention is how the present models behave in the presence of more complex variation in substitution rates. Some studies have suggested that these rates may vary in a non-equilibrium fashion (e.g., [START_REF] Magallón | Fossils to break long branches in molecular dating: a comparison of relaxed clocks applied to the origin of angiosperms[END_REF], and others suggest a relationship with life-history traits [START_REF] Lartillot | A phylogenetic model for investigating correlated evolution of substitution rates and continuous phenotypic characters[END_REF][START_REF] Mayrose | A likelihood method for detecting trait-dependent shifts in the rate of molecular evolution[END_REF]. Since the penalised models estimate a different substitution rate for each branch, they seem appropriate to handle more complex cases than considered here. This clearly requires further investigation.

The present study clealy calls for a wider assessment of the statistical properties of the methods developed here in comparison with other methods. This will require a much larger simulation survey (e.g., [START_REF] Wertheim | Relaxed molecular clocks, the bias-variance trade-off, and the quality of phylogenetic inference[END_REF]. An issue that also needs more attention is the characterisation of the relative merits of different methods. For instance, the mean path lengths method [START_REF] Britton | Phylogenetic dating with confidence intervals using mean path lengths[END_REF] provides a test of the molecular clock for each node of the tree. Whether this is a valuable tool for preliminary explorative analyses before running more complex model fitting, with likelihood or Bayesian models, is not yet clear.

A point of the present study that will require further study is how to select the best model in a penalised likelihood framework. The proposed information criterion, ΦIC, gave some equivocal results. The criterion correctly selected the "relaxed" model in the majority of cases which appears as a positive result. On the other hand, when the data were simulated with the correlated model, the ΦIC selected the clock model in most cases. This may be understood easily because fitting the correlated model tends to minimise the difference in substitution rates among contiguous branches, thus tending to a clock model. In other words, if the rate of autocorrelation is strong, the variation among branches will be low so the observed variation will tend to a clock model. On the other hand, if the autocorrelation is weak, then the rates will not be serially correlated among branches and they will tend to vary in a "relaxed" way. This may explain why the correlated model is likely to not be selected by the proposed criterion. A possible solution to this problem would be to use the ΦIC criterion only to compare the "relaxed" and the clock models. A related issue, which was left aside for the present study, is that the amount of smoothing is actually controlled by a parameter, denoted as λ by [START_REF] Sanderson | Estimating absolute rates of molecular evolution and divergence times: a penalized likelihood approach[END_REF], which was fixed to one in all the results reported here. This clearly needs further study. Furthermore, [START_REF] Kim | Penalized likelihood phylogenetic inference: bridging the parsimony-likelihood gap[END_REF] showed some difficulties in finding the best value of λ by cross-validation and called for the development of a better scheme to this end.

The present study aimed to investigate how Sanderson's penalised likelihood approach could be improved by extending it to other models of substitution rate variation. Some computational improvements have also been done, and a new information criterion for model selection has been developed. These improvements provide a framework for molecular dating of phylogenies which seems an attractive alternative to Bayesian methods, particularly considering the possibility to analyse large trees (n = 200) in short times (a few minutes). Four perspectives may be drawn from this work. First, the interplay between the clock and correlated models may lead to improved estimation. Since the clock model provided the best estimates of branch lengths, and the correlated model provided the best estimates of substitution rates, both models could be combined. Second, the derivation of standard-errors of the estimates was left aside for the present study. As done elsewhere [START_REF] Xia | A distance-based least-square method for dating speciation events[END_REF], bootstrap or other resampling methods can be used. Another approach might be to use profile likelihood-based methods to infer uncertainty of the parameter estimates. Because of the penalty term, this does not seem straightforward and is currently under study. Third, it will be necessary to study the issue of the smoothing parameter λ and its role in parameter estimation. Finally, the statistical properties and scope of the proposed information criterion ΦIC will need to be further investigated.

We note that t i , the length of branch i, is not a parameter and is given by t i = d a i -d b i where the d's are dates (measured from present back to the past) and a i is the node ancestor of node b i along branch i. The dates are parameters unless they are fixed calibration points.

We can now derive the first partial derivatives of ln L which is simple for the rates:

∂ ln L ∂r i = x i r i -t i ,
and slightly more complicated for the dates. We first rewrite the log-likelihood as:

ln L = ∑ i x i ln r i + x i ln(d a i -d b i ) -r i (d a i -d b i ) -ln x i !,
which makes clearer how we obtain the derivates with respect to the unknown dates (this uses ln(u) = u /u):

∂ ln L ∂d k = ∑ i x i t i -r i A(k, i) - x i t i -r i D(k, i),
where the index k is for all unknown dates, A(k, i) is an indicator function taking the value 1 if k is a node ancestor (basal) of branch i or 0 otherwise, and similarly for D(k, i) if k is a node descendant (terminal). The advantage of this formulation, from a computational point of view, is that A and D can be built once as incidence matrices and the above equation is then calculated with two matrix products after calculating the terms within parentheses for all branches.

We now consider the penalty term Φ which differs depending on the model. Only the rates are concerned here since the dates are not involved in the calculation of Φ.

For the auto-correlated model, we rewrite Φ as:

Φ = ∑ k, j r 2 k -2r k r j + r 2 j + var(r basal ),
reminding that (k, j) is for all pairs of contiguous branches. So for branch i we have (ignoring the variance term for the moment):

∂Φ ∂r i
= ∑ i, j

(2r i -2r j ) + ∑ k,i (-2r k + 2r i ).
The first sum is for all pairs (i, j) where i is an ancestral branch connected to branch j, and the second sum has only one term for the pair (k, i) with k being the ancestral branch of i. So, this can be simplified as:

∂Φ ∂r i = 2 η i r i - η i ∑ j r j ,
where η i is the number of branches connected to branch i (η i = 1 for the terminal branches), and the sum is now over all branches connected to i (ancestral or not, but not sister-ones). Like discussed above, this can be computed in a straightforward way since the η i 's can be built once, and also an incidence matrix indicating whether two branches are connected.

For the basal branches (i.e., connected to the root), a term related to the variance of their rates must be added to the above, that is for branch i:

2 r i 1 - 1 n basal - 1 n basal n basal ∑ j =i r j n basal -1 ,
where i and j is for the basal branches and n basal is their number. This simplifies to r i -r j if

n basal = 2.
For the "relaxed" model, we use the fact that by definition the derivative of the empirical cumulative distribution function, which is a step function, is zero, and the derivative of the theoretical cumulative distribution function F(x) is, also by definition, the density function f (x). So, using

(u 2 ) = 2u u: ∂Φ ∂r i = -2 f (r i ) F (r i ) -F(r i )

Hessian

The matrix of second partial derivatives is obtained by deriving the first partial derivatives with respect to all parameters:

∂ 2 Ψ ∂θ i ∂θ j .
The Hessian is a square matrix with its diagonal elements being ∂ 2 Ψ/∂θ 2 i . As above, me way decompose for the two components of Ψ, so:

∂ 2 ln L ∂r 2 i = - x i r 2 i ∂ 2 ln L ∂r i ∂r j = 0 i = j.
Rewriting the above first derivative as ∂ ln L/∂r i = x i /r i -(d a i -d b i ), we easily find:

∂ 2 ln L ∂r i ∂d k = -A(k, i) + D(k, i).
We now turn to the first partial derivatives with respect to d and derive them a second time with respect to r:

∂ 2 ln L ∂d k ∂r i = -A(k, i) + D(k, i).
The second partial derivatives with respect to d are more complicated to derive (this uses (1/u) = -u /u 2 ):

∂ 2 ln L ∂d 2 k = ∑ i x i t 2 i A(k, i) + x i t 2 i D(k, i), ∂ 2 ln L ∂d k ∂d l = ∑ i - x i t 2 i A(k, i)D(l, i) - x i t 2 i A(l, i)D(k, i) k = l.
However, this last equation is simpler, since this is equivalent to consider only the branch i defined by nodes k and l and thus simplifies to

-x i /t 2 i . This is because, either A(k, i)D(l, i) is equal to 1 only if k is the ancestor of l along i, or A(l, i)D(k, i) is equal to 1 only if l is the ancestor of k along i.
We now consider Φ. For the auto-correlated model, we have:

∂ 2 Φ ∂r 2 i = 2η i ∂ 2 Φ ∂r i ∂r j = -2,
where j are the branches connected to i (see above). Clearly, these depend only on the η i 's and so do not need to be updated during optimisation.

For the "relaxed" model, the generic form may be written as:

∂ 2 Φ ∂r 2 i = -2 ∂ f ∂r i F (r i ) + 2 ∂ f ∂r i F(r i ) + 2 f (r i ) 2 , = -2 ∂ f ∂r i F (r i ) -F(r i ) .
Because of the difficulty in deriving the density function f , this is not considered further. Wrong model Correct model q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 1

Wrong model Correct model q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 5

-10 -5 0 5 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 10

Fitted model

Error on branch lengths -10 -5 0 5 "Relaxed" Correlated Clock q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Simulated model: "Relaxed" "Relaxed" Correlated Clock q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Simulated model: Correlated "Relaxed" Correlated Clock q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Simulated model: Clock

Error on branching times -20 -10 0 10 Wrong model Correct model q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 1

Wrong model Correct model q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 3 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 5

-20 -10 0 10 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Nb of calibration points: 10

  ln T ∝ ln n + Ncal + Model + Sim + n : Model with T the computing time, Ncal the number of calibration points, Model a categorical variable specifying the fitted model, Sim a categorical variable specifying the simulated model, and n : Model the first-order interaction term between these two variables. The effect of ln n was positive (estimated coefficient: 2.18, SE = 0.09) while the effect of Ncal was negative (-0.061, SE = 0.007). Interestingly, n was a better predictor than the number of estimated parameters (k which is highly correlated with n). The explained variance was high (adjusted R 2 = 0.79), and the residual standard-deviation was σ = 1.059 (as a comparison the standard-deviation of ln T was 2.32). An examination of the residuals of this model showed a relatively homogeneous distribution suggesting a satisfactory overall fit (see Supplementary Information). Giving values of n, Ncal, Model, and Sim it is possible to predict the mean computing time, E(T ), taking the exponential of the above linear model. Because of the logarithmic transformation of T , a 95% prediction interval has to be computed with E(T )/e 1.96 σ and E(T ) × e 1.96 σ instead of the usual ±1.96 σ. An easy-to-use R program is provided in the Supplementary Information to perform these calculations. For instance, for n = 1000 and Ncal = 1, E(T ) = 7 min 38 sec for fitting a clock model on data following the same model (95% prediction interval: [57 sec, 1 h 1 min]). On the other hand, for the same values of n and Ncal, if the data come from a correlated model and one wishes to fit the same model, then E(T ) = 34 h ([4 h 16 min, 271 h 6 min]).

Fig. 1 .

 1 Fig. 1. (a) One of the 400 random ultrametic trees (n = 100). The diamonds indicate the selected calibration points. The other trees were generated from this one as explained in the text: (b) under the "relaxed" model with rates following a Γ distribution, (c) under the correlated model, and (d) under the clock model.

Fig. 2 .

 2 Fig. 2. Distribution of the mean error on branch length estimation with respect to whether the correct model was used for fitting the chronogram (x-axis) and the number of calibration points (panels). The mean error was calculated as the difference of the estimated and the true values of the parameters averaged over the whole tree. The box limits give the first and third quartiles and the whiskers extend to 1.5 times the box limits. The extreme values going beyond the whiskers are shown with open circles. The filled circle indicates the median.

Fig. 3 .

 3 Fig. 3. Distribution of the mean error on branch length estimation with respect to the fitted model of the chronogram (x-axis) and the model used to simulate the data (panels). See Fig. 2 for details.

Fig. 4 .

 4 Fig. 4. Distribution of the mean error on substitution rate estimation with respect to whether the correct model was used for fitting the chronogram (x-axis) and the number of calibration points (panels). See Fig. 2 for details.

Fig. 5 .

 5 Fig. 5. Distribution of the mean error on substitution rate estimation with respect to the fitted model of the chronogram (x-axis) and the model used to simulate the data (panels). See Fig. 2 for details.

Table 1 :

 1 Mean, median, standard-deviation (SD), minimum, and maximum from a sample of 1920 computing times (in seconds). n: number of taxa.

	n Mean Median	SD Min	Max
	20	0.67	0.14	2.04 0.02	16.04
	50	3.71	1.10	8.06 0.10	69.84
	100 15.65	5.41 22.73 0.35 170.13
	200 62.84	26.49 71.65 1.39 325.78

Table 2 :

 2 Number of times when the fitted model was selected with ΦIC with respect to the simulated model.

	Simulated model	Selected model
		"Relaxed" Correlated Clock
	"Relaxed"	1276	3	263
	Correlated	292	0 1250
	Clock	0	0 1542

Acknowledgements

I am grateful to two anonymous reviewers and the Associate Editor for their constructive comments on a previous version of my manuscript. The simulations benefited from the ISEM computing cluster platform. This is publication IRD-ISEM 201x-xxx.

Appendix A. Derivation of gradients and Hessian

Gradients

The objective function is made of the difference of two functions, so the derivatives can be decomposed in two parts (with θ i being any free parameter of the model):

So we consider first the log-likelihood part:

where the index i is for all branches of the tree. Since ζ i = r i t i , this has to be written as:

Fitted model

Error on branching times -20 -10 0 10 "Relaxed" Correlated Clock q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Simulated model: "Relaxed" "Relaxed" Correlated Clock q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Simulated model: Correlated "Relaxed" Correlated Clock q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Simulated model: Clock