Nonlocal Low-Rank Tensor Factor Analysis for Image Restoration

TL;DR

This paper introduces NLR-TFA, a novel low-rank tensor factor analysis method for image restoration that effectively captures non-local correlations, improving reconstruction especially at low sampling rates.

Contribution

It proposes a new low-rank tensor factorization approach combined with ADMM and deep learning to enhance image reconstruction in compressive sensing.

Findings

Outperforms existing methods at low CS sampling rates

Effective in both noiseless and noisy measurements

Demonstrates superior image quality in experiments

Abstract

Low-rank signal modeling has been widely leveraged to capture non-local correlation in image processing applications. We propose a new method that employs low-rank tensor factor analysis for tensors generated by grouped image patches. The low-rank tensors are fed into the alternative direction multiplier method (ADMM) to further improve image reconstruction. The motivating application is compressive sensing (CS), and a deep convolutional architecture is adopted to approximate the expensive matrix inversion in CS applications. An iterative algorithm based on this low-rank tensor factorization strategy, called NLR-TFA, is presented in detail. Experimental results on noiseless and noisy CS measurements demonstrate the superiority of the proposed approach, especially at low CS sampling rates.

Figures

Tables

Table 1. Table 1: Reconstruction PSNR (dB) and SSIM of 8 8 8 images using different CS recovery methods at different CSr.

Image	Method	CSr
		0.02	0.04	0.06	0.08	0.10
Barbara	BM3D-CS	16.44, 0.389	18.03, 0.422	19.70, 0.470	21.85, 0.561	23.70, 0.635
	TVAL3	21.32, 0.588	22.58, 0.620	23.28, 0.639	23.78, 0.654	24.25, 0.667
	NLR-CS	16.79, 0.400	18.90, 0.433	20.56, 0.501	24.00, 0.659	25.13, 0.725
	D-AMP	15.08, 0.354	17.76, 0.417	20.44, 0.496	22.66, 0.623	24.05, 0.660
	Ours	26.15, 0.812	27.29, 0.835	27.40, 0.837	28.05, 0.848	28.30, 0.852
Boats	BM3D-CS	16.76, 0.398	17.98, 0.440	20.11, 0.507	23.21, 0.624	24.65, 0.687
	TVAL3	21.58, 0.634	23.66, 0.681	24.57, 0.699	25.22, 0.713	25.75, 0.725
	NLR-CS	17.20, 0.409	19.41, 0.489	21.76, 0.559	23.76, 0.638	25.52, 0.711
	D-AMP	15.52, 0.365	18.41, 0.456	20.76, 0.528	22.85, 0.600	25.22, 0.703
	Ours	30.79, 0.898	31.73, 0.909	31.84, 0.910	32.35, 0.916	32.33, 0.912
Cameraman	BM3D-CS	16.84, 0.432	18.34, 0.496	20.62, 0.585	22.32, 0.654	23.85, 0.711
	TVAL3	20.33, 0.557	21.53, 0.579	22.02, 0.580	22.37, 0.578	23.19, 0.688
	NLR-CS	17.00, 0.436	19.10, 0.531	21.60, 0.620	23.99, 0.709	26.25, 0.789
	D-AMP	16.13, 0.401	17.51, 0.461	19.84, 0.562	21.74, 0.629	23.61, 0.705
	Ours	25.88, 0.828	26.47, 0.829	26.66, 0.823	27.16, 0.824	27.64, 0.819
Foreman	BM3D-CS	18.96, 0.610	21.06, 0.680	25.48, 0.726	29.80, 0.831	31.69, 0.866
	TVAL3	23.35, 0.777	26.81, 0.808	28.47, 0.821	29.45, 0.828	30.02, 0.833
	NLR-CS	19.64, 0.645	22.33, 0.695	28.09, 0.809	31.68, 0.865	34.15, 0.900
	D-AMP	19.42, 0.639	22.07, 0.687	26.31, 0.750	29.12, 0.813	31.90, 0.872
	Ours	34.01, 0.934	34.61, 0.936	35.18, 0.937	35.33, 0.939	36.02, 0.939
House	BM3D-CS	18.12, 0.518	20.35, 0.590	25.79, 0.755	29.78, 0.816	30.50, 0.824
	TVAL3	22.57, 0.706	24.58, 0.724	26.36, 0.740	27.63, 0.766	28.68, 0.814
	NLR-CS	18.80, 0.554	21.77, 0.620	27.16, 0.761	31.23, 0.835	33.66, 0.869
	D-AMP	17.74, 0.492	19.93, 0.577	24.00, 0.693	27.39, 0.758	29.84, 0.818
	Ours	30.42, 0.873	31.76, 0.884	33.07, 0.898	34.32, 0.910	34.80, 0.913
