# On the Maximal Independent Sets of $k$-mers with the Edit Distance

Leran Ma  Schreyer Honors College, The Pennsylvania State University, United
States Department of Computer Science and Engineering, The Pennsylvania State
University, United States  Ke Chen  Department of Computer Science and
Engineering, The Pennsylvania State University, United States  Mingfu Shao
Department of Computer Science and Engineering, The Pennsylvania State
University, United States  Huck Institutes of the Life Sciences, The
Pennsylvania State University, United States

###### Abstract

In computational biology, $k$-mers and edit distance are fundamental concepts.
However, little is known about the metric space of all $k$-mers equipped with
the edit distance. In this work, we explore the structure of the $k$-mer space
by studying its maximal independent sets (MISs). An MIS is a sparse sketch of
all $k$-mers with nice theoretical properties, and therefore admits critical
applications in clustering, indexing, hashing, and sketching large-scale
sequencing data, particularly those with high error-rates. Finding an MIS is a
challenging problem, as the size of a $k$-mer space grows geometrically with
respect to $k$. We propose three algorithms for this problem. The first and
the most intuitive one uses a greedy strategy. The second method implements
two techniques to avoid redundant comparisons by taking advantage of the
locality-property of the $k$-mer space and the estimated bounds on the edit
distance. The last algorithm avoids expensive calculations of the edit
distance by translating the edit distance into the shortest path in a
specifically designed graph. These algorithms are implemented and the
calculated MISs of $k$-mer spaces and their statistical properties are
reported and analyzed for $k$ up to 15. Source code is freely available at
https://github.com/Shao-Group/kmerspace.

##  1 Introduction

Given an alphabet $\Sigma$, we denote by $S_{k}$ the set of all sequences of
length $k$ over $\Sigma$. Clearly, $|S_{k}|=|\Sigma|^{k}$. In this work, we
call a sequence of length $k$ a $k$-mer, and call $S_{k}$ the $k$-mer space,
for the sake of convenience. The Levenshtein distance [9], also known as the
edit distance, between two sequences $u$ and $v$, denoted as
$\textrm{edit}(u,v)$, is defined as the minimum number of insertions,
deletions, and substitutions needed to transform $u$ into $v$. The edit
distance is a metric; the $k$-mer space together with the edit distance
$(S_{k},\textrm{edit})$ forms a metric space.

Despite the ubiquitous use of $k$-mers and the edit distance in computational
biology, the properties and structures of the metric space
$(S_{k},\textrm{edit})$ remain largely unknown. It is worth noting that
$(S_{k},\textrm{edit})$ is intrinsically different from the well-studied
metric spaces on sequences such as the Hamming distance space or other normed
spaces (when a $k$-mer is represented as a vector of dimension $k$), as
evidenced by the fact that embedding the edit distance on $\\{0,1\\}^{k}$ into
the $L_{1}$ space requires a distortion of $\Omega(\log k)$ [14, 8].

In this work, we study the _maximal independent set_ (MIS) of a $k$-mer space
parameterized with an integer $d\geq 0$. Formally, we say a subset $M\subset
S_{k}$ is _independent_ if $\textrm{edit}(u,v)>d$ for every two $k$-mers $u$
and $v$ in $M$, and we say an independent subset $M$ is _maximal_ if there
does not exist another independent subset $M^{\prime}$ such that $M\subsetneq
M^{\prime}$. These definitions can be equivalently stated using the language
of graph. Given an integer $d\geq 0$, we define an undirected graph
$G^{d}_{k}=(V_{k}=S_{k},E^{d}_{k})$ where there is an edge $(u,v)\in
E_{k}^{d}$ if $\textrm{edit}(u,v)\leq d$. The above defined MISs of the
$k$-mer space with respect to an integer $d$ is equivalent to MISs in the
graph $G^{d}_{k}$.

We address the problem of finding an MIS of a $k$-mer space. Solving this
problem helps understand the structure of the metric space
$(S_{k},\textrm{edit})$ as an MIS indicates how dense a $k$-mer space is under
the edit distance. More importantly, finding MISs of $G^{d}_{k}$ has critical
applications in large-scale sequence analysis thanks to its nice properties.
We elaborate on a few below.

First, an MIS $M$ naturally can be used as a set of “centers” for clustering
$k$-mers. The centers ($k$-mers in an $M$) are distant from each other which
means that clusters are not crammed together. In addition, every $k$-mer $v$
can find a “nearby” center, as there exists a $k$-mer $u\in M$ such that
$\textrm{edit}(u,v)\leq d$ (otherwise $M\cup\\{v\\}$ is an independent set
that strictly contains $M$). These properties make an MIS an ideal set of
centers for clustering $k$-mers.

Second, $k$-mers of an MIS can be used to design a new seeding scheme for
aligning error-prone sequencing reads. Modern fast aligners often use exact
$k$-mer matches as anchors (or seeds) [3, 10]. However, such methods exhibits
low sensitivity when aligning reads with high error-rate (i.e., long-reads
data generated by PacBio [16] and Oxford Nanopore [5] protocols), because
under high error-rate, two biologically related sequences hardly share any
identical $k$-mers (for a reasonable choice of $k$ such as $k=15$ and an error
rate of 15%). One can choose to use smaller $k$ but this will lead to high
false positives as unrelated sequences can, by chance, share many short
$k$-mers. A new seeding scheme can be designed by first mapping a $k$-mer into
its nearest $k$-mer $u$ in an MIS $M$: two $k$-mers form an anchor if they are
mapped to the same $k$-mer in $M$. The advantage of this MIS-based seeding
scheme is that it tolerates edits in $k$-mers and is therefore more sensitive
for error-prone sequencing data. Comparing to other $k$-mer-alternative
seeding methods such as spaced seeds [2, 11] and indel seeds [12] where only
restricted types of edits are allowed (for example, spaced seeds cannot model
indels), this MIS-based scheme can recognize similar $k$-mers with both indels
and substitutions. Moreover, the intrinsic properties of an MIS also control
the false-positives: by triangle inequality, two distant $k$-mers that are
$2d$ edits apart will never be mapped to the same $k$-mer.

Third, an MIS of the $k$-mer space can be used to improve _sketching_
approaches. Sketching enables scaling to large datasets by only selecting a
subset of representative $k$-mers in a sequence. Existing $k$-mer sketching
methods such as FracMinHash [1, 4] often use a random permutation of $k$-mers
and pick the top fraction as representative $k$-mers. By using an MIS as a
representative subset instead, the selected $k$-mers are guaranteed to be at
least of $d$ edits apart of each other, and are hence more efficient.

All above applications require finding an MIS in the first place. However,
computing an MIS of a $k$-mer space is difficult, simply because the number of
$k$-mers $|S_{k}|=|\Sigma|^{k}$ (equivalently, the number of vertices in the
graph $G_{k}^{d}$) grows exponentially with $k$ and the number of edges in the
graph $|E_{k}^{d}|=O(|\Sigma|^{2k})$ grows even faster. Existing algorithms
for finding an MIS in a general graph, for example, the greedy algorithm that
takes $O(m)$ time for a graph with $m$ edges, are therefore not suitable for
our problem. In fact, even constructing $E_{k}^{d}$ explicitly is often
unaffordable. In this work, we design three algorithms for finding an MIS of a
$k$-mer space. Our algorithms take into account the special properties of
$k$-mers and the edit distance so that the underlying graph $G_{k}^{d}$ does
not need to be explicitly constructed. Our algorithms are able to scale to
instances with $k=15$. We implemented these algorithms and calculated an MIS
for all combinations of $k$ and $d$ for $k\leq 15$ and reported their
statistics. We also analyzed and concluded when to use which algorithm for
different combinations of $k$ and $d$.

##  2 Algorithms

We design three algorithms for finding a MIS of all $k$-mers, given $k$ and
$d$. The first one is a greedy algorithm that is similar to the greedy
algorithm for computing an MIS on general graphs but without explicitly
building all edges of the underlying graph. The second algorithm improves the
first one by reducing redundant comparisons through recognizing the locality
properties of $k$-mers and estimating and incorporating the bounds of the edit
distance into the algorithm. The last algorithm transforms calculating the
edit distance into finding the shortest path in a specifically constructed
graph, and finding an MIS is then transformed into efficient graph traversing
together with accompanied data structures to speed up.

###  2.1 A Simple Greedy Algorithm

This algorithm maintains a (dynamic) array $M$ that stores the current MIS,
initialized as an empty array. The algorithm examines each $k$-mer in $S_{k}$:
for the current $k$-mer $v$, it compares with all $k$-mers in $M$, and adds
$v$ to $M$ if $\textrm{edit}(u,v)>d$ for every $u\in M$. See the pseudocode
given below.

We now show that this algorithm is correct, i.e., the returned $M$ is indeed
an MIS. According to the algorithm, $k$-mer $v$ is added to $M$ only if
$\textrm{edit}(u,v)>d$ for all $u\in M$. This means that any two $k$-mers in
$M$ have an edit distance larger than $d$, i.e., $M$ is independent. On the
other hand, for any $k$-mer $v$ not included in $M$, the algorithm guarantees
that there exists a $k$-mer $u$ in $M$ such that $\textrm{edit}(u,v)\leq d$;
this implies that $M$ is maximal.

Algorithm 1 Simple Greedy Algorithm

initialize an empty array $M$

for each $k$-mer $v\in S_{k}$ do

