Manifold Based Low-rank Regularization for Image Restoration and Semi-supervised Learning

TL;DR

This paper introduces a manifold-based low-rank regularization technique that effectively handles nonlinear data structures, improving performance in various image restoration and semi-supervised learning tasks.

Contribution

It proposes a novel manifold-based low-rank regularization method that offers greater flexibility for nonlinear data, outperforming existing approaches in image and data science applications.

Findings

Effective in image inpainting and super-resolution

Improves X-ray CT image reconstruction quality

Enhances semi-supervised learning accuracy on handwritten digits

Abstract

Low-rank structures play important role in recent advances of many problems in image science and data science. As a natural extension of low-rank structures for data with nonlinear structures, the concept of the low-dimensional manifold structure has been considered in many data processing problems. Inspired by this concept, we consider a manifold based low-rank regularization as a linear approximation of manifold dimension. This regularization is less restricted than the global low-rank regularization, and thus enjoy more flexibility to handle data with nonlinear structures. As applications, we demonstrate the proposed regularization to classical inverse problems in image sciences and data sciences including image inpainting, image super-resolution, X-ray computer tomography (CT) image reconstruction and semi-supervised learning. We conduct intensive numerical experiments in several image restoration problems and a semi-supervised learning problem of classifying handwritten digits using the MINST data. Our numerical tests demonstrate the effectiveness of the proposed methods and illustrate that the new regularization methods produce outstanding results by comparing with many existing methods.

Equations

$$\mathcal{P}:{\mathbb{R}}^{m\times n}$$

$$f$$