Lena	BM3D-CS	16.09, 0.427	18.22, 0.511	20.47, 0.603	23.55, 0.694	24.63, 0.729
	TVAL3	21.60, 0.671	23.23, 0.693	23.88, 0.697	24.41, 0.701	24.81, 0.705
	NLR-CS	16.71, 0.461	19.07, 0.545	21.36, 0.649	25.47, 0.746	25.82, 0.771
	D-AMP	15.53, 0.399	17.21, 0.485	19.93, 0.567	22.29, 0.681	24.60, 0.733
	Ours	28.51, 0.881	29.70, 0.894	29.85, 0.889	30.42, 0.898	30.65, 0.896
Monarch	BM3D-CS	14.55, 0.352	15.24, 0.390	17.65, 0.500	19.75, 0.581	21.64, 0.666
	TVAL3	17.31, 0.537	18.91, 0.577	19.82, 0.595	20.56, 0.611	21.24, 0.624
	NLR-CS	14.69, 0.357	16.03, 0.413	18.53, 0.552	21.84, 0.686	25.12, 0.796
	D-AMP	14.01, 0.334	14.86, 0.366	17.55, 0.497	19.84, 0.587	21.77, 0.675
	Ours	26.14, 0.884	27.28, 0.895	28.21, 0.901	28.97, 0.904	29.25, 0.904
Parrots	BM3D-CS	16.62, 0.506	19.32, 0.613	23.07 0.693	25.04, 0.761	26.10, 0.787
	TVAL3	21.11, 0.714	22.90, 0.722	23.47, 0.714	24.38, 0.770	25.84, 0.804
	NLR-CS	17.39, 0.550	20.84, 0.681	24.95, 0.782	28.18, 0.835	30.35, 0.875
	D-AMP	16.89, 0.520	19.55, 0.619	22.62, 0.684	24.10, 0.711	26.31, 0.792
	Ours	28.28, 0.899	28.88, 0.893	29.22, 0.890	29.55, 0.887	29.84, 0.885
Average	BM3D-CS	16.80, 0.454	18.57, 0.518	21.61, 0.605	24.42, 0.690	25.85, 0.738
	TVAL3	21.15, 0.648	23.02, 0.676	23.98, 0.685	24.72, 0.703	25.47, 0.732
	NLR-CS	17.28, 0.477	19.68, 0.551	23.00, 0.654	26.27, 0.747	28.25, 0.805
	D-AMP	16.29, 0.438	18.41, 0.509	21.43, 0.610	23.75, 0.675	25.91, 0.745
	Ours	28.77, 0.876	29.71, 0.884	30.18, 0.886	30.77, 0.891	31.10, 0.890

Table 2. Table 2: PSNR(dB) of reconstructions from measurements generated by random Gaussian sampling of 8 images at different CSr.

Image	Method	CSr
		0.02	0.04	0.06	0.08
Barbara	D-AMP	15.79	17.56	19.20	21.96
	Ours	24.13	25.52	26.38	27.22
Boats	D-AMP	16.06	17.93	19.84	22.47
	Ours	26.04	27.44	28.59	29.51
Cameraman	D-AMP	16.66	18.32	19.79	22.42
	Ours	23.59	24.66	25.41	26.15
Foreman	D-AMP	20.50	22.78	25.75	28.40
	Ours	27.72	29.23	30.66	31.47
House	D-AMP	18.61	20.85	24.42	26.97
	Ours	26.68	28.20	29.48	30.17
Lena	D-AMP	16.34	18.15	19.77	22.15
	Ours	24.63	26.04	26.99	27.78
Monarch	D-AMP	14.50	15.68	17.94	19.93
	Ours	21.57	23.01	24.19	24.89
Parrots	D-AMP	17.49	19.38	22.56	25.24
	Ours	25.92	27.19	27.96	28.54
Average	D-AMP	16.99	18.83	21.15	23.69
	Ours	25.04	26.41	27.46	28.22

Table 3. Table 3: PSNR(dB) of reconstructions with variant patch sizes at different CSr for foreman image.

Patch Size	CSr
	0.02	0.04	0.06	0.08	0.10
$4\times 4$	34.01	34.61	35.18	35.33	36.02
$6\times 6$	31.16	31.95	32.72	32.88	33.59
$8\times 8$	29.63	30.41	31.37	31.68	32.24
$10\times 10$	28.57	29.45	30.44	30.86	31.33
$12\times 12$	27.39	28.61	29.60	30.11	30.54

Equations

$$\mathbf{y}=\mathbf{\Phi }\mathbf{x}+\mathbf{e},$$

$$\mathcal{T}$$

$$\mathcal{T}(i,j,k)$$

$$\mathcal{T}$$

$$\mathcal{T}(i,j,k)$$

$$\underset{\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{\lambda }}{\min }\Vert \mathcal{T}-\sum _{r=1}^R{\lambda }_r{\mathbf{a}}_r\otimes {\mathbf{b}}_r\otimes {\mathbf{c}}_r{\Vert }_F^2$$