$\mathit{isMapped}\leftarrow\text{false}$

for each $k$-mer $u\in M$ do

if $\textrm{edit}(u,v)\leq d$ then

$\mathit{isMapped}\leftarrow\text{true}$

break

end if

end for

if $\mathit{isMapped}==\text{false}$ then

add $v$ into $M$

end if

end for

return $M$

The above algorithm runs in $O(|M|\cdot|\Sigma|^{k}\cdot d\cdot k)$ time, as
in the worst case it compares each $k$-mer in $S_{k}$ with all $k$-mers in
$M$, and determining if $\textrm{edit}(u,v)\leq d$ for two $k$-mers $u$ and
$v$ takes $O(d\cdot k)$ time. The running time is output-sensitive. It is in
favor of instances with small MIS.

###  2.2 An Improved Greedy Algorithm

We say a $k$-mer $v$ is _mapped_ to a $k$-mer $u$ in the (partially)
constructed MIS $M$ if $\textrm{edit}(u,v)\leq d$. Such a mapping explains why
$v$ is not selected into the MIS. Algorithm 1 essentially finds a mapping for
each $k$-mer $v\in S_{k}\setminus M$ using an iterative, exhaustive search. As
the size of an MIS increases, this can significantly impair the performance of
the algorithm. Observe that the search order matters: if a mapping of a
$k$-mer is found early, all the following comparisons can be avoided.
Furthermore, based on the results of previous comparisons, it can be inferred
that some $k$-mers in the MIS are close to $v$ while some others are too far
away to serve as mappings. We design two techniques that allow us to quickly
determine if $v$ can be mapped to some $k$-mer in $M$, by utilizing the
locality-properties, and to quickly filter out those $k$-mers in $M$ that $v$
cannot be mapped to, by estimating the bounds on the edit distance. These two
techniques serve as extensions to the first algorithm to avoid unnecessary
comparisons.

The first technique is based on the observation that a $k$-mer is likely to
have a shared mapping with its “neighboring” $k$-mers. For a $k$-mer $v$, we
define its _nearest neighbors_ to be the subset of $S_{k}$ in which each
$k$-mer has edit distance $1$ from $v$, i.e., a substitution. For any given
$k$-mer, clearly, its nearest neighbor set contains $k(|\Sigma|-1)$ $k$-mers.
For example, the set of nearest neighbors of the 3-mer $AAA$ is
$\\{AAC,AAG,AAT,ACA,AGA,ATA,CAA,GAA,TAA\\}$. If $AAA$ is selected into the
MIS, all the above 9 $k$-mers will have it as a shared mapping. This locality-
property suggests a modification to Algorithm 1: store the mapping information
of $k$-mers (once found) so that potential shared mappings can be examined
first. Note that a $k$-mer can have multiple valid mappings in the MIS, but we
only keep record of one to save both time and space. When checking whether a
$k$-mer $v$ should be added to the MIS or not, we first compute the nearest
neighbors of $v$ and check if $v$ shares a mapping with them. If a shared
mapping is found, we store this mapping information for $v$ and directly
conclude that $v$ cannot be included in the MIS.

Algorithm 2 Improved Greedy Algorithm

initialize empty arrays $M$, $mapping$, and $d\sigma$ //$d\sigma[u]$ is later
used to store $\textrm{edit}(u,\sigma^{k})$, for each $\sigma\in\Sigma$

for each $k$-mer $v\in S_{k}$ do

$\mathit{isMapped}\leftarrow\text{false}$

for each nearest neighbor $u$ of $v$ do

if $\textrm{edit}\left(mapping[u],v\right)\leq d$ then

$\mathit{isMapped}\leftarrow\text{true}$

$mapping[v]\leftarrow mapping[u]$

break

end if

end for

if $\mathit{isMapped}==\text{true}$ then

continue

end if

$v_{\sigma}\leftarrow\textrm{edit}\left(v,\sigma^{k}\right)$,
$\forall\sigma\in\Sigma$

for each $k$-mer $u\in M$ do

if $\max_{\sigma\in\Sigma}\\{|v_{\sigma}-d\sigma[u]|\\}>d$ then

continue

else if $\min_{\sigma\in\Sigma}\\{v_{\sigma}+d\sigma[u]\\}\leq d$ or
$\textrm{edit}(u,v)\leq d$ then

$mapping[v]\leftarrow u$

$\mathit{isMapped}\leftarrow\text{true}$

break

end if

end for

if $\mathit{isMapped}==\text{false}$ then

add $v$ into $M$

$mapping[v]\leftarrow v$

$d\sigma[v]\leftarrow v_{\sigma}$, $\forall\sigma\in\Sigma$

end if

end for

output $M$

The second technique uses the estimated bounds of the edit distance. Recall
that $(S_{k},\textrm{edit})$ is a metric space, in particular, the triangle
inequality holds. Let $\sigma^{k}$ be the $k$-mer that is purely composed of
the character $\sigma$. For any two $k$-mers $u,v\in S_{k}$ and
$\sigma\in\Sigma$, we have

$\displaystyle\textrm{edit}(u,v)$
$\displaystyle\geq\left|\textrm{edit}\left(u,\sigma^{k}\right)-\textrm{edit}\left(\sigma^{k},v\right)\right|,$
(1) $\displaystyle\textrm{edit}(u,v)$
$\displaystyle\leq\textrm{edit}\left(u,\sigma^{k}\right)+\textrm{edit}\left(\sigma^{k},v\right).$
(2)

Note that the calculation of $\textrm{edit}(u,\sigma^{k})$ is simple: all
characters in $u$ that are not $\sigma$ should be substituted by $\sigma$
(insertions and deletions cannot reduce the number of edits needed). Hence,
the upper and lower bounds provided by inequalities (1) and (2) can be used as
a filter before performing the expensive calculation of $\textrm{edit}(u,v)$.
Using this idea, we modify the algorithm as follows. For each $k$-mer $u$ in
the constructed MIS, we store additional values $\textrm{edit}(u,\sigma^{k})$
for each $\sigma\in\Sigma$. When searching for a mapping of a $k$-mer $v$, for
each $k$-mer $u$ in the MIS, we check if the parameter $d$ is within the range
given by inequalities (1) and (2). If so, we calculate the exact distance
$\textrm{edit}(u,v)$ and compare it with $d$; otherwise, if $d$ is less than
$\max_{\sigma\in\Sigma}\left\\{\left|\textrm{edit}\left(u,\sigma^{k}\right)-\textrm{edit}\left(\sigma^{k},v\right)\right|\right\\}$,
we directly conclude that $\textrm{edit}(u,v)>d$ and therefore $u$ cannot be a
valid mapping of $v$; if $d$ is greater than
$\min_{\sigma\in\Sigma}\left\\{\textrm{edit}\left(u,\sigma^{k}\right)+\textrm{edit}\left(\sigma^{k},v\right)\right\\}$,
then we know that $\textrm{edit}(u,v)<d$ and $v$ cannot be included in the
MIS.

Algorithm 2 incorporates above two techniques. Note that Algorithm 2 is also
correct, i.e., the returned $M$ is guaranteed an MIS, following the
correctness of Algorithm 1 and above analysis. The worst-case running time of
Algorithm 2 is the same with Algorithm 1 but it runs faster in practice (see
Section 3).

###  2.3 A BFS-based Algorithm

In the third algorithm, we build a graph in which vertices correspond to all
$k$-mers and ($k-1$)-mers, and edges correspond to pairs of vertices with edit
distance being exactly 1 except that there is no edge between any two
$(k-1)$-mers. For example, if we take $k=4$, the graph contains all the 4-mers
and 3-mers. There is no edge between any pair of 3-mers. Two 4-mers are
connected if one can be obtained from the other by a single substitution. A
4-mer and a 3-mer are connected if the 4-mer can be obtained by inserting a
character to the 3-mer.

The key property of this graph is that two $k$-mers $u$ and $v$ have
$\textrm{edit}(u,v)=d$ if and only if the distance (the length of the shortest
path) between $u$ and $v$ in this graph is $d$. To formally see this, by the
construction of the graph, every path between two $k$-mers gives a valid
sequence of edits transforming one to the other. Hence, we only need to show
that given a sequence of $d$ edits that transforms a $k$-mer $u$ to a $k$-mer
$v$, the edits can be rearranged such that it corresponds to a path of length
$d$ from $u$ to $v$ in the graph. The only issue comes from the fact that the
given sequence of edits may contain several insertions (or deletions) in a row
so the intermediate strings may not be $k$-mers or $(k-1)$-mers. However,
because both $u$ and $v$ are $k$-mers, the sequence of edits must contain the
same number of insertions and deletions. Thus the sequence can be rearranged
so that each deletion is followed immediately by an insertion. Such a sequence
has a direct representation in the graph. For example, the transformation from
the 5-mer $TGATT$ to the 5-mer $ATTGA$ can be represented by the following
shortest path of length 4: $TGATT\rightarrow GATT\rightarrow GATTG\rightarrow
ATTG\rightarrow ATTGA$.

Following the above property, the calculation of all $k$-mers within an edit
distance of at most $d$ from a $k$-mer $u$ can be (efficiently) achieved by
_exploring_ $u$ using breadth-first search, i.e., traversing all vertices
reachable from $u$ within a distance of $d$. Again, if $u$ is already added to
the MIS, then all vertices found during the exploration can be marked as
mapped to $u$.

