Riemannian Trust Region Method for Haplotype Assembly

TL;DR

This paper introduces a Riemannian trust region optimization method to accurately solve the haplotype assembly problem by modeling it as a sphere-constrained maximization, effectively escaping local maxima and saddle points.

Contribution

It presents a novel manifold optimization approach using trust region methods for haplotype assembly, addressing nonconvexity challenges.

Findings

High accuracy in haplotype estimation demonstrated

Effective escape from local maxima and saddle points

Outperforms existing methods in simulation results

Abstract

In this letter we model the Haplotype assembly problem (HAP) as a maximization problem over an $(n-1)$-dimensional sphere. Due to nonconvexity of the feasible set, we propose a manifold optimization approach to solve the mentioned maximization problem. To escape local maxima as well as saddle points we utilize trust region method. Simulation results show that our proposed method is with high accuracy in estimation of Haplotype.

Equations

$$\begin{array}[]{l}\mathop{\min}\limits_{{\bf{X}}}\,\,\,\,\left\|{{{\rm{P}}_{\Omega}}\left({\bf{M}}\right)-{{\rm{P}}_{\Omega}}\left(\bf{X}\right)}\right\|_{F}^{2}\\ \,\,\,\,\,{\rm{subject\,\,to}}\,\,\,\,\,\,rank({\bf{X}})=1,\end{array}$$