$$\mathbf{x}=\arg \underset{x}{\min }\Vert \mathbf{x}{\Vert }_1,s.t.\Vert \mathbf{y}-\mathbf{\Phi }\mathbf{x}\Vert \le ϵ.$$

$$\mathbf{x}=\arg \underset{x}{\min }\Vert \mathbf{y}-\mathbf{\Phi }\mathbf{x}{\Vert }_2^2+\lambda \Vert \mathbf{x}{\Vert }_1,$$

$${\mathcal{T}}_p=[{\mathbf{x}}_{p,1},{\mathbf{x}}_{p,2},...,{\mathbf{x}}_{p,k}].$$

$${\mathcal{L}}_p=\sum _{r\in S_{\ell }}{\lambda }_r{\mathbf{a}}_r\otimes {\mathbf{b}}_r\otimes {\mathbf{c}}_r,$$

$$\hat{\mathbf{x}}=\arg \underset{\mathbf{x}}{\min }\Vert \mathbf{y}-\mathbf{\Phi }\mathbf{x}{\Vert }_2^2+\eta \sum _p\Vert {\tilde{\mathcal{T}}}_p\mathbf{x}-{\mathcal{L}}_p{\Vert }_F^2,$$

$$\mathbf{x}=({\mathbf{\Phi }}^H\mathbf{\Phi }+\eta \sum _p{\tilde{\mathcal{T}}}_p^{\top }{\tilde{\mathcal{T}}}_p{)}^{-1}({\mathbf{\Phi }}^H\mathbf{y}+\eta \sum _p{\tilde{\mathcal{T}}}_p^{\top }{\mathcal{L}}_p),$$

$$(\hat{\mathbf{x}},\hat{\mathbf{z}})=\arg \underset{\mathbf{x},\mathbf{z}}{\min }$$

$$+\eta \sum _p\Vert {\tilde{\mathcal{T}}}_p\mathbf{z}-{\mathcal{L}}_p{\Vert }_F^2,$$

$${\mathbf{z}}^{t+1}$$

$${\mathbf{x}}^{t+1}$$

$${\mathbf{\mu }}^{t+1}$$

$${\mathbf{z}}^{t+1}=(\eta \sum _p{\tilde{\mathcal{T}}}_p^{\top }{\tilde{\mathcal{T}}}_p+\beta \mathbf{I}{)}^{-1}(\beta {\mathbf{x}}^t+\frac{{\mathbf{\mu }}^t}{2}+\eta \sum _p{\tilde{\mathcal{T}}}_p^{\top }{\mathcal{L}}_p).$$

$${\mathbf{x}}^{t+1}={\mathbf{U}}^{-1}({\mathbf{\Phi }}^H\mathbf{y}+\beta {\mathbf{z}}^{t+1}-\frac{{\mathbf{\mu }}^t}{2}),$$

$${\mathbf{U}}^{-1}={\beta }^{-1}(\mathbf{I}-{\mathbf{\Phi }}^H{\mathbf{V}}^{-1}\mathbf{\Phi }),$$

$$\arg \underset{\theta }{\min }{\mathbb{E}}_{ϵ}[\Vert ϵ-{\mathcal{C}}_{\theta }\mathbf{V}ϵ{\Vert }_2^2+\Vert ϵ-\mathbf{V}{\mathcal{C}}_{\theta }ϵ\Vert ],$$

$${\mathbf{x}}^{t+1}={\beta }^{-1}(\mathbf{I}-{\mathbf{\Phi }}^H{\mathcal{C}}_{\theta }\mathbf{\Phi })({\mathbf{\Phi }}^H\mathbf{y}+\beta {\mathbf{z}}^{t+1}-\frac{{\mathbf{\mu }}^t}{2}).$$

$${\mathbf{x}}^{t+1}=({\mathbf{D}}^H\mathbf{D}+\beta \mathbf{I}{)}^{-1}({\mathbf{F}}^H{\mathbf{D}}^H\mathbf{y}+(\beta {\mathbf{z}}^{t+1}-\frac{{\mathbf{\mu }}^t}{2})),$$

$${\mathbf{x}}^{t+1}={\mathbf{F}}^H\{({\mathbf{D}}^H\mathbf{D}+\beta \mathbf{I}{)}^{-1}({\mathbf{D}}^H\mathbf{y}+\mathbf{F}(\beta {\mathbf{z}}^{t+1}-\frac{{\mathbf{\mu }}^t}{2}))\}.$$

Full text

1 Introduction