We do not need to fully explore every single $k$-mer. Instead, we can reuse
the information stored in exploring previous vertices to stop early.
Specifically, for each vertex $u$ (not in the current MIS), we store the
distance from $u$ to any vertex in the MIS, and update it in the exploring.
See Algorithm 3 for the complete pseudocode.

Algorithm 3 BFS-based Algorithm

initialize empty arrays $\mathit{distance}$ and $M$

$\mathit{distance}[v]\leftarrow\infty$ for all $k$-mers $v$

$\mathit{distance}[u]\leftarrow\infty$ for all ($k-1$)-mers $u$

for each $k$-mer $v\in S_{k}$ do

if $\mathit{distance}[v]==\infty$ then //explore $v$

add $v$ to $M$

$\mathit{distance}[v]\leftarrow 0$

initialize an empty queue $\mathit{frontier}$

$\mathit{frontier}.$add$(v)$

while $\mathit{frontier}$ is not empty do

$u\leftarrow\mathit{frontier}.$pop$()$

for each adjacent vertex $w$ of $u$ in the graph do

if $\mathit{distance}[w]>\mathit{distance}[u]+1$ then

$\mathit{distance}[w]\leftarrow\mathit{distance}[u]+1$

if $\mathit{distance}[w]<d$ then

$\mathit{frontier}.$add$(w)$

end if

end if

end for

end while

end if

end for

return $M$

The exploring step (variant of BFS) guarantees that all the unexplored
$k$-mers must have an edit distance greater than $d$ from all vertices in the
MIS, and all the explored $k$-mers must have an edit distance less than or
equal to $d$ from a vertex in the MIS. Hence, $k$-mers in the returned MIS are
at least $d+1$ edits apart (i.e., independent) and no other $k$-mers can be
added to the resulting MIS (i.e., maximal). This shows that Algorithm 3 is
correct.

To see its time complexity, note that the initial distance of a vertex is at
most $d$, hence it can be updated (decreased) at most $d-1$ times by the
subsequent explorings. In other words, each vertex can be explored at most $d$
times. Each time a $k$-mer is explored, we compute its neighbors with distance
$1$. The number of neighbors increases linearly with respect to $k$ (for a
constant-size alphabet). Thus, the running time for this algorithm is
$O(d\cdot k\cdot|V|)$, where $|V|=|\Sigma|^{k}+|\Sigma|^{k-1}$ is the total
number of $k$-mers and ($k-1$)-mers. The space complexity of this algorithm is
$O(|V|)$. For large values of $k$, the graph is dense and too large to be
stored in memory, we compute edges of the graph on the fly. The space is
mainly used to store the distance array, which grows linearly with respect to
the number of vertices.

##  3 Results

We implemented all three algorithms described above and conducted experiments
for $2\leq k\leq 15$ and $d<k$ (with the DNA alphabet $\Sigma=\\{A,C,G,T\\}$).
Table 1 reports the sizes of the resulting MISs which are depicted in Fig. 1.
For each row (i.e., with the same $d$ value), the size of the MIS increases
(approximately) geometrically with respect to $k$, which is consistent with
the geometric growth of the size of the $k$-mer space $S_{k}$. For each column
(i.e., with the same $k$ value), the size of the MIS decreases geometrically
with respect to the growth of $d$. As discussed before, $d$ specifies the
sparsity of the resulting MIS. With larger $d$ value, fewer vertices can be
selected into the MIS.

Table 2 reports the fastest algorithm to calculate the MIS for each
combination of $k$ and $d$. Table 2 together with Table 1 suggests a general
strategy to choose the best algorithm for a specific combination of $k$ and
$d$. Algorithm 1 is efficient for $d$ values within the interval $[k-4,k-1]$
since these $d$ values result in an MIS with less than 15 vertices
approximately. The BFS-based Algorithm (Algorithm 3) is the best for $d$
values between 1 and 4 because such a small $d$ leads to an MIS that is too
large to perform pairwise comparisons efficiently. The rest of the $d$ values
are better be handled by Algorithm 2.

Table 3 and Table 4 record the corresponding time and memory usage,
respectively. The rows of Table 3 show that the running time generally
increases geometrically with respect to the parameter $k$. The current
computation bottleneck is the case with $k=15$ and $d=5$, which takes
approximately 15 hours using Algorithm 2. Algorithm 2 always has the largest
memory overhead comparing with other two. The extra memory is mainly used to
store the mapping information of all $k$-mers. The peak memory usage is
approximately 2 GB for $k=15$ and $d\in\\{5,6,\ldots,10\\}$.

Table 1: The size of calculated MIS (by the fastest algorithm) for different
$k$ and $d$.

| $k$
---|---
|  | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
$d$ | 1 | 4 | 16 | 64 | 256 | 1024 | 4096 | 16384 | 65536 | 262144 | 1048576 | 4194304 | 16777216 | 67108864 | 268435456
2 |  | 4 | 12 | 36 | 96 | 311 | 1025 | 3451 | 11743 | 40604 | 141943 | 500882 | 1782677 | 6388106
3 |  |  | 4 | 8 | 20 | 57 | 164 | 481 | 1463 | 4574 | 14522 | 46908 | 153767 | 510118
4 |  |  |  | 4 | 4 | 14 | 34 | 90 | 242 | 668 | 1894 | 5517 | 16440 | 49992
5 |  |  |  |  | 4 | 4 | 12 | 25 | 57 | 133 | 338 | 879 | 2346 | 6486
6 |  |  |  |  |  | 4 | 4 | 10 | 17 | 38 | 79 | 188 | 448 | 1107
7 |  |  |  |  |  |  | 4 | 4 | 9 | 13 | 28 | 54 | 112 | 251
8 |  |  |  |  |  |  |  | 4 | 4 | 4 | 12 | 20 | 37 | 75
9 |  |  |  |  |  |  |  |  | 4 | 4 | 4 | 11 | 14 | 30
10 |  |  |  |  |  |  |  |  |  | 4 | 4 | 4 | 10 | 13
11 |  |  |  |  |  |  |  |  |  |  | 4 | 4 | 4 | 8
12 |  |  |  |  |  |  |  |  |  |  |  | 4 | 4 | 4
13 |  |  |  |  |  |  |  |  |  |  |  |  | 4 | 4
14 |  |  |  |  |  |  |  |  |  |  |  |  |  | 4

Figure 1: The size of an MIS (found by the fastest algorithm) with respect to parameters $k$ and $d$. Table 2: The fastest algorithm for each combination of $k$ and $d$. Numbers 1, 2, and 3 correspond to Algorithms 1, 2, and 3.  | $k$
---|---
|  | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
$d$ | 1 | 1 | 1 | 1 | 1 | 1 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3
2 |  | 1 | 1 | 1 | 1 | 1 | 3 | 3 | 3 | 3 | 3 | 3 | 3 | 3
3 |  |  | 1 | 1 | 1 | 1 | 2 | 3 | 3 | 3 | 3 | 3 | 3 | 3
4 |  |  |  | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 3 | 3 | 3
5 |  |  |  |  | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 2 | 2 | 2
6 |  |  |  |  |  | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2 | 2
7 |  |  |  |  |  |  | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
8 |  |  |  |  |  |  |  | 1 | 1 | 1 | 1 | 1 | 2 | 2
9 |  |  |  |  |  |  |  |  | 1 | 1 | 1 | 1 | 2 | 2
10 |  |  |  |  |  |  |  |  |  | 1 | 1 | 1 | 1 | 2
11 |  |  |  |  |  |  |  |  |  |  | 1 | 1 | 1 | 1
12 |  |  |  |  |  |  |  |  |  |  |  | 1 | 1 | 1
13 |  |  |  |  |  |  |  |  |  |  |  |  | 1 | 1
14 |  |  |  |  |  |  |  |  |  |  |  |  |  | 1
Table 3: Running time (seconds) of the fastest algorithm for each $(k,d)$. | $k$
---|---
|  | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
$d$ | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 5 | 25 | 114 | 528 | 3382 | 14028
2 |  | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 11 | 46 | 202 | 977 | 5444 | 22380
3 |  |  | 0 | 0 | 0 | 0 | 0 | 3 | 18 | 78 | 344 | 2161 | 10317 | 46701
4 |  |  |  | 0 | 0 | 0 | 0 | 1 | 8 | 55 | 460 | 4191 | 21801 | 36614
5 |  |  |  |  | 0 | 0 | 0 | 0 | 4 | 25 | 152 | 1409 | 11887 | 53478
6 |  |  |  |  |  | 0 | 0 | 0 | 3 | 17 | 94 | 659 | 4423 | 23922
7 |  |  |  |  |  |  | 0 | 0 | 1 | 13 | 77 | 504 | 3008 | 14864
8 |  |  |  |  |  |  |  | 0 | 1 | 7 | 52 | 438 | 2422 | 11308
9 |  |  |  |  |  |  |  |  | 0 | 5 | 31 | 258 | 2128 | 10003
10 |  |  |  |  |  |  |  |  |  | 4 | 23 | 210 | 1809 | 9197
11 |  |  |  |  |  |  |  |  |  |  | 21 | 187 | 1225 | 7322
12 |  |  |  |  |  |  |  |  |  |  |  | 184 | 988 | 5548
13 |  |  |  |  |  |  |  |  |  |  |  |  | 931 | 4743
14 |  |  |  |  |  |  |  |  |  |  |  |  |  | 4332
Table 4: Memory usage (kB) of the fastest algorithm for each $(k,d)$.