$$\underset{\mathcal{M}\subset {\mathbb{R}}^{{\tau }^2},f}{\min }\sum _{x\in \mathcal{I}}\mathrm{Rank}(R_{\mathcal{M},x}(\mathcal{P}(f)),\text{s.t.}\mathcal{P}(f)\subset \mathcal{M},\mathcal{D}(f)=g,$$

$$({\nabla }_{\mathcal{M}}f)(x,y):=(f(y)-f(x))\sqrt{\omega (x,y)},x,y\in \mathcal{I}.$$

$$\underset{\mathcal{M}\subset {\mathbb{R}}^{{\tau }^2},f}{\min }\sum _{x\in \mathcal{I}}\Vert R_{\mathcal{M},x}(\mathcal{P}(f)){\Vert }_{\ast }+\frac{\lambda }{2}\Vert {\nabla }_{\mathcal{M}}f{\Vert }_2^2,\text{s.t.}\mathcal{P}(f)\subset \mathcal{M},\mathcal{D}(f)=g,$$

$$\underset{f,\mathcal{M}\subset {\mathbb{R}}^{{\tau }^2}}{\min }\sum _{x\in \mathcal{I}}\Vert R_{\mathcal{M},x}(\mathcal{P}(f)){\Vert }_{\ast }+\frac{\lambda }{2}\Vert {\nabla }_{\mathcal{M}}f{\Vert }_2^2,\text{s.t.}\mathcal{P}(f)\subset \mathcal{M},f{|}_{\Omega }=h{|}_{\Omega }.$$

$$\{\begin{array}{l}{f}^{k+1}=\arg {\min }_f{\sum }_{x\in \mathcal{I}}\Vert R_{{\mathcal{M}}^k,x}(\mathcal{P}(f){\Vert }_{\ast }+\frac{\lambda }{2}\Vert {\nabla }_{\mathcal{M}}^{k}f{\Vert }_2^2,\\ \text{s.t.}\mathcal{P}(f)\subset {\mathcal{M}}^k,f{|}_{\Omega }=h{|}_{\Omega },\\ {\mathcal{M}}^{k+1}=P(f^{k+1}).\end{array}$$

$$\underset{f,\mathbf{\alpha }}{\min }\sum _{x\in \mathcal{I}}\Vert R_{{\mathcal{M}}^k,x}\mathbf{\alpha }{\Vert }_{\ast }+\frac{\lambda }{2}\Vert {\nabla }_{{\mathcal{M}}^k}f{\Vert }_2^2,\text{s.t.}\mathcal{P}(f)=\mathbf{\alpha },f{|}_{\Omega }=h{|}_{\Omega }.$$

$$\mathcal{Q}(\mathbf{\alpha })=\{{\mathcal{Q}}_x(\mathbf{\alpha })=R_{{\mathcal{M}}^k,x}\mathbf{\alpha },x\in \mathcal{I}\}$$

$$\underset{f,\{{\beta }_x\}}{\min }\sum _{x\in \mathcal{I}}\Vert {\beta }_x{\Vert }_{\ast }+\frac{\lambda }{2}\Vert {\nabla }_{{\mathcal{M}}^k}f{\Vert }_2^2,\text{s.t.}{\mathcal{Q}}_x(\mathcal{P}(f))={\beta }_x,f{|}_{\Omega }=h{|}_{\Omega }.$$

$$\begin{array}{r}\underset{f,\{{\beta }_x\}}{\min }\underset{\{D_x\}}{\max }\sum _{x\in \mathcal{I}}\Vert {\beta }_x{\Vert }_{\ast }+\frac{\lambda }{2}\Vert {\nabla }_{{\mathcal{M}}^k}f{\Vert }_2^2+\sum _{x\in \mathcal{I}}\frac{\mu }{2}\Vert {\mathcal{Q}}_x(\mathcal{P}(f))-{\beta }_x+D_x{\Vert }_2^2,\\ \text{s.t.}f{|}_{\Omega }=h{|}_{\Omega }.\end{array}$$

$$\{\begin{array}{l}{\beta }_x^{l+1}=\arg \underset{{\beta }_x}{\min }\Vert {\beta }_x{\Vert }_{\ast }+\frac{\mu }{2}\Vert {\beta }_x-{\mathcal{Q}}_x(\mathcal{P}(f^l))-D_x^l{\Vert }_2^2,\forall x\in \mathcal{I},\\ f^{l+1}=\arg \underset{f}{\min }\frac{\lambda }{2}\Vert {\nabla }_{{\mathcal{M}}^k}f{\Vert }_2^2+\sum _{x\in \mathcal{I}}\frac{\mu }{2}\Vert {\mathcal{Q}}_x(\mathcal{P}(f))-{\beta }_x^{l+1}+D_x^l{\Vert }_2^2,\text{ s.t.}f{|}_{\Omega }=h{|}_{\Omega },\\ D_x^{l+1}=D_x^l+({\mathcal{Q}}_x(\mathcal{P}(f^{l+1}))-{\beta }_x^{l+1}),\forall x\in \mathcal{I}.\end{array}$$

$${\beta }_x^{l+1}={\mathcal{T}}_{1/\mu }({\mathcal{Q}}_x(\mathcal{P}(f^l))+D_x^l).$$

$${\mathcal{T}}_t(X)=U{S}_{\mathcal{T}}V,S_{\mathcal{T}}=\max (S-t,0).$$

$$\{\begin{array}{l}\left(-\lambda {\Delta }_{{\mathcal{M}}^k}+\sum _{x\in \mathcal{I}}\mu {\mathcal{P}}^{\top }{\mathcal{Q}}_x^{\top }{\mathcal{Q}}_x\mathcal{P}\right)f=\mu {\mathcal{P}}^{\top }\left(\sum _{x\in \mathcal{I}}{\mathcal{Q}}_x^{\top }({\beta }_x^{l+1}-D_x^l)\right).\text{vspace}0.2cm\\ f{|}_{\Omega }=h{|}_{\Omega }.\end{array}$$

$$\{\begin{array}{l}\left(-\lambda {\Delta }_{{\mathcal{M}}^k}+\mu W\right)(f)=\mu {\mathcal{P}}^{\top }\left(\sum _{x\in \mathcal{I}}{\mathcal{Q}}_x^{\top }({\beta }_x^{l+1}-D_x^l)\right),\\ f{|}_{\Omega }=h{|}_{\Omega },\end{array}$$

$$f^{l+1}{|}_{{\Omega }^c}=(A{|}_{{\Omega }^c}{)}^{-1}(\mu {\mathcal{P}}^{\top }(\sum _{x\in \mathcal{I}}{\mathcal{Q}}_x^{\top }({\beta }_x^{l+1}-D_x^l)-A{|}_{\Omega }h{|}_{\Omega }).$$

$$\underset{f,\mathcal{M}\subset {\mathbb{R}}^{{\tau }^2}}{\min }\sum _{x\in \mathcal{I}}\Vert R_{\mathcal{M},x}(\mathcal{P}(f)){\Vert }_{\ast },\text{s.t.}\mathcal{P}(f)\subset \mathcal{M},\mathcal{A}f=g.$$

$$\begin{array}{r}\underset{f,\mathbf{\beta }}{\min }\underset{\{D_{1,x}\},D_2}{\max }\sum _{x\in \mathcal{I}}\Vert {\beta }_x{\Vert }_{\ast }+\frac{{\mu }_1}{2}\sum _{x\in \mathcal{I}}\Vert {\mathcal{Q}}_x(\mathcal{P}(f))-{\beta }_x+D_{1,x}{\Vert }_2^2+\frac{{\mu }_2}{2}\Vert \mathcal{A}f-g+D_2{\Vert }_2^2.\end{array}$$

$${\mathbf{\varphi }}_i(x)=\{\begin{array}{l}1,L(x)=i.\\ 0,\text{otherwise}.\end{array},x\in S,i=0,1,2,\dots ,l.$$

$$L(x)=\arg \underset{i}{\max }{\mathbf{\varphi }}_i(x),\forall x\in P\backslash S.$$

$$\underset{\Phi }{\min }\sum _{x\in \mathcal{I}}\Vert (R_{\mathcal{M},x})\Phi {\Vert }_{\ast },\text{s.t.}P\subset \mathcal{M},\Phi (x,i){|}_{x\in S}=\{\begin{array}{l}1,L(x)=i.\\ 0,\text{otherwise}.\end{array}$$

$$\begin{array}{r}\underset{\Phi ,\{{\psi }_x\}}{\min }\underset{\{D_x\}}{\max }\sum _x\left(\Vert {\psi }_x{\Vert }_{\ast }+\frac{\mu }{2}\Vert {\psi }_x-{\mathcal{Q}}_x(\Phi )-D_x{\Vert }_2^2\right)\\ \text{s.t.}\Phi (x,i){|}_{x\in S}=\{\begin{array}{l}1,L(x)=i.\\ 0,\text{otherwise}.\end{array}\end{array}$$

$$\{\begin{array}{l}{\psi }_x^{k+1}=\arg {\min }_{{\psi }_x}\Vert {\psi }_x{\Vert }_{\ast }+\frac{\mu }{2}\Vert {\psi }_x-{\mathcal{Q}}_x({\Phi }^k)-D_x^k{\Vert }_2^2,\forall x\in P,\\ {\Phi }^{k+1}=\arg {\min }_{\Phi }{\sum }_x\frac{\mu }{2}\Vert {\psi }_x^{k+1}-{\mathcal{Q}}_x(\Phi )-D_x^k{\Vert }_2^2,\text{s.t.}\Phi (x,i){|}_{x\in S}=\{\begin{array}{l}1,L(x)=i.\\ 0,\text{otherwise}.\end{array}\\ D_x^{k+1}=D_x^k+{\mathcal{Q}}_x({\Phi }^{k+1})-{\psi }_x^{k+1},\forall x\in P.\end{array}$$

$$\{\begin{array}{l}{\psi }_x^{k+1}={\mathcal{T}}_{1/\mu }({\mathcal{Q}}_x({\Phi }^k)-D_x^l),\forall x\in P,\\ {\Phi }^{k+1}=\tilde{\mathcal{Q}}(\{{\psi }_x^{k+1}-D_x^k{\}}_x){\chi }_{S^c}+{\Phi }^0{\chi }_S,\\ D_x^{k+1}=D_x^k+{\mathcal{Q}}_x({\Phi }^{k+1})-{\psi }_x^{k+1},\forall x\in P.\end{array}$$

$$\text{mbox}PSNR(f,\tilde{f})=10{\log }_{10}\frac{MN(f_{\max }-f_{\min }{)}^2}{\Vert f-\tilde{f}{\Vert }_2^2},$$

Full text

1 Introduction

Regularization methods play important roles in many ill-posed inverse problems arising in science and engineering. Examples include inverse problems considered in signal processing and image sciences such as image denoising, image impainting, image deconvolution [13, 1], just to name a few. Mathematically, a image restoration problem can be viewed as reconstructing a clean image $f$ from a degraded image $g$ based on the degradation relationship $\mathcal{D}(f)=g$. It is challenging to reconstruct $f$ from $g$ as the problem is usually ill-posed due to the highly underdetermined constraints and possible noise. Observations of natural image with prior information such as piecewise smoothness, shape edges, textures, repetitive patterns and sparse representations under certain transformations make regularization methods quite effective to handle image processing problems. Successful methods include the total variation (TV) methods, nonlocal methods and wavelet tight frame methods [40, 4, 22, 19] and many others. Moreover, regularization methods can also be considered in problems arising from data science. A typical example is semi-supervised learning, where tasks aim at labeling data from a small amount of labeled training data set. Regularization methods such as the harmonic extension method [47] have been considered to this type of ill-posed problem. In this paper, we consider a different regularization, called manifold based low-rank (MLR) regularization as a linearization of manifold dimension, which generalizes the global low-rank prior knowledge for linear objects to manifold-region-based locally low-rank for nonlinear objects.

The idea of the MLR proposed in this paper is inspired by a recent method called the low-dimensional manifold model (LDMM) discussed in [36]. Using the image patches discussed in nonlocal methods [4, 37], the LDMM interprets image patches as a point cloud sampled in a low-dimensional manifold embedded in a high dimensional ambient space, which provides a new way of regularization by minimizing the dimension of the corresponding image patch manifold. This can be explained as a natural extension of the idea of low-rank regularization for linear objects to data with more complicated structures. Moreover, the authors in [36] elegantly find that the point-wisely defined manifold dimension can be computed as a Dirichlet energy of the coordinate functions on the manifold, whose corresponding boundary value problem can be further solved by a point integral method proposed in [34]. The LDMM performs very well in image inpainting and super-resolution. This model is later considered in collaborative ranking problems [31]. Based on weighted graph laplacian (WGL), an improvement of LDMM called LDMM+WGL is proposed more recently in [42].

In this paper, instead of representing the manifold dimension as a manifold-derivative involved quantity [36], we propose a linear approximation of the manifold dimension. Note that the quantity of the dimension at each point $x\in\mathcal{M}$ is the same as the dimension of the tangent space at $x$. This quantity only depends on a local neighborhood of $x$ on $\mathcal{M}$, which can be approximated as the rank of the covariance matrix generated by the set of $K$-nearest-neighbourhood (KNN) points of $x$ on $\mathcal{M}\subset\mathbb{R}^{n}$ in the discretized sense. In other words, the low-dimensional property of $\mathcal{M}$ at $x$ can be linearly approximated as the low-rank property of the this corresponding covariance matrix, which is essentially the same as the low-rank property of the matrix $R_{\mathcal{M},x}$ formed by those KNN points near $x$. As an example illustrated in Figure 1, we construct a patch manifold of the Barbara image using patch size $11\times 11$. This leads to a set of image patches represented as a point clouds in $\mathbb{R}^{121}$. The rank of $R_{\mathcal{M},x}$ for the Barbara image is color-coded in the right image of Figure 1, which clearly illustrates that $\mathrm{Rank}(R_{\mathcal{M},x})$ has low value for this natural image. As a linear approximation of the $Dim_{\mathcal{M}}(x)$ proposed in [36], the manifold based quantity $\mathrm{Rank}(R_{\mathcal{M},x})$ does not involve with any manifold differential operators, which has potential to apply this concept to more general data processing problem such as a preliminary example demonstrate in section 3. On the other hand, this consideration is reasonable as the globally defined “Rank” can only handle linear objects, while this manifold based locally defined $\mathrm{Rank}$ has advantages to regularize data with nonlinear structures.

Based on the MLR prior knowledge, we use the matrix nuclear norm relaxation for matrix rank as the method considered in low-rank matrix completion theory [9] and apply MLR to the image patch manifold for image restoration problems including image inpainting, image super-resolution and X-ray computer tomography (CT) image reconstruction. It is clear the definition MLR relies on the construction of KNN which is essentially dependent on the manifold structure. Therefore, a split-Bregman method [26] is considered to solve the proposed model by iteratively updating the manifold structure and the objective image. Moreover, we also apply the proposed regularization for a semi-supervised learning problem, where MLR is applied to a labeling matrix with a fixed manifold structure provided by the input data. Our numerical results tested for a benchmark data set of handwritten digits illustrate the effectiveness of the proposed method.

The rest part of this paper is organized as follows. In Section 2, we discuss our manifold based low-rank regularization for the image restoration problems including image inpainting, image super-resolution and X-ray CT image reconstruction. Detailed models and numerical algorithms for various image processing problems are discussed. In Section 3, we consider the manifold based low-rank model to a semi-supervised learning problem. Intensive numerical experiments and comparisons with existing methods are conducted in Section 4. We conclude our work in Section 5.

2 Manifold based low-rank regularization for image restoration

In this section, we consider the MLR method for image restoration problems including image inpainting, image super-resolution and X-ray CT image reconstruction. The idea of MLR is applied to a image patch manifold with a fixed patch size similar as the way proposed in [36]. We further relax the problem of matrix rank minimization as a problem of matrix nuclear norm optimization and solve the proposed optimization problem based on the split Bregman iteration [26] and the singular value thresholding algorithm [6].

The classical image restoration models mainly focus on local properties of the objective image such as smoothes and jumps. Image features can be further enhanced due to its possible repetitive patterns non-locally. The nonlocal based image restoration methods [4, 15, 22] extract and match non-local repetitive structures of images using image patches. Given a discrete image $f\in\mathbb{R}^{m\times n}$ defined on a domain $\mathscr{I}=\{1,2,\ldots,m\}\times\{1,2,\ldots,n\}$, a size $\tau=2\eta+1$ patch transform $\mathcal{P}$ can be defined by:

[TABLE]

where $x$ is the center of each patch, $\mathscr{P}=\{-\eta,-\eta+1,\ldots,0,1,\ldots,\eta-1,\eta\}^{2}$ represents the patch index set and $\tilde{f}\in\mathbb{R}^{(m+2\eta)\times(n+2\eta)}$ is a proper extension (symmetric extension in this paper) of $f$ such that $\tilde{f}(x)=f(x),\forall x\in\mathscr{I}$. An essential observation of nonlocal methods is that images can be restored by enhancing similar patterns which may not lie in nearby regions of $\mathscr{I}$ domain. Therefore, comparing with the direct regularization methods on the image domain of $f$, the quality of image restoration can be usually improved using nonlocal methods. For instance, nonlocal based variational methods [4, 22, 46] and nonlocal based wavelet frame based methods [38] demonstrate outstanding image restoration results.

Given a patch matrix $\mathcal{P}(f)$, one can regard each patch $\mathcal{P}(f)(\cdot,x)$ as a $\tau^{2}$ dimensional column vector. Consequently, $\mathcal{P}(f)$ can be viewed as a set of points in $\mathbb{R}^{\tau^{2}}$. To conduct further analysis of this point cloud, we model $\mathcal{P}(f)$ as a set of points sampled on a manifold $\mathcal{M}\in\mathbb{R}^{\tau^{2}}$. Thereafter, we also abuse the notation $x$ as the corresponding point $\mathcal{P}(f)(\cdot,x)$ on $\mathcal{M}$. This manifold interpretation has been proposed in existing work [37, 36]. More recently, [36] proposes a low dimensional manifold model (LDMM) for image restoration. This work observes that the dimension of patch manifold $\mathcal{M}$ should intrinsically have a low-dimensional structure and proposes to regularize the dimension of the patch manifold $\mathcal{M}$ for image restoration. Moreover, the authors elegantly show that the dimension function $Dim(\mathcal{M})$ at $x\in\mathcal{M}$ can be represented by $Dim_{\mathcal{M}}(x)=\sum_{1\leq s\leq\tau^{2}}\|\nabla_{\mathcal{M}}(\mathcal{P}(f)(s,\cdot))(x)\|_{2}^{2}$, which transforms the dimension regularization problem to be a variational partial differential equation model that is proposed to solve using a point integral method discussed in [34]. Later on, [31] generalized the LDMM model into matrix completions with better performance than traditional low-rank regularized model in completing the Netflix matrix [2] which does not have exactly global low rank.

2.1 Manifold based low-rank regularization for the patch manifold

Inspired by the regularization of the manifold dimension represented as a manifold derivative involved quantity [36], we propose a linear approximation of the manifold dimension in the following way. Note that the quantity of dimension at each point $x\in\mathcal{M}$ is the same as the dimension of the tangent space $\mathcal{T}_{x}\mathcal{M}$ at $x$ which only relies on a local neighborhood of $x$ on $\mathcal{M}$. In the discrete sense of $\mathcal{M}$ sampled as the patch matrix $\mathcal{P}(f)$, the quantity $\mathrm{dim}(\mathcal{T}_{x}\mathcal{M})$ can be approximated as the rank of the covariance matrix generated by the set of $K$-nearest-neighbourhood (KNN) points of $x$ in $\mathcal{P}(f)$. In other words, the low-dimensional property of $\mathcal{M}$ at $x$ can be linearly approximated as the low-rank property of the corresponding covariance matrix, which is essentially the same as the low-rank property of the matrix formed by those KNN points near $x$. More precisely, if we define the restriction operator $R_{\mathcal{M},x}$ as the KNN points near $x$, then the low-dimensional prior knowledge of the patch manifold $\mathcal{M}$ at $x$ can be linearly approximated as the low-rank prior knowledge of the matrix formed by points in $R_{\mathcal{M},x}$ denoted as $R_{\mathcal{M},x}(\mathcal{P}(f))$. Namely, we define the manifold based rank at $x$ as $\text{Rank}_{\mathcal{M}}(x)=\text{Rank}(R_{\mathcal{M},x}(\mathcal{P}(f)))$.

For image restoration problems, if the fidelity information is $\mathcal{D}(f)=g$ as a constraint where $\mathcal{D}$ is a degradation operator, by regularizing $\text{Rank}_{\mathcal{M}}(x)$ for all the point $x$, we consider the following the manifold based low-rank regularization for image restoration:

[TABLE]

On the one hand, the minimization of the rank, or the $\ell_{0}$ norm of the singular value, is NP hard to be optimized generally. Therefore, $\ell_{1}$ norm of the singular value, or the nuclear norm of the localized matrix, is an appropriate way to relax the local rank as the pioneer work of low-rank matrix completion theory developed in [9]. The minimization of nuclear norm can be solved by applying the singular value thresholding (SVT) algorithm [6]. On the other hand, we observe that it is necessary to smoothen the images, or enhance the features and textures in practice. Therefore, one can apply some positive/negative diffusion based regularization [24] to $f$ to guarantee the smoothness of the object image. For example, we choose the diffusion term as the non-local gradient operator defined in (3).

[TABLE]

Therefore, a MLR image restoration model can be stated as:

[TABLE]

when $\lambda>0$ the regularization term $\frac{\lambda}{2}\|\nabla_{\mathcal{M}}f\|_{2}^{2}$ represents a diffusion term which can smoothen the regions. When $\lambda<0$, the regularization represents inverse diffusion which can enhance the patterns [24, 23, 5]. Otherwise, $\lambda=0$ leads the model (4) identical to model (2) as a pure MLR regularized image restoration model.

We remark that a close related work [20] imposed the low-rank regularization in a non-local transform domain of images, which is applicable to recover images from missing Fourier coefficients. In particular, to improve the robustness of the algorithm, the low-rank regularization is considered to the grouped “similar patches”, which can be regarded as a type of “locally low-rank regularization” although [20] did not explicitly view the “low-rank” in manifold sense. In addition, this method considers to group patches without sufficient overlapping, thus it only includes a rough sampling on the patch manifold which may not be able to accurately reflect the low-dimensional structure of the patch manifold.

2.2 MLR for image inpainting

Image inpainting [3] is a process to restore images whose pixels are missing, over-written or corrupted. More precisely, the inpainting problem aims at reconstructing an image $f$ only based on its partial information on a given set $\Omega\subset\mathscr{I}$. Such ill-posed problem is generally based on some assumptions such that the object image $f$ is piecewise smooth, or has repetitive textures. With these assumptions, regularization inpainting methods, such as variational PDE based models [12, 41, 43], wavelet based models [14, 10, 7, 8, 17] and low dimensional manifold model [36] have been proposed.

We would like to demonstrate that MLR model can restore images and preserve both the piecewise smooth regions and textures from a small random portion of information. As a special case of (4), the low-rank regularized image inpainting model can be stated as:

[TABLE]

In particular, if the index set $\Omega$ is picked as $\{1,s+1,2s+1,\ldots\}\times\{1,s+1,2s+1,\ldots\}$, the problem is called sub-sampled super-resolution problem. As the nuclear norm in the first term of the above problem depends on the manifold structure, we consider to solve this problem by alternatively updating the manifold $\mathcal{M}$ and solving $f$ similar as the method considered in [36]. The outline of solving (5) can be stated as follows:

[TABLE]

To solve $f^{k+1}$ from the first step in (6) with a fixed manifold structure $\mathcal{M}^{k}$, we use the split Bregman iteration [26]. After introducing an auxiliary variable $\bm{\alpha}=\mathcal{P}(f)\in\mathbb{R}^{\tau^{2}\times mn}$, this problem can be reinterpreted as:

[TABLE]

Since each column of $\bm{\alpha}$ may occur multiply times in different $\|R_{\mathcal{M}^{k},x}\bm{\alpha}\|_{*}$, it is difficult to simultaneously optimize several nuclear norms together. Therefore, denote the image size as $m\times n$, patch size as $\tau\times\tau$, and the KNN size $K$, we introduce the duplicate operator $\mathcal{Q}:\mathbb{R}^{\tau^{2}\times mn}\rightarrow\mathbb{R}^{K\tau^{2}\times mn}$ can be defined as:

[TABLE]

Then, we denote $\mathcal{Q}_{x}(\bm{\alpha})=\beta_{x},\forall x\in\mathscr{I}$ such that $\|(R_{\mathcal{M}^{k},x})\bm{\alpha}\|_{*}=\|{\beta}_{x}\|_{*}$. As a result, $\sum_{x\in\mathscr{I}}\|(R_{\mathcal{M}^{k},x})\bm{\alpha}\|_{*}=\sum_{x\in\mathscr{I}}\|{\beta}_{x}\|_{*}$ becomes a separable formula. Thus, Step 1 in (6) can be reinterpreted as:

[TABLE]

Therefore, the above the equality constraint $\mathcal{Q}_{x}(\mathcal{P}(f))={\beta}_{x}$ can be solved by considering the following saddle point problem using a augmented Lagrangian formula with the dual variable $\{D_{x}\}$:

[TABLE]

where $\mu$ is the parameter to control the augmented Lagrangian. Similar to the one-step iterative method in the alternating direction method of multipliers (ADMM) and split Bregman iteration [25, 26], The optimization problem (10) can be iteratively solved as:

[TABLE]

The first sub-optimization problem has a closed-form solution provided by the singular value thresholding [6]. Namely,

[TABLE]

where for any matrix $X$ with a singular value decomposition $X=USV$, the singular value thresholding operator $\mathcal{T}$ is provided as

[TABLE]

Next, we solve $f^{l+1}$ in (11), The solution of the linear constrained minimization problem satisfies the following Dirichlet boundary value problem:

[TABLE]

In (14), since the duplication operators $\{\mathcal{Q}_{x}\}$ have only one non-zero element in each row, we have that for all $x$, $(\mathcal{Q}_{x}^{\top}\mathcal{Q}_{x})_{ij}=\sum_{p}(\mathcal{Q}_{x}^{\top})_{ip}(\mathcal{Q}_{x})_{pj}=\sum_{p}(\mathcal{Q}_{x})_{pi}(\mathcal{Q}_{x})_{pj}$ which is always [math] if $i\neq j$. Therefore, $\sum\limits_{x\in\mathscr{I}}\mathcal{Q}_{x}^{\top}\mathcal{Q}_{x}=W_{\mathcal{Q}}$ becomes a diagonal weight matrix. Similarly, the patch manifold transform operator $\mathcal{P}$ also has only one non-zero element in each row. After left multiplied by a diagonal matrix, $\left(\sum\limits_{x\in\mathscr{I}}\mathcal{Q}_{x}^{\top}\mathcal{Q}_{x}\right)\mathcal{P}=W_{\mathcal{Q}}\mathcal{P}$ is still a matrix with only one non-zero element in each row. Therefore, $\sum\limits_{x\in\mathscr{I}}\mathcal{P}^{\top}\mathcal{Q}_{x}^{\top}\mathcal{Q}_{x}\mathcal{P}=W$ is a diagonal weight matrix for the input image whose entries is the occurrence of each pixel in all local regions of patch manifold $\{\beta_{x}\}$. We can consequently rewrite (14) as:

[TABLE]

Denote the left hand side of the linear system as $A=-\lambda\Delta_{\mathcal{M}^{k}}+\mu W$, plugging the boundary condition $f_{\Omega}=h_{\Omega}$ into the first equation, we can solve $f^{l+1}$ restricted in $\Omega^{c}$ as follows:

[TABLE]

Therefore, combining (6), (11), (12) and (16), we can solve the MLR based image inpainting model (5) as Algorithm 1. Note that the max number of inner iterations can be chosen as $1$ to reduce the computational time.

2.3 MLR for X-ray CT reconstruction

As a special case of image restoration, medical imaging plays important role in different clinical applications. Here, we consider an application of our method to X-ray Computed Tomography (CT), which aims at reconstructing images from their Radon transform. Mathematically, the X-ray CT reconstruction problem can be essentially represented as a linear inverse problem: $\mathcal{A}f=g$, where $\mathcal{A}\in\mathbb{R}^{m\times n}$ is a measurement matrix representing the collection of discrete line integrations with different projection angles and along different beamlets, $f\in\mathbb{R}^{n}$ is vectorized 2 dimensional image and $g\in\mathbb{R}^{m}$ is the corresponding measurement. Given the geometry matrix $\mathcal{A}$ and $g$, the task of X-ray CT reconstruction is to find an appropriate value of $f$ [39, 28]. In literature, there are some classical methods available, such as the filtered back projection (FBP) type methods [21, 16, 35, 33], the algebraic reconstruction techniques (ART) [27]. In practice, however, to minimize the radiation dose by reducing the number of projection angles and beamlets, the amount of measurement $m$ becomes much less than the dimension of the object image $n$, which makes the reconstruction becoming an under-determined problem with infinitely many solutions. As a result, previously mentioned FBP and ART methods usually suffer from artifacts because of the insufficient measurements. Regularization methods such as TV based medical imaging models [30] and wavelet regularization based medical imaging models [29, 18] makes it possible to reconstruct piecewise smooth or piecewise constant object images. However, it is still a big challenge to preserve tiny features due to possible over-smoothing, which motivate us to propose a MLR CT imaging model to preserve both smooth pieces and tiny features. This model is a special case of (4) with linear degradation operator $\mathcal{D}=\mathcal{A}$ and $\lambda=0$ as follows:

[TABLE]

Note that this model is also applicable for average filter based super resolution problem, Fourier domain inpainting problem, and image deconvolution problems.

To solve (17), similar as (10), after defining the duplication operator $\{\mathcal{Q}_{x}\}$ and localized patch manifold $\{\beta_{x}\}$, by splitting the linear constraints $\mathcal{Q}_{x}(\mathcal{P}(f))=\beta_{x},\forall x\in\mathscr{I}$ and $\mathcal{A}f=g$, we obtain the saddle point problem of model (17) using the augmented Lagrangian:

[TABLE]

Similar as Algorithm 1, applying the ADMM we can design algorithm 2 for solving CT reconstruction model (17).

3 Semi-supervised learning using MLR

As another advantage of the proposed MLR, this idea can be adapted to handle various data processing problem. Here, we propose the extension of this approach to a semi-supervised learning problem. Many other potential applications in data science will be investigated in our future work.

Semi-supervised learning is a learning paradigm aiming at labeling data from a small amount of labeled training data set [48]. Mathematically speaking, given a data set $P=\{x_{1},x_{2},\ldots,x_{n}\}\subset\mathbb{R}^{d}$, the semi-supervised learning problem is to find a label function $L:P\rightarrow\{0,1,2,\ldots,l\}$ representing the label index of the each $x_{i}$ with given prior knowledge of $L$ in a labeled subset set $S\subset P$. The challenge of a semi-supervised learning problem is to estimate an accurate assignment of $L$ based on a vey small portion information $L(S)$. The general idea of semi-supervised learning is to explore the manifold structure of the data based on an assumption that similar unlabeled samples should be assigned the same classification. Based on this, diffusion based models [47, 44, 43] has been considered to tackle this problem. In this section, we would like to formulate a different way of estimating $L$ from highly insufficient labeled samples based on the MLR method.

Similar as notations discussed in [47, 44, 43], to solve the semi-supervised learning problem, we define the cluster functions $\{\bm{\phi}_{i}(x)\}$ which is partially assigned from the training data $S$.

[TABLE]

By viewing $\bm{\phi}_{i}(x)$ a column vector with length $n$, we obtain a cluster matrix $\Phi=(\phi_{0},\cdots,\phi_{l})\in\mathbb{R}^{n\times(l+1)}$. Therefore, if we can estimate all the components of $\Phi$, or all $\{\phi_{i}(x)\}$, the value of all unknown $L(x)$ for $x\in P\backslash S$ can be estimated by:

[TABLE]

Assume the point matrix $P$ is sampled on a manifold $\mathcal{M}$ and define the local restriction operator $R_{\mathcal{M},x}$ as the restriction of a matrix to $x$-th point and its $K$-nearest neighbourhood (KNN). Then by definition of $\Phi$ and $\bm{\phi}_{i}(x)$, the rank of $R_{\mathcal{M},x}\Phi$ equals to the number of different labels occurred in the KNN. Based on the assumption that similar data samples or nearby points should have similar classification, localization of $\Phi$ should only include a few different labels, i.e., $R_{\mathcal{M},x}\Phi$ has low-rank structure although $\Phi$ might be a full-rank matrix. As an example, we consider the public available MINST data set [32] which includes $70,000$ handwritten digits images. We simply view each image as a point in $\mathbb{R}^{d}$ and pick the KNNs of each point (image) in terms of Euclidean distance. Left image in Figure 2 shows that majority part of $\{R_{\mathcal{M},x}\Phi\}$ has low-rank structure from the ground truth of cluster matrix $\Phi$. Interestingly, right image in Figure 2 shows that the $20$-nearest neighborhood of the first image, in which two digits $5$ and $3$ appear because of their similar distribution in terms of Euclidean distance. Therefore, the rank of $R_{\mathcal{M},1}\Phi$ equals to $2$.

Based on the observation that $R_{\mathcal{M},x}\Phi$ has low-rank structure, the corresponding MLR model for cluster matrix estimation can be stated as follows:

[TABLE]

Different from the previous image restoration models, the geometric of manifold $\mathcal{M}$ is only determined by information from the data set $P$ which is fixed and irrelevant to the evolution of $\Phi$. Correspondingly, with fixed localization of $\Phi$, the model (19) is convex and can be solved via standard ADMM. Since it is difficult to simultaneously minimize all the restrictions of $\Phi$, similar as the image restoration cases, we define a duplication operator $\mathcal{Q}=\{\mathcal{Q}_{x}\}_{x}$ such that $\mathcal{Q}_{x}\Phi=R_{\mathcal{M},x}\Phi=\psi_{x}$ and $\|R_{\mathcal{M},x}\Phi\|_{*}=\|\psi_{x}\|_{*}$. With the auxiliary variables $\{\psi_{x}\}$ and linear constraint $\mathcal{Q}_{x}\Phi=\psi_{x}$, we introduce a group of dual variables $\{D_{x}\}$ and obtain the following saddle point problem with the augmented Lagrangian:

[TABLE]

Similar as the image restoration case, with the definition of the duplication operator $\mathcal{Q}$, because $\mathcal{Q}^{\top}\mathcal{Q}=\sum_{x}\mathcal{Q}_{x}^{\top}\mathcal{Q}_{x}=W_{\mathcal{Q}}$ which is a diagonal matrix, we can define the left inverse operator as $\tilde{\mathcal{Q}}=W_{\mathcal{Q}}^{-1}\mathcal{Q}^{\top}$ such that $\tilde{\mathcal{Q}}\mathcal{Q}=I$. Standard ADMM brings the outline of the iteration as follows:

[TABLE]

In (21), the first step can be solved by singular value thresholding operator defined in (13) as $\psi_{x}^{k+1}=\mathcal{T}_{1/\mu}(\mathcal{Q}_{x}(\Phi^{k})-D_{x}^{k})$. The equality constraint $\Phi(x,i)|_{x\in S}=\begin{cases}1,\ \ \ L(x)=i.\\ 0,\ \ \ \text{otherwise}.\end{cases}$ in the second step is an orthogonal projection operator. Therefore, ${\Phi}^{k+1}=\tilde{\Phi}^{k+1}\chi_{S^{c}}+\Phi^{0}\chi_{S}$, where $\tilde{\Phi}^{k+1}=(\sum_{x}\mathcal{Q}_{x}^{\top}\mathcal{Q}_{x})^{-1}(\sum_{x}\mathcal{Q}_{x}^{\top}({\psi_{x}^{k+1}}-D_{x}^{k}))=W_{{\mathcal{Q}}}^{-1}\sum_{x}\mathcal{Q}_{x}^{\top}({\psi_{x}^{k+1}}-D_{x}^{l}))=\tilde{\mathcal{Q}}(\{\psi_{x}^{k+1}-D_{x}^{k}\}_{x}).$ Then the iteration can be re-sketched as:

[TABLE]

Given an appropriate initialization and sufficient iterations, we obtain the solution of $\Phi$ and the corresponding columns $\phi_{i}(x)$. Therefore, the index set $L(x)$ for $x\notin S$ can be estimated by $L(x)=\max_{i\in\{0,1,2,\ldots,l\}}\phi_{i}(x),x\notin S$, which completes the full estimation of $L(x)$.

It is clear that a better initial guess of $\Phi^{0}$ can further improve the index completion result. Therefore, we propose to recursively update the initial guess $\Phi^{0}$ based on the result from (22), the ultimate algorithm for semi-supervised learning can be summarized in Algorithm 3.

4 Numerical Experiments

In this section, we conduct numerical experiments for the proposed MLR models to various image restoration problems, X-ray CT imaging and semi-supervised learning. Our results validate that the proposed method can successfully reduce the reconstruction error and preserve both edges and repetitive patterns. For all image restoration results, besides the visual quality, we also quantitatively evaluate the results of image restoration using the peak signal-to-noise ratios (PSNR) value:

[TABLE]

with the ground truth image $\tilde{f}$, where $f_{\max}$ and $f_{\min}$ are its maximal and minimal pixel values respectively and $M$, $N$ are the size of the image. All the numerical simulations are implemented by MATLAB in a PC with 32GB RAM and 2.7 GHz CPUs.

4.1 Image inpainting and super-resolution

In the first experiment, we test Algorithm 1 to inpaint images from random missing pixels, in which the index set $\Omega$ is uniformly randomly chosen with fixed rate. Figure 3 shows the restoration results of Barbara image from same $10\%$ random available pixels using different methods. It can be seen that the traditional wavelet based method [7], the classical harmonic extension method and TV based method [11] cannot preserve the textures in this low rate of available information because given information in the texture part is recognized as some noise in these two restored images. Both purely manifold based low-rank model and the LDMM method [36] have much better estimation and preservation of the textures, while the low-rank regularization of the patch manifold may generate some artifacts which breaks some smooth regions. The proposed method include both manifold based low-rank and inverse diffusion ($\lambda=-20$ for image inpainting) can enhance the recovered image to obtain a better texture and smooth region representation. Our method provides comparable results with the most recent proposed LDMM + Weighted graph laplacian (LDMM+WGL) method [42].

Due to non-convexity of the model, we also numerically verify the convergence of the algorithm 1. For the numerical simulations shown as above, the convergence curves of the object function $\sum_{x\in\mathscr{I}}\|{\beta}_{x}\|_{*}+\frac{\lambda}{2}\|\nabla_{\mathcal{M}^{k}}f\|_{2}^{2}$ and the relative error of linear constraints $\sum_{x\in\mathscr{I}}\|{\mathcal{Q}}_{x}(\mathcal{P}(f))-{\beta}_{x}\|_{2}$ are shown in Figure 4, which validate that for the proposed Algorithm 1, the object function converges to a stable value and the relative error of linear constraint converges to zero.

We further test the proposed model for different level of available information and conduct comparisons with the LDMM method. Figure 5 shows other Barbara image inpainting results from $5\%,20\%$ and $40\%$ random available information. In the case of using $5\%$ available information, the MLR model produces a qualitatively and quantitatively better result than the one obtained from LDMM. However, the image from LDMM+WGL method has the highest PSNR, although it qualitatively produces more artifacts near the mouth region. In the case of using $10\%$ available information, although the proposed MLR model produces an image with the highest PSNR value, it is hard visually distinct results from MLR and LDMM+WGL. Thus, MLR and LDMM+WGL are comparable and better than LDMM in this case. MLR and LDMM+WGL methods produce similar high quality results when the sampling rate increases to $20\%$ available information although this rate of information may also be quite challenging to other existing methods. All three methods produces very good results with $40\%$ information. Moreover, we also apply the proposed image inpainting model to other images to test the capability of the MLR for handling texture and carton parts. For images with more textures such as the fingerprint image, the baboon image and the boat image, Figure 6 shows that the proposed MLR method can still preserve more features. In particular, at the bottom part of the fingerprint image highlighted by the red box, the LDMM method generates some vertical artifacts while the MLR method produce more accurate estimation. The LDMM+WGL method successfully improves the inpainting results from the LDMM method, but some vertical artifacts still remain. For the boat image in Figure 6, we observe that the proposed MLR method can restore more isolated line structures on the top of the boat as highlighted by red boxes while the LDMM method tends to remove the thin lines. The LDMM+WGL method produces a comparable result with the one from MLR method. For the baboon image, since the texture is too tiny and not repeated frequently, all methods do not provide a result with clear skin and beard structure. The LDMM+WGL method seems to enlarge the artifacts in this case. On the other hand, for images with less textures such as the peppers image, Figure 6 shows that the proposed method can reduce the possibility of generating artifacts which should not exist. For example, at the center of the green pepper (highlighted by the red box), and at the center of the camera support (highlighted by the red box), the artifacts from the LDMM method and the LDMM+WGL method break the smooth regions while the proposed MLR method preserves the smooth parts because the smooth regions also include repetitive patterns and formulate the low-rank structure.

Additionally, we also implement the MLR method for image inpainting from manual scratches. Figure 7 shows that compared to the wavelet based image inpainting model [7], the proposed model has much better quality of recovering the fingerprint structure in terms of both the visualization and the PSNR value. Moreover, for the second row with wider scratches, the proposed MLR model has better estimation of the fingerprint pattern other than simply smoothen the scratched regions.

In the second experiment, we show the results of super-resolution. In [36], the authors conduct the super resolution as a special type of image inpainting problem with highly coherent fixed index set $\Omega=\{1,s+1,2s+1,\ldots\}\times\{1,s+1,2s+1,\ldots\}$. Using the same model and algorithm as the image inpainting problem, the results of this super-resolution problem from sub-sampled pixel are shown as follows in 8. It can be seen that the super-resolution result is better than results from traditional bi-cubic interpolation and comparable to results from the LDMM method and the LDMM+WGL method.

As another case of super resolution, the problem is assumed as image restoration from filtered low resolution version of images. Define an average operator $\mathcal{A}$, the input low resolution image $f_{L}=\mathcal{A}(f)$, which provide a linear constraint fidelity condition and similar as the medical imaging model (17). Using the formula (18) and applying Algorithm 2, the super resolution results from $4\times 4$ and $8\times 8$ average filtered low resolution images are shown in Figure 9. The proposed MLR method produces more detailed information and sharper images than bi-cubic interpolation and LDMM in [36].

4.2 X-ray CT Reconstruction

It is quite challenging to reconstruct satisfactory image for the X-ray CT problem with a small amount of radiation dose. In this section, we apply the model (18) and Algorithm 2 to the fan-beam projection measurement of images with reduced number of projection views. We consider the CT imaging for a human chest slice (See Figure 10) from the data of ”Low Dose CT Grand Challenge” provided by Dr. Cynthia McCollough, the Mayo Clinic, the American Association of Physicists in Medicine, and supported by grants EB017095 and EB017185 from the National Institute of Biomedical Imaging and Bioengineering. Regarding to the linear fidelity $\mathcal{A}f=g$, the ground truth image and the object image $f$ has resolution $256\times 256$ and the Radon transform measurement $g$ in this section always includes $512$ projection lines in each projection view. Therefore, $\#PROJ$ projection views represents the measurements has $\frac{\text{Card}(g)}{\text{Card}(f)}=\frac{\#PROJ\times 512}{256^{2}}=\frac{\#PROJ}{128}$ portion of the object function. The huge sparse geometric matrix $\mathcal{A}$ is generated by Siddon’s method [45] as pre-process.

For CT imaging from $15$, $30$ and $60$ views, the CT reconstruction results from the proposed MLR method are shown in Figure 10. It can be seen that the proposed model performs better than the wavelet based method [18] in term of both the visual quality and the PSNR value. For wavelet based method, stronger regularization as in Figure 10 would remove the small features since they would be recognized as artifacts or noise, while weaker regularization cannot remove the artifacts caused by insufficient projection angles. In particular, in the case with 15 projections, the wavelet based method cannot recover the main vessels at the right side while our method still produce very good results. Moreover, for 60 projections, the zoom-in part shows that the proposed model can successfully reconstruct these tiny features, which is important for futher clinical diagnosis and therapy.

Additionally, to further illustrate the effectiveness of the MLR method for CT image reconstruction, we test our method by applying the geometric matrix to a natural image. Figure 11 show that for inverse Radon transform of natural image with apparent textures from all $15$, $20$ and $30$ projection views, the proposed method has even greater advantage comparing to the wavelet based CT reconstruction method since the traditional wavelet based method cannot distinguish the texture from the artifacts caused by the low-dose projection.

4.3 Semi-supervised Learning

Our final experiment is conducted to test the proposed MLR method for handwritten digits recognition based on the MINST data which is initially provided and processed in [32], as shown in Figure 12, including totally 70, 000 different $28\times 28$ “handwritten digits” images. As a special case of semi-supervised learning problem, we regard each image as a $784$ dimensional vector, and view all the images as a set of 70, 000 points in $\mathbb{R}^{784}$. Therefore, the vectorized images can formulate a point matrix $P\in\mathbb{R}^{784\times 70000}$. The labels $\{L(x)\}$ can possibly take the value from ${0,1,2,\ldots,9}$.

For initial purpose of MINST data, the given indices set $S$ has size $60,000$ and one need to estimate the rest $10,000$ index with lowest error. Recently, the full $70,000$ indices set can be roughly reconstructed from $50-100$ given indices and some diffusion based methods. For example, [47] proposed an initial graph Laplacian based method. Later on, [43] proposed a weighted graph Laplacian method, from which the inpainting accuracy can exceed $80\%$ from merely $70$ of given indices.

In this experiment, we apply the MLR based Algorithm 3 to this semi-supervised learning problem. In particular, we attempt to reconstruct all the $70,000$ labels of the MINST data [32] from uniformly random sampled $35,50,70,100$, and $700$ labels. For each sampling rate, we take 10 different random samples for comparisons. Figure 13 shows the success rate of label estimation by graph Laplacian (GL) [47], weighted graph Laplacian (WGL) [43], and the proposed manifold based locally low-rank approximation based model (MLR). The first five images in Figure 13 shows the success rate for each individual random sample with a fixed number of sample indices. The last image in Figure 13 shows the average success rate which is naturally monotone increasing with respect to the number of sample indices. It can be clearly observed that the proposed method has the highest accuracy of estimation for almost all the random samples. In terms of average success rate, the proposed model outperforms the previously proposed graph Laplacian and weighted graph Laplacian based methods. We remark that further improvement can be expected if special treatments for shape recognition and similarity can be conducted which will be our future work.

5 Conclusions

In this paper, we propose a manifold based low-rank regularization method for image restoration and semi-supervised learning. The proposed regularization can be viewed as a point-wise linearization of the manifold dimension, which generalize the concept of low-rank regularization for linear objects as a concept of manifold based low-rank for nonlinear objects. Using the proposed regularization, we investigate new methods of image inpaining, image super-resolution and X-ray CT image reconstruction. We further extend this method to a general data analysis problem, semi-supervised learning. Intensive numerical experiments demonstrate that the proposed MLR method is comparable to or even outperforms the existing wavelet based models [7, 18] and PDE based models [47, 43, 36].

Several directions will be investigated in our future work. For instance, the current method can be adapted to handle images with noisy input. It is also an important problem to explore a better method to pick the “local regions” or manifold representation. For example, for semi-supervised learnings, the left image in Figure 2 shows that the KNN obtained by Euclidean distance may still include some ambiguity. In particular, some KNNs may have local rank as high as 7 or 8, which reduces the reliability of local low rank regularization. Therefore, developing a data-driven approach to non-Euclidean geometry for MLR will be a very interesting direction to investigate in our future work.

Acknowledgement

We thank Prof. Stanley Osher, Prof. Zuoqiang Shi and Mr. Wei Zhu kindly share their valuable comments and codes of both LDMM and LDMM+WGL for comparisons.