Inspired by the nonlocal self-similarity of image patches [3], various image processing algorithms have been proposed to investigate the low-rank property of image patch groups [9, 11, 17, 31, 44]. In general, these methods first select a reference patch and then search for similar patches across the image to form a group. Following this, the patches in this group are vectorized and stacked to a matrix. Since these patches are similar, the constructed matrix has a low-rank property. Via performing this low-rank model on every (overlapping) patch in the image, state-of-the-art image restoration results have been achieved.

One common issue in the above algorithms is that the original two-dimensional (2D) patches are vectorized to construct the group matrix, which loses the spatial structure within the image patch. We propose a tensor based algorithm to retain this structure while still leveraging the advantages of low-rank patch models. Tensor factorization methods [19] offer a useful way to learn latent structure from complex multiway data, and have been used in image processing tasks [43]. These methods decompose the tensor data into a set of factor matrices (one for each mode or “way” of the tensor), that can be used as a latent feature representation for the objects in each of the tensor modes [33]. Recently tensor approaches have been applied in computer vision, such as image denoising [32], and video denoising [42]. We consider the general image restoration problem with specific applications to compressive sensing [12] using the tensor factor analysis (TFA) approach. In particular, instead of vectorizing the image patches, we impose low-rank TFA to the 3D image patch groups.

To be concrete, the image restoration problem aims to estimate the latent clean image $\boldsymbol{x}$ given the observation $\boldsymbol{y}$ (the corrupted/compressive measurement) and the measurement matrix $\boldsymbol{\Phi}$, which can be formulated as

[TABLE]

where ${\boldsymbol{e}}$ denotes the measurement noise and is usually modeled by ${\boldsymbol{e}}\sim\mathcal{N}(0,\sigma^{2}{\bf I})$. The problem has different variants for different $\boldsymbol{\Phi}$: when $\boldsymbol{\Phi}$ is the identity matrix, this is a denoising problem [16]; when $\boldsymbol{\Phi}$ is a diagonal matrix whose diagonal elements are either 1 or 0, keeping or removing corresponding pixels, this becomes image inpainting [63]; and when $\boldsymbol{\Phi}\in\mathbb{C}^{M\times N}$ and $M\ll N$, this is a compressive sensing (CS) problem.

We focus on the CS problem [14, 28], and consider the image reconstruction with a very limited number of compressive measurements [50, 51, 52, 58]. A framework of the proposed NLR-TFA method is shown in Fig. 1. There has been recent interest in using deep learning technique for CS problems [30]. A reconstruction network is developed in [20, 58], and an ADMM-net is proposed in [47]; while the former paper focused on block-wise CS, the latter one only considered CS-MRI [28]. By contrast, we develop a tensor-based framework for general CS applications. Further, to overcome the large-scale matrix inversion in CS, we further propose to pre-train a deep convolutional neural network to approximate the expensive matrix inversion operator [41].

The remainder of this paper is organized as follows. Section 2 presents background knowledge on tensors and CS. Section 3 derives our algorithm. Extensive results are provided in Section 4 to demonstrate the superiority of the proposed algorithm relative to other leading approaches, and Section 5 concludes the paper.

2 Background

2.1 Tensor Decomposition

Tensors (of order higher than two) are arrays indexed by three or more indices. Specifically, matrices are two-way arrays, and tensors are three- or higher-way arrays. In the sequel, we mainly focus on $3\text{-}$way tensors and everything naturally generalizes to higher-order tensors.

A $rank\text{-}1$ $3\text{-}$way tensor $\mathcal{T}\in\mathbb{R}^{I\times J\times K}$ is an outer product of $3$ unit vectors ${\boldsymbol{a}}\in\mathbb{R}^{I}$, ${\boldsymbol{b}}\in\mathbb{R}^{J}$, and ${\boldsymbol{c}}\in\mathbb{R}^{K}$:

[TABLE]

where $\otimes$ is the vector outer product and $\lambda$ is a scalar.

The $rank$ of a tensor $\mathcal{T}$ is the minimum number of $rank\text{-}1$ tensors that sum to $\mathcal{T}$. Therefore, a $rank\text{-}R$ $3\text{-}$way tensor $\mathcal{T}$ can be written as

[TABLE]

where unit vectors ${\boldsymbol{a}}_{r}$, ${\boldsymbol{b}}_{r}$, and ${\boldsymbol{c}}_{r}$ are tensor factors and $\lambda_{r}$ is a scalar that evaluates the significance of tensor factors ${\boldsymbol{a}}_{r},{\boldsymbol{b}}_{r},{\boldsymbol{c}}_{r}$. For brevity, $\mathcal{T}=\sum_{r=1}^{R}\lambda_{r}{\boldsymbol{a}}_{r}\otimes{\boldsymbol{b}}_{r}\otimes{\boldsymbol{c}}_{r}$ can be denoted as $\mathcal{T}=\langle{\boldsymbol{\lambda}},{{\bf A}},{\bf B},{\bf C}\rangle$, where ${\boldsymbol{\lambda}}:=[\lambda_{1},...,\lambda_{R}]$, ${{\bf A}}:=[{\boldsymbol{a}}_{1},...,{\boldsymbol{a}}_{R}]$, ${\bf B}:=[{\boldsymbol{b}}_{1},...,{\boldsymbol{b}}_{R}]$, and ${\bf C}:=[{\boldsymbol{c}}_{1},...,{\boldsymbol{c}}_{R}]$. When the rank $R$ is low, the number of elements needed to represent $\mathcal{T}$ is dramatically decreased after tensor decomposition, i.e., from $I\times K\times J$ to $(I+K+J)\times R$.