| $k$
---|---
|  | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15
$d$ | 1 | 1076 | 1076 | 1076 | 1076 | 1336 | 1332 | 1328 | 1324 | 3648 | 4568 | 8524 | 23916 | 85316 | 331164
2 |  | 1076 | 1076 | 1076 | 1076 | 1076 | 1324 | 1320 | 3544 | 4596 | 8364 | 24012 | 85396 | 331220
3 |  |  | 1076 | 1076 | 1076 | 1072 | 1476 | 1532 | 4008 | 5148 | 8888 | 24336 | 85984 | 331908
4 |  |  |  | 1076 | 1076 | 1076 | 1076 | 1624 | 5460 | 11508 | 36224 | 29512 | 92524 | 341752
5 |  |  |  |  | 1076 | 1072 | 1076 | 1648 | 5352 | 11468 | 36016 | 134416 | 527728 | 2100852
6 |  |  |  |  |  | 1072 | 1076 | 1072 | 1952 | 11584 | 36272 | 134484 | 527672 | 2100688
7 |  |  |  |  |  |  | 1076 | 1076 | 1880 | 2012 | 36164 | 134480 | 527688 | 2100568
8 |  |  |  |  |  |  |  | 1072 | 1920 | 1868 | 1952 | 1920 | 527772 | 2100568
9 |  |  |  |  |  |  |  |  | 1980 | 2012 | 3464 | 1920 | 527760 | 2100548
10 |  |  |  |  |  |  |  |  |  | 1980 | 1868 | 3456 | 1964 | 2100636
11 |  |  |  |  |  |  |  |  |  |  | 2032 | 3396 | 3184 | 3308
12 |  |  |  |  |  |  |  |  |  |  |  | 1976 | 3380 | 1992
13 |  |  |  |  |  |  |  |  |  |  |  |  | 3404 | 1964
14 |  |  |  |  |  |  |  |  |  |  |  |  |  | 3384

##  4 Conclusion and Discussion

We studied the problem of extracting an MIS as a representative substructure
from a $k$-mer space. Three algorithms are designed to efficiently solve the
problem for different ranges of parameters $k$ and $d$. The first one is a
simple greedy algorithm similar to the greedy graph for general graphs. The
second algorithm extends the first by implementing two techniques to avoid
some redundant comparisons. The third algorithm represents the edit distance
as a shortest path in an extended graph and uses a variant of BFS. Experiments
are done for $k$ up to 15. The computation bottleneck occurs at $k=15$ and
$d=5$ where the second algorithm performs the best. The corresponding peak
running time is approximately 15 hours, and the peak memory usage is about 2
GB.

For future work, we would like to extend this study to larger $k$-mer spaces.
Considering the current computation bottleneck, one potential improvement is
to design more heuristics to partition the $k$-mer space so that only a small
subset of $k$-mers will be involved in the expensive pairwise edit distance
computation.

Another idea is to revise the BFS-based algorithm. The extended graph with
both $k$-mers and $(k-1)$-mers is highly symmetric with respect to
permutations of the alphabet. For example, the search tree for the 5-mers
$AGAAC$ and $TATTG$ are isomorphic where $A$, $G$, and $C$ are replaced with
$T$, $A$, and $G$, respectively. This algorithm may take advantage of such
symmetries to explore the graph more efficiently. Identifying and utilizing
symmetric substructures of the $k$-mer space is of independent interest and
has applications in other related fields [17].

Regarding the practical applications of the constructed MISs, one potentiality
is their usage in $k$-mer-based reads clustering problems [6, 7], where the
choice of centers is crucial. Because an MIS contains a group of
representatives for the $k$-mer space that are guaranteed to be a certain edit
distance apart, it is ideal for generating even-sized clusters. Moreover, the
MISs can also be used to generate keys for $k$-mer-based reads hashing
problems [13, 15] where nearby $k$-mers are desired to share a hash value
(i.e., locality-sensitive hashing), which is particularly useful in aligning
error-prone long reads—a challenging yet unsolved problem. Last, it is
interesting to apply the constructed MISs as the set of representative
$k$-mers in sketching large-scale sequences.

## References

  * [1] C. Titus Brown and Luiz Irber.  sourmash: a library for minhash sketching of dna.  Journal of Open Source Software, 1(5):27, 2016.  doi:10.21105/joss.00027.
  * [2] Andrea Califano and Isidore Rigoutsos.  FLASH: A fast look-up algorithm for string homology.  In Proceedings of IEEE Conference on Computer Vision and Pattern Recognition, pages 353–359. IEEE, 1993.
  * [3] Alexander Dobin, Carrie A Davis, Felix Schlesinger, Jorg Drenkow, Chris Zaleski, Sonali Jha, Philippe Batut, Mark Chaisson, and Thomas R Gingeras.  Star: ultrafast universal rna-seq aligner.  Bioinformatics, 29(1):15–21, 2013.
  * [4] Luiz Irber, Phillip T. Brooks, Taylor Reiter, N. Tessa Pierce-Ward, Mahmudur Rahman Hera, David Koslicki, and C. Titus Brown.  Lightweight compositional analysis of metagenomes with fracminhash and minimum metagenome covers.  bioRxiv, 2022.  doi:10.1101/2022.01.11.475838.
  * [5] Miten Jain, Sergey Koren, Karen H Miga, Josh Quick, Arthur C Rand, Thomas A Sasani, John R Tyson, Andrew D Beggs, Alexander T Dilthey, Ian T Fiddes, et al.  Nanopore sequencing and assembly of a human genome with ultra-long reads.  Nature Biotechnology, 36(4):338–345, 2018.
  * [6] Chang Sik Kim, Martyn D Winn, Vipin Sachdeva, and Kirk E Jordan.  $k$-mer clustering algorithm using a mapreduce framework: application to the parallelization of the inchworm module of trinity.  BMC Bioinformatics, 18(1):1–15, 2017.  doi:10.1186/s12859-017-1881-8.
  * [7] Jessime M Kirk, Susan O Kim, Kaoru Inoue, Matthew J Smola, David M Lee, Megan D Schertzer, Joshua S Wooten, Allison R Baker, Daniel Sprague, David W Collins, et al.  Functional classification of long non-coding rnas by $k$-mer content.  Nature Genetics, 50(10):1474–1482, 2018.  doi:10.1038/s41588-018-0207-8.
  * [8] Robert Krauthgamer and Yuval Rabani.  Improved lower bounds for embeddings into $l_{1}$.  SIAM Journal on Computing, 38(6):2487–2498, 2009.
  * [9] Vladimir I Levenshtein et al.  Binary codes capable of correcting deletions, insertions, and reversals.  In Soviet Physics Doklady, volume 10, no. 8, pages 707–710. Soviet Union, 1966.
  * [10] Heng Li.  Minimap2: pairwise alignment for nucleotide sequences.  Bioinformatics, 34(18):3094–3100, 2018.
  * [11] Bin Ma, John Tromp, and Ming Li.  Patternhunter: faster and more sensitive homology search.  Bioinformatics, 18(3):440–445, 2002.
  * [12] Denise Mak, Yevgeniy Gelfand, and Gary Benson.  Indel seeds for homology search.  Bioinformatics, 22(14):e341–e349, 2006.
  * [13] Hamid Mohamadi, Justin Chu, Benjamin P Vandervalk, and Inanc Birol.  nthash: recursive nucleotide hashing.  Bioinformatics, 32(22):3492–3494, 2016.  doi:10.1093/bioinformatics/btw397.
  * [14] Rafail Ostrovsky and Yuval Rabani.  Low distortion embeddings for edit distance.  Journal of the ACM (JACM), 54(5):23–es, 2007.
  * [15] Enrico Petrucci, Laurent Noé, Cinzia Pizzi, and Matteo Comin.  Iterative spaced seed hashing: closing the gap between spaced seed hashing and $k$-mer hashing.  Journal of Computational Biology, 27(2):223–233, 2020.  doi:10.1089/cmb.2019.0298.
  * [16] Anthony Rhoads and Kin Fai Au.  PacBio sequencing and its applications.  Genomics, Proteomics & Bioinformatics, 13(5):278–289, 2015.
  * [17] Klaus U Schulz and Stoyan Mihov.  Fast string correction with levenshtein automata.  International Journal on Document Analysis and Recognition, 5(1):67–85, 2002.  doi:10.1007/s10032-002-0082-8.


11institutetext:  Department of Astronomy, Beijing Normal University, Beijing
100875, P. R. China<EMAIL_ADDRESS>

# Four New Observational $H(z)$ Data From Luminous Red Galaxies of Sloan
Digital Sky Survey Data Release Seven

Cong Zhang 11 Han Zhang 11 Shuo Yuan 11 Siqi Liu 11 Tong-Jie Zhang and Yan-
Chun Sun 1111

###### Abstract

By adopting the differential age method, we utilize selected 17832 luminous
red galaxies (LRGs) from Sloan Digital Sky Survey Data Release Seven (SDSS
DR7) covering redshift $0<z<0.4$ to measure Hubble parameters. Using a full
spectrum fitting package UlySS, these spectra are reduced with single stellar
population (SSP) models and optimal age information of our selected sample are
derived. With the decreasing age-redshift relation, four new observational
$H(z)$ data (OHD) points are obtained, which are $H(z)=69.0\pm 19.6$ km s-1
Mpc-1 at $z=0.07$, $H(z)=68.6\pm 26.2$ km s-1 Mpc-1 at $z=0.12$,
$H(z)$=$72.9\pm 29.6$ km s-1 Mpc-1 at $z=0.2$ and $H(z)$=$88.8\pm 36.6$ km s-1
Mpc-1 at $z=0.28$, respectively. Combined with other 21 available OHD data
points, a performance of constraint on both flat and non-flat $\Lambda$CDM
model is presented.

