Haplotype Assembly Using Manifold Optimization and Error Correction Mechanism

TL;DR

This paper introduces a novel manifold optimization approach with error correction for haplotype assembly, effectively handling noisy data and outperforming existing methods in accuracy.

Contribution

It proposes a new MEC-based matrix completion model on rank-one manifolds and proves its convergence, advancing haplotype assembly techniques.

Findings

Outperforms existing matrix completion methods in haplotype assembly accuracy

Provides a convergent iterative algorithm for the proposed optimization problem

Achieves more accurate haplotype estimation compared to recent MEC algorithms

Abstract

Recent matrix completion based methods have not been able to properly model the Haplotype Assembly Problem (HAP) for noisy observations. To cope with such a case, in this letter we propose a new Minimum Error Correction (MEC) based matrix completion optimization problem over the manifold of rank-one matrices. The convergence of a specific iterative algorithm for solving this problem is proved. Simulation results illustrate that the proposed method not only outperforms some well-known matrix completion based methods, but also presents a more accurate result compared to a most recent MEC based algorithm for haplotype estimation.

Equations

$$\left\{{\begin{array}[]{*{20}{c}}{{r_{ij}}=\pm 1\,\,\,\,\,\,\left({i,j}\right)\in\Omega}\\ {{r_{ij}}=\times\,\,\,\,\,\,\left({i,j}\right)\notin\Omega,}\end{array}}\right.$$

$$\overline{\mathbf{R}}={\mathbf{c}}_{m\times 1}{\mathbf{h}}_{n\times 1}^T,$$

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

$$hd\left({\mathbf{r}}_i,\mathbf{h}\right)=\sum _{j|\left(i,j\right)\in \Omega }d\left(r_{ij},h_j\right),$$

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

$${{\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.$$

$$\underset{\mathbf{X}}{\min }{\Vert {\mathrm{P}}_{\Omega }\left(\mathbf{R}\right)-{\mathrm{P}}_{\Omega }\left(\mathbf{X}\right)\Vert }_F^2+\lambda {\Vert \mathbf{X}\Vert }_{\ast },$$

$$\begin{array}[]{l}{T_{\bf{X}}}{{\mathcal{M}}^{(1)}}=\left\{{\left[{{\bf{U}}\,\,{{\bf{U}}_{\bot}}}\right]\left[{\begin{array}[]{*{20}{c}}{{\mathbb{R}}}&{{\mathbb{R}^{1\times\left({n-r}\right)}}}\\ {{\mathbb{R}^{\left({m-1}\right)\times 1}}}&{{{\bf{0}}_{\left({m-1}\right)\times\left({n-1}\right)}}}\end{array}}\right]{{\left[{{\bf{V}}\,\,{{\bf{V}}_{\bot}}}\right]}^{T}}}\right\}\\ \,\,\,\,\,\,\,\,\,\,\,\,\,=\left\{{{\bf{UM}}{{\bf{V}}^{T}}+{{\bf{U}}_{p}}{{\bf{V}}^{T}}+{\bf{UV}}_{p}^{T}:}\right.{\bf{M}}\in{\mathbb{R}},\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,{{\bf{U}}_{p}}\in{\mathbb{R}^{m\times 1}},{\bf{U}}_{p}^{T}{\bf{U}}={\bf{0}},{{\bf{V}}_{p}}\in{\mathbb{R}^{n\times 1}},{\bf{V}}_{p}^{T}{\bf{V}}={\bf{0}},\end{array}$$

$${⟨\mathrm{grad}f(\mathbf{X}),\xi ⟩}_{\mathbf{X}}=\mathrm{D}f(\mathbf{X})[\xi ]\forall \xi \in T_{\mathbf{X}}\mathcal{M},$$

$$\overline{\bf{R}}=\left[{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{1}\\ {1}\\ {1}\end{array}}&{\begin{array}[]{*{20}{c}}{-1}\\ {-1}\\ {-1}\end{array}}&{\begin{array}[]{*{20}{c}}1\\ {1}\\ 1\end{array}}&{\begin{array}[]{*{20}{c}}{-1}\\ {-1}\\ {-1}\end{array}}&{\begin{array}[]{*{20}{c}}{-1}\\ {-1}\\ {-1}\end{array}}\end{array}}\right].$$

$${{\rm{P}}_{\Omega}}\left({{{\overline{\bf{R}}}_{E}}}\right)=\left[{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{1}\\ {-1}\\ {1}\end{array}}&{\begin{array}[]{*{20}{c}}{0}\\ {-1}\\ {-1}\end{array}}&{\begin{array}[]{*{20}{c}}1\\ {1}\\ {-1}\end{array}}&{\begin{array}[]{*{20}{c}}{-1}\\ {-1}\\ {-1}\end{array}}&{\begin{array}[]{*{20}{c}}{1}\\ {-1}\\ {-1}\end{array}}\end{array}}\right],$$

$${\overline{\bf{R}}_{1}}=\left[{\begin{array}[]{*{20}{c}}{0.2255}&{-0.763}&{0.2255}&{-0.8165}&{-0.3655}\\ {0.2955}&{-1}&{0.2955}&{-1.07}&{-0.{\rm{4790}}}\\ {0.2955}&{-1}&{0.2955}&{-1.07}&{-0.{\rm{4790}}}\end{array}}\right]$$

$${\overline{\bf{R}}_{2}}=\left[{\begin{array}[]{*{20}{c}}{-\varepsilon\gamma}&{{\rm{0}}{\rm{.9975}}\varepsilon}&{\varepsilon\gamma}&{{\rm{0}}{\rm{.9975}}\varepsilon}&{{\rm{0}}{\rm{.9975}}\varepsilon}\\ \gamma&{{\rm{-0}}{\rm{.9975}}}&{-\gamma}&{{\rm{-0}}{\rm{.9975}}}&{{\rm{-0}}{\rm{.9975}}}\\ \gamma&{{\rm{-0}}{\rm{.9975}}}&{-\gamma}&{{\rm{-0}}{\rm{.9975}}}&{{\rm{-0}}{\rm{.9975}}}\end{array}}\right],$$

$$\begin{array}[]{l}hd\left({vec\left({{{\rm{P}}_{\Omega}}\left({sign\left({{{\overline{\bf{R}}}_{1}}}\right)}\right)}\right),vec\left({\left({{{\rm{P}}_{\Omega}}\left({{{\overline{\bf{R}}}_{E}}}\right)}\right)}\right)}\right)<\\ hd\left({vec\left({{{\rm{P}}_{\Omega}}\left({sign\left({{{\overline{\bf{R}}}_{2}}}\right)}\right)}\right),vec\left({\left({{{\rm{P}}_{\Omega}}\left({{{\overline{\bf{R}}}_{E}}}\right)}\right)}\right)}\right),\end{array}$$

$$\begin{array}[]{l}{\left\|{{{{\rm{P}}_{\Omega}}\left({sign\left({{{\overline{\bf{R}}}_{1}}}\right)}\right)}-{{\rm{P}}_{\Omega}}\left({{{\overline{\bf{R}}}_{E}}}\right)}\right\|_{0}<}\\ {\left\|{{{{\rm{P}}_{\Omega}}\left({sign\left({{{\overline{\bf{R}}}_{2}}}\right)}\right)}-{{\rm{P}}_{\Omega}}\left({{{\overline{\bf{R}}}_{E}}}\right)}\right\|_{0}},\end{array}$$

$$\underset{\mathbf{X}\in {\mathcal{M}}^{(1)}}{\min }{\Vert {\mathrm{P}}_{\Omega }\left({\overline{\mathbf{R}}}_E\right)-{\mathrm{P}}_{\Omega }\left(\mathrm{sign}\left(\mathbf{X}\right)\right)\Vert }_0.$$

$${\left({\mathrm{P}}_{\Omega }\left({\overline{\mathbf{R}}}_E\right)-{\mathrm{P}}_{\Omega }\left({\gamma }_1{\tan }^{-1}\left({\gamma }_2\mathbf{X}\right)\right)\right)}_{ij}=0$$

$$\underset{\mathbf{X}\in {\mathcal{M}}^{(1)}}{\min }f(\mathbf{X})={\Vert {\mathrm{P}}_{\Omega }\left({\overline{\mathbf{R}}}_E\right)-{\mathrm{P}}_{\Omega }\left({\gamma }_1{\tan }^{-1}\left({\gamma }_2\mathbf{X}\right)\right)\Vert }_p^p,$$

$$\begin{array}[]{l}{\rm{NMSE}}=\frac{{\left\|{{\overline{\bf{R}}}-\widehat{{\overline{\bf{R}}}}}\right\|_{F}^{2}}}{{\left\|{\overline{\bf{R}}}\right\|_{F}^{2}}}\end{array},$$

Full text

I Introduction

Study of genetic variations is a critical task to disease recognition and also drug development. Such variations can be recognized using haplotypes which are known as strings of Single Nucleotide Polymorphisms (SNPs) of chromosomes [1]. However, due to hardware limitations, haplotypes cannot be entirely measured. Instead, they are supposed to be reconstructed from the short reads obtained by high-throughput sequencing systems. For diploid organisms, each SNP site is modeled by either $+1$ or $-1$. Hence, a haplotype sequence $\bf{h}$ corresponding to one of chromosomes is formed by a vector of $\pm 1$ entries, while the haplotype of the other one is given by $-\bf{h}$. Accordingly, the short reads are obtained from either $\bf{h}$ or $-\bf{h}$. Such a sampling procedure can be modeled by a matrix, say ${\bf{R}}$, whose entries are either $\pm 1$ or unobserved. Specifically, if we denote the set of observed entries by $\Omega$, the entries of ${\bf{R}}\in{\left\{-1,+1,\times\right\}}^{m\times n}$ are given by

[TABLE]

where $\times$ shows unobserved entries. The completed matrix of $\bf{R}$, defined by $\overline{\bf{R}}$, is supposed to be a rank-one matrix whose entries are $\pm 1$ and can be factorized as

[TABLE]

where ${\bf{h}}$ is the intended haplotype sequence and c is a sequence of $+1$ and $-1$ entries corresponding to ${\bf{h}}$ and $-{\bf{h}}$, respectively. When the number of reads (observed entries) are sufficient and there is no error in reads, it is easy to obtain $\overline{\bf{R}}$ and consequently $\bf{h}$.

In practice, however, measurements are corrupted by noise which usually leads to erroneous signs in some observed entries. Many haplotype assembly methods deal with erroneous measurements mostly based on the MEC defined as [2]

[TABLE]

where the Hamming distance $hd(\cdot,\cdot)$ calculates the number of mismatch entries using

[TABLE]

with $d(\cdot,\cdot)$ being [math] for equal inputs, and $1$ otherwise.

For optimal solution of (3), which is an NP-hard problem, some heuristic methods have already been addressed [2, 3, 4]. Apart from MEC approaches, some other methods are developed based on the fact that the rank of $\overline{\bf{R}}$ ought to be 1. In [5], matrix factorization is used in the minimization problem,

[TABLE]

where ${{\rm{P}}_{\Omega}}(\cdot)$ is a sampling operator defined as

[TABLE]

and $\|\cdot\|_{F}$ shows the Frobenius norm. Then, $\bf{h}$ is calculated by applying the sign function over $\bf{v}$. The nonconvex problem of (5) can be solved using the alternating minimization approach. For fast matrix completion, QR decomposition may also be applied [6].

From the low-rank matrix completion point of view, haplotype estimation can be performed using some other approaches as follows. Based on [7], it is possible to complete ${{\rm{P}}_{\Omega}}\left({\bf{R}}\right)$ using the following optimization problem,

[TABLE]

where ${\left\|\cdot\right\|_{*}}$ shows the nuclear norm and $\lambda$ is a regularization factor. Afterwards, we can estimate $\bf{h}$ by applying the sign function over the right singular vector corresponding to the largest singular value. In addition, minimizing the cost function $\left\|{{{\rm{P}}_{\Omega}}\left({\bf{R}}\right)-{{\rm{P}}_{\Omega}}\left({\bf{X}}\right)}\right\|_{F}^{2}$ over the manifold of rank-one matrices [8] or Grassmann manifold [9] will result in completion of ${{\rm{P}}_{\Omega}}\left({\bf{R}}\right)$ and estimating $\bf{h}$.

In this paper, we propose a new error correction mechanism over the manifold of rank one matrices to solve the noisy HAP. This approach benefits the underlying structure of the HAP, which seeks for a rank one matrix with the entries of $\pm 1$. As a result, unlike common matrix completion methods, we propose a Hamming distance cost function over the manifold of rank one matrices to complete the desired matrix. We next present a surrogate for the nondifferentiable structure of this cost function and analyze its convergence behaviour. In this way, as opposed to the existing MEC based algorithms, we will be able to derive the performance guarantees for our optimization problem. Simulation results confirm that this method is more reliable for solving noisy HAPs.

The paper is organized as follows. In Section II some required concepts of manifold optimization are presented. Section III is devoted to our main result. Simulation results are illustrated in Section IV and Section V concludes the paper.

II Preliminaries

As mentioned, a HAP can be modeled as a rank-one matrix completion problem. On the other hand, the set of real valued rank-one matrices can be considered as a smooth manifold [8]. To optimize a differentiable function over such a manifold, we first remind some concepts for optimization over manifolds.

Theorem 1**.**

[8]** The set of all $m\times n$ real valued matrices of rank $1$ is a smooth manifold, $\mathcal{M}^{(1)}$, whose tangent space at point ${\bf{X}}\in\mathcal{M}^{(1)}$ is defined as

[TABLE]

where ${\bf{U\Sigma V}}^{T}$ is the singular value decomposition (SVD) of ${\bf{X}}$.

Convergence of manifold optimization schemes is evaluated by a metric on the manifold. Inspired from the Euclidean space $\mathbb{R}^{m\times n}$, the metric on the manifold of rank-one matrices is defined by the inner product $\left<\xi,\eta\right>_{\bf{X}}={\rm{tr}}\left({{\xi^{T}}\eta}\right)$, where the subscript $\bf{X}$ shows the restriction to the tangent space $T_{\bf{X}}\mathcal{M}^{(1)}$ and the tangent vectors $\xi,\eta\in T_{\bf{X}}\mathcal{M}^{(1)}$. A smoothly varying inner product is Riemannian metric and the manifold endowed with this metric is called Riemannian. Also, the gradient of a function on such a manifold is given as follows.

Definition 1**.**

[10]** Let $f$ be a scalar valued function over the Riemannian manifold $\mathcal{M}$ endowed with the inner product $\left<\cdot,\cdot\right>_{\bf{X}}$. Then, the gradient of $f$ at point ${\bf{X}}\in\mathcal{M}$ denoted by ${\rm{grad}}f({\bf{X}})$ is the unique element of the tangent space $T_{\bf{X}}\mathcal{M}$ satisfying

[TABLE]

where ${\rm{D}}f$ denotes the directional derivative acting on the tangent vector $\xi$.

By iteratively solving an optimization problem over the manifold $\mathcal{M}$, we may find a point on $T_{\bf{X}}\mathcal{M}$ which would not necessarily belong to the manifold. To bring this point back to the manifold, the retraction function can be utilized.

Definition 2**.**

[10*]** Let $T\mathcal{M}:=\mathop{\cup}\limits_{{\bf{X}}\in{{\mathcal{M}}}}\left\{{\bf{X}}\right\}\times{T_{\bf{x}}}{{\mathcal{M}}}$ be the tangent bundle and $R:T\mathcal{M}\to\mathcal{M}$ be a smooth mapping whose restriction to $T_{\bf{X}}\mathcal{M}$ is $R_{\bf{X}}$. Then, $R$ is a retraction on the manifold $\mathcal{M}$, if,

$(1)$ $R_{\bf{X}}({\bf{0}}_{\bf{X}})={\bf{X}}$, where, ${\bf{0}}_{\bf{X}}$ is the zero element of $T_{\bf{X}}\mathcal{M}$,

$(2)$ ${\rm{D}}R_{\bf{X}}({\bf{0}}_{\bf{X}})={\rm{id}}_{T_{\bf{X}}\mathcal{M}}$, where, ${\rm{id}}_{T_{\bf{X}}\mathcal{M}}$ is the identity mapping on ${T_{\bf{X}}\mathcal{M}}$.*

More clearly, we find the result of $R_{\bf{X}}(\xi)$ by first calculating ${\bf{Y}}={\bf{X}}+\xi$ where ${\bf{X}}\in\mathcal{M}$ and $\xi\in T_{\bf{X}}\mathcal{M}$, and next applying the retraction mapping on ${\bf{Y}}$. It has been shown that the retraction $R_{\bf{X}}(\xi)$ to the manifold of rank-one matrices is equivalent to taking the SVD of ${\bf{Y}}$ and then making all the singular values equal to zero except the largest one [11].

Another concept to note is the gradient descent algorithm on Riemannian manifolds which is given in Alg. 1 [10]. As seen, Steps 1 and 2 concern with the search direction and convergence evaluation. Step 3, known as Armijo backtracking, guarantees a sequence of points which constitutes a descent direction [10]. Step 4 performs retraction on the manifold.

[TABLE]

III Main result

First, using an example for the noisy HAP, we show that minimizing the cost function $\left\|{{{\rm{P}}_{\Omega}}\left({\bf{R}}\right)-{{\rm{P}}_{\Omega}}\left({\bf{X}}\right)}\right\|_{F}^{2}$ and calculating the haplotype by applying the sign function over the right singular vector corresponding to the largest singular value would not lead to the desired result. Then, we propose an optimization problem which can properly model the HAP.

Example 1**.**

Substitute {\bf{h}}=\left[{\begin{array}[]{*{20}{c}}{1}&{-1}&1&{-1}&{-1}\end{array}}\right]^{T} and {\bf{c}}=\left[{\begin{array}[]{*{20}{c}}{1}&{1}&1\end{array}}\right]^{T} in $\rm{(}$2$\rm{)}$ to build up $\overline{\bf{R}}$ as

[TABLE]

Next, let an erroneous sampling from $\overline{\bf{R}}$ be given as

[TABLE]

where subscript $E$ stands for erroneous observation. As seen, the only unobserved entry is ${\left\{{\left(1,2\right)}\right\}}$ and the set of erroneous observations is ${\left\{{\left(1,5\right)},{\left(2,1\right)},{\left(3,3\right)}\right\}}$. Based on our simulations, minimization of $\left\|{{{\rm{P}}_{\Omega}}\left(\overline{\bf{R}}_{E}\right)-{{\rm{P}}_{\Omega}}\left({\bf{X}}\right)}\right\|_{F}^{2}$ over the manifold of rank-one matrices could lead to either

[TABLE]

or

[TABLE]

*where $\varepsilon$ and $\gamma$ are infinitesimal positive numbers. However, one can easily check that although for this example we get

$\left\|{{{\rm{P}}_{\Omega}}\left(\overline{\bf{R}}_{E}\right)-{{\rm{P}}_{\Omega}}\left({\overline{\bf{R}}_{1}}\right)}\right\|_{F}\simeq\left\|{{{\rm{P}}_{\Omega}}\left(\overline{\bf{R}}_{E}\right)-{{\rm{P}}_{\Omega}}\left({\overline{\bf{R}}_{2}}\right)}\right\|_{F}\simeq 2.8284$, their haplotypes are estimated differently as {\bf{h}}_{1}=\left[{\begin{array}[]{*{20}{c}}{1}&{-1}&{1}&{-1}&{-1}\end{array}}\right]^{T} and {\bf{h}}_{2}=\left[{\begin{array}[]{*{20}{c}}{1}&{-1}&{-1}&{-1}&{-1}\end{array}}\right]^{T}, while only ${\bf{h}}_{1}$ is the correct answer. Note that haplotypes are estimated by applying the sign function over the right singular vector of the completed rank-one matrix.*

For the above example, having the Hamming distance in mind, it is easy to verify that the following inequality holds,

[TABLE]

where $vec\left(\cdot\right)$ vectorizes the input matrix. Equivalently, we have

[TABLE]

where $\|\cdot\|_{0}$ shows the $l_{0}$-norm which counts the number of nonzero entries of a matrix.

Therefore, we propose the following minimization cost function for noisy HAPs as

[TABLE]

Using (16), we can find a rank-one matrix whose sign is as similar as possible to the observations. However, since (16) is nondifferentiable, we replace ${\rm{sign}}(x_{ij})$ by the differentiable function $\gamma_{1}{\tan^{-1}}(\gamma_{2}x_{ij})$ where the positive values $\gamma_{1}$ and $\gamma_{2}$ are selected to approximate the sign function. Another difficulty with nondifferentiability of (16) is the discontinuity of $l_{0}$-norm which is solved by replacing the $l_{p}$-norm $\|\cdot\|_{p}$ that is applied to the vector form of the matrix for $p>0$. However, since $l_{p}$-norm is yet nondifferentiable in the case of

[TABLE]

for any $\left({i,j}\right)\in\Omega$, we limit the value of $\gamma_{1}$ to \left(0,{\raise 3.01385pt\hbox{2}\!\mathord{\left/{\vphantom{2\pi}}\right.\kern-1.2pt}\!\lower 3.01385pt\hbox{\pi}}\right) . In this way, (16) is reformulated in a new problem as

[TABLE]

where ${\tan^{-1}}(\cdot)$ acts elementwise. Although, for $0<p\leq 1$ the sparsity is promoted in (18), here we choose $p$ equal or slightly greater than 1 to easily use the triangle inequality to prove the convergence of Alg. 1. Simulation results will show that such a choice will lead to more accurate estimates of haplotypes. Now, we prove the convergence of Alg. 1 for (18).

Theorem 2**.**

Let positive values \gamma_{1}<{\raise 3.01385pt\hbox{2}\!\mathord{\left/{\vphantom{2\pi}}\right.\kern-1.2pt}\!\lower 3.01385pt\hbox{\pi}} and $\gamma_{2}$ be given. Then, by choosing the initial matrix ${\bf{X}}_{0}$ in such a way that $\|{\bf{X}}_{0}\|_{F}^{2}$ is bounded and $f({\bf{X}}_{0})<|\Omega|$, Alg. 1 will converge.

Proof.

Using Theorem 4.3.1 in [10], we know that every limit point of the infinite sequence $\{{\bf{X}}_{i}\}$ generated by Alg. 1 is a critical point of $f(\bf{X})$ in (18). Accordingly, we need to prove that under the mentioned conditions of Theorem 2, $\{{\bf{X}}_{i}\}$ owns a limit point and also belongs to the manifold of rank-one matrices, meaning that it converges. To do so, we show that this sequence lies in a compact set in which any sequence owns a convergent subsequence [12]. Consequently, due to the fact that $\{{\bf{X}}_{i}\}$ is a decreasing sequence, having a convergent subsequence guarantees its convergence.

To show the compactness of the set containing $\{{\bf{X}}_{i}\}$, we need to prove that the set is bounded and closed [12]. First, we prove the closedness of the set. Suppose that in our problem there exists a limit point for $\{{\bf{X}}_{i}\}$, which can be either a zero matrix, or a rank-one matrix belonging to the manifold of rank-one matrices. However, we ought to show that under the conditions of Theorem 2 the limit point of $\{{\bf{X}}_{i}\}$ is a rank-one matrix. To see this, it is enough to prove that the zero matrix is not a limit point of $\{{\bf{X}}_{i}\}$. By Step 3 of Alg. 1, we know that $\{{\bf{X}}_{i}\}$ is a decreasing sequence for $f({\bf{X}})$. Suppose that $\{{\bf{Y}}_{j}\}$ is a subsequence of $\{{\bf{X}}_{i}\}$ converging to the zero matrix. Then, there is an integer $K$ such that $\forall j>K$, $\|{\bf{Y}}_{j}\|_{p}^{p}\leq\|{\bf{Y}}_{j}\|_{F}^{2}<2^{-j}$. From the triangle inequality, we get $f({\bf{Y}}_{j})\geq\mathop{\left\|{{{\rm{P}}_{\Omega}}\left({{{\overline{\bf{R}}}_{E}}}\right)}\right\|_{p}^{p}}-{\left\|{{\rm{P}}_{\Omega}}\left({\gamma_{1}{\tan^{-1}}\left(\gamma_{2}{\bf{Y}}_{j}\right)}\right)\right\|_{p}^{p}}$, $\forall j>K$, and thus, $\mathop{\lim}\limits_{j\to\infty}f({{\bf{Y}}_{j}})\geq\left|\Omega\right|$. Now, we easily enforce Alg. 1 to have an initial matrix $\{{\bf{X}}_{0}\}$ for which $f({\bf{X}}_{0})<\left|\Omega\right|$. Hence, due to the fact that $\{{\bf{X}}_{i}\}$ is a decreasing sequence for $f({\bf{X}})$, the zero matrix will not be the limit point. This ensures that $\{{\bf{X}}_{i}\}$ is in a closed set.

Now, we show that $\{{\bf{X}}_{i}\}$ stays in a bounded set. For this to happen, we modify (18) to define $f_{\rm{m}}(\bf{X})$ by adding the regularization term $\mu\left\|{\bf{X}}\right\|_{F}^{2}$, where $\mu$ is a positive value. Since we enforced Alg. 1 to have the initial matrix ${\bf{X}}_{0}$ satisfying $f({\bf{X}}_{0})<\left|\Omega\right|$, we get $f_{\rm{m}}({\bf{X}}_{0})<\left|\Omega\right|+\mu\left\|{\bf{X}}_{0}\right\|_{F}^{2}$. Also, let ${\bf{X}}_{0}$ be a norm bounded initial matrix. Consequently, $f_{\rm{m}}({\bf{X}}_{0})$ will be bounded to a value, say $z$. Therefore, the decreasing sequence $\{{\bf{X}}_{i}\}$ generated by Alg. 1 for $f_{\rm{m}}$ satisfies $f({\bf{X}}_{i})+\mu\left\|{\bf{X}}_{i}\right\|_{F}^{2}<z$ showing that $\left\|{\bf{X}}_{i}\right\|_{F}^{2}$ is bounded. Considering the fact that by choosing an infinitesimal value for $\mu$, the sequence generated by Alg. 1 for either $f({\bf{X}})$ or $f_{\rm{m}}({\bf{X}})$ is the same, $\{{\bf{X}}_{i}\}$ stays in a bounded set.

From the above reasoning, $\{{\bf{X}}_{i}\}$ stays in a compact set, meaning that the proposed problem in (18) will converge. ∎

Note that although replacing ${\rm{sign}}\left({\bf{X}}\right)$ by ${\gamma_{1}}{\tan^{-1}}\left({{\gamma_{2}}{\bf{X}}}\right)$ in general leads to slightly reducing the haplotype estimation accuracy, in effect, it makes the optimization problem mathematically tractable and also practical.

IV Simulation Results

For comparison purposes two criteria are considered. First, we evaluate the completion performance of different methods using the Normalized Mean Square Error (NMSE) defined as

[TABLE]

where $\widehat{{\overline{\bf{R}}}}$ is the estimation of the original matrix ${\overline{\bf{R}}}$. The other criterion is the Hamming distance ($\rm{hd}$) of the original haplotype and its estimate. Note that since the original haplotype can be either $\bf{h}$ or $-\bf{h}$, we calculate the Hamming distance of the estimated haplotype for both cases and choose the minimum one. Our results are compared with the matrix completion based methods [5],[8],[9],[13] and the most recent MEC based algorithm [4].

Simulations are performed with synthetic data. To generate ${\overline{\bf{R}}}_{m\times n}$ in (2), we randomly generate ${\bf{h}}_{n\times 1}$ and ${\bf{c}}_{m\times 1}$. The random set of observations $\Omega$ is produced while the probability of observation is $\rm{pd}$. Moreover, erroneous observations $\Omega_{E}$ are defined by changing the sign of some entries observed by $\Omega$, where $\Omega_{E}\subset\Omega$. Also, we have set $\gamma_{1}=0.5$, $\gamma_{2}=2$, and $p=1.2$. For the initial matrix ${\bf{X}}_{0}$, we use the rank-one approximation of ${\rm{P}}_{\Omega}({\overline{\bf{R}}}_{E})$. Simulation results show that this initial matrix satisfies the conditions of Theorem 2. To obtain smooth curves, the results are averaged over 50 independent trials, each one with different sets of $\Omega$ and $\Omega_{E}$. Figs. 2 and 2 are depicted for $m=250$, $n=300$, and $0.05\leq$pd$\leq 0.2$. 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.25$, where $\left|\cdot\right|$ denotes the cardinality of a set. One can observe that the proposed method outperforms the other ones in estimating haplotypes by generating lower NMSEs and $\rm{hd}$s. As seen in Fig. 2, the MEC based algorithm in [4] can only compete with our proposed algorithm for the observation probabilities less than 0.07. This is in accordance with Theorem 2 of [4] which states that for a nonzero observation error probability, the increment of observation probability would increase the haplotype estimation error.

To generate Figs. 4 and 4, we change ${{\left|\Omega_{E}\right|}\mathord{\left/{\vphantom{{\left|\Omega\right|}{\left|{{\Omega}}\right|}}}\right.\kern-1.2pt}{\left|{{\Omega}}\right|}}$ from $0.14$ to $0.28$ and proportionally $p$ from $1.05$ to $1.2$, and set $\rm{pd}=0.07$. Once more, these results demonstrate the outperformance of the proposed method. As mentioned in Theorem 2, ${\gamma_{1}}<{2\mathord{\left/{\vphantom{2\pi}}\right.\kern-1.2pt}\pi}$ is required for our convergence criterion. Moreover, both ${\gamma_{1}}$ and ${\gamma_{2}}$ can be selected arbitrarily, so that the function ${\gamma_{1}}{\tan^{-1}}\left({{\gamma_{2}}{\bf{X}}}\right)$ can satisfyingly approximate ${\rm{sign}}\left({\bf{X}}\right)$. However, based on our simulation results, for a large value of ${\gamma_{2}}$ the runtime increases significantly, and for a very small value of ${\gamma_{2}}$, we will lose the haplotype estimation accuracy.

V Conclusion

The haplotype assembly problem was investigated. A new matrix completion minimization problem was proposed which benefits from an error correction mechanism over the manifold of rank-one matrices. The convergence of an iterative algorithm for this problem was theoretically proved. Simulation results demonstrated the validation of the proposed optimization problem by generating more accurate haplotype estimates compared to some recent related methods.