Practically, most tensor decomposition problems are NP-hard. However, in most real applications, as long as the tensor does not have too many components, and the components are not adversarially chosen, tensor decomposition can be computed in polynomial time [34]. The tensor-decomposition problem seeks to estimate ${{\bf A}}$, ${\bf B}$, ${\bf C}$ and coefficients ${\boldsymbol{\lambda}}$ from a tensor $\mathcal{T}$. Adopting a least-squares fitness criterion, the problem is

[TABLE]

In this work, we employ Jenrich’s algorithm [21] to solve this problem. In Algorithm 1, “$+$” denotes pseudo-inverse of a matrix.

2.2 Compressive Sensing

Compressive sensing is a signal acquisition technique that enables sampling a signal at sub-Nyquist rates [12]. In CS, a reconstruction algorithm is used to recover the original high-dimensional signal from a small number of random linear measurements. Taking compressive measurements can be viewed as a process of linearly mapping an $N$-dimensional signal vector $\boldsymbol{x}$ to an $M$-dimensional measurement vector $\boldsymbol{y}$, $M\ll N$, using a measurement matrix $\boldsymbol{\Phi}\in\mathbb{C}^{M\times N}$, i.e., $\boldsymbol{y}=\boldsymbol{\Phi}\boldsymbol{x}$. Since the matrix $\boldsymbol{\Phi}$ is rank-deficient, there exists more than one $\boldsymbol{x}$ that yields the same measurement $\boldsymbol{y}$. In this paper, the compressive sensing rate (CSr) is defined as CSr$=M/N$.

To recover $\boldsymbol{x}$, one searches for the vector that possesses a certain structure among all the vectors $\boldsymbol{x}$ that satisfy $\boldsymbol{y}\approx\boldsymbol{\Phi}\boldsymbol{x}$. In the case of a sparse $\boldsymbol{x}$, a popular method is to solve the optimization problem

[TABLE]

This problem is convex and known formally as basis pursuit denoising (BPDN). It has been shown in [6][13] that if $\boldsymbol{x}$ is sufficiently sparse and $\boldsymbol{\Phi}$ satisfies certain properties, then the $s$-sparse signal can be accurately recovered from $m=\mathcal{O}(s\log(n/s))$ random linear measurements [5]. Equation (5) can be equivalently translated to the following unconstrained optimization problem

[TABLE]

where $\lambda$ is a regularization parameter. Various methods can be used to solve the above minimization problem [49, 59, 57], and in this work we adopt the alternative direction multiplier method (ADMM) framework [2, 23]. Specifically, we consider the application of image CS [14, 51, 52], and beyond sparsity, we propose a new tensor factorization approach to exploit the high-order structure in the image patches, seeking high reconstruction performance at extremely low CSr, e.g., CSr<0.05. Refer to Fig. 3 for one example of reconstructed image using our proposed algorithm, compared with other leading algorithms at CSr$=0.02$ (with image size $256\times 256$).

3 Method

We propose a new model that recovers compressively sensed images using low-rank tensor factor analysis and ADMM. First, we generate a tensor from the estimated image based on patch grouping. Then we impose low rankness on the tensor after tensor decomposition. This new low-rank tensor is fed into a global objective function which is solved by ADMM. These two steps are iteratively performed until satisfying some criterion. The complete algorithm is shown in Algorithm 2.

3.1 Patch-Based Low-Rank Tensor Factorization

Given the observation $\boldsymbol{y}$, we first obtain the estimated image $\hat{\boldsymbol{x}}$ using some fast algorithm, e.g., the DCT or wavelet based algorithm [18, 59]. In denoising, $\hat{\boldsymbol{x}}$ can be simply set equal to $\boldsymbol{y}$. Then the estimated image is divided into $P$ overlapping patches $\{\boldsymbol{x}_{1},...,\boldsymbol{x}_{P}\}$. The basic assumption underlying the proposed approach is that a sufficient number of similar patches can be found for any reference patch, a.k.a, the nonlocal self-similarity (NSS) prior [3]. For each reference patch $\boldsymbol{x}_{p}$, we perform a $k$-nearest-neighbor search for the nonlocal similar patches across the image to form a group $\{\boldsymbol{x}_{p,1},\boldsymbol{x}_{p,2},...,\boldsymbol{x}_{p,k}\}$, where $k$ is the number of similar patches (including the reference patch itself). Here, the Euclidean distance of pixel intensity is used as the metric to group patches. By concatenating the grouped patches on the third dimension in ascending order of Euclidean distances, we generate a tensor $\mathcal{T}_{p}$ for reference patch $\boldsymbol{x}_{p}$