###### keywords:

Cosmology:cosmological parameters -Cosmology: observations -Galaxies:
evolution

##  1 Introduction

A variety of cosmological observations are used for a better understanding of
the expansion of the Universe quantitatively, for example the mapping of the
cosmic microwave background (CMB) anisotropies (Spergel et al., 2007; Komatsu
et al., 2011), the measurement of baryon acoustic oscillation (BAO) peaks
(Eisenstein et al., 2005; Percival et al., 2010), measurements of ‘standard
candles’ such as the redshift-distance relationship of type Ia supernovae
(SNIa) (SNIa; Riess et al., 1998; Hicken et al., 2009) and gamma-ray bursts
(GRBs) (GRBs; Ghirlanda et al., 2004; Li et al., 2008) . Hubble parameter
$H(z)$, which is defined as : $H(z)=\dot{a}/a$, where $a$ denotes the cosmic
scale factor and $\dot{a}$ is its rate of change with respect to the cosmic
time, is directly related to the expansion history of the Universe. The method
based on the observational $H(z)$ data (OHD) has been used to test
cosmological models (e.g., Yi & Zhang, 2007; Chen & Ratra, 2011). Besides
parameters constraints, OHD can also be used as an auxiliary model selection
criterion (Li et al., 2009).

In practice, the Hubble parameter $H(z)$ is usually evaluated as a function of
the redshift $z$, with $a(t)/a(t_{0})=1/(1+z)$, where $t_{0}$ is the current
cosmic time:

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

$z$ is the cosmological redshift and $t$ is the age of the Universe when the
observed photon is emitted. Derivative of redshift with respect to cosmic
time, $dz/dt$ has a direct determination on $H(z)$. $H(z)$ has been measured
through the differential method according to Eq. (1), which was first put
forward by Jimenez & Loeb (2002). This differential method has been
demonstrated in Jimenez et al. (2003). However, it may be difficult to select
galaxies as ‘cosmic chronometers’ and determine the accurate age of a galaxy
considering stars in a galaxy born continuously, and a young stellar
population may dominate their emission spectra (Zhang et al., 2010). Luminous
red galaxies (LRGs) which have photometric properties consistent with an old,
passively evolving stellar population (Roseboom et al., 2006) are regarded as
a good candidate of this ‘cosmic chronometers’ (Crawford et al., 2010).

We employ ULySS111ULySS is available at: http://ulyss.univ-lyon1.fr/, an
available package on-line to fit a full-length spectrum. As the relatively
homogeneous stellar component of LRGs, we use single stellar population (SSP)
fitting and gain age information, which will be discussed in detail in Section
3. With the age-redshift relation, four OHD points are deduced accordingly.

There are already 21 OHD points got from both differential age method (Jimenez
et al., 2003; Simon et al., 2005; Stern et al., 2010; Moresco et al., 2012)
and baryon acoustic oscillation (BAO) method (Gaztañaga et al., 2009) .
Currently the number of OHD points are still scarce compared with SNIa
luminosity distance data. The potential power of OHD in constraining
cosmological parameters is explored in Ma & Zhang (2011) in detail. It has
achieved a conclusion that the constraining power of OHD can be as strong as
that of SNIa when its quantity reaches a certain value which, depending on the
error model used in that paper, is 64, thus it is significant to gain new
independent Hubble parameters.

This paper is organized as follows. We briefly give our LRGs selection sample
in Section 2. In Section 3, we describe explicitly how we gain the age
information from the LRG spectra using ULySS and in Section 4 we present our
method to get the OHD from the galaxy ages. A cosmology constraint using all
available OHD including our 4 new ones is given in Section 5. Finally, in the
last section, we discuss the limitation and prospect of our results.

##  2 Source Selection of LRG Sample

It is important and necessary to select a large homogeneous passively evolving
sample of LRGs to obtain the age-redshift relation. The SDSS (York et al.,
2000; Stoughton et al., 2002; Abazajian et al., 2003, SDSS) is currently the
largest photometric and spectroscopic sky survey, which include five-band
imagine over 104 deg2 with accurate photometric calibration and spectroscopy
of 106 galaxies(Abazajian et al., 2009). The SDSS spectroscopic survey
consists of two samples of galaxies selected with different criteria, which
are named the MAIN sample (Strauss et al., 2002) and the LRG sample
(Eisenstein et al., 2003) respectively. The wavelength of these galaxy spectra
cover the range from 3800Åto 9200Åwith spectral resolution
$\lambda/(\Delta\lambda)$ = 1850-2200, and with dedicated software the spectra
are automatically reduced, which flux calibrates the spectra and references
them to the heliocentric frame and converts to vacuum wavelengths.

We chose the sample from LRG selection criterion of SDSS DR7, which is:
i)Selecting galaxies from the Catalog Archive Server (CAS) database with the
TARGET_GALAXY_RED flag. ii)The S/N of the $r$-waveband should be greater than
10. iii)The restrictions, which are that specClass = ‘SPEC_GALAXY’, that zStar
= ‘XCORR_HIC’, that zWarning = 0, that eClass $<$ 0, that z $<$ 0.4 and that
fracDev_r $>$ 0.8, should also be satisfied. The LRG selection criterion of
SDSS (Eisenstein et al., 2001) is on the basis of color and magnitude to yield
a sample of luminous, intrinsically red galaxies. However, the sample selected
from SDSS according to Eisenstein’s selection criterion are not very
homogeneous. Furthermore, we use the sample from Carson & Nichol’s sample
(Carson & Nichol, 2010) with a restraint on signal to noise ratio (SNR). As
SNR has an impact on our next fitting step, it is necessary to demand SNR of
R-band to be greater than 10. The 17832 selected quiescent, luminous red
galaxies basing on ‘Carson & Nichol sample’ method from SDSS DR7 cover the
redshift range from 0-0.4. Redshift distribution of our sample is shown in
Fig. 1.

Figure 1: Redshift distribution of 17832 LRGs.

Compared with other spectra selection work, ‘Carson & Nichol sample’
calibrates, for the first time, the SDSS spectra on Lick/Image Dissector
Scanner (IDS) system. Their general selected criterion steps are as follows.
First step is to get LRGs from Catalog Archive Server (CAS) as outlined in
Eisenstein et al. (2001). Further, restraining the spectrum using standard
emission lines such as H$\alpha$, H$\beta$ and OIII 5007 to increase the
number of truly quiescent galaxies in the previous sample. To fully picture
physical properties such as velocity dispersion and absolute luminosity of
these LRGs, through correcting velocity dispersions for aperture effects and
performing $K+e$ corrections to the magnitude, four subsamples are produced
with different absolute magnitude and velocity dispersion. For a explicit
description of this method, please refer to Carson & Nichol (2010).

##  3 Age-Redshift Relation

We proceed to describe our way of obtaining age information of LRGs. There are
many methods to learn the age and metallicity of stellar systems from a
spectrum, such as SED fitting, spectrophotometric indices (e.g., Lick, Rose
indices) and full spectrum fitting (Koleva et al., 2008). In this paper, we
adopt the full spectrum fitting method to analyse the physical properties of
stellar populations. Full spectrum fitting, which makes use of all the
information contained in the signal, is insensitive to extinction or flux
calibrating errors and independently from the shape of continuum. We adopt an
open-source software, ULySS, to explore the history of stellar populations.

ULySS is a full spectrum fitting package developed by a group at
Universit$\acute{e}$ de Lyon (Koleva et al., 2009b). Its principle is to seek
the minimum $\chi^{2}$ in the process of fitting an observed spectrum with a
model spectrum in the pixel space with MPFIT function. When fitting the
observed spectrum ($F_{obs}(\lambda)$), the package uses a linear combination
of $k$ non-linear components ($CMP_{i}$) with weight $W_{i}$ respectively to
approximate it. In this process, the composite model is possibly convolved
with a line-of-sight velocity distribution (LOSVD), multiplied by an $n^{th}$
order polynomial of $P_{n}(\lambda)$ and summed to another polynomial of
$Q_{m}(\lambda)$ (for more details please refer to Koleva et al., 2009b):

$F_{obs}(\lambda)=P_{n}(\lambda)\times\\{LOSVD(v_{sys},\sigma,h3,h4)\otimes\sum_{i=0}^{i=k}W_{i}CMP_{i}(a_{1},a_{2},a_{3},...,\lambda)\\}+Q_{m}(\lambda).$
(2)

For the study of stellar population, the CMPi is characterized by age and
[Fe/H]. It uses Levenberg-Marquardt routine to evaluate parameters in
individual CMPi and the coefficients of $P_{n}(\lambda)$ and $Q_{m}(\lambda)$
(Koleva et al., 2009b). This method has already been successfully tested in Wu
et al. (2011). The reliability and robustness of using ULySS to the study of
the history of stellar population has been verified (e.g., Koleva et al.,
2009c, a).

###  3.1 Model Selection and Matching resolution