$${{\rm{P}}_{\Omega}}({\bf{Q}})=\left\{{\begin{array}[]{*{20}{c}}{{{\rm{P}}_{\Omega}}({q_{ij})}={q_{ij}}\,\,\,\,\,\,\left({i,j}\right)\in\Omega}\\ {{{\rm{P}}_{\Omega}}({q_{ij})}=0\,\,\,\,\,\,\left({i,j}\right)\notin\Omega}\end{array}}\right.,$$

$$\begin{array}[]{l}\mathop{\min}\limits_{{\bf{X}}}\,\,\,\,\,rank({\bf{X}})\\ \,\,\,{\rm{subject\,\,to}}\,\,\,\,\,\,{\left\|{{{\rm{P}}_{\Omega}}\left({\bf{M}}\right)-{{\rm{P}}_{\Omega}}\left(\bf{X}\right)}\right\|_{F}^{2}}<\delta,\end{array}$$

$$\underset{\mathbf{u},\mathbf{v}}{\min }{\Vert {\mathrm{P}}_{\Omega }\left(\mathbf{M}\right)-{\mathrm{P}}_{\Omega }\left({\mathbf{u}}_{m\times 1}{\mathbf{v}}_{n\times 1}^T\right)\Vert }_F^2,$$

$$\mathrm{MEC}\left(\mathbf{M},\mathbf{z}\right)=\sum _{i=1}^m\min \left(hd\left({\mathbf{m}}_i,\mathbf{z}\right),hd\left({\mathbf{m}}_i,-\mathbf{z}\right)\right),$$

$$hd\left({\mathbf{m}}_i,\mathbf{z}\right)=\sum _{j|\left(i,j\right)\in \Omega }d\left(m_{ij},z_j\right),$$

$$\begin{array}[]{l}\mathop{\max}\limits_{\bf{x}}\,\,\,\,\,{{\bf{x}}^{T}}{\bf{Wx}}\\ {\rm{subject\,\,to}}\,\,{\rm{}}x_{i}^{2}=1\end{array},$$

$${\xi }_{\mathbf{x}}=\mathrm{grad}f\left(\mathbf{x}\right)=\left(I-\mathbf{x}{\mathbf{x}}^T\right)\mathrm{Grad}f\left(\mathbf{x}\right),\forall \mathbf{x}\in S^{n-1},$$

$$\alpha \left(t\right)=\alpha \left(0\right)\cos \left(\Vert \dot{\alpha }\left(0\right)\Vert \upsilon t\right)+\dot{\alpha }\left(0\right)\frac{1}{\upsilon \Vert \dot{\alpha }\left(0\right)\Vert }\sin \left(\Vert \dot{\alpha }\left(0\right)\Vert \upsilon t\right),$$

$$i\left(\mathcal{M}\right):=\underset{\mathbf{x}\in \mathcal{M}}{\inf }\sup \left\{\epsilon >0:{{\mathrm{Exp}}_{\mathbf{x}}|}_{B_{\epsilon }\left(0_{\mathbf{x}}\right)}\mathrm{is}\mathrm{diffeomorphism}\right\},$$

$$\Vert P_{\alpha }^{0\leftarrow 1}\mathrm{grad}f\left(\mathbf{y}\right)-\mathrm{grad}f\left(\mathbf{x}\right)\Vert \le \beta \mathrm{dist}\left(\mathbf{x},\mathbf{y}\right),$$

$$\text{dist}\left(\mathbf{x},\mathbf{y}\right)={\int }_0^1{⟨\dot{\alpha }\left(t\right),\dot{\alpha }\left(t\right)⟩}^{\frac{1}{2}}dt,$$

$$P_{\alpha }^{0\leftarrow 1}{\zeta }_{\mathbf{y}}={\zeta }_{\mathbf{x}}+{\nabla }_{\xi }\zeta +{\int }_0^1\left(P_{\alpha }^{0\leftarrow \tau }{\nabla }_{\dot{\alpha }\left(\tau \right)}\zeta -{\nabla }_{\xi }\zeta \right)\text{d}\tau ,$$

$${\nabla }_{\xi }{\zeta }_{\mathbf{x}}=\left(I-\mathbf{x}{\mathbf{x}}^T\right)\text{D}{\zeta }_{\mathbf{x}}\left(\mathbf{x}\right)\left[\xi \right],$$

$$\mathrm{Hess}f(\mathbf{x}):T_{\mathbf{x}}\mathcal{M}\to T_{\mathbf{x}}\mathcal{M}:\xi \to {\nabla }_{\xi }\mathrm{grad}f,\forall \xi \in T_{\mathbf{x}}\mathcal{M}.$$

$$\begin{array}[]{*{20}{l}}{\mathop{\min}\limits_{\eta\in{T_{{{\bf{x}}}}}{\cal M}}{{\hat{m}}_{{{\bf{x}}}}}\left(\eta\right)=f\left({{{\bf{x}}}}\right)+\left\langle{{\rm{grad}}f\left({{{\bf{x}}}}\right),\eta}\right\rangle_{{{\bf{x}}}}+0.5\left\langle{{H}\left[\eta\right],\eta}\right\rangle_{{{\bf{x}}}}}\\ {{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\rm{s}}.{\rm{t}}.{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{\mkern 1.0mu}{{\left\langle{\eta,\eta}\right\rangle}_{{{\bf{x}}}}}\leq\Delta^{2}}\end{array},$$

$$\Vert \xi \Vert \ge \mu \mathrm{dist}\left(\mathbf{x},R_{\mathbf{x}}\xi \right)\forall \mathbf{x}\in \mathcal{M},\forall \xi \in T_{\mathbf{x}}\mathcal{M},\Vert \xi \Vert \le {\delta }_{\mu },$$

$$\left|\frac{d}{d\tau }{\hat{f}}_{\mathbf{x}}\left(\tau \xi \right){|}_{\tau =t}-\frac{d}{d\tau }{\hat{f}}_{\mathbf{x}}\left(\tau \xi \right){|}_{\tau =0}\right|\le {\beta }_{RL}t,$$

$$\begin{array}[]{l}{{\hat{m}}_{{{\bf{x}}_{k}}}}\left(0\right)-{{\hat{m}}_{{{\bf{x}}_{k}}}}\left({{\eta_{k}}}\right)\geq\\ \,\,\,\,\,\,\,\,\,\,\,\,\,{c_{1}}\left\|{{\rm{grad}}f\left({{{\bf{x}}_{k}}}\right)}\right\|\min\left({{\Delta_{k}},\frac{{\left\|{{\rm{grad}}f\left({{{\bf{x}}_{k}}}\right)}\right\|}}{{\left\|{{H_{k}}}\right\|}}}\right)\end{array},$$

$$\underset{k\to \infty }{\lim }\mathrm{grad}f\left({\mathbf{x}}_k\right)=0.$$

$$\underset{\genfrac{}{}{0ex}{}{\mathbf{x}\in {\mathbb{R}}^n,\text{hfill}}{{\Vert \mathbf{x}\Vert }_2\le 1\text{hfill}}}{\max }{\Vert \overline{\mathbf{M}}\mathbf{x}\Vert }_1$$

$$\underset{\mathbf{x}\in S^{n-1}\text{hfill}}{\max }{\Vert \overline{\mathbf{M}}\mathbf{x}\Vert }_1.$$

$$\underset{\mathbf{x}\in S^{n-1}\text{hfill}}{\max }{\Vert \left({\mathrm{P}}_{\Omega }\left(\mathbf{M}\right)\right)\mathbf{x}\Vert }_1,$$

$$\underset{\mathbf{x}\in S^{n-1}\text{hfill}}{\min }-{\Vert \left({\mathrm{P}}_{\Omega }\left(\mathbf{M}\right)\right)\mathbf{x}\Vert }_1.$$

$$\underset{\mathbf{x}\in S^{n-1}}{\min }f(\mathbf{x})=-\sum _{i=1}^m{\left({\left({\mathbf{M}}_{{\Omega }_i}\mathbf{x}\right)}^2+\epsilon \right)}^{\frac{1}{2}},$$

$$\begin{array}{c}\text{dist}\left(\mathbf{x},\mathbf{y}\right)={\Vert \text{grad}f\left(\mathbf{x}\right)\Vert }_2.\end{array}$$

$$\begin{array}{c}\Vert P_{\alpha }^{0\leftarrow 1}{\zeta }_{\mathbf{y}}-{\zeta }_{\mathbf{x}}\Vert =\Vert {\int }_0^1{P}_{\alpha }^{0\leftarrow \tau }{\nabla }_{\dot{\alpha }\left(\tau \right)}\zeta \text{d}\tau \Vert \text{hfill}\\ ⩽{\int }_0^1{\Vert {\nabla }_{\dot{\alpha }\left(\tau \right)}\zeta \Vert }_2\text{d}\tau ,\text{hfill}\end{array}$$

$$\begin{array}{c}{\nabla }_{\dot{\alpha }\left(\tau \right)}\zeta =\left(I-\alpha \left(\tau \right)\alpha {\left(\tau \right)}^T\right)\underset{t\to 0}{\lim }\left(\frac{{\zeta }_{\alpha \left(\tau \right)+t\dot{\alpha }\left(\tau \right)}-{\zeta }_{\alpha \left(\tau \right)}}{t}\right).\text{hfill}\end{array}$$

$$\text{Grad}f\left(\alpha \left(\tau \right)\right)=-\sum _{i=1}^m\left({\mathbf{M}}_{{\Omega }_i}^T{\mathbf{M}}_{{\Omega }_i}\alpha \left(\tau \right)\right){\left({\left({\mathbf{M}}_{{\Omega }_i}\alpha \left(\tau \right)\right)}^2+\epsilon \right)}^{-\frac{1}{2}},$$

$$\underset{t\to 0}{\lim }\left(\frac{{\zeta }_{\alpha \left(\tau \right)+t\dot{\alpha }\left(\tau \right)}-{\zeta }_{\alpha \left(\tau \right)}}{t}\right)=\sum _{i=1}^m\left(\frac{-{\mathbf{M}}_{{\Omega }_i}^T{\mathbf{M}}_{{\Omega }_i}\dot{\alpha }\left(\tau \right)}{{\left({\left({\mathbf{M}}_{{\Omega }_i}\alpha \left(\tau \right)\right)}^2+\epsilon \right)}^{1/2}}\right).$$

Full text

I Introduction

Haplotype is a string of single nucleotide polymorphisms (SNPs) of chromosomes [1]. Haplotypes are useful in studying human evolutionary history and drug discovery and development. Haplotypes can be assembled utilizing sequenced reads, where each read is a fragment of the two chromosomes [2]. Such assembly can be performed by solving mathematical models of sequenced reads. For diploid organisms, each SNP site is considered to be either $+1$ or $-1$. This means that the result of haplotype assembly from the sequenced reads ought to be a vector of $\pm 1$ elements. Moreover, the obtained haplotype, say ${\bf{h}}_{n\times 1}$, is corresponding to one of the two chromosomes and $-\bf{h}$ is corresponding to the other one.

One mathematical approach to find the haplotype is matrix completion. In this approach a read matrix, say ${\bf{M}}_{m\times n}$, containing sequenced reads is created. The elements of this matrix are either $\pm 1$ corresponding to the reads or $\times$ where there is no read. Then, it is required to replace the $\times$ elements of the matrix ${\bf{M}}$ with $+1$ or $-1$ so that the rank of the newly generated matrix, say $\overline{\bf{M}}$, be one. Then the factorization $\overline{\bf{M}}=\overline{\bf{c}}_{m\times 1}\overline{\bf{h}}_{n\times 1}^{T}$ gives the haplotype ${\bf{h}}_{n\times 1}$ through applying sign function over $\overline{\bf{h}}_{n\times 1}$. When the reads are affected by noise, the sign of some of the entries of the matrix ${\bf{M}}_{m\times n}$ is changed. In this case using the following optimization problem the completed matrix $\overline{\bf{M}}$ is obtained. Please note that in HAP, whether to estimate ${\bf{h}}$ or $-{\bf{h}}$, the estimation is accurate.

[TABLE]

where ${{\rm{P}}_{\Omega}}$ is the sampling operator acting as follows

[TABLE]

and $\Omega$ is the set of the observations. There exist some approaches to solve the problem (1). Based on [3], the problem (1) can be directly solved over the manifold of rank-one matrices. Also, a variation of the formulation of the problem (1) has been solved in [4] over the Cartesian product of Grassmann manifolds. In fact, instead optimizing over the variable $\bf{X}$, the problem is solved over the variables $\bf{U}$, $\bf{\Sigma}$ and $\bf{V}$ where ${\bf{X}}={\bf{U}\bf{\Sigma}\bf{V}}^{T}$ is the singular value decomposition of $\bf{X}$. Consider the following equivalent optimization problem for the problem (1).

[TABLE]

where $\delta$ is a given value. In [5] the nonconvex objective function $rank({\bf{X}})$ has been replaced by the nuclear norm and the minimization problem turns into a convex one [6]. Moreover, the problem (1) can be modeled as the following optimization problem.

[TABLE]

Alternating minimization is utilized to solve the above nonconvex problem [7]. In [8], the authors have given an example to show that the proposed objective function of the problem (1) is not a reliable cost function in haplotype assembly. They instead have proposed an objective function to truly model the haplotype assembly problem. Using simulation results they have shown their approach is more accurate than previously proposed algorithms in haplotype assembly .

Apart from the aforementioned optimization approaches, it has been proposed to minimize the following minimum error correction (MEC) function to estimate the haplotype.

[TABLE]

where ${\bf{m}}_{i}$ is the $i$-th row of $\bf{M}$ and

[TABLE]

where $d(\cdot,\cdot)$ equals [math] for the same inputs and $1$ otherwise. It is NP-Hard to obtain the optimal solution of the MEC problem (5). Therefore, heuristic methods have been proposed to solve the problem (5). For example, in [9], a heuristic combinatorial approach is proposed to find the haplotype. However, no provable guarantee for a good solution is proposed therein. Another heuristic approaches can be found in [10] and [11]. In [12], the authors have firstly offered the following optimization problem to find the haplotype,

[TABLE]

where ${\bf{W}}$ is the adjacency matrix defined to evaluate the similarity of each pair of rows of the read matrix and $x_{i}$ is the $i$th entry of $\bf{x}$. Then, they proposed a semidefinite program and solved it rapidly and accurately.

In this paper, we talk about some facts about the haplotype assembly problem (HAP). Then, based on such facts we propose a maximization problem to estimate haplotypes. Note that unlike the maximization problem 7, our formulation directly uses the reads and moreover is not a quadratic form. The proposed problem is defined over the $(n-1)$-dimensional sphere. We use an algorithm to solve the mentioned problem and prove the convergence of the algorithm. Simulation results illustrate that for a large observation error probability, our method outweighs some of the previously proposed methods in the sense of haplotype estimation accuracy.

The remainder of the paper is organized as follows. In Section II some preliminaries are introduced to state our main result. Section III talks about the main result of the paper. Simulation results are presented in Section IV. Conclusion is presented in Section V.

II Preliminaries

In Section III we formulate HAP as a manifold optimization problem. To do so, we require some related concepts as discussed in the following.

Definition 1**.**

[13]** Let $\mathcal{M}$ be a manifold. The set of all tangent vectors at point ${{\bf{x}}}\in{\mathcal{M}}$ is called the tangent space at point ${\bf{x}}$ and denoted by $T_{{\bf{x}}}{\mathcal{M}}$. Moreover, the disjoint collection of the tangent spaces is called tangent bundle and denoted by $T{\mathcal{M}}$.

Definition 2**.**

A given subset $\mathcal{U}$ of the manifold $\mathcal{M}$ along with a bijective mapping $\varphi$ between $\mathcal{U}$ and an open subset of $\mathbb{R}^{d}$ consist a pair $\left(\mathcal{U},\varphi\right)$ which is called a $d$-dimensional chart of $\mathcal{M}$.

Definition 3**.**

[13]** A differentiable bijective mapping between two manifolds is called diffeomorphism provided that its inverse is differentiable too.

Definition 4**.**

[13]** A tangent vector field $\xi$ on $\mathcal{M}$ is a smooth function which assigns to each point of the manifold a tangent vector belonging to tangent bundle. The gradient of a real valued smooth function $f$ over a manifold is an example of vector field which is denoted by ${\rm{grad}}f$. For example, tangent vector field $\xi$ over sphere $S^{n-1}$ is specified as:

[TABLE]

where ${\rm{Grad}}f\left({\bf{x}}\right)$ is the gradient over the Euclidean space and $I$ is the identity matrix [13].

Definition 5**.**

[13]** Let the manifold $\mathcal{N}$ be a subset of the manifold $\mathcal{M}$. When the manifold topology of $\mathcal{N}$ coincides the induced topology of $\mathcal{M}$, $\mathcal{N}$ is called embedded submanifold.

Definition 6**.**

[13]** A manifold endowed with a smoothly varying inner product, is called Riemannian manifold. The inner product over the manifold is denoted by either $g$ or $\left<\cdot,\cdot\right>$. Moreover, $\left<\cdot,\cdot\right>_{\bf{x}}$ illustrates the restriction of the inner product to the tangent space $T_{\bf{x}}\mathcal{M}$. Also, the metric induced by this norm is called Riemannian distance and denoted by $\rm{dist}(\cdot,\cdot)$.

Definition 7**.**

[13]** A locally distance minimizing curve over a manifold is called geodesic. Specifically, straight lines over Euclidean spaces are geodesics.

Definition 8**.**

[13]** A globally distance minimizing curve over a manifold is called minimizing geodesic. For example, $\alpha$ is the minimizing geodesic between points ${\bf{x}},{\bf{y}}\in S^{n-1}$ over sphere $S^{n-1}$ as:

[TABLE]

where ${\dot{\alpha}\left(0\right)}\in T_{\bf{x}}{S^{n-1}}$, $\alpha(0)={\bf{x}}$, and $\alpha(1)={\bf{y}}$. For simplicity we consider $\upsilon=1$.

Definition 9**.**

[13]** Let $\mathcal{M}$ be a Riemannian manifold whose tangent bundle is $T\mathcal{M}$. Exponential map $\rm{Exp}$ is a mapping from $T\mathcal{M}$ to $\mathcal{M}$ so that for $v\in T\mathcal{M}$, $\rm{Exp(v)}$ is equal to $h$ at time $1$; where $h$ is the unique geodesic starting from the base point of $v$ with velocity $v$ at time [math].

Definition 10**.**

[13]** Injectivity radius of the manifold $\mathcal{M}$ is defined as follows,

[TABLE]

where ${{\left.{{\rm{Exp}}_{\bf{x}}}\right|}_{{B_{\varepsilon}}\left({{0_{\bf{x}}}}\right)}}$ shows the restriction of the mapping ${{\rm{Exp}}_{\bf{x}}}$ to the ball ${{B_{\varepsilon}}\left({{0_{\bf{x}}}}\right)}$. Moreover, ${{B_{\varepsilon}}\left({{0_{\bf{x}}}}\right)}$ is a normal neighborhood.

Definition 11**.**

[13]** Let $(\mathcal{M},g)$ own a positive injectivity radius of $i(\mathcal{M})$. Then the real valued function $f$ on $\mathcal{M}$ is $\rm{Lipschitz\,continuously\,differentiable}$ ($L-C^{1}$) provided that

$f$* is differentiable,* 2. 2.

$\forall{\bf{x}},{\bf{y}}\in\mathcal{M}$* with ${\rm{dist}}({\bf{x}},{\bf{y}})<i(\mathcal{M})$, there exists $\beta$ for which*

[TABLE]

where $\alpha$ is the unique minimizing geodesic satisfying $\alpha(0)={\bf{x}}$ and $\alpha(1)={\bf{y}}$. Moreover, $P_{\alpha}^{0\leftarrow 1}$ is an isometry operator; which translates the tangent vector ${\rm{grad}}f\left(y\right)\in T_{\bf{y}}{\mathcal{M}}$ to the tangent space $T_{\bf{x}}{\mathcal{M}}$ making it possible to differentiate the tangent vectors ${\rm{grad}}f\left({\bf{y}}\right)$ and ${\rm{grad}}f\left({\bf{x}}\right)$, and is called parallel translation. Also, ${\rm{dist}}({\bf{x}},{\bf{y}})$ is obtained by taking infimum over the length of all curves joining ${\bf{x}}$ to ${\bf{y}}$ (see formula (3-30) of **[13]**). For example, for the compact embedded submanifold $S^{n-1}$ we have

[TABLE]

where $\left\langle{a,b}\right\rangle={a^{T}}b$.

Proposition 1**.**

(Lemma 7.4.7 of [13]) Consider ${{B_{\varepsilon}}\left({{{\bf{x}}}}\right)}$ as a normal neighborhood of ${\bf{x}}\in\mathcal{M}$ and $\zeta$ as a continuously differentiable tangent vector field over $\mathcal{M}$. Also, let $\alpha$ be the unique minimizing geodesic with $\alpha(0)={\bf{x}}$, $\alpha(1)={\bf{y}}$, and $\dot{\alpha}(1)=\xi$. Then $\forall{\bf{y}}\in{{B_{\varepsilon}}\left({{{\bf{x}}}}\right)}$:

[TABLE]

where $\nabla$ is Riemannian connection which generalizes the concept of directional derivative of a vector field. Specifically, $\nabla$ for Sphere at point ${\bf{x}}$ is as follows:

[TABLE]

where ${\text{D}}{\zeta_{\bf{x}}}\left({\bf{x}}\right)\left[\xi\right]$ is the conventional directional derivative on a Euclidean space.

Definition 12**.**

[13]** The Riemannian Hessian of the real valued function $g$ on the Riemannian manifold $\mathcal{M}$ at ${\bf{x}}\in\mathcal{M}$ is defined as follows:

[TABLE]

To solve a manifold optimization problem, Riemannian Trust Region (RTR) method can be utilized. RTR methods utilize second order geometry of the cost function which lets them escape saddle points and obtain better results compared to first order line search methods [13]. RTR algorithm to minimize given cost function $f$ over the manifold $\mathcal{M}$ is presented in Table 1. Before presenting RTR algorithm, let us define the trust region subproblem.

Definition 13**.**

[13]** Let $f({\bf{x}})$ be a real valued cost function over the manifold $\mathcal{M}$, the following quadratic optimization problem is called trust region subproblem,

[TABLE]

where $H$ is some symmetric operator on ${T_{{{\bf{x}}}}}\mathcal{M}$.

[TABLE]

Proposition 2**.**

[13]** Let $\{{\bf{x}}_{k}\}$ be a sequence generated by Alg. 1. Also, the following conditions are satisfied,

Mapping $f$ is $C^{1}$ and bounded below at the level set $\{{\bf{x}}\in{\mathcal{M}}:f({\bf{x}})\leq f({\bf{x}}_{0})\}$, 2. 2.

Mapping $f$ is $L-C^{1}$, 3. 3.

There exist $\mu>0$ and $\delta_{\mu}>0$ such that the retraction function $R:T\mathcal{M}\to\mathcal{M}$ (see Definition 4.1.1 of **[13]**) satisfies

[TABLE] 4. 4.

Mapping $\widehat{f}:T\mathcal{M}\to{\mathbb{R}}:\xi\to f(R\xi)$ is radially $L-C^{1}$, i.e., $\exists{\beta_{RL}},{\delta_{RL}}>0$ such that

[TABLE]

for all ${\bf{x}}\in\mathcal{M}$, $t<\delta_{RL}$ and $\xi\in T_{\bf{x}}\mathcal{M}$ with $\|\xi\|=1$. 5. 5.

There is a constant $\beta_{H}$ such that $\|H_{k}\|\leq\beta_{H}$ for all $k$, where $H_{k}$ is the symmetric operator defined in Definition 16 at iteration $k$. 6. 6.

Any $\eta_{k}$ obtained in Step 1 of Alg. 1 satisfies inequality

[TABLE]

for some constant $c_{1}>0$, where ${{\left\|{{H_{k}}}\right\|}}$ is the operator norm of ${{H_{k}}}$. Then, the following holds:

[TABLE]

III Main result

In this section we fistly present some facts in regard of HAP. Then, we propose an optimization problem to estimate haplotypes. Finally, we discuss the convergence of an algorithm for solving the proposed optimization problem.

Let $\bf{M}$ be the noiseless read matrix and $\overline{\bf{M}}$ be the completion of $\bf{M}$. Then, the following statements hold.

$\overline{\bf{M}}$ is a rank one matrix for which there exists the factorization $\overline{\bf{M}}=\overline{\bf{c}}_{m\times 1}\overline{\bf{h}}_{n\times 1}^{T}$. 2. 2.

There is no difference for HAP to estimate ${\bf{h}}$ as $sign(\overline{\bf{h}})$ or $sign(-\overline{\bf{h}})$. In other words, $\bf{h}$ is equivalent to $-\bf{h}$. 3. 3.

There exists a unique maximizer for the problem,

[TABLE]

up to the equivalence of $\bf{h}$ and $-\bf{h}$. 4. 4.

It can be verified that (21) is equivalent to the following optimization problem over Sphere $S^{n-1}$,

[TABLE]

We therefore propose the next optimization problem to estimate the haplotype for the noisy HAP,

[TABLE]

Please note that, even though objective function of Problem 21 is convex, we intend to maximize the objective function. Therefore, the solution of the optimization problem cannot be trivially obtained through convex optimization approaches. Moreover, due to Alg. 1 is a descent algorithm, we rewrite Problem (23) in the following form:

[TABLE]

Let us make the objective function of (24) differentiable to easily use smooth optimization approaches for finding the solution. For this purpose, we propose the subsequent problem,

[TABLE]

where ${{\bf{M}}_{{\Omega_{i}}}}$ is the $i$th row of the matrix ${{\rm{P}}_{\Omega}}\left({\bf{M}}\right)$ and $\varepsilon$ is a very small positive value. Please note that for $\varepsilon=0$, the objective function of (25) turns into the objective function of (24).

Theorem 1**.**

For the optimization problem (25), Conditions of Proposition 2 are satisfied and consequently Alg. 1 converges.

Proof.

We require to prove that the conditions of Proposition 2 are satisfied for $f({\bf{x}})$ in problem (25).

Condition 1; It is easy to see that $f({\bf{x}})$ is continuous. Let $(\mathcal{U},\varphi_{i})$ be a chart for $S^{n-1}$, then $f$ is differentiable over $S^{n-1}$ if $f\circ\varphi_{i}^{-1}$ is differentiable [13]. Let ${\bf{x}}=(x_{1},\cdots,x_{n})\in{S}^{n}$, then, a straightforward choice for $\varphi_{i}$ is to choose $\varphi_{i}$ so that $\varphi_{i}({\bf{x}})=(x_{1},\cdots,x_{i-1},x_{i+1},\cdots,x_{n})$, when $x_{i}\neq 0$. Accordingly, it is easy to see that $f\circ\varphi_{i}^{-1}$ is differentiable and subsequently $f$ is differentiable. Therefore, $f({\bf{x}})$ is $C^{1}$. Moreover, for ${\bf{x}}\in{S^{n-1}}$ it is obvious that $f({\bf{x}})$ is lower bounded.

Condition 2; Let $\alpha(t)$ be the geodesic defined in Definition (8) with $\upsilon=1$, and distance function be as Equation (12). Now, let $\xi(0)=\dot{\alpha}(0)={\rm{grad}}f({\bf{x}})$. Being aware of the fact that ${{\bf{x}}^{T}}\rho=0,\,\,\,\,\forall\rho\in{T_{\bf{x}}}{S^{n-1}}$, specifically for $\rho={\rm{grad}}f({\bf{x}})$, it is easy to verify that,

[TABLE]

Then, by considering $\zeta$ as the gradient vector field in equality (13), we have:

[TABLE]

where we used the isometry property of ${P_{\alpha}^{0\leftarrow\tau}}$. Using Equation (14), we have:

[TABLE]

By calculating the Euclidean gradient for the cost function of Problem (25) as

[TABLE]

it is not difficult to verify that

[TABLE]

Then, using Equality (8), Cauchy-Shwarz inequality and considering that ${\left\|{\dot{\alpha}\left(\tau\right)}\right\|_{2}}={\left\|{{\text{grad}}f\left({\bf{x}}\right)}\right\|_{2}}$ and ${\left\|{\left({I-\alpha\left(\tau\right)\alpha{{\left(\tau\right)}^{T}}}\right)}\right\|_{F}}\leqslant n$, we will have:

[TABLE]

The aforementioned argument along with Equality (26) show that $f$ is $L-C^{1}$ with $\beta=n\sum\limits_{i=1}^{m}{{{\left\|{\frac{{\left({{\mathbf{M}}_{{\Omega_{i}}}^{T}{{\mathbf{M}}_{{\Omega_{i}}}}}\right)}}{{\sqrt{\varepsilon}}}}\right\|}_{F}}}$.

Condition 3; Let $R_{\bf{x}}$ be the restriction of retraction function $R$ to the tangent space $T_{\bf{x}}S^{n-1}$ which can be defined as ${R_{\bf{x}}}(\xi)=\frac{{{\bf{x}}+\xi}}{{{{\left\|{{\bf{x}}+\xi}\right\|}_{2}}}},\forall\xi\in T_{\bf{x}}S^{n-1}$ [13]. It is not difficult to verify that for a given $\xi\in T_{\bf{x}}S^{n-1}$ the geodesic $\alpha(t)$, defined in Formula (9), satisfies $\alpha(0)={\bf{x}}$, $\alpha(1)=R_{\bf{x}}(\xi)$ when $\dot{\alpha}\left(0\right)=\frac{\xi}{{{{\left\|\xi\right\|}_{2}}}}{\cos^{-1}}\left({\left\|{{\bf{x}}+\xi}\right\|_{2}^{-1}}\right)$. Then, using Equation (12) we have:

[TABLE]

Now, let us find $\delta_{\mu}$ and $\mu$ for sake of inequality (17). It is easy to numerically verify that for $\mu=1$ and any arbitrary $\delta_{\mu}\geq 0$, inequality (17) holds.

Condition 4; Let us evaluate radially $L-C^{1}$ property of $\widehat{f}$. We have:

[TABLE]

where $\gamma=\frac{{\left({{\bf{x}}+\tau\xi}\right)}}{{{{\left\|{{\bf{x}}+\tau\xi}\right\|}_{2}}}}$. Now, using the fact that ${\bf{x}}$ is orthogonal to $\xi$, $\|{\bf{x}}\|_{2}=\|\xi\|_{2}=1$ and some simple calculus, for any arbitrary $t>0$ we obtain that:

[TABLE]

Meaning that radially $L-C^{1}$ condition has been satisfied with ${\beta_{RL}}=t\left({\sum\limits_{i=1}^{m}{\left\|{{{\mathbf{M}}_{{\Omega_{i}}}}}\right\|_{2}^{3}+\left\|{{{\mathbf{M}}_{{\Omega_{i}}}}}\right\|_{2}^{2}}}\right)$.

**Condition 5;**Let us consider the symmetric operator $H_{k}$ be Hessian of function $f$ at point ${\bf{x}}_{k}$. Then, based on Definition 12, the fact that $\|H_{k}\|\leq\|H_{k}\|_{F}$ and also an argument discussed in verifying Condition 2, we just require to prove the boundedness of $\|{\rm{grad}}f({\bf{x}}_{k})\|$ (see Inequality (31)). Based on Equation (8) and using the Eulidean gradient of function $f$, as defined in Equation (29), we will have:

[TABLE]

This results in

[TABLE]

Condition 6; It is shown in [13] that using the truncated conjugate gradient method (see Alg. 11 of [13]), Inequality (19) is satisfied with $c_{1}=1/2$.

The proof is complete. ∎

IV Simulation Results

Hamming distance of the original haplotype and its estimation, denoted by $\rm{hd}$, is our criterion to evaluate the performance of different methods. Please note that in HAP, it does not matter that the original haplotype is either $\bf{h}$ or $-\bf{h}$. Accordingly, to obtain $\rm{hd}$ we calculate the Hamming distance of the estimated haplotype with both $\bf{h}$ and $-\bf{h}$ and choose the minimum one. Simulations are performed using synthetic data created by random generation of bipolar vectors ${\bf{h}}_{n\times 1}$ and ${\bf{c}}_{m\times 1}$ to construct ${\bf{M}}={\bf{c}}_{m\times 1}{\bf{h}}_{n\times 1}^{T}$. The observation set $\Omega$ is randomly produced with the probability of observation $\rm{pd}$. Also, some observed samples are erroneous after changing their original sign; we show the set of erroneous samples by $\Omega_{E}$, where $\Omega_{E}\in\Omega$. Fig. 1 is depicted for $m=250$, $n=300$, and $0.25\leq$pd$\leq 0.7$. Moreover, ${{\left|\Omega_{E}\right|}\mathord{\left/{\vphantom{{\left|\Omega\right|}{\left|{{\Omega_{E}}}\right|}}}\right.\kern-1.2pt}{\left|{{\Omega}}\right|}}$ is set to $0.35$, where $\left|\cdot\right|$ denotes the cardinality of a set. As seen, our method outperforms the other approaches by generating lower $\rm{hd}$.

V Conclusion

In this letter we proposed a new method for haplotype estimation. We properly modeled HAP over an (n-1)-dimensional Sphere. We also discussed the convergence of a Riemannian trust region method. Simulation results confirmed our method outperforms some of the other haplotype assembly methods.