[TABLE]

As $\boldsymbol{x}_{p}$ has zeros distance to itself, it is always found as the leading patch in $\mathcal{T}_{p}$, i.e., $\mathcal{T}_{p}(:,:,1)=\boldsymbol{x}_{p}$. Eventually we have $P$ tensors and each tensor corresponds to a reference patch. The coordinates of the grouped patches are also recorded for later image aggregation. Suppose the size of each patch is $m\times n$, then the size of generated tensor is $\mathcal{T}_{p}\in{\mathbb{R}}^{m\times n\times k}$. It has been shown that the grouped patches can be denoised by low-rank approximation [4]. In this work, the low-rankness is imposed on the tensor by taking the most significant tensor factors after tensor decomposition

[TABLE]

where $S_{\ell}$ selects the most significant $\ell$ tensor factors of $\mathcal{T}_{p}$. Therefore, $\mathcal{L}_{p}$ has rank $\ell$. The significance of tensor factors is evaluated by $\lambda_{r}$’s in (3).

Under the assumption that the grouped patches have similar structure, $\mathcal{T}_{p}$ has a low-rank property, which ensures that $\mathcal{T}_{p}$ can be represented by a relatively low-rank tensor $\mathcal{L}_{p}$. This low-rank imposition shares the same spirit with the hard thresholding algorithm [1]. In Fig. 2, we impose low rankness on tensors with size $8\times 8\times 20$ generated from a clean image using the above approach, and then aggregate the tensors back to images. As can be seen, reconstructed images with low-rank tensors can accurately approximate the original image. The images aggregated from $rank$-$3$ to $rank$-$8$ tensors are highly similar (with PSNR merely increased by $0.55$ dB), which indicates that in this case $rank$-$3$ tensors are adequate to represent the original image.

3.2 Image Recovery via ADMM

With the reconstructed low-rank tensors $\mathcal{L}_{p}$, the following optimization problem is proposed for CS recovery:

[TABLE]

where $\tilde{\mathcal{T}}_{p}\boldsymbol{x}=\mathcal{T}_{p}$ denotes the tensor formed for each reference patch $\boldsymbol{x}_{p}$ and $\eta$ is a regularization parameter. The closed-form solution for this quadratic optimization problem is

[TABLE]

where $H$ is the Hermitian transpose operation, $\sum_{p}\tilde{\mathcal{T}}_{p}^{\top}\mathcal{L}_{p}$ denotes the results of averaging all of the similar patches for each reference patch, and each entry of $\sum_{p}\tilde{\mathcal{T}}_{p}^{\top}\tilde{\mathcal{T}}_{p}$ corresponds to an image pixel location whose value is the number of overlapping patches that cover this pixel location. In (10), the matrix to be inverted can be large, for which conjugate gradient descent (CG) is usually applied [8]. We adopt ADMM to solve this problem, introducing auxiliary variable $\boldsymbol{z}$. Applying ADMM to (9), we obtain the global objective function

[TABLE]

where ${\boldsymbol{\mu}}$ is the Lagrange multiplier, and $\beta>0$ is the penalty parameter. Instead of minimizing $\boldsymbol{x}$ and $\boldsymbol{z}$ simultaneously, ADMM decomposes the problem into two subproblems that minimizes w.r.t $\boldsymbol{x}$ and $\boldsymbol{z}$, respectively. More specifically, the optimization problem in (3.2) consists of the following iterations:

[TABLE]

Both $\boldsymbol{x}$ and $\boldsymbol{z}$ admit closed-form solutions. For fixed $\boldsymbol{x}^{t}$ and ${\boldsymbol{\mu}}^{t}$,

[TABLE]

Then we use the updated $\boldsymbol{z}^{t+1}$ to update $\boldsymbol{x}$,

[TABLE]

where ${{\bf U}}^{-1}=(\boldsymbol{\Phi}^{H}\boldsymbol{\Phi}+\beta{\bf I})^{-1}$.

3.2.1 Pretrained Deep Convolutional Architectures

In (16), the term ${{\bf U}}^{-1}$ involves an expensive matrix inversion, which makes direct computation of $\boldsymbol{x}^{t+1}$ impractical. An inversion-free approach proposed by [41] addresses this problem by learning a convolutional neural network to approximate the matrix inversion. Note that training this neural network is data-independent since ${{\bf U}}^{-1}$ is only dependent on $\boldsymbol{\Phi}$. Applying the Sherman-Morrison-Woodbury formula, we reduce this matrix inversion to a smaller scale,