There exist several population models and in our paper, three of them provided
by the package ULySS, which are Pegase-HR/ELODIE3.1, Galaxev/STELIB (hereafter
BC03) and Vazdekis/Miles, are tested. The information of these models is
listed in Tab. 1. Koleva et al. (2008)also test these three models and
verifies their reliability.

The full spectrum fitting uses the redundancy of the spectrum and the
multiplicative polynomial could decrease the influence of flux calibration and
Galactic extinction. On the other aspect, the fitting method is more sensitive
to the wavelength range of the spectrum. The spectra from SDSS cover the
wavelength from 3800Åto 9200Å, as we could see, only the wavelength range of
BC03 could cover the whole wavelength of spectra from SDSS. After comparison
of the three different models, we choose BC03 as the reliable model for use.

Table 1: Information in Three Models

Model Library Resolution /Å Wavelength /Å Age/Gyr [Fe/H]/dex IMF Pegase
ELODIE3.1 0.55 4000-6800 0.1-20 -3.21-1.62 Salpeter Galaxev STELIB 3 3200-9500
0.1-20 -2.3-0.4 Chabrier Vazdekis MILES 2.3 3525-7500 0.1-17.5 -1.7-0.2
Salpeter

The first step in fitting is to match resolutions between the observed spectra
with the model we chose. There are two ways for matching, either by
transforming the resolution of the model or the observed spectrum. In our
paper, we choose to transform the model provided by ULySS by injecting
relative line spread function (LSF) between our spectrum and the model. When
determining the LSF, we make use of the available velocity dispersion template
stars as a standard star from http://www.sdss.org/dr7/algorithms/veldisp.html.
We make a five times linear interpolation in wavelength and convolve with the
model. In this way, a new matching model is generated, and all this process
could be accomplished through ULySS function (Koleva et al., 2009b).

###  3.2 Single Stellar Population Fit

It is essential to study the stellar population of galaxies if we want to
reconstruct the star formation history (SFH) of galaxies. We fit the spectrum
with a single SSP for the following reasons. Firstly, LRG is believed to be
drawn from a same parent population with the most percent of their stars are
formed from a single-burst. Secondly, SSP-equivalent properties correspond to
the ‘luminosity-weighted’ average over the distributions. Thus by a single SSP
fitting, a general view of a galaxy could be gotten (Du et al., 2010). We
carry out an initial study of age-dating using the stellar population models
of (Bruzual & Charlot, 2003; Stern et al., 2010) to synthesize spectra, which
is also provided by ULySS.

But when interpreting the galaxy spectra, we find both of age and metallicity
guess value tend to affect the fitting result. It is difficult to remove such
influence especially at the low resolving spectra (Du et al., 2010). The
fitting of ULySS starts from a point (initial value) in the parameter
space(age and metallicity), therefore we stress the importance of that point
as the returning fitting value will be initial value dependent and the
minimization will be local. In order to identify and understand the age-
metallicity degeneracy and local minimizations, we analyse that age-
metallicity degeneracy and construct $\chi^{2}$ map which is a function
provided by ULySS. The $\chi^{2}$ map is a visualization of parameters space,
and the principle of this map is to give a grid of nodes in a 2-D projection
of the parameter space. $\chi^{2}$ map returns the minimization of each node
(Koleva et al., 2009b). A subsample of 164 LRG spectra from redshift 0.03 to
0.179 is dealt with by building their $\chi^{2}$ maps. It has been discovered
that the metallicity of LRGs are consistent with the assumption that selected
galaxies are thought to be have similar metallicities (Jimenez et al., 2003)
and it floats between about 0.1 dex and 0.2 dex.

Figure 2: Best fit with BC03 model and the residual spectrum for a galaxy. The
top panel shows the spectrum in black and the best fit in blue. The red
regions are rejected from the fit, the indigo line is the multiplicative
polynomial. The bottom panel shows the residuals from the best fit. The
continuous green lines represent the 1$\sigma$ deviation.

ULySS could provide initial value vectors as initial value, that is to say, if
giving points intensive enough, it is possible to find the global minimization
and thus break the degeneracy. We set our metallicity range guess points from
0.1-0.2 dex to every galaxy in our sample, while age guess point span from
5000 Myr-11200 Myr. Fig. 2 shows the best fit with BC03 model and the residual
spectrum for a galaxy when a specific group of initial values is given.

For each galaxy, groups of age and metallicity values can be gained as
results. We ignore the metallicity and just focus on the age. To identify the
final accurate age of each galaxy, we base our selections on the minimal
$\chi^{2}$ criterion. In particular, if all fitting ages for a single galaxy
surpass the age of our universe, it is not convincing as there are several
reasons causing this situation such as the low resolution of spectra and model
dependence. We set a limit of that $a+\sigma<16Gyr$, where the $a$ is the
fitting age and the $\sigma$ is the error, in the fitting process to exclude
the situation that the age of a galaxy excess the age of the Universe. The
fitting result is shown in Fig. 3. For clarity, the ages less than 7Gyr have
not been plotted. In this figure, a clear tend that galaxies age decrease with
the redshift is displayed.

Figure 3: Fitting results of the 17832 LRGs. For clarity, the ages less than
7Gyr have not been plotted. The red line shows the theoretical age of the
universe $t(z)$ for a $\Lambda$CDM model with $\Omega_{m}=0.29$ and $H_{0}=69$
km s-1 Mpc-1 . The blue line indicates $t(z)+3$Gyr, whose derivation of the
3Gyr comes from the systematic errors. A clear trend is present: the oldest
ages of the galaxies decrease with redshift.

##  4 The determination of OHD