[TABLE]

where ${{\bf V}}^{-1}=(\beta{\bf I}+\boldsymbol{\Phi}\boldsymbol{\Phi}^{H})^{-1}$ has the dimension $m\times m$. To approximate ${{\bf V}}^{-1}$, a trainable deep convolutional neural network $\mathcal{C}_{\theta}$ parameterized by $\theta$ is employed, i.e., $\mathcal{C}_{\theta}\approx{{\bf V}}^{-1}$. $\mathcal{C}_{\theta}$ is learned by minimizing the sum of two reconstruction losses of two auto-encoders with shared weights

[TABLE]

where ${{\bf V}}=\beta{\bf I}+\boldsymbol{\Phi}\boldsymbol{\Phi}^{H}$ can be computed directly, and $\epsilon$ is sampled from publicly available image datasets [10]. Note that $\mathcal{C}_{\theta}$ is pretrained for different $\boldsymbol{\Phi}$ and $\beta$. By plugging the learned $\mathcal{C}_{\theta}$ into (17), we obtain a reusable term ${{\bf U}}^{-1}=\beta^{-1}({\bf I}-\boldsymbol{\Phi}^{H}\mathcal{C}_{\theta}\boldsymbol{\Phi})$ as the replacement for the cumbersome inversion matrix. Hence $\boldsymbol{x}^{t+1}$ is updated by

[TABLE]

3.2.2 Fourier space solution

Equation (16) can be solved by transforming the problem from the image space into the Fourier space when $\boldsymbol{\Phi}$ is a partial Fourier transform matrix [11]. For a down-sampling matrix ${\bf D}$ and a Fourier transform matrix ${{\bf F}}$, $\boldsymbol{\Phi}={\bf D}{{\bf F}}$ is substituted into (16)

[TABLE]

where the inverse matrix $({\bf D}^{H}{\bf D}+\beta{\bf I})^{-1}$ is a diagonal matrix and thus can be computed easily. Equation (20) is equivalent to

[TABLE]

Therefore, $\boldsymbol{x}^{t+1}$ can be obtained by applying inverse Fourier transform to terms in the brackets of the right hand side of (21).

4 Experiments

We conduct experiments on both noisy and noiseless CS mearurements of $8$ test images: Barbara, Boats, Cameraman, Foreman, House, Lena, Monarch, and Parrots. All images are resized to $256\times 256$. Since excellent results have been obtained when the CS rate is large, i.e. ${\rm CSr}>0.1$ [11, 29], we here focus on testing cases with limited number of measurements (${\rm CSr}<0.1$). The CS measurements are generated by pseudo-radial sampling of the test images in the Fourier domain. Unlike traditional random sampling schemes, pseudo-radial sampling produces streaking artifacts which are more difficult to remove. We also perform experiments using an $M\times N$ random sensing matrix generated by a standard Gaussian distribution. The deep convolutional architecture $C_{\theta}$ is pretrained by the sensing matrix then used in ADMM to solve the expensive matrix inversion. A DCT based CS recovery algorithm [52] is used for the initial image estimate. The main parameters of the proposed method are set as follows: patch size $4\times 4$, number of similar patches for each reference patch $50$, and the imposed low rank is $20$. The regularization parameter $\eta$ and the ADMM parameter $\beta$ are tuned separately for each sensing rate. In practice, we have found that ADMM converges fast after a few iterations and the performance gain is mainly from our low-rank TFA. Thus we set the number of outer-loop iterations to $I=50$ and the number of inner loop iterations $J=2$. Experimental results for noiseless CS measurements and noisy CS measurements are evaluated by peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) [40].

4.1 Baselines

We compare our results with four competitive CS image restoration algorithms: BM3D-CS [15], TVAL3 [22], NLR-CS [11] and D-AMP [29]. BM3D-CS is a nonparametric method that applies block matching and 3D filtering algorithm [9] on compressive sensing. TVAL3 restores images by combining the classic augmented Lagrangian multiplier method with total-variation regularization. NLR-CS is currently a state-of-the-art method that uses a nonlocal low-rank regularization method along with ADMM to solve image CS problems. D-AMP employs a denoiser in an approximate message passing framework, achieving state-of-the-art performance on noisy measurements, especially at high CS rates. The source codes of these baseline methods are downloaded from the respective author’s website and parameters for these algorithms are set to their default values.

4.2 Experiments with Noiseless CS Measurements