According to the Eq. 1, the slope of the linear fit of data of $t(z_{i})$,
which is the oldest age at redshift $z_{i}$, relates to the Hubble parameter
directly. That is, the Hubble parameter at $z_{eff}$, $H(z_{eff})$, can be
calculated by formula $H(z_{eff})=-[1/(1+z_{eff}](\Delta t/\Delta z)^{-1}$,
where $z_{eff}=(z_{max}+z_{min})/2$.

Mathematically, the value of $\Delta z=z_{max}-z_{min}$ should not be too
large, so we divide our data into four redshift regions, $0.03\leq z\leq
0.11$, $\ 0.08\leq z\leq 0.16$, $0.16\leq z\leq 0.24$, and $0.24\leq z\leq
0.32$.

To calculate the $t(z_{i})$ of each subsample, we adopt a common bin-dividing
method. Take the first subsample for example. We divide it into several bins
from 3 to 20, select the oldest galaxy of each bin as a set of $t(z_{i})$ and
fit the corresponding points of $t(z_{i})$ with a straight line to get a
candidate of OHD of this subsample. Here we get $18$ candidates of OHD of this
subsample. The number of bins corresponding to the selected candidate is named
as $n_{best}$. Three parameters are used here for selecting the most suitable
number of bins that we choose: SSE$/n$ for the average of the Sum of Squares
for Error, $\sigma_{slope}/slope$ for the relative error of the fitting result
of the slope and $P_{\geq\chi^{2}}$ for the goodness of fit. The left
subsamples follow the same method.

Here we explain our reason of choosing the above three parameters as
criterion. Imagine that an envelop line can be obtained naturally with ideal
data, therefore, we should take the bin-dividing number $n$ as large as
possible. Unfortunately, for the real case, the envelop fluctuates strongly.
This is due to the existence of ‘fake-oldest galaxies’, which formed too late
to be considered as ‘cosmic chronometer’ compared with the LRGs at other
redshifts. Therefore, large $n$ would increase the risk of selecting the
‘fake-oldest galaxies’ as ‘cosmic chronometer’. That is what the parameter of
SSE$/n$ can indicate, with value $0$ representing the ideal case. Conversely,
too small $n$ represents large statistic errors, which is evaluated by
parameter $\sigma_{slope}/slope$. Moreover, $P_{\geq\chi^{2}}$ is calculated
to represent the goodness of fit. Taking the three parameters into account, we
finally obtain the value of the $n_{best}$ of each subsample.

In every subsample, we expect the $n_{best}$ satisfies the the smallest
SSE$/n$, the largest $P_{\geq\chi^{2}}$ and the smallest absolute value of
$\sigma_{slope}/slope$. In the first subsample, $n=6$ case meet the criterion
perfectly. So does $n=7$ in the second subsample and $n=6$ in the third
subsample. For the fourth case, judging the $\sigma_{slope}/slope$ first, we
can see for both $n=10$ and $n=12$ share the least value. Then, considering
the second parameter SSE$/n$, we chose $n=12$ for smaller value. Besides, its
$P_{\geq\chi^{2}}$ is also satisfied.

Fig.4 shows the oldest ages in each bin when we divide the subsamples into
their corresponding $n_{best}$ bins and their optimal fits are also plotted.

Then we will introduce how we get the best parameters of linear fitting and
their error bars.

When we fit $n$ set of data ($z_{i}$, $t_{i}$) with a straight line $t=kz+b$
and the error of $t_{i}$ is $\sigma_{i}$ ($i=1,2,...,n$), the best parameters
of the linear fitting are obtained by minimizing the $\chi^{2}$:

$\chi^{2}=\sum_{i=1}^{n}\frac{(t_{i}-kz_{i}-b)^{2}}{\sigma_{i}^{2}},$ (3)

To get the error of $\sigma_{k}$, we must rewrite the Eq. 3 as follows so that
we can use some formulae of the regressive parameter errors:

$\chi^{2}=\sum_{i=1}^{n}({\frac{t_{i}}{\sigma_{i}}-k\frac{z_{i}}{\sigma_{i}}-b\frac{1}{\sigma_{i}}})^{2}.$
(4)

This form of $\chi^{2}$ is same as that when we fit the data of
($t_{i}/\sigma_{i}$,$z_{i}/\sigma_{i},1/\sigma_{i}$) with linear function of
$\frac{t}{\sigma}=k\frac{z}{\sigma}+b\frac{1}{\sigma}$ with the same method of
minimum $\chi^{2}$.

Then, using the well known formula on confidence intervals of regression
coefficients at a confidence level of $1-\alpha$ (seeing any textbook on
regression such as He & Liu, 2011, for detail), we can get the error of the
$k$ which can be regarded as $\sigma_{k}$ (i.e. the $\sigma_{slope}$ in Tab.
2):

$\sigma_{k}=\sigma_{slope}=|t_{(\alpha/2,\
n-2)}|\sqrt{\frac{\textrm{S(1,1)}}{n-2}\times\textrm{SSE}},$ (5)

where $|t_{(\alpha/2,\ n-2)}|$ is the absolute value of the inverse of
Student’s t-CDF (Cumulative Distribution Function) with degrees of freedom
$n-2$ for the corresponding probabilities in $\alpha/2$, which can be
calculated with function of tinv in MATLAB, SSE is the Sum of Squares for
Error and the $S(1,1)$ is the element in the first row and the first column in
the inverse of the following matrix:

$\left(\begin{array}[]{cc}\sum\limits_{i=1}^{n}\frac{z_{i}^{2}}{\sigma_{i}^{2}},&\sum\limits_{i=1}^{n}\frac{z_{i}}{\sigma_{i}^{2}}\\\
\\\
\sum\limits_{i=1}^{n}\frac{z_{i}}{\sigma_{i}^{2}},&\sum\limits_{i=1}^{n}\frac{1}{\sigma_{i}^{2}}\end{array}\right).$
(6)

In addition, the steps above can be completed by MATLAB Toolboxes easily.

With the Eq. 1, we can get the relation between the error of $H(z)$ and the
$\sigma_{slope}$:

$\sigma_{H}=\frac{1}{1+z_{eff}}\frac{1}{\sigma_{slope}^{2}},$ (7)

with which the error of $H(z)$ is calculated finally.

Table 2: Fitting results

$n$ 4 5 6 7 8 5 6 7 8 9 $-\sigma_{slope}/slope$ 0.64 0.44 0.28 0.31 0.28 0.49
1.45 0.38 0.54 0.74 SSE$/n$ 0.053 0.076 0.049 0.106 0.094 0.251 0.964 0.319
0.802 1.575 $P_{\geq\chi^{2}}$† 0.900 0.945 0.990 0.980 0.993 0.969 0.915
0.997 0.992 0.980 $n$ 4 5 6 7 8 10 11 12 13 14 $-\sigma_{slope}/slope$ 1.08
0.43 0.41 0.50 0.44 0.41 0.43 0.41 0.43 0.43 SSE$/n$ 0.028 0.019 0.017 0.034
0.031 0.208 0.113 0.175 0.163 0.211 $P_{\geq\chi^{2}}$† 0.945 0.992 0.999
0.999 0.999 1.000 1.000 1.000 1.000 1.000

†Too big $P_{\geq\chi^{2}}$ in these two subsamples result from too great
errors. The two other parameters
are taken account of mainly because of the very close $P_{\geq\chi^{2}}$ in
this subsample.

Figure 4: The oldest ages in each bin when we divide the subsamples into their
corresponding $n_{best}$ bins and their optimal fitted results. The solid line
represents the best fitting for each subsample. In the first subsample
$z_{min}=0.033$, $z_{max}=0.109$, and $z_{eff}=0.07$. In the second subsample,
$z_{min}=0.090$, $z_{max}=0.156$, $z_{eff}=0.12$. In the third subsample,
$z_{min}=0.170$, $z_{max}=0.236$, $z_{eff}=0.20$, and in the forth subsample,
$z_{min}=0.243$, $z_{max}=0.315$, $z_{eff}=0.28$. Figure 5: All available OHD
points. The solid line plots the theoretical Hubble parameter $H_{fid}$ as a
function of $z$ from the spatially flat $\Lambda$CDM model with
$\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}$ = 72 km s-1Mpc-1. The
OHD points are listed in Tab. 3 Table 3: All available OHD

$z$ $H(z)^{\dagger}$ $\sigma_{H(z)^{\dagger}}$ Ref. 0.090 69 12 Jimenez et al.
(2003) 0.170 83 8 Simon et al. (2005) 0.270 77 14 Simon et al. (2005) 0.400 95
17 Simon et al. (2005) 0.900 117 23 Simon et al. (2005) 1.300 168 17 Simon et
al. (2005) 1.430 177 18 Simon et al. (2005) 1.530 140 14 Simon et al. (2005)
1.750 202 40 Simon et al. (2005) 0.480 97 62 Stern et al. (2010) 0.880 90 40
Stern et al. (2010) 0.179 75 4 Moresco et al. (2012) 0.199 75 5 Moresco et al.
(2012) 0.352 83 14 Moresco et al. (2012) 0.593 104 13 Moresco et al. (2012)
0.680 92 8 Moresco et al. (2012) 0.781 105 12 Moresco et al. (2012) 0.875 125
17 Moresco et al. (2012) 1.037 154 20 Moresco et al. (2012) 0.24 79.69 3.32
Gaztañaga et al. (2009) 0.43 86.45 3.27 Gaztañaga et al. (2009) 0.07 69.0 19.6
†† 0.12 68.6 26.2 †† 0.20 72.9 29.6 †† 0.28 88.8 36.6 ††

† The unit is km s-1Mpc-1.
†† Work in this paper.

##  5 Cosmological constraints from OHD

Now we have totally 25 available OHD which are listed in Tab. 3. Using these
25 OHD, we constrain the cosmological parameters. The best fit parameters of
the model via OHD are determined by minimizing the $\chi^{2}$:

$\chi^{2}_{OHD}(\mathbf{p_{model}})=\sum\limits_{i=1}^{N}{\frac{(H_{obs}(z_{i})-H_{th}(z_{i}))^{2}}{\sigma_{obs,i}^{2}}},$
(8)

where $\mathbf{p_{model}}$ is a free parameter. Based on the basic equations
of $\Lambda$CDM model, we have the following two models: for a flat
$\Lambda$CDM($\Omega_{k}=0$),
$H_{th}(z)=H_{0}\sqrt{\Omega_{m}(1+z)^{3}+(1-\Omega_{m})}$ with
$\mathbf{p_{flat}}=(\Omega_{m},H_{0})$; for a non-flat $\Lambda$CDM,
$H_{th}(z)=H_{0}\sqrt{\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}+(1-\Omega_{\Lambda}-\Omega_{m})(1+z)^{2}}$,
with $\mathbf{p_{non-flat}}=(\Omega_{\Lambda},\Omega_{m},H_{0})$. For both
models, the likelihood function can be written as $L\propto\exp$
$(-\chi_{OHD}^{2}/2)$.

We use MCMC method to calculate the likelihood in the independent parameter
space. The Markov chains are generated and analysed via the Python MCMC code –
Pymc. The parameter $H_{0}$ and the density parameter
$\mathbf{\Omega}=(\Omega_{m},\ \Omega_{\Lambda})$ are treated as independent
parameters in these two models (Wei, 2010), and all of the prior distributions
of them are set as the uniform distribution (Tab. 4).

Table 4: The Prior of Model Parameters

Model Parameters Prior Distribution $\Omega_{m}$ Uniform(0.0,1.5)†
$\Omega_{\Lambda}$ Uniform(0.0,2.5)† $H_{0}$ Uniform(50,100) †

† _Uniform(lower limit,upper limit)_ stand for an uniform distribution in the
interval [lower limit,upper limit], and the same below.

Fig. 6 and Fig. 7 show the one-dimensional marginalized probability
distribution for each parameter at the diagonal positions and the two-
dimensional marginalized confidence regions at the non-diagonal positions of
the flat $\Lambda$CDM model and non-flat $\Lambda$CDM model respectively. The
$3\sigma$-s confidence regions ($68.7{\%},\ 95.45{\%},\ 99.73{\%}$ ) for each
parameter are calculated carefully. The best fit parameters are pointed out in
the contour plots with a long vertical line and a long horizontal line across.

To compare the constraining results of the ‘new version OHD’ with the new
points we obtained in this paper added and the ‘old version OHD’ (see Tab. 3
for detail), we display both results of the same model in Fig. 6 and Fig. 7
with different line colors to distinguish: the red lines referring to the ‘new
version OHD’ and the blue ones means the ‘old version OHD’.

Figure 6: OHD constraint region of flat $\Lambda$CDM. Red lines represent the
‘new version’ of OHD, with the blue ones for the ‘old version’. The diagonal
plots show the 1-D PDF and the vertical lines from center to edge match the
$3\sigma$-s intervals; The non-diagonal plots show the 2-D confidence regions
with contours for inner to outer match the $3\sigma$-s levels, and the open
blue circle points out the best fit point by ‘old version’ while the red cross
shows the ones via ‘new version’.  Figure 7:  OHD constraint region of non-
flat $\Lambda$CDM. Red lines represents the ‘new version’ of OHD, with blue
ones for the ‘old version’. The diagonal plots show the 1-D PDF and the
vertical lines from center to edge match the $3\sigma$-s intervals; The non-
diagonal plots show the 2-D confidence regions with contours for inner to
outer match the $3\sigma$-s levels, and the open blue circle points out the
best fit point by ‘old version’ while the red cross shows the ones via ‘new
version’.

The constraints in Fig. 6 and Fig. 7 show the well performances on
constraining $\Lambda$CDM. In Fig 6 and Fig 7 the 1-D marginalized constraints
are more stringent than the old one and for the same model the contour plots
against the same parameters pair have more shrunk confidence region at the
same level. Since the number of new added-on points are so small, the degree
of shrinking is weak as well.

##  6 Conclusion and Discussion

In this paper, we present our measurements of four new OHD data points from
the ages of passively-evolving galaxies at redshift $0<$z$<0.4$. A large
sample of LRGs spectra has been fitted by us with SSP models and a age-
redshift relation is obtained and displayed. By computing the relative ages of
these LRGs, we gain four new OHD data points. Combining them with other
available 21 OHD points, we constrain cosmological parameters using the
updated dataset of OHD. It should be mentioned the similar work Liu et al.
(2012), in which they used the SDSS DR7 to constrain the $H_{0}$, a particular
Hubble parameter at $z=0$. We hope to give a tighter constraint on cosmology
parameters with our new OHD data points. Unfortunately, these four points do
limited improvements for constraining cosmological parameters more accurately,
because of the relatively large error bar of these points.

Here, we explain the possible reasons. Firstly, the low SNR of spectra from
SDSS leads to the large uncertainty when calculating the age of the galaxies.
Combining spectra from various observation project may be a solution, as
discussed in Jimenez et al. (2003); Stern et al. (2010). Comparing our points
with the previous, especially those in Simon et al. (2005) and Moresco et al.
(2012), although there really exists a big difference in terms of the
accuracy, we would like to illustrate as followings. One is a theoretically
good method to obtain the age of LRG in Moresco et al. (2012). Still, there
are many problems in the method to be solved. The other is the results in
Simon et al. (2005) with more accurate than ours because of the high quality
of their spectra. OHD data points in Stern et al. (2010), whose quality of
spectrum and method of fitting spectra are similar to ours, share the
comparable degree of accuracy with us, which also proves the objectivity of
our OHD data points.

Second, as our selected LRGs are incomplete and unable to cover all old-age
galaxies in the universe, there would be some difficulty in tracing the
‘cosmic chronometers’. That is to say, oldest age in our sample may not
represent cosmic age at that redshift. It would be improved through future
redshift surveys of observation and lead to successful realization of
differential age method. In addition, as high SNR is essential for a precise
fitting which determines the age of galaxies, we suggest that the accuracy of
OHD would be improved if the SNR of spectra would increase. We employ ULySS to
reconstruct the stellar population of galaxies. Though the robustness of ULySS
has been illustrated, the fitting minimal may still be local because of the
limitation of groups of initial values. There has been significant advancement
in modelling stellar populations of LRG galaxies and we expect these to
improve the accuracy in this process in the future.

The number of OHD is still scarce compared with SNIa data sets. Advantages of
constraining cosmological models with the OHD (Jimenez & Loeb, 2002; Maor et
al., 2001) are proved, therefore increasing the number of OHD is imperative.
OHD plays almost the same role as that of SNIa for the joint constraints on
the $\Lambda$CDM model. The number of OHD points would extended in further
decades with more and deeper observation of galaxies and at that time OHD set
lone is capable to be used in place of current SNIa data sets (Ma & Zhang,
2011). Fortunately, we have seen that both project and proposed Sandage-Loeb
observational plan (Corasaniti et al., 2007) can be used to extend our
knowledge of cosmic expansion into the even deeper redshift. Finally, we think
it is reasonable to expect that OHD will complement SNIa, BAO and weak lensing
and help us detect more information of the evolution history of our universe.

ACKNOWLEDGMENTS. Special thanks to Robert C. Nichol and Dan P. Carson for
providing us the LRG sample data. We are grateful to G. Bruzual, M. Koleva,
Wei Du, Gao-Chao Liu, Xian-Min Meng, Xue-Lei Chen, You-Jun Lu and Cong Ma for
their helpful suggestions. This work was supported by the National Science
Foundation of China (Grants No. 11173006), the Ministry of Science and
Technology National Basic Science program (project 973) under grant No.
2012CB821804, and the Fundamental Research Funds for the Central Universities.

## References

  * Abazajian et al. (2003) Abazajian, K., et al. 2003, AJ, 126, 2081
  * Abazajian et al. (2009) Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJS, 182, 543
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000
  * Carson & Nichol (2010) Carson, D. P., & Nichol, R. C. 2010, MNRAS, 408, 213
  * Chen & Ratra (2011) Chen, Y., & Ratra, B. 2011, Physics Letters B, 703, 406
  * Corasaniti et al. (2007) Corasaniti, P.-S., Huterer, D., & Melchiorri, A. 2007, Phys. Rev. D, 75, 062001
  * Crawford et al. (2010) Crawford, S. M., et al. 2010, MNRAS, 406, 2569
  * Du et al. (2010) Du, W., Luo, A. L., Prugniel, P., Liang, Y. C., & Zhao, Y. H. 2010, MNRAS, 409, 567
  * Eisenstein et al. (2003) Eisenstein, D. J., Hogg, D. W., et al. 2003, ApJ, 585, 694
  * Eisenstein et al. (2005) Eisenstein, D. J., Zehavi, I., et al. 2005, ApJ, 633, 560
  * Eisenstein et al. (2001) Eisenstein, D. J., et al. 2001, AJ, 122, 2267
  * Gaztañaga et al. (2009) Gaztañaga, E., Cabré, A., & Hui, L. 2009, MNRAS, 399, 1663
  * Ghirlanda et al. (2004) Ghirlanda, G., Ghisellini, G., Lazzati, D., & Firmani, C. 2004, ApJ, 613, L13
  * He & Liu (2011) He, X., & Liu, W. 2011, Applied regression analysis (China Renmin University Press)
  * Hicken et al. (2009) Hicken, M., et al. 2009, ApJ, 700, 1097
  * Jimenez & Loeb (2002) Jimenez, R., & Loeb, A. 2002, ApJ, 573, 37
  * Jimenez et al. (2003) Jimenez, R., Verde, L., Treu, T., & Stern, D. 2003, ApJ, 593, 622
  * Koleva et al. (2009a) Koleva, M., de Rijcke, S., Prugniel, P., Zeilinger, W. W., & Michielsen, D. 2009a, MNRAS, 396, 2133
  * Koleva et al. (2009b) Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009b, A&A, 501, 1269
  * Koleva et al. (2009c) Koleva, M., Prugniel, P., De Rijcke, S., Zeilinger, W. W., & Michielsen, D. 2009c, Astronomische Nachrichten, 330, 960
  * Koleva et al. (2008) Koleva, M., Prugniel, P., Ocvirk, P., Le Borgne, D., & Soubiran, C. 2008, MNRAS, 385, 1998
  * Komatsu et al. (2011) Komatsu, E., et al. 2011, ApJS, 192, 18
  * Li et al. (2008) Li, H., et al. 2008, ApJ, 680, 92
  * Li et al. (2009) Li, M., Li, X.-D., Wang, S., & Zhang, X. 2009, J. Cosmology Astropart. Phys, 6, 036
  * Liu et al. (2012) Liu, G., Lu, Y., Chen, X., et al. 2012, ApJ, 758, 107
  * Ma & Zhang (2011) Ma, C., & Zhang, T.-J. 2011, ApJ, 730, 74
  * Maor et al. (2001) Maor, I., Brustein, R., & Steinhardt, P. J. 2001, Physical Review Letters, 86, 6
  * Moresco et al. (2012) Moresco, M., et al. 2012, J. Cosmology Astropart. Phys, 8, 006
  * Percival et al. (2010) Percival, W. J., et al. 2010, MNRAS, 401, 2148
  * Riess et al. (1998) Riess, A. G., et al. 1998, AJ, 116, 1009
  * Roseboom et al. (2006) Roseboom, I. G., et al. 2006, MNRAS, 373, 349
  * Simon et al. (2005) Simon, J., Verde, L., & Jimenez, R. 2005, Phys. Rev. D, 71, 123001
  * Spergel et al. (2007) Spergel, D. N., et al. 2007, ApJS, 170, 377
  * Stern et al. (2010) Stern, D., Jimenez, R., Verde, L., Kamionkowski, M., & Stanford, S. A. 2010, J. Cosmology Astropart. Phys, 2, 008
  * Stoughton et al. (2002) Stoughton, C., et al. 2002, AJ, 123, 485
  * Strauss et al. (2002) Strauss, M. A., et al. 2002, AJ, 124, 1810
  * Wei (2010) Wei, H. 2010, Physics Letters B, 687, 286
  * Wu et al. (2011) Wu, Y., Singh, H. P., Prugniel, P., Gupta, R., & Koleva, M. 2011, A&A, 525, A71
  * Yi & Zhang (2007) Yi, Z.-L., & Zhang, T.-J. 2007, Modern Physics Letters A, 22, 41
  * York et al. (2000) York, D. G., et al. 2000, AJ, 120, 1579
  * Zhang et al. (2010) Zhang, T.-J., Ma, C., & Lan, T. 2010, Advances in Astronomy, 2010, 184284