We first perform CS image restoration experiments from noiseless measurements using pseudo-radial sampling scheme at five different CS rates, i.e., CSr = {$0.02$, $0.04$, $0.06$, $0.08$, $0.10$}. Table 1 summarizes the results of our proposed algorithm compared with various CS inversion algorithms at different CSr values. Both BM3D-CS and the state-of-the-art methods NLR-CS and D-AMP suffer at extremely low CSr. As CSr increases, NLR-CS yields significant performance improvements, which leads one of its results at higher CSr ($0.10$) to surpass other methods. Our proposed method achieves the best performance in all cases when CSr$<0.1$. On average, our proposed NLR-TFA algorithm outperforms all other competing methods at CSr = {$0.02$, $0.04$, $0.06$, $0.08$, and $0.10$}. The PSNR gains of our NLR-TFA over BM3D-CS, TVAL3, NLR-CS and D-AMP can be as much as $15.04$ dB, $10.65$ dB, $14.36$ dB, and $14.58$ dB, respectively. Furthermore, it can be observed that our average reconstruction PSNR only decreases by $2.33$ dB as CSr decreases from $0.10$ to $0.02$, while this number for NLR-CS and D-AMP are $10.67$ dB and $9.62$ dB, respectively, indicating that the proposed method is very stable at low CSr, i.e., with limited number of measurements.

To evaluate the reconstruction visually, two examples of restored Foreman and Monarch images at CSr of $0.02$ and $0.06$ are shown in Figs. 3 and 4. It is evident that our method recovers the best visual quality among all competiting methods. Large-scale sharp edges and small-scale fine structures are both reconstructed in two images. In particular, at extremely low CSr of $0.02$, our NLR-TFA method (Fig. 3(d)) can effectively approximate the original image while some other methods, such as NLR-CS, can only reconstruct scratches (Fig. 3(c)).

We next generate CS measurements by sampling the test images with a random Gaussian sensing matrix. A deep convolutional neural network $\mathcal{C}_{\theta}$ is pretrained using the sensing matrix to solve the large-scale matrix inversion problem in ADMM. We compare our results with D-AMP which provides random Gaussian sensing implementations in their sourcing code. The reconstruction results at CSr of $0.02$, $0.04$, $0.06$ and $0.08$ are shown in Table 2. Our proposed method outperforms D-AMP on all test images and sensing rates. In addition, the runtime of the algorithm using this inversion-free approach is also much faster than using approaches with inner-loop updates, such as conjugate gradient (CG). This demonstrates the efficiency and accuracy of the deep convolutional architecture on approximating the inversion matrix.

Patch Size Selection

We have conducted experiments using different patch sizes. As can be seen from Table 3 (in the pseudo-radial sampling scheme), smaller patch sizes lead to better reconstruction results. When we adopt patch size of $4\times 4$ on the foreman image, the PSNR is significantly increased compared to results with larger patch sizes. Similar conclusions are found on other images as well. In addition, the running time of our model is largely influenced by the tensor size. For tensor size of $4\times 4\times 50$ (image patch size $4\times 4$ and $50$ similar patches for each reference patch), the average running time at CSr=$0.06$ is about 10 minutes on a desktop with i7 CPU @3.4GHz and 24G RAM.

4.3 Experiments with Noisy CS Measurements

Similar experiments are conducted with noisy CS measurements to show the robustness of our algorithm to noise. The noisy CS measurements are obtained by adding random zero-mean Gaussian noise to the noiseless CS measurements. We test the algorithms at three levels of noise, with standard deviation $\sigma$={$10$, $20$, $30$}. The PSNR comparison of all methods for Boats and Monarch images at sensing rates of $0.10$ and $0.04$, respectively, are shown in Fig. 6. Our proposed NLR-TFA outperforms other competing methods at all levels of noise. Furthermore, D-AMP shows great robustness to noise while BM3D-CS and NLR-CS are relatively sensitive to high-level noise. As the noise level increases, the reconstruction performance of our algorithm decreases slowly. This shows that the proposed method is robust to noise. In Fig. 5, we show the reconstructions of Boats images using various algorithms at CSr=$0.10$ and $\sigma=30$. The proposed NLR-TFA clearly yields the best reconstruction and is robust to noise.

5 Conclusion

We have presented a low-rank tensor-factor-analysis-based approach to solve image-restoration problems. A tensor is generated by concatenating similar patches from the estimated image for each exemplar patch. Low-rank is imposed on the tensors to exploit non-local correlation and the high-order structure information via tensor factorization. ADMM is employed on low-rank tensors, where we use either a pretrained convolutional architecture or Fourier approaches to solve a matrix inversion. Experimental results demonstrate that the proposed NLR-TFA method outperforms state-of-the-art algorithms on CS image reconstruction at low CS sampling rates.

We have demonstrated the superiority of our model in image CS with three-way tensors, which can also be used in depth CS [25, 26, 53], polarization CS  [39] and dynamic range CS [48]. Future work includes extending the model to four-way tensors for video CS [24, 27, 35, 36, 45, 46, 54, 55, 56, 60, 61], hyperspectral image CS [7, 62] and to five-way tensors for joint spectral-temporal CS [37